Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Updates for handling CrIS SWIR data #587

Closed
wants to merge 17 commits into from
Closed

Updates for handling CrIS SWIR data #587

Show file tree
Hide file tree
Changes from 15 commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
d94604a
updates for CrIS SW and LW-SW QC and CrIS SW scene dependent obs erro…
erinjones2 Feb 11, 2021
c8013a1
merge with develop and add updates to setuprad
Mar 2, 2023
60ebc53
adding some comments/documentation to setuprad, qcmod, crtm_interface…
Mar 3, 2023
21eb45d
fix a line in crtm_interface.f90
Mar 7, 2023
488cf6e
Issue #546 - modify read_cris to add a surface sensitive MW channel f…
Mar 10, 2023
8478f9d
Issue #546 - fix sfc_channel variable declaration syntax in read_cris
Mar 10, 2023
7dc9f4e
modifying MW surface channel in read_cris to be one that is actively …
Apr 24, 2023
4d4d7d2
Merge remote-tracking branch 'upstream/develop' into master
Apr 25, 2023
40aad5c
keeping up with develop and adding a small change to setuprad for GSI…
May 15, 2023
759f80f
resolving conflicts with develop; GSI issue #546
May 16, 2023
cac4bc3
Merge remote-tracking branch 'upstream/develop'
erinjones2 May 24, 2023
ff91f56
bug fix for setuprad plus keeping up with develop - GSI Issue #546
erinjones2 May 26, 2023
a88ff51
Merge remote-tracking branch 'upstream/develop' into master
Jun 26, 2023
10fe207
Merge branch 'master' of https://github.com/erinjones2/GSI into master
Jun 26, 2023
08af125
Merge remote-tracking branch 'upstream/develop' into master
Jun 27, 2023
2a4c143
Merge remote-tracking branch 'upstream/develop' into master
Jul 5, 2023
0131a0b
generalizing sun glint check for hyperspectral IR SW channels in qcmod
Jul 5, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
66 changes: 66 additions & 0 deletions src/gsi/crtm_interface.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,12 +40,16 @@ module crtm_interface
! 2019-03-13 eliu - add quality control to identify areas with cold-air outbreak
! 2019-03-13 eliu - add calculation of GFDL cloud fraction
! 2019-03-22 Wei/Martin - Added VIIRS AOD capability alongside MODIS AOD
! 2019-07-24 ejones - add get_crtm_temp_tl subroutine to get tangent linear for obs Tb using
! crtm outputs

!
!
! subroutines included:
! sub init_crtm
! sub call_crtm
! sub destroy_crtm
! sub get_crtm_temp_tl
!
! attributes:
! language: f90
Expand Down Expand Up @@ -83,6 +87,7 @@ module crtm_interface
public init_crtm ! Subroutine initializes crtm for specified instrument
public call_crtm ! Subroutine creates profile for crtm, calls crtm, then adjoint of create
public destroy_crtm ! Subroutine destroys initialization for crtm
public get_crtm_temp_tl ! Subroutine computes tangent linear for obs Tb
public sensorindex
public surface
public isatid ! = 1 index of satellite id
Expand Down Expand Up @@ -3199,4 +3204,65 @@ subroutine get_lai(data_s,nchanl,nreal,itime,ilate,lai_type,lai)
return
end subroutine get_lai

subroutine get_crtm_temp_tl(nchanl,sc_index,tb_obs,error0,tref,tl_tbobs)
!$$$ subprogram documentation block
! . . . .
! subprogram: get_crtm_temp_tl
!
! prgmmr: ejones date: 2019-07-24
!
! abstract: get tangent linear temperature from the CRTM to appropriately
! inflate observation error when applicable (e.g. for shortwave
! hyperspectral IR channels, which have higher instrument noise at
! cooler temperatures).
!
! program history log:
! 2019-07-24 ejones
!
! input argument list:
! nchanl - number of channels for sensor
! sc_index - index of channel numbers in the full 2211 channel set
! tref - reference temperature to compute reference radiance_tl at
! tb_obs - brightness temperature array
! error0 - array of original obs errors from satinfo file (specified
! for reference temperature)
!
! output argument list:
! tl_tbobs - array of scene/temperature dependent obs errors
!
! attributes:
! language: f90
!
!$$$
!----------
use crtm_planck_functions, only: crtm_planck_radiance_tl, crtm_planck_radiance, &
crtm_planck_temperature_tl

implicit none

