Merge pull request NOAA-EMC#374 from ADCollard/IRBugFix

GitHub Issue NOAA-EMC#371. Minor bug fixes to IR processing
RussTreadon-NOAA · May 25, 2022 · 8c22b33 · 8c22b33
2 parents a90b632 + cf4b323
commit 8c22b33
Show file tree

Hide file tree

Showing 5 changed files with 65 additions and 23 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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,13 +218,15 @@ 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
 
 
 ! 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
@@ -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)