Merge remote-tracking branch 'remotes/upstream/release/gfsda.v16.3.0'…

… into feature/hafs_sync_202210

modulefiles/gsi_common.lua

@@ -14,7 +14,7 @@ local sfcio_ver=os.getenv("sfcio_ver") or "1.4.1"
 local nemsio_ver=os.getenv("nemsio_ver") or "2.5.2"
 local wrf_io_ver=os.getenv("wrf_io_ver") or "1.2.0"
 local ncio_ver=os.getenv("ncio_ver") or "1.1.2"
-local crtm_ver=os.getenv("crtm_ver") or "2.3.0"
+local crtm_ver=os.getenv("crtm_ver") or "2.4.0"
 local ncdiag_ver=os.getenv("ncdiag_ver") or "1.0.0"
 
 load(pathJoin("netcdf", netcdf_ver))

regression/global_4denvar_T126.sh

@@ -5,11 +5,6 @@ set -x
 
 exp=$jobname
 
-# Set path/file for gsi executable
-#basedir=/scratch1/portfolios/NCEPDEV/da/save/Daryl.Kleist
-#gsipath=$basedir/gsi/
-#gsiexec=$gsipath/trunk/src/global_gsi.x
-
 # Set the JCAP resolution which you want.
 # All resolutions use LEVS=64
 export JCAP=126

regression/global_T62.sh

@@ -5,12 +5,6 @@ set -x
 
 exp=$jobname
 
-# Set path/file for gsi executable
-#basedir=/scratch1/portfolios/NCEPDEV/da/save/Michael.Lueken
-#gsiexec=$gsiexec
-#gsiexec=$basedir/EXP-port/src/global_gsi.x
-
-
 # Set the JCAP resolution which you want.
 # All resolutions use LEVS=64
 #export JCAP=62

regression/global_T62_ozonly.sh

@@ -5,12 +5,6 @@ set -x
 
 exp=$jobname
 
-# Set path/file for gsi executable
-#basedir=/scratch1/portfolios/NCEPDEV/da/save/Michael.Lueken
-#gsiexec=$gsiexec
-#gsiexec=$basedir/EXP-port/src/global_gsi.x
-
-
 # Set the JCAP resolution which you want.
 # All resolutions use LEVS=64
 #export JCAP=62

regression/global_fv3_4denvar_C192.sh

@@ -5,11 +5,6 @@ set -x
 
 exp=$jobname
 
-# Set path/file for gsi executable
-#basedir=/scratch1/portfolios/NCEPDEV/da/save/Daryl.Kleist
-#gsipath=$basedir/gsi/
-#gsiexec=$gsipath/trunk/src/global_gsi.x
-
 # Set the JCAP resolution which you want.
 export JCAP=190
 export LEVS=127

regression/global_fv3_4denvar_T126.sh

@@ -5,11 +5,6 @@ set -x
 
 exp=$jobname
 
-# Set path/file for gsi executable
-#basedir=/scratch1/portfolios/NCEPDEV/da/save/Daryl.Kleist
-#gsipath=$basedir/gsi/
-#gsiexec=$gsipath/trunk/src/global_gsi.x
-
 # Set the JCAP resolution which you want.
 # All resolutions use LEVS=64
 export JCAP=62

regression/global_hybrid_T126.sh

@@ -5,11 +5,6 @@ set -x
 
 exp=$jobname
 
-# Set path/file for gsi executable
-#basedir=/scratch1/portfolios/NCEPDEV/da/save/Daryl.Kleist
-#gsipath=$basedir/gsi/
-#gsiexec=$gsipath/trunk/src/global_gsi.x
-
 # Set the JCAP resolution which you want.
 # All resolutions use LEVS=64
 export JCAP=126

regression/netcdf_fv3_regional.sh

@@ -7,15 +7,6 @@ set -x
 # Set experiment name
 exp=$jobname
 
-#TM=00
-#TM2=03
-#tmmark=tm${TM}
-
-
-# Set path/file for gsi executable
-#gsiexec=/meso/save/Wanshu.Wu/Code/trunk/trunk_40320/src/global_gsi_org
-#gsiexec=/da/save/Michael.Lueken/trunk/src/global_gsi.x
-
 # Set runtime and save directories
 tmpdir=$tmpdir/tmpreg_netcdf_fv3_regional/${exp}
 savdir=$savdir/outreg_netcdf_fv3_regional/${exp}

regression/nmmb_nems_4denvar.sh

@@ -7,15 +7,6 @@ set -x
 # Set experiment name
 exp=$jobname
 