! passed variables
integer(i_kind) ,intent(in) :: nchanl
real(r_kind) ,intent(in) :: tref ! reference temperature (K)
integer(i_kind),dimension(nchanl) ,intent(in) :: sc_index ! channel number in 2211 set
real(r_kind),dimension(nchanl) ,intent(in) :: tb_obs ! tbs
real(r_kind),dimension(nchanl) ,intent(in) :: error0 ! error from satinfo
real(r_kind),dimension(nchanl) ,intent(out) :: tl_tbobs ! new error to use

! local variables
integer(i_kind) :: i
integer(i_kind) :: crtm_sensorindex ! CRTM Sensor_Index
integer(i_kind) :: crtm_channelindex ! CRTM Channel_Index
real(r_kind) :: tl_refrad,obsrad ! computed radiances and tl_radiances

i=0
crtm_sensorindex=channelinfo(1)%sensor_index
do i=1,nchanl
crtm_channelindex=channelinfo(1)%channel_index(sc_index(i))
call crtm_planck_radiance_tl(crtm_sensorindex,crtm_channelindex,tref,error0(i),tl_refrad)
call crtm_planck_radiance(crtm_sensorindex,crtm_channelindex,tb_obs(i),obsrad)
call crtm_planck_temperature_tl(crtm_sensorindex,crtm_channelindex,obsrad,tl_refrad,tl_tbobs(i))
end do
return
end subroutine get_crtm_temp_tl

end module crtm_interface
113 changes: 94 additions & 19 deletions src/gsi/qcmod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,7 @@ module qcmod
! for all variables
! 2019-03-27 h. liu - add ABI QC
! 2019-06-10 h. liu - add Geostationary satellites CSR data QC to replace qc_abi,qc_seviri
! 2019-07-20 ejones - mods to qc_irsnd for CrIS, add sun glint QC flag for CrIS
! 2019-09-29 X.Su - add troflg and lat_c for hilbert curve tunning
! 2019-04-19 eliu - add QC flag for cold-air outbreak
! 2021-04-29 Jung/Collard - Fix numerics for emissivity check
Expand Down Expand Up @@ -305,6 +306,9 @@ module qcmod
integer(i_kind),parameter:: ifail_satzen_qc=52
! Reject because of surface emissivity/temperature influence in subroutine qc_irsnd
integer(i_kind),parameter:: ifail_sfcir_qc=53
! Reject because of sun glint in subroutine qc_irsnd
integer(i_kind),parameter:: ifail_glint_irqc=54


! QC_AMSUA
! Reject because factch6 > limit in subroutine qc_amsua
Expand Down Expand Up @@ -2065,10 +2069,10 @@ subroutine qc_saphir(nchanl,sfchgt,luse,sea, &
return
end subroutine qc_saphir

subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &
cris, hirs, zsges,cenlat,frac_sea,pangs,trop5,zasat,tzbgr,tsavg5,tbc,tb_obs,tbcnob,tnoise, &
wavenumber,ptau5,prsltmp,tvp,temp,wmix,emissivity_k,ts, &
id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole)
subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &
cris,hirs,zsges,cenlat,frac_sea,pangs,trop5,zasat,solazi,satazi,tzbgr,tsavg5, &
tbc,tb_obs,tbcnob,tnoise,wavenumber,ptau5,prsltmp,tvp,temp,wmix,emissivity_k,ts, &
id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole,cris_sw)
! id_qc,aivals,errf,varinv,varinv_use,cld,cldp,kmax,zero_irjaco3_pole,radmod) ! all-sky

!$$$ subprogram documentation block
Expand All @@ -2084,8 +2088,10 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &
! 2010-08-10 derber transfered from setuprad
! 2011-08-20 zhu add cloud qc for passive channels based on the cloud
! level determined by channels with irad_use=1 and 0
! 2015-03-26 mkim add extra qc for sfc sensitive channels
! 2015-03-26 mkim add extra qc for sfc sensitive channels
! These qc are optional.(On/off by depending on the number in satinfo table)
! 2019-07-20 ejones Add sun glint check for CrIS SW obs
! 2019-07-24 ejones Add cris_sw logic so observations aren't counted twice in aivals stats
!
! input argument list:
! nchanl - number of channels per obs
Expand All @@ -2102,9 +2108,11 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &
! zsges - elevation of guess
! cenlat - latitude of observation
! frac_sea - fraction of grid box covered with water
! pangs - solar zenith angle
! pangs - solar zenith angle (deg)
! trop5 - tropopause pressure
! zasat - satellite zenith angle
! zasat - satellite zenith angle (rad)
! solazi - solar azimuth angle (deg)
! satazi - satellite azimuth angle (deg)
! tzbgr - Tz over water
! tsavg5 - surface skin temperature
! tbc - simulated - observed BT with bias correction
Expand All @@ -2123,6 +2131,7 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &
! errf - criteria of gross error
! varinv - observation weight (modified obs var error inverse)
! varinv_use - observation weight used(modified obs var error inverse)
! cris_sw - logical for cris sw case
!
! output argument list:
! id_qc - qc index - see qcmod definition
Expand All @@ -2146,14 +2155,15 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &

! Declare passed variables

logical, intent(in ) :: sea,land,ice,snow,luse,goessndr, cris, hirs
logical, intent(inout) :: zero_irjaco3_pole
logical, intent(in ) :: sea,land,ice,snow,luse,goessndr,cris,hirs
logical, intent(inout) :: zero_irjaco3_pole,cris_sw
integer(i_kind), intent(in ) :: nsig,nchanl,ndat,is
integer(i_kind),dimension(nchanl), intent(in ) :: ich
integer(i_kind),dimension(nchanl), intent(inout) :: id_qc
integer(i_kind),dimension(nchanl), intent(in ) :: kmax
real(r_kind), intent(in ) :: zsges,cenlat,frac_sea,pangs,trop5
real(r_kind), intent(in ) :: tzbgr,tsavg5,zasat
real(r_kind), intent(in ) :: solazi,satazi
real(r_kind), intent( out) :: cld,cldp
real(r_kind),dimension(40,ndat), intent(inout) :: aivals
real(r_kind),dimension(nchanl), intent(in ) :: tbc,emissivity_k,ts,wavenumber,tb_obs,tbcnob
Expand All @@ -2168,6 +2178,7 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &


real(r_kind) :: demisf,dtempf,efact,dtbf,term,cenlatx,sfchgtfact
real(r_kind) :: pangs_rad,solazi_rad,satazi_rad,relazi_rad,glint_rad,glint
real(r_kind) :: sum,sum2,sum3,cloudp,tmp,dts,delta
real(r_kind),dimension(nchanl) :: dtb
integer(i_kind) :: i,j,k,kk,lcloud,m
Expand All @@ -2181,10 +2192,31 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &
! solar zenith angle tiny_r_kind
irday = 1
if (pangs <= 89.0_r_kind .and. frac_sea > zero) then
if(.not. cris_sw) then
! QC2 in statsrad
if(luse)aivals(9,is) = aivals(9,is) + one
if(luse)aivals(9,is) = aivals(9,is) + one
endif
do i=1,nchanl
if(wavenumber(i) > r2000)then
! check for sun glint for CrIS at wavenumbers shortward of 2386.88
if(cris)then
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there a reason for this if (cris) logic? Wouldn't the sun-glint test be applicable to all ir sounders?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It should be. We were focused only on looking at and making modifications for CrIS SW, hence the if (cris) logic. I can modify the code to make this generic for all SW, if that's preferred.

if(wavenumber(i) > 2386.88)then
! calculate sun glint
! pangs, solazi, satazi need conversion to radians
! zasat passed from reader in radians, take abs avalue
pangs_rad=pangs*deg2rad
solazi_rad=solazi*deg2rad
satazi_rad=satazi*deg2rad
relazi_rad=(180.0+solazi-satazi)*deg2rad
glint_rad=acos(cos(abs(zasat))*cos(pangs_rad)+sin(abs(zasat))*sin(pangs_rad)*cos(relazi_rad))
glint=glint_rad*rad2deg
! QC low peaking CrIS SW obs with sun glint less than or equal to 15deg
if (glint <= 15.0) then
varinv(i)=zero
varinv_use(i)=zero
if(id_qc(i) == igood_qc)id_qc(i)=ifail_glint_irqc
endif
endif
else if(wavenumber(i) > r2000)then
if(wavenumber(i) > r2400)then
varinv(i)=zero
varinv_use(i)=zero
Expand Down Expand Up @@ -2223,7 +2255,7 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &
if (qc_noirjaco3_pole .and. (abs(cenlat)>r60)) zero_irjaco3_pole=.true.

! If GOES and lza > 60. do not use
if( goessndr .and. zasat*rad2deg > r60) then
if( goessndr .and. zasat*rad2deg > r60 .and. (.not. cris_sw)) then
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why would cris_sw be true for GOES? Or do you want to generalise this for all IR instruments?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cris_sw shouldn't be true for GOES; this is a redundant check because I'm paranoid. It can be removed if preferred.

! QC5 in statsrad
if(luse)aivals(12,is) = aivals(12,is) + one
do i=1,nchanl
Expand All @@ -2237,7 +2269,9 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &
sfchgtfact=one
if (zsges > r2000) then
! QC1 in statsrad
if(luse)aivals(8,is) = aivals(8,is) + one
if (.not. cris_sw) then
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why does the SW changes affect the QC decisions over high topography?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It doesn't; this is included so obs aren't counted more than once, since qc_irsnd is called twice in setuprad (calling qc_irsnd twice pertains to the LW vs SW cloud detection... see cloud detection comment)

