Assimilate GMI in GSI (#689)

regression/global_4denvar.sh

@@ -55,14 +55,15 @@ cycg=`echo $gdate | cut -c9-10`
 dumpobs=gdas
 prefix_obs=${dumpobs}.t${cyca}z
 prefix_ges=gdas.t${cycg}z
-prefix_ens=gdas.t${cycg}z
+prefix_ens=enkfgdas.t${cycg}z
 suffix=tm00.bufr_d
 
 dumpges=gdas
 COMROOTgfs=$casesdir/gfs/prod
-datobs=$COMROOTgfs/$dumpobs.$PDYa/${cyca}/atmos
-datges=$COMROOTgfs/$dumpges.$PDYg/${cycg}/atmos
-datens=$COMROOTgfs/enkfgdas.$PDYg/${cycg}/atmos
+datobs=$COMROOTgfs/$dumpobs.$PDYa/${cyca}/obs
+dathis=$COMROOTgfs/$dumpges.$PDYg/${cycg}/model_data/atmos/history
+datanl=$COMROOTgfs/gdas.$PDYg/${cycg}/analysis/atmos
+datens=$COMROOTgfs/enkfgdas.$PDYg/${cycg}
 
 
 # Set up $tmpdir
@@ -265,28 +266,28 @@ $nln $datobs/${prefix_obs}.esatms.${suffix}        ./atmsbufrears
 ## $nln $datobs/${prefix_obs}.amsr2.tm00.bufr_d    ./amsr2bufr
 
 # Copy bias correction, atmospheric and surface files
-$nln $datges/${prefix_ges}.abias                      ./satbias_in
-$nln $datges/${prefix_ges}.abias_pc                   ./satbias_pc
-$nln $datges/${prefix_ges}.abias_air                  ./aircftbias_in
-$nln $datges/${prefix_ges}.radstat                    ./radstat.gdas
-
-$nln $datges/${prefix_ges}.sfcf003.nc         ./sfcf03
-$nln $datges/${prefix_ges}.sfcf004.nc         ./sfcf04
-$nln $datges/${prefix_ges}.sfcf005.nc         ./sfcf05
-$nln $datges/${prefix_ges}.sfcf006.nc         ./sfcf06
-$nln $datges/${prefix_ges}.sfcf007.nc         ./sfcf07
-$nln $datges/${prefix_ges}.sfcf008.nc         ./sfcf08
-$nln $datges/${prefix_ges}.sfcf009.nc         ./sfcf09
-
-$nln $datges/${prefix_ges}.atmf003.nc         ./sigf03
-$nln $datges/${prefix_ges}.atmf004.nc         ./sigf04
-$nln $datges/${prefix_ges}.atmf005.nc         ./sigf05
-$nln $datges/${prefix_ges}.atmf006.nc         ./sigf06
-$nln $datges/${prefix_ges}.atmf007.nc         ./sigf07
-$nln $datges/${prefix_ges}.atmf008.nc         ./sigf08
-$nln $datges/${prefix_ges}.atmf009.nc         ./sigf09
-
-$nln $datens/${prefix_ens}.sfcf006.ensmean.nc         ./sfcf06_anlgrid
+$nln $datanl/${prefix_ges}.abias                      ./satbias_in
+$nln $datanl/${prefix_ges}.abias_pc                   ./satbias_pc
+$nln $datanl/${prefix_ges}.abias_air                  ./aircftbias_in
+$nln $datanl/${prefix_ges}.radstat                    ./radstat.gdas
+
+$nln $dathis/${prefix_ges}.sfcf003.nc         ./sfcf03
+$nln $dathis/${prefix_ges}.sfcf004.nc         ./sfcf04
+$nln $dathis/${prefix_ges}.sfcf005.nc         ./sfcf05
+$nln $dathis/${prefix_ges}.sfcf006.nc         ./sfcf06
+$nln $dathis/${prefix_ges}.sfcf007.nc         ./sfcf07
+$nln $dathis/${prefix_ges}.sfcf008.nc         ./sfcf08
+$nln $dathis/${prefix_ges}.sfcf009.nc         ./sfcf09
+
+$nln $dathis/${prefix_ges}.atmf003.nc         ./sigf03
+$nln $dathis/${prefix_ges}.atmf004.nc         ./sigf04
+$nln $dathis/${prefix_ges}.atmf005.nc         ./sigf05
+$nln $dathis/${prefix_ges}.atmf006.nc         ./sigf06
+$nln $dathis/${prefix_ges}.atmf007.nc         ./sigf07
+$nln $dathis/${prefix_ges}.atmf008.nc         ./sigf08
+$nln $dathis/${prefix_ges}.atmf009.nc         ./sigf09
+
+$nln $datens/ensstat/model_data/atmos/history/${prefix_ens}.sfcf006.ensmean.nc         ./sfcf06_anlgrid
 
 export ENS_PATH='./ensemble_data/'
 mkdir -p ${ENS_PATH}
@@ -296,7 +297,7 @@ for fh in $flist; do
     imem=1
     while [[ $imem -le $NMEM_ENKF ]]; do
 	member="mem"`printf %03i $imem`
-	$nln $datens/$member/$sigens ${ENS_PATH}sigf${fh}_ens_${member}
+	$nln $datens/$member/model_data/atmos/history/$sigens ${ENS_PATH}sigf${fh}_ens_${member}
 	(( imem = $imem + 1 ))
     done
 done

regression/global_enkf.sh

@@ -51,17 +51,14 @@ cyca=`echo $global_adate | cut -c9-10`
 PDYg=`echo $gdate | cut -c1-8`
 cycg=`echo $gdate | cut -c9-10`
 
-dumpobs=gdas
-prefix_obs=${dumpobs}.t${cyca}z
-prefix_ges=gdas.t${cycg}z
-prefix_ens=gdas.t${cycg}z
+prefix_obs=enkfgdas.t${cyca}z
+prefix_ens=enkfgdas.t${cycg}z
 suffix=tm00.bufr_d
 
 dumpges=gdas
 COMROOTgfs=$casesdir/gfs/prod
-datobs=$COMROOTgfs/enkfgdas.$PDYa/${cyca}/atmos
-datges=$COMROOTgfs/$dumpges.$PDYg/${cycg}/atmos
-datens=$COMROOTgfs/enkfgdas.$PDYg/${cycg}/atmos
+datobs=$COMROOTgfs/enkfgdas.$PDYa/${cyca}/ensstat/analysis/atmos
+datens=$COMROOTgfs/enkfgdas.$PDYg/${cycg}
 
 
 # Set up $tmpdir
@@ -166,19 +163,19 @@ nfhrs=`echo $IAUFHRS_ENKF | sed 's/,/ /g'`
 for fhr in $nfhrs; do
     for imem in $(seq 1 $NMEM_ENKF); do
 	memchar="mem"$(printf %03i $imem)
-	$nln $datens/$memchar/${prefix_ens}.atmf00${fhr}.nc sfg_${global_adate}_fhr0${fhr}_${memchar}
+	$nln $datens/$memchar/model_data/atmos/history/${prefix_ens}.atmf00${fhr}.nc sfg_${global_adate}_fhr0${fhr}_${memchar}
         if [ $cnvw_option = ".true." ]; then
-            $nln $datens/$memchar/${prefix_ens}sfcf00${fhr}.nc sfgsfc_${global_adate}_fhr0${fhr}_${memchar}
+            $nln $datens/$memchar/model_data/atmos/history/${prefix_ens}sfcf00${fhr}.nc sfgsfc_${global_adate}_fhr0${fhr}_${memchar}
         fi
 	(( imem = $imem + 1 ))
     done
-    $nln $datens/${prefix_ens}.atmf00${fhr}.ensmean.nc sfg_${global_adate}_fhr0${fhr}_ensmean
+    $nln $datens/ensstat/model_data/atmos/history/${prefix_ens}.atmf00${fhr}.ensmean.nc sfg_${global_adate}_fhr0${fhr}_ensmean
     if [ $cnvw_option = ".true." ]; then
         $nln $datens/${prefix_ens}.sfcf00${fhr}.ensmean.nc sfgsfc_${global_adate}_fhr0${fhr}_ensmean
     fi
 done
 
-$nln $datobs/${prefix_obs}.abias_int ./satbias_in
+$nln $datobs/${prefix_obs}.abias_int.ensmean ./satbias_in
 
 
 cd $tmpdir

regression/regression_namelists.sh

@@ -181,7 +181,7 @@ OBS_INPUT::
    sstviirs       viirs-m     j1          viirs-m_j1          0.0     4     0
    abibufr        abi         g18         abi_g18             0.0     1     0
    ahibufr        ahi         himawari9   ahi_himawari9       0.0     1     0
-   atmsbufr       atms        n21         atms_n21            0.0     1     1
+   atmsbufr       atms        n21         atms_n21            0.0     1     0
    crisfsbufr     cris-fsr    n21         cris-fsr_n21        0.0     1     0
    sstviirs       viirs-m     j2          viirs-m_j2          0.0     4     0
    ompsnpbufr     ompsnp      n21         ompsnp_n21          0.0     0     0

regression/regression_var.sh

@@ -190,7 +190,7 @@ export savdir="$ptmp"
 export JCAP="62"
 
 # Case Study analysis dates
-export global_adate="2022110900"
+export global_adate="2024022300"
 export rtma_adate="2020022420"
 export fv3_netcdf_adate="2017030100"
 export rrfs_3denvar_glbens_adate="2021072518"

src/gsi/clw_mod.f90

@@ -2019,7 +2019,7 @@ subroutine gmi_37pol_diff(tb37v,tb37h,tsim37v,tsim37h,clw,ierrret)
 
      clw = one - (tb37v-tb37h)/(tsim37v-tsim37h)
      clw=max(zero,clw)
-     if(tb37h > tb37v)  then
+     if ((tb37h > tb37v) .or. (tb37h > 500_r_kind )) then
         ierrret = 1
         clw= r1000
      endif

src/gsi/deter_sfc_mod.f90

@@ -33,7 +33,7 @@ module deter_sfc_mod
   use satthin, only: sno_full,isli_full,sst_full,soil_moi_full, &
       soil_temp_full,soil_type_full,veg_frac_full,veg_type_full, &
       fact10_full,zs_full,sfc_rough_full,zs_full_gfs
-  use constants, only: zero,one,two,one_tenth,deg2rad,rad2deg
+  use constants, only: zero,one,two,one_tenth,deg2rad,rad2deg, rearth
   use gridmod, only: nlat,nlon,regional,tll2xy,nlat_sfc,nlon_sfc,rlats_sfc,rlons_sfc, &
       rlats,rlons,dx_gfs,txy2ll,lpl_gfs
   use guess_grids, only: nfldsfc,hrdifsfc,ntguessfc
@@ -1331,7 +1331,7 @@ subroutine deter_sfc_amsre_low(dlat_earth,dlon_earth,isflg,sfcpct)
    end subroutine deter_sfc_amsre_low
 
 
-subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg,sfcpct)
+subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
 ! subprogram:    deter_sfc_gmi           determine land surface type