-#TM=00
-#TM2=03
-#tmmark=tm${TM}
-
-
-# Set path/file for gsi executable
-#gsiexec=/meso/save/Wanshu.Wu/Code/trunk/trunk_40320/src/global_gsi_org
-#gsiexec=/da/save/Michael.Lueken/trunk/src/global_gsi.x
-
 # Set runtime and save directories
 tmpdir=$tmpdir/tmpreg_nmmb_nems_4denvar/${exp}
 savdir=$savdir/outreg_nmmb_nems_4denvar/${exp}

regression/regression_var.sh

@@ -131,18 +131,6 @@ case $machine in
   ;;
 esac
 
-if [[ "$cmaketest" = "false" ]]; then
-  export builddir=$noscrub/build
-  export gsisrc="$basedir/$updat/src"
-  export gsiexec_updat="$gsisrc/global_gsi.x"
-  export gsiexec_contrl="$basedir/$contrl/src/global_gsi.x"
-  export enkfexec_updat="$gsisrc/enkf/global_enkf.x"
-  export enkfexec_contrl="$basedir/$contrl/src/enkf/global_enkf.x"
-  export fixgsi="$basedir/$updat/fix"
-  export scripts="$basedir/$updat/regression"
-  export ush="$basedir/$updat/ush"
-fi
-
 # We are dealing with *which* endian files
 export endianness="Big_Endian"
 

src/gsi/aircraftinfo.f90

@@ -57,7 +57,7 @@ module aircraftinfo
   logical :: cleanup_tail ! logical to remove tail number no longer used
   logical :: upd_aircraft ! indicator if update bias at 06Z & 18Z
 
-  integer(i_kind), parameter :: max_tail=10000  ! max tail numbers
+  integer(i_kind), parameter :: max_tail=100000  ! max tail numbers
   integer(i_kind) npredt          ! predictor number
   integer(i_kind) ntail           ! total tail number
   integer(i_kind) ntail_update    ! new total tail number

src/gsi/cplr_gfs_ensmod.f90

@@ -385,7 +385,7 @@ subroutine move2bundle_(grd3d,en_loc3,atm_bundle,m_cvars2d,m_cvars3d,iret)
     character(len=70) :: filename
 
     integer(i_kind) :: ierr
-    integer(i_kind) :: km,m
+    integer(i_kind) :: km1,m
     integer(i_kind) :: icw,iql,iqi,iqr,iqs,iqg  
     real(r_kind),pointer,dimension(:,:) :: ps
     !real(r_kind),pointer,dimension(:,:) :: sst
@@ -440,31 +440,31 @@ subroutine move2bundle_(grd3d,en_loc3,atm_bundle,m_cvars2d,m_cvars3d,iret)
 !      if(trim(cvars2d(m))=='sst') sst=en_loc3(:,:,m_cvars2d(m)) !no sst for now
     enddo
 
-    km = en_perts(1,1)%grid%km
+    km1 = en_perts(1,1)%grid%km - 1
 !$omp parallel do  schedule(dynamic,1) private(m) 
     do m=1,nc3d
        if(trim(cvars3d(m))=='sf')then
-          u    = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km)
+          u    = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1)
        else if(trim(cvars3d(m))=='vp') then
-          v    = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km)
+          v    = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1)
        else if(trim(cvars3d(m))=='t')  then
-          tv   = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km)
+          tv   = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1)
        else if(trim(cvars3d(m))=='q')  then
-          q    = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km)
+          q    = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1)
        else if(trim(cvars3d(m))=='oz') then
-          oz   = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km)
+          oz   = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1)
        else if(trim(cvars3d(m))=='cw') then
-          cwmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km)
+          cwmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1)
        else if(trim(cvars3d(m))=='ql') then
-          qlmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km)
+          qlmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1)
        else if(trim(cvars3d(m))=='qi') then
-          qimr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km)
+          qimr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1)
        else if(trim(cvars3d(m))=='qr') then
-          qrmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km)
+          qrmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1)
        else if(trim(cvars3d(m))=='qs') then
-          qsmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km)
+          qsmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1)
        else if(trim(cvars3d(m))=='qg') then
-          qgmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km)
+          qgmr = en_loc3(:,:,m_cvars3d(m):m_cvars3d(m)+km1)
        end if
     enddo
 
