From cf4b3232efc8856f62733fcf2e6642d128817999 Mon Sep 17 00:00:00 2001 From: Andrew Collard Date: Fri, 29 Apr 2022 21:34:22 +0000 Subject: [PATCH] GitHub Issue NOAA-EMC/GSI#371: Bug fixes to IR processing --- src/gsi/qcmod.f90 | 7 +++++-- src/gsi/read_airs.f90 | 14 ++++++++++--- src/gsi/read_cris.f90 | 14 ++++++++----- src/gsi/read_iasi.f90 | 48 ++++++++++++++++++++++++++++++++----------- src/gsi/setuprad.f90 | 5 ++++- 5 files changed, 65 insertions(+), 23 deletions(-) diff --git a/src/gsi/qcmod.f90 b/src/gsi/qcmod.f90 index 54ffac432a..f4afdbae9d 100644 --- a/src/gsi/qcmod.f90 +++ b/src/gsi/qcmod.f90 @@ -78,6 +78,7 @@ module qcmod ! 2019-06-10 h. liu - add Geostationary satellites CSR data QC to replace qc_abi,qc_seviri ! 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 ! ! subroutines included: ! sub init_qcvars @@ -2335,8 +2336,10 @@ subroutine qc_irsnd(nchanl,is,ndat,nsig,ich,sea,land,ice,snow,luse,goessndr, & sum=zero sum2=zero do i=1,nchanl - sum=sum+tbc(i)*ts(i)*varinv_use(i) - sum2=sum2+ts(i)*ts(i)*varinv_use(i) + if ( varinv_use(i) > tiny_r_kind .and. ts(i) > 0.0001_r_kind) then + sum=sum+tbc(i)*ts(i)*varinv_use(i) + sum2=sum2+ts(i)*ts(i)*varinv_use(i) + endif end do if (abs(sum2) < tiny_r_kind) sum2 = sign(tiny_r_kind,sum2) dts=abs(sum/sum2) diff --git a/src/gsi/read_airs.f90 b/src/gsi/read_airs.f90 index ab6f770a91..16e890abd1 100644 --- a/src/gsi/read_airs.f90 +++ b/src/gsi/read_airs.f90 @@ -561,7 +561,7 @@ subroutine read_airs(mype,val_airs,ithin,isfcalc,rmesh,jsatid,gstime,& endif crit1 = crit1 + rlndsea(isflg) - call checkob(dist1,crit1,itx,iuse) + call checkob(one,crit1,itx,iuse) if(.not. iuse)cycle read_loop ! Set common predictor parameters @@ -747,7 +747,11 @@ subroutine read_airs(mype,val_airs,ithin,isfcalc,rmesh,jsatid,gstime,& ! Compute "score" for observation. All scores>=0.0. Lowest score is "best" crit1 = crit1+pred - call checkob(dist1,crit1,itx,iuse) + if (pred == zero) then + call checkob(dist1,crit1,itx,iuse) + else + call checkob(one,crit1,itx,iuse) + endif if(.not. iuse)cycle read_loop ! check for missing channels (if key channel reject) @@ -771,7 +775,11 @@ subroutine read_airs(mype,val_airs,ithin,isfcalc,rmesh,jsatid,gstime,& if( iskip >= satinfo_nchan )cycle read_loop ! Map obs to grids - call finalcheck(dist1,crit1,itx,iuse) + if (chsst > tsavg) then + call finalcheck(dist1,crit1,itx,iuse) + else + call finalcheck(one,crit1,itx,iuse) + endif if(.not. iuse) cycle read_loop ! Replace popped AIRS channel Tb with zero diff --git a/src/gsi/read_cris.f90 b/src/gsi/read_cris.f90 index 29067b3f69..a257480c9a 100644 --- a/src/gsi/read_cris.f90 +++ b/src/gsi/read_cris.f90 @@ -163,7 +163,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& real(r_kind) :: dlon, dlat real(r_kind) :: dlon_earth,dlat_earth,dlon_earth_deg,dlat_earth_deg real(r_kind) :: rsat - real(r_kind) :: pred, crit1, dist1 + real(r_kind) :: pred, pred1, pred2, crit1, dist1 real(r_kind) :: sat_zenang, sat_look_angle, look_angle_est real(crtm_kind) :: radiance real(r_kind) :: tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr @@ -696,12 +696,16 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& ! Compute "score" for observation. All scores>=0.0. Lowest score is "best" if ( cloud_properties(1) < one ) then !Assume clear clear = .true. - else ! Assume a lapse rate to convert hgt to delta TB. - pred = cloud_properties(2) *7.0_r_kind / r1000 + else + pred1 = cloud_properties(2) *7.0_r_kind / r1000 ! Assume a lapse rate to convert hgt to delta TB. + radiance = allchan(2,sfc_channel_index) * r1000 ! Conversion from W to mW + call crtm_planck_temperature(sensorindex,sfc_channel,radiance,temperature(sfc_channel_index)) ! radiance to BT calculation + pred2 = tsavg *0.98_r_kind - temperature(sfc_channel_index) + pred = max(pred1,pred2) ! use the largest of lapse rate (pred1) or sfc channel-surface difference (pred2) endif else -! If cloud_properties is missing from BUFR, use proxy of warmest fov. +! If cloud_properties are missing from BUFR, use proxy of warmest fov. ! the surface channel is fixed and set earlier in the code (501). radiance = allchan(2,sfc_channel_index) * r1000 ! Conversion from W to mW @@ -710,7 +714,7 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,& if ( tsavg*0.98_r_kind <= temperature(sfc_channel_index)) then ! 0.98 is a crude estimate of the surface emissivity clear = .true. else - pred = (tsavg * 0.98_r_kind - temperature(sfc_channel_index)) + pred = tsavg * 0.98_r_kind - temperature(sfc_channel_index) endif else cycle read_loop diff --git a/src/gsi/read_iasi.f90 b/src/gsi/read_iasi.f90 index b26de400d8..367c224508 100644 --- a/src/gsi/read_iasi.f90 +++ b/src/gsi/read_iasi.f90 @@ -67,6 +67,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& ! 2015-10-22 Jung - added logic to allow subset changes based on the satinfo file ! 2016-04-28 jung - added logic for RARS and direct broadcast from NESDIS/UW ! 2018-05-21 j.jin - added time-thinning. Moved the checking of thin4d into satthin.F90. +! 2022-04-29 Jung/Collard - allow thinning based on clear sky if AVHRR is missing ! ! input argument list: ! mype - mpi task id @@ -187,7 +188,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& ! Other work variables - real(r_kind) :: clr_amt,piece + real(r_kind) :: piece real(r_kind) :: rsat, dlon, dlat real(r_kind) :: dlon_earth,dlat_earth,dlon_earth_deg,dlat_earth_deg real(r_kind) :: lza, lzaest,sat_height_ratio @@ -205,7 +206,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& real(r_kind) cdist,disterr,disterrmax,dlon00,dlat00 logical :: outside,iuse,assim,valid - logical :: iasi,quiet + logical :: iasi,quiet,cloud_info integer(i_kind) :: ifov, instr, iscn, ioff, sensorindex integer(i_kind) :: i, j, l, iskip, ifovn, bad_line, ksatid, kidsat, llll @@ -217,6 +218,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& integer(i_kind):: error_status, irecx,ierr integer(i_kind):: radedge_min, radedge_max integer(i_kind) :: subset_start, subset_end, satinfo_nchan, sc_chan, bufr_chan + integer(i_kind) :: sfc_channel_index integer(i_kind),allocatable, dimension(:) :: channel_number, sc_index, bufr_index integer(i_kind),allocatable, dimension(:) :: bufr_chan_test character(len=20),dimension(1):: sensorlist @@ -224,6 +226,7 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& ! Set standard parameters character(8),parameter:: fov_flag="crosstrk" + integer(i_kind),parameter:: sfc_channel=1271 integer(i_kind),parameter:: ichan=-999 ! fov-based surface code is not channel specific for iasi real(r_kind),parameter:: expansion=one ! exansion factor for fov-based surface code. ! use one for ir sensors. @@ -633,22 +636,23 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& ! Set common predictor parameters crit1 = crit1 + rlndsea(isflg) - call checkob(dist1,crit1,itx,iuse) + call checkob(one,crit1,itx,iuse) if(.not. iuse)cycle read_loop ! Clear Amount (percent clear) - call ufbrep(lnbufr,cloud_frac,1,7,iret,'FCPH') - clr_amt = cloud_frac(1) - clr_amt=max(clr_amt,zero) - clr_amt=min(clr_amt,100.0_r_kind) - ! Compute "score" for observation. All scores>=0.0. Lowest score is "best" - pred = 100.0_r_kind - clr_amt + pred = r100 + cloud_info = .false. + call ufbrep(lnbufr,cloud_frac,1,7,iret,'FCPH') + if (iret == 7 .and. cloud_frac(1) <= r100 .and. cloud_frac(1) >= zero) then + pred = r100 - cloud_frac(1) + cloud_info = .true. + endif crit1 = crit1 + pred - - call checkob(dist1,crit1,itx,iuse) + call checkob(one,crit1,itx,iuse) if(.not. iuse)cycle read_loop + call ufbseq(lnbufr,cscale,3,10,iret,'IASIL1CB') if(iret /= 10) then write(6,*) 'READ_IASI read scale error ',iret @@ -680,18 +684,25 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& ! Coordinate bufr channels with satinfo file channels ! If this is the first time or a change in the bufr channels is detected, sync with satinfo file if (ANY(int(allchan(1,:)) /= bufr_chan_test(:))) then + sfc_channel_index = 0 bufr_index(:) = 0 bufr_chans: do l=1,bufr_nchan bufr_chan_test(l) = int(allchan(1,l)) ! Copy this bufr channel selection into array for comparison to next profile satinfo_chans: do i=1,satinfo_nchan ! Loop through sensor (iasi) channels in the satinfo file if ( channel_number(i) == int(allchan(1,l)) ) then ! Channel found in both bufr and satinfo file bufr_index(i) = l + if ( channel_number(i) == sfc_channel) sfc_channel_index = l exit satinfo_chans ! go to next bufr channel endif end do satinfo_chans end do bufr_chans endif + if (sfc_channel_index == 0) then + write(6,*)'READ_IASI: ***ERROR*** SURFACE CHANNEL USED FOR QC WAS NOT FOUND' + cycle read_loop + endif + iskip = 0 jstart=1 channel_loop: do i=1,satinfo_nchan @@ -729,8 +740,21 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,& crit1=crit1 + ten*float(iskip) +! If the surface channel exists (~960.0 cm-1) and the AVHRR cloud information is missing, use an +! estimate of the surface temperature to determine if the profile may be clear. + if (.not. cloud_info) then + pred = tsavg*0.98_r_kind - temperature(sfc_channel_index) + pred = max(pred,zero) + endif + + crit1=crit1 + pred + ! Map obs to grids - call finalcheck(dist1,crit1,itx,iuse) + if (pred == zero) then + call finalcheck(dist1,crit1,itx,iuse) + else + call finalcheck(one,crit1,itx,iuse) + endif if(.not. iuse)cycle read_loop ! diff --git a/src/gsi/setuprad.f90 b/src/gsi/setuprad.f90 index 491afb402d..92767d4391 100644 --- a/src/gsi/setuprad.f90 +++ b/src/gsi/setuprad.f90 @@ -1165,9 +1165,12 @@ subroutine setuprad(obsLL,odiagLL,lunin,mype,aivals,stats,nchanl,nreal,nobs,& end if end if - do j=1, npred-angord +! Zero out air mass terms if required + do j=2, npred-angord pred(j,i)=pred(j,i)*air_rad(mm) end do + +! Zero out angle terms if required if (adp_anglebc) then do j=npred-angord+1, npred pred(j,i)=pred(j,i)*ang_rad(mm)