@@ -1354,11 +1354,6 @@ subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg,sfcpct)
 !                2 sea ice
 !                3 snow
 !                4 mixed
-!      sfcpct(0:3)- percentage of 4 surface types
-!                 (0) - sea percentage
-!                 (1) - land percentage
-!                 (2) - sea ice percentage
-!                 (3) - snow percentage
 !
 ! attributes:
 !   language: f90
@@ -1370,15 +1365,11 @@ subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg,sfcpct)
 
    real(r_kind)               ,intent(in   ) :: dlat_earth,dlon_earth
    integer(i_kind)            ,intent(  out) :: isflg
-   real(r_kind),dimension(0:3),intent(  out) :: sfcpct
-
-   integer(i_kind) jsli,it
-   integer(i_kind):: klat1,klon1,klatp1,klonp1
-   real(r_kind):: dx,dy,dx1,dy1,w00,w10,w01,w11
-   real(r_kind) :: dlat,dlon
+   integer(i_kind) jsli,it, i, j
+   integer(i_kind):: klat1,klon1,klatp1,klonp1, ksmall, klarge, n_grid
+   real(r_kind) :: dlat,dlon, grid_dist
+   integer(i_kind):: klatn,klonn,klatpn,klonpn
    logical :: outside
-   integer(i_kind):: klat2,klon2,klatp2,klonp2
-
 !
 !  For interpolation, we usually use o points (4points for land sea decision)
 !  In case of lowfreq channel (Large FOV), add the check of x points(8 points)