@@ -930,7 +930,7 @@ subroutine parallel_read_gfsnc_state_(en_full,m_cvars2d,m_cvars3d,nlon,nlat,nsig
 
    atmges = open_dataset(filename,errcode=ierror)
    if (ierror /=0) then
-      write(6,*)' PARALLEL_READ_GFSNC_STATE:  ***ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS'
+      write(6,*)' PARALLEL_READ_GFSNC_STATE:  ***FATAL ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS'
       call stop2(999)
    endif
    ! get dimension sizes

src/gsi/general_read_gfsatm.f90

@@ -2825,6 +2825,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
    call gsi_bundlegetpointer(gfs_bundle,'q' ,g_q   ,ier);istatus=istatus+ier
    call gsi_bundlegetpointer(gfs_bundle,'oz',g_oz  ,ier);istatus=istatus+ier
 !  call gsi_bundlegetpointer(gfs_bundle,'cw',g_cwmr,ier);istatus=istatus+ier
+   istatus1=0   
    call gsi_bundlegetpointer(gfs_bundle,'ql',g_ql  ,ier);istatus1=istatus1+ier
    call gsi_bundlegetpointer(gfs_bundle,'qi',g_qi  ,ier);istatus1=istatus1+ier
    call gsi_bundlegetpointer(gfs_bundle,'qr',g_qr  ,ier);istatus1=istatus1+ier

src/gsi/gesinfo.F90

@@ -148,7 +148,7 @@ subroutine gesinfo
      write(filename,'("sigf",i2.2)')nhr_assimilation
      inquire(file=filename,exist=fexist)
      if(.not.fexist) then
-        write(6,*)' GESINFO:  ***ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS'
+        write(6,*)' GESINFO:  ***FATAL ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS'
         call stop2(99)
         stop
      end if
@@ -339,12 +339,12 @@ subroutine gesinfo
         ! open the netCDF file
         atmges = open_dataset(filename,errcode=iret)
         if (iret /=0) then
-           write(6,*)'GESINFO:  ***ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS'
+           write(6,*)'GESINFO:  ***FATAL ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS'
            call stop2(99)
         endif
         sfcges = open_dataset(sfilename,errcode=iret)
         if (iret /=0) then
-           write(6,*)'GESINFO:  ***ERROR*** ',trim(sfilename),' NOT AVAILABLE: PROGRAM STOPS'
+           write(6,*)'GESINFO:  ***FATAL ERROR*** ',trim(sfilename),' NOT AVAILABLE: PROGRAM STOPS'
            call stop2(99)
         endif
         ! get dimension sizes
@@ -451,7 +451,7 @@ subroutine gesinfo
 !    Check for consistency with namelist settings
      if (gfshead%jcap/=jcap_b.and..not.regional .or. gfshead%levs/=nsig) then
         if (gfshead%levs/=nsig) then
-           write(6,*)'GESINFO:  ***ERROR*** guess levels inconsistent with namelist'
+           write(6,*)'GESINFO:  ***FATAL ERROR*** guess levels inconsistent with namelist'
            write(6,*)'      guess nsig=',gfshead%levs
            write(6,*)'   namelist nsig=',nsig
            fatal = .true.
@@ -466,7 +466,7 @@ subroutine gesinfo
               fatal = .false.
            else
               if ( mype == mype_out ) &
-                 write(6,*)'GESINFO:  ***ERROR*** guess jcap inconsistent with namelist'
+                 write(6,*)'GESINFO:  ***FATAL ERROR*** guess jcap inconsistent with namelist'
               fatal = .true.
            endif
            if ( mype == mype_out ) &

src/gsi/netcdfgfs_io.f90

@@ -1300,7 +1300,7 @@ subroutine read_sfc_anl_(isli_anl)
     ! open the netCDF file
     sfcges = open_dataset(filename,errcode=iret)
     if (iret/=0) then
-       write(6,*) trim(my_name),':  ***ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS'
+       write(6,*) trim(my_name),':  ***FATAL ERROR*** ',trim(filename),' NOT AVAILABLE: PROGRAM STOPS'
        call stop2(999)
     endif
 

src/gsi/radinfo.f90

@@ -615,7 +615,7 @@ subroutine radinfo_read
 ! !USES:
 
     use obsmod, only: iout_rad
-    use constants, only: zero,one,zero_quad
+    use constants, only: zero,one,zero_quad, r10
     use mpimod, only: mype
     use mpeu_util, only: perr,die
     implicit none
@@ -855,7 +855,8 @@ subroutine radinfo_read
                          varA(i,j)=varx(i)
                       end do
                       ostats(j)=ostatsx
-                      if ((any(varx/=zero) .and. iuse_rad(j)>-2) .or. iuse_rad(j)==4) & 
+                      if ((all(varx==zero) .and. iuse_rad(j)>-2) .or. iuse_rad(j)==4) cycle read3 
+                      if ((any(varx/=r10) .and. iuse_rad(j)>-2) .or. iuse_rad(j)==4) &
                          inew_rad(j)=.false.
                       cycle read3
                    end if
@@ -1867,6 +1868,7 @@ subroutine init_predx
          end do 
       end do loop_a
 
+      write(6,*) 'INIT_PREDX: inst_sat  new_chan = ', trim(fdiag_rad), new_chan
       if (.not. update .and. new_chan==0) then 
          call close_radiag(fdiag_rad,lndiag)
          cycle loopf

src/gsi/read_diag.f90

@@ -1165,7 +1165,9 @@ subroutine read_radiag_data_nc(diag_status,header_fix,data_fix,data_chan,data_ex
 
   data_fix        = diag_status%all_data_fix(diag_status%cur_ob_idx)
   data_chan(:)    = diag_status%all_data_chan(diag_status%cur_ob_idx,:)
-  data_extra(:,:) = diag_status%all_data_extra(diag_status%cur_ob_idx,:,:) 
+  if (header_fix%iextra > 0) then
+     data_extra(:,:) = diag_status%all_data_extra(diag_status%cur_ob_idx,:,:)
+  endif
 
   diag_status%cur_ob_idx = diag_status%cur_ob_idx + 1
 

src/gsi/read_files.f90

@@ -620,7 +620,7 @@ subroutine read_files(mype)
   endif
   if (l4densvar .and. nfldsig/=ntlevs_ens) then
      if (mype==0) then
-        write(6,*)'READ_FILES: ***ERROR*** insufficient atm fcst for 4densvar:  PROGRAM STOPS'
+        write(6,*)'READ_FILES: ***FATAL ERROR*** insufficient atm fcst for 4densvar:  PROGRAM STOPS'
         do i=1,ntlevs_ens
            ihr=nhr_obsbin*(i-1)+nhr_half
            present=.false.
@@ -629,7 +629,7 @@ subroutine read_files(mype)
            end do
            if (.not.present) then
               write(filename,'(''sigf'',i2.2)')ihr
-              write(6,*)'READ_FILES: ***ERROR*** file ',trim(filename),' missing:  PROGRAM STOPS'
+              write(6,*)'READ_FILES: ***FATAL ERROR*** file ',trim(filename),' missing:  PROGRAM STOPS'
            endif
         end do
      endif
@@ -652,7 +652,7 @@ subroutine read_files(mype)
   endif
   if (l4densvar .and. nfldsfc/=ntlevs_ens) then
      if (mype==0) then
-        write(6,*)'READ_FILES: ***ERROR*** insufficient sfc fcst for 4densvar:  PROGRAM STOPS'
+        write(6,*)'READ_FILES: ***FATAL ERROR*** insufficient sfc fcst for 4densvar:  PROGRAM STOPS'
         do i=1,ntlevs_ens
            ihr=nhr_obsbin*(i-1)+nhr_half
            present=.false.
@@ -661,7 +661,7 @@ subroutine read_files(mype)
            end do
            if (.not.present) then
               write(filename,'(''sfcf'',i2.2)')ihr
-              write(6,*)'READ_FILES: ***ERROR*** file ',trim(filename),' missing:  PROGRAM STOPS'
+              write(6,*)'READ_FILES: ***FATAL ERROR*** file ',trim(filename),' missing:  PROGRAM STOPS'
            endif
         end do
      endif

src/gsi/read_prepbufr.f90

@@ -475,7 +475,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
   if(tob)then
      nreal=25
   else if(uvob) then 
-     nreal=27
+     nreal=26
   else if(spdob) then
      nreal=24
   else if(psob) then
@@ -2239,8 +2239,8 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
                  cdata_all(25,iout)=var_jb(5,k)            ! non linear qc parameter
                  cdata_all(26,iout)=one                    ! hilbert curve weight, modified later 
                  if(perturb_obs)then
-                    cdata_all(28,iout)=ran01dom()*perturb_fact ! u perturbation
-                    cdata_all(29,iout)=ran01dom()*perturb_fact ! v perturbation
+                    cdata_all(27,iout)=ran01dom()*perturb_fact ! u perturbation
+                    cdata_all(28,iout)=ran01dom()*perturb_fact ! v perturbation
                  endif
 
               else if(spdob) then 
@@ -3057,7 +3057,11 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
 
         do k=1,ndata
            ikx=nint(cdata_out(10,k))
-           itype=ictype(ikx)
+           if (ikx>0) then
+              itype=ictype(ikx)
+           else
+              itype=0
+           endif
            if( itype ==230 .or. itype ==231 .or. itype ==233) then
               prest=r10*exp(cdata_out(4,k))
               if (prest <100.0_r_kind) cycle