if(luse)aivals(8,is) = aivals(8,is) + one
endif
sfchgtfact = (r2000/zsges)**4
endif

Expand Down Expand Up @@ -2309,7 +2343,15 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &
j=ich(i)
if (passive_bc .and. iuse_rad(j)==-1) then
if (lcloud .ge. kmax(i)) then
if(luse)aivals(11,is) = aivals(11,is) + one
if(.not.cris) then
if(luse)aivals(11,is) = aivals(11,is) + one
else if(cris) then
if(.not. cris_sw .and. wavenumber(i) < 2000.0) then
if(luse)aivals(11,is) = aivals(11,is) + one
else if(cris_sw .and. wavenumber(i) >= 2000.0) then
if(luse)aivals(11,is) = aivals(11,is) + one
end if
end if
varinv(i) = zero
if(id_qc(i) == igood_qc)id_qc(i)=ifail_cloud_qc
cycle
Expand All @@ -2321,7 +2363,15 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &

if ( ptau5(lcloud,i) > 0.02_r_kind) then
! QC4 in statsrad
if(luse)aivals(11,is) = aivals(11,is) + one
if(.not.cris) then
if(luse)aivals(11,is) = aivals(11,is) + one
else if(cris) then
if(.not. cris_sw .and. wavenumber(i) < 2000.0) then
if(luse)aivals(11,is) = aivals(11,is) + one
else if(cris_sw .and. wavenumber(i) >= 2000.0) then
if(luse)aivals(11,is) = aivals(11,is) + one
end if
end if
varinv(i) = zero
if(id_qc(i) == igood_qc)id_qc(i)=ifail_cloud_qc
end if
Expand Down Expand Up @@ -2350,7 +2400,18 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &
delta=max(r0_05*tnoise(i),r0_02)
if(abs(dts*ts(i)) > delta)then
! QC3 in statsrad
if(luse .and. varinv(i) > zero) aivals(10,is) = aivals(10,is) + one
if(.not.cris) then
if(luse .and. varinv(i) > zero) &
aivals(10,is) = aivals(10,is) + one
else if(cris) then
if(.not. cris_sw .and. wavenumber(i) < 2000.0) then
if(luse .and. varinv(i) > zero) &
aivals(10,is) = aivals(10,is) + one
else if(cris_sw .and. wavenumber(i) >= 2000.0) then
if(luse .and. varinv(i) > zero) &
aivals(10,is) = aivals(10,is) + one
end if
end if
varinv(i) = zero
if(id_qc(i) == igood_qc)id_qc(i)=ifail_sfcir_qc
end if
Expand All @@ -2365,7 +2426,11 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &
do i=1,nchanl
if (ts(i) > 0.2_r_kind) then
! QC3 in statsrad
if(luse .and. varinv(i) > zero) aivals(10,is) = aivals(10,is) + one
if(.not. cris_sw .and. wavenumber(i) < 2000.0) then
if(luse .and. varinv(i) > zero) aivals(10,is) = aivals(10,is) + one
else if(cris_sw .and. wavenumber(i) >= 2000.0) then
if(luse .and. varinv(i) > zero) aivals(10,is) = aivals(10,is) + one
end if
varinv(i) = zero
if(id_qc(i) == igood_qc)id_qc(i)=ifail_sfcir_qc
end if
Expand Down Expand Up @@ -2393,7 +2458,15 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &
if ( abs(dtz) > tzchks ) then
varinv(i) = zero
if ( id_qc(i) == igood_qc ) id_qc(i) = ifail_tzr_qc
if(luse)aivals(13,is) = aivals(13,is) + one
if(.not.cris) then
if(luse)aivals(13,is) = aivals(13,is) + one
else if(cris) then
if(.not. cris_sw .and. wavenumber(i) < 2000.0) then
if(luse) aivals(13,is) = aivals(13,is) + one
else if(cris_sw .and. wavenumber(i) >= 2000.0) then
if(luse) aivals(13,is) = aivals(13,is) + one
end if
endif
endif
endif
enddo
Expand All @@ -2402,7 +2475,9 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, &

cenlatx=abs(cenlat)*r0_04
if (cenlatx < one) then
if(luse)aivals(6,is) = aivals(6,is) + one
if(.not. cris_sw) then
if(luse)aivals(6,is) = aivals(6,is) + one
endif
efact = half*(cenlatx+one)
do i=1,nchanl
if(varinv(i) > tiny_r_kind) errf(i)=efact*errf(i)
Expand Down
Loading