@@ -1407,90 +1398,55 @@ subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg,sfcpct)
      end if
 
      klon1=int(dlon); klat1=int(dlat)
-     dx  =dlon-klon1; dy  =dlat-klat1
-     dx1 =one-dx;    dy1 =one-dy
-     w00=dx1*dy1; w10=dx1*dy; w01=dx*dy1; w11=dx*dy
 
      klat1=min(max(1,klat1),nlat_sfc); klon1=min(max(0,klon1),nlon_sfc)
      if(klon1==0) klon1=nlon_sfc
      klatp1=min(nlat_sfc,klat1+1); klonp1=klon1+1
-     if(klonp1==nlon_sfc+1) klonp1=1
-     klonp2 = klonp1+1
-     if(klonp2==nlon_sfc+1) klonp2=1
-     klon2=klon1-1
-     if(klon2==0)klon2=nlon_sfc
-     klat2=max(1,klat1-1)
-     klatp2=min(nlat_sfc,klatp1+1)
 
 !    Set surface type flag.  Begin by assuming obs over ice-free water
 
-     sfcpct = zero
-
-     jsli = isli_full(klat1 ,klon1 )
-     if(sno_full(klat1 ,klon1 ,it) > one .and. jsli == 1)jsli=3
-     sfcpct(jsli)=sfcpct(jsli)+one
-
-     jsli = isli_full(klatp1,klon1 )
-     if(sno_full(klatp1 ,klon1 ,it) > one .and. jsli == 1)jsli=3
-     sfcpct(jsli)=sfcpct(jsli)+one
-
-     jsli = isli_full(klat1 ,klonp1)
-     if(sno_full(klat1 ,klonp1 ,it) > one .and. jsli == 1)jsli=3
-     sfcpct(jsli)=sfcpct(jsli)+one
-
-     jsli = isli_full(klatp1,klonp1)
-     if(sno_full(klatp1 ,klonp1 ,it) > one .and. jsli == 1)jsli=3
-     sfcpct(jsli)=sfcpct(jsli)+one
-
-     jsli = isli_full(klatp2,klon1)
-     if(sno_full(klatp2 ,klon1 ,it) > one .and. jsli == 1)jsli=3
-     sfcpct(jsli)=sfcpct(jsli)+one
-
-     jsli = isli_full(klatp2,klonp1)
-     if(sno_full(klatp2 ,klonp1 ,it) > one .and. jsli == 1)jsli=3
-     sfcpct(jsli)=sfcpct(jsli)+one
-
-     jsli = isli_full(klatp1,klon2)
-     if(sno_full(klatp1 ,klon2 ,it) > one .and. jsli == 1)jsli=3
-     sfcpct(jsli)=sfcpct(jsli)+one
-
-     jsli = isli_full(klatp1,klonp2)
-     if(sno_full(klatp1 ,klonp2 ,it) > one .and. jsli == 1)jsli=3
-     sfcpct(jsli)=sfcpct(jsli)+one
-
-     jsli = isli_full(klat1,klon2)
-     if(sno_full(klat1 ,klon2 ,it) > one .and. jsli == 1)jsli=3
-     sfcpct(jsli)=sfcpct(jsli)+one
-
-     jsli = isli_full(klat1,klonp2)
-     if(sno_full(klat1 ,klonp2 ,it) > one .and. jsli == 1)jsli=3
-     sfcpct(jsli)=sfcpct(jsli)+one
-
-     jsli = isli_full(klat2,klon1)
-     if(sno_full(klat2 ,klon1 ,it) > one .and. jsli == 1)jsli=3
-     sfcpct(jsli)=sfcpct(jsli)+one
-
-     jsli = isli_full(klat2,klonp1)
-     if(sno_full(klat2 ,klonp1 ,it) > one .and. jsli == 1)jsli=3
-     sfcpct(jsli)=sfcpct(jsli)+one
-
-     sfcpct=sfcpct/12.0_r_kind
-
-!     sfcpct(3)=min(sfcpct(3),sfcpct(1))
-!     sfcpct(1)=max(zero,sfcpct(1)-sfcpct(3))
-
-     if(sfcpct(0) > 0.99_r_kind)then
-        isflg = 0
-     else if(sfcpct(1) > 0.99_r_kind)then
-        isflg = 1
-     else if(sfcpct(2) > 0.99_r_kind)then
-        isflg = 2
-     else if(sfcpct(3) > 0.99_r_kind)then
-        isflg = 3
-     else
-        isflg = 4
-     end if
 
+     grid_dist=rearth * (rlats_sfc(klatp1) - rlats_sfc(klat1))
+     n_grid=int(40000 / grid_dist) + 1
+     klatn = max(klat1 - n_grid, 1)
+     klonn = klon1 - n_grid
+     if (klonn < 0)  klonn = nlon_sfc - klonn
+     klatpn = min((klat1 + n_grid), nlat_sfc)
+     klonpn = klon1 + n_grid
+     if (klonpn > nlon_sfc)  klonpn = klonpn - nlon_sfc
+
+     isflg=0
+     outer: do i = klatn, klatpn
+       ! assume n_grid > 2
+       if (0 < klonpn - klonn .and. klonpn - klonn < nlon_sfc / 2) then
+         do j = klonn, klonpn
+           if (isli_full(i, j) /= 0) then
+             isflg = 1
+             exit outer
+           end if
+         end do
+       else
+         if (klonpn < klonn) then
+           ksmall = klonpn
+           klarge = klonn
+         else
+           ksmall = klonn
+           klarge = klonpn
+         end if
+         do j = 1, ksmall
+           if (isli_full(i, j) /= 0) then
+             isflg = 1
+             exit outer
+           endif
+         end do
+         do j = klarge, nlon_sfc
+           if (isli_full(i, j) /= 0) then
+             isflg = 1
+             exit outer
+           end if
+         end do
+       end if
+     end do outer
      return
 
    end subroutine deter_sfc_gmi

src/gsi/radiance_mod.f90

@@ -1326,11 +1326,7 @@ subroutine radiance_ex_biascor_gmi(radmod,clw_obs,clw_guess_retrieval,nchanl,cld
 
     do i=1,nchanl
        if (radmod%lcloud4crtm(i)<0) cycle
-       if (clw_obs <= cclr(i) .and. clw_guess_retrieval <= cclr(i) .and. abs(clw_obs-clw_guess_retrieval) < 0.001_r_kind) then
-           cld_rbc_idx(i)=one   !clear/clear
-       else
-           cld_rbc_idx(i)=zero
-       endif
+       if ((clw_obs-cclr(i))*(clw_guess_retrieval-cclr(i))<zero .and. abs(clw_obs-clw_guess_retrieval)>=0.005_r_kind) cld_rbc_idx(i)=zero
     end do
     return
 

src/gsi/read_gmi.f90

@@ -184,7 +184,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,&
   integer(i_kind) :: jc,bufsat,n    
   integer(i_kind),dimension(5):: iobsdate
   integer(i_kind):: method,iobs,num_obs
-  integer(i_kind),parameter   :: maxobs=4000000
+  integer(i_kind),parameter   :: maxobs=6000000
   !-- integer(i_kind),parameter   :: nscan=74       ! after binning ifov, 221/3 + 1
   integer(i_kind),parameter   :: nscan=221
 
@@ -414,7 +414,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,&
           call ufbrep(lnbufr,var_check1,1,nchanl,iret,'GMICHQ')
           !call ufbrep(lnbufr,gmirfi,1,nchanl,iret,'GMIRFI')
           call ufbrep(lnbufr,pixelsaza,1,ngs,iret,'SAZA')
-          call ufbrep(lnbufr,val_angls,n_angls,ngs,iret,'SAMA SZA SMA SGA')
+          call ufbrep(lnbufr,val_angls,n_angls,ngs,iret,'BEARAZ SOZA SOLAZI SSGA')
           call ufbint(lnbufr,pixelloc,2, 1,iret,'CLATH CLONH')
 
         if (any(var_check1 < 99999999999_r_double)) then   ! 100000000000 seems to be the missing value
@@ -696,7 +696,7 @@ subroutine read_gmi(mype,val_gmi,ithin,rmesh,jsatid,gstime,&
 
     call deter_sfc(dlat,dlon,dlat_earth,dlon_earth,t4dv,isflg,idomsfc,sfcpct, &
          ts,tsavg,vty,vfr,sty,stp,sm,sn,zz,ff10,sfcr)
-    call deter_sfc_gmi(dlat_earth,dlon_earth,isflg,sfcpct)
+    call deter_sfc_gmi(dlat_earth,dlon_earth,isflg)
 
 
 !   Only keep obs over ocean    - ej