Merge branch 'NOAA-EMC:master' into feature/rm_sfcanl

scripts/exgdas_efsoi_update.sh

@@ -361,6 +361,10 @@ cat > enkf.nml << EOFnml
    sattypes_rad(68)= 'ahi_himawari8', dsis(68)= 'ahi_himawari8',
    sattypes_rad(69)= 'abi_g16',       dsis(69)= 'abi_g16',
    sattypes_rad(70)= 'abi_g17',       dsis(70)= 'abi_g17',
+   sattypes_rad(71)= 'iasi_metop-c',  dsis(71)= 'iasi_metop-c',
+   sattypes_rad(72)= 'viirs-m_npp',   dsis(72)= 'viirs-m_npp',
+   sattypes_rad(73)= 'viirs-m_j1',    dsis(73)= 'viirs-m_j1',
+
    $SATOBS_ENKF
 /
 &ozobs_enkf

scripts/exgdas_enkf_update.sh

@@ -350,6 +350,8 @@ cat > enkf.nml << EOFnml
    sattypes_rad(69)= 'abi_g16',       dsis(69)= 'abi_g16',
    sattypes_rad(70)= 'abi_g17',       dsis(70)= 'abi_g17',
    sattypes_rad(71)= 'iasi_metop-c',  dsis(71)= 'iasi_metop-c',
+   sattypes_rad(72)= 'viirs-m_npp',   dsis(72)= 'viirs-m_npp',
+   sattypes_rad(73)= 'viirs-m_j1',    dsis(73)= 'viirs-m_j1',
    $SATOBS_ENKF
 /
 &ozobs_enkf

scripts/exglobal_atmos_analysis.sh

@@ -127,6 +127,7 @@ GMI1CRBF=${GMI1CRBF:-${COMIN_OBS}/${OPREFIX}gmi1cr.tm00.bufr_d${OSUFFIX}}
 SAPHIRBF=${SAPHIRBF:-${COMIN_OBS}/${OPREFIX}saphir.tm00.bufr_d${OSUFFIX}}
 SEVIRIBF=${SEVIRIBF:-${COMIN_OBS}/${OPREFIX}sevcsr.tm00.bufr_d${OSUFFIX}}
 AHIBF=${AHIBF:-${COMIN_OBS}/${OPREFIX}ahicsr.tm00.bufr_d${OSUFFIX}}
+SSTVIIRS=${SSTVIIRS:-${COMIN_OBS}/${OPREFIX}sstvcw.tm00.bufr_d${OSUFFIX}}
 ABIBF=${ABIBF:-${COMIN_OBS}/${OPREFIX}gsrcsr.tm00.bufr_d${OSUFFIX}}
 CRISBF=${CRISBF:-${COMIN_OBS}/${OPREFIX}cris.tm00.bufr_d${OSUFFIX}}
 ESCRIS=${ESCRIS:-${COMIN_OBS}/${OPREFIX}escris.tm00.bufr_d${OSUFFIX}}
@@ -510,6 +511,7 @@ $NLN $B1AVHPM          avhpmbufr
 $NLN $AHIBF            ahibufr
 $NLN $ABIBF            abibufr
 $NLN $HDOB             hdobbufr
+$NLN $SSTVIIRS         sstviirs
 
 [[ $DONST = "YES" ]] && $NLN $NSSTBF nsstbufr
 
@@ -778,7 +780,7 @@ cat > gsiparm.anl << EOF
   $OBSQC
 /
 &OBS_INPUT
-  dmesh(1)=145.0,dmesh(2)=150.0,dmesh(3)=100.0,time_window_max=3.0,
+  dmesh(1)=145.0,dmesh(2)=150.0,dmesh(3)=100.0,dmesh(4)=70.0,time_window_max=3.0,
   $OBSINPUT
 /
 OBS_INPUT::
@@ -869,10 +871,10 @@ OBS_INPUT::
    gsnd1bufr      sndrd4      g15         sndrD4_g15          0.0     1     0
    oscatbufr      uv          null        uv                  0.0     0     0
    mlsbufr        mls30       aura        mls30_aura          0.0     0     0
-   avhambufr      avhrr       metop-a     avhrr3_metop-a      0.0     1     0
-   avhpmbufr      avhrr       n18         avhrr3_n18          0.0     1     0
-   avhambufr      avhrr       metop-b     avhrr3_metop-b      0.0     1     0
-   avhpmbufr      avhrr       n19         avhrr3_n19          0.0     1     0
+   avhambufr      avhrr       metop-a     avhrr3_metop-a      0.0     4     0
+   avhpmbufr      avhrr       n18         avhrr3_n18          0.0     4     0
+   avhambufr      avhrr       metop-b     avhrr3_metop-b      0.0     4     0
+   avhpmbufr      avhrr       n19         avhrr3_n19          0.0     4     0
    amsr2bufr      amsr2       gcom-w1     amsr2_gcom-w1       0.0     3     0
    gmibufr        gmi         gpm         gmi_gpm             0.0     3     0
    saphirbufr     saphir      meghat      saphir_meghat       0.0     3     0
@@ -888,6 +890,8 @@ OBS_INPUT::
    amsuabufr      amsua       metop-c     amsua_metop-c       0.0     1     1
    mhsbufr        mhs         metop-c     mhs_metop-c         0.0     1     1
    iasibufr       iasi        metop-c     iasi_metop-c        0.0     1     1
+   sstviirs       viirs-m     npp         viirs-m_npp         0.0     4     0
+   sstviirs       viirs-m     j1          viirs-m_j1          0.0     4     0
 ::
 &SUPEROB_RADAR
   $SUPERRAD

src/gsi/gsi_obOperTypeManager.F90

@@ -345,6 +345,7 @@ function dtype2index_(dtype) result(index_)
         !
     case("avhrr_navy"); index_= iobOper_rad
     case("avhrr"  ); index_= iobOper_rad
+    case("viirs-m"  ); index_= iobOper_rad
 
   case("tcp"    ,"[tcpoper]"    ); index_= iobOper_tcp
 

src/gsi/gsi_rfv3io_mod.f90

@@ -1848,6 +1848,7 @@ subroutine wrfv3_netcdf(fv3filenamegin)
     real(r_kind),pointer,dimension(:,:,:):: ges_w   =>NULL()
     real(r_kind),pointer,dimension(:,:,:):: ges_delzinc   =>NULL()
     real(r_kind),pointer,dimension(:,:,:):: ges_delp  =>NULL()
+    real(r_kind),dimension(:,:  ),allocatable:: ges_ps_write
 
 
     real(r_kind), dimension(lat2,lon2,nsig) :: io_arr_qr, io_arr_qs
@@ -1918,8 +1919,10 @@ subroutine wrfv3_netcdf(fv3filenamegin)
       deallocate(g_prsi,ges_ps_inc)
 
     else
-      ges_ps=ges_ps*1000.0_r_kind
-      call gsi_bundleputvar (gsibundle_fv3lam_dynvar_nouv,'ps',ges_ps,istatus)
+      allocate(ges_ps_write(lat2,lon2))
+      ges_ps_write=ges_ps*1000.0_r_kind
+      call gsi_bundleputvar (gsibundle_fv3lam_dynvar_nouv,'ps',ges_ps_write,istatus)
+      deallocate(ges_ps_write)
     endif
 !   write out
     if (ier/=0) call die('get ges','cannot get pointers for fv3 met-fields, ier =',ier)

src/gsi/ncepnems_io.f90

@@ -132,6 +132,7 @@ module ncepnems_io
   public read_nemsatm
   public read_nemssfc
   public read_nemssfc_anl
+  public berror_read_hsst
   public write_fv3atm_nems  
   public write_nemsatm
   public write_nemssfc
@@ -168,6 +169,10 @@ module ncepnems_io
      module procedure read_nemssfc_anl_
   end interface
 
+  interface berror_read_hsst
+     module procedure berror_read_hsst_
+  end interface
+
   interface read_nemsnst
      module procedure read_nemsnst_
   end interface
@@ -1788,6 +1793,152 @@ subroutine read_nemssfc_anl_(iope,isli_anl)
 
   end subroutine read_nemssfc_anl_
 
+  subroutine read_hsst_(hsst,cwoption)
+!
+! abstract : read sst bg error correlation length from berror fix file
+!
+!     03/13/2020 - Xu Li (based on read_wgt)
+!
+   use kinds,     only : i_kind,r_single,r_kind
+   use gridmod,   only : nlat,nlon,nsig
+   use mpeu_util, only : check_iostat
+
+   implicit none
+
+   integer(i_kind)                    , intent(in   ) :: cwoption
+   real(r_single),dimension(nlat,nlon), intent(out  ) :: hsst
+
+   real(r_single),dimension(nlat,nlon) :: corsst
+   character(len=*),parameter :: myname_=myname//'::read_hsst'
+
+   real(r_single),dimension(nlat,nsig,nsig):: agvin
+   real(r_single),dimension(nlat,nsig) :: wgvin,bvin
+
+   integer(i_kind) :: inerr,istat,ier
+   integer(i_kind) :: nsigstat,nlatstat
+   integer(i_kind) :: isig
+   character(len=5) :: var
+
+   real(r_single), allocatable, dimension(:,:) :: hwllin
+   real(r_single), allocatable, dimension(:,:) :: corzin
+   real(r_single), allocatable, dimension(:,:) :: corq2
+   real(r_single), allocatable, dimension(:,:) :: vscalesin
+
+   character(len=256) :: berror_stats = "berror_stats"      ! filename
+
+   ! Open background error statistics file
+   inerr = 22
+   open(inerr,file=berror_stats,form='unformatted',status='old',iostat=ier)
+   call check_iostat(ier,myname_,'open('//trim(berror_stats)//')')
+
+   ! Read header.  Ensure that vertical resolution is consistent
+   ! with that specified via the user namelist
+
+   rewind inerr
+   read(inerr,iostat=ier)nsigstat,nlatstat
+   call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (nsigstat,nlatstat)')
+
+   if ( nsigstat/=nsig .or. nlatstat/=nlat ) then
+      write(6,*)'PREBAL: **ERROR** resolution of berror_stats incompatiable with GSI'
+
+      write(6,*)'PREBAL:  berror nsigstat,nlatstat=', nsigstat,nlatstat, &
+         ' -vs- GSI nsig,nlat=',nsig,nlat
+      call stop2(101)
+   endif
+
+   write(6,*) myname_,'(read_hsst):  read error amplitudes ', &
+     '"',trim(berror_stats),'".  ', &
+     'nsigstat,nlatstat =',nsigstat,nlatstat
+
+   read(inerr,iostat=ier) agvin,bvin,wgvin
+   call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (agvin,bvin,wgvin)')
+
+   write(*,*) 'in read_hsst_, cwoption = ',cwoption
+
+   readloop: do
+      read(inerr,iostat=istat) var, isig
+      if ( istat/=0 ) exit
+      write(*,*) 'var, isig, istat : ',var, isig, istat
+      allocate(corzin(nlat,isig))
+      if ( var=='q' .or. var=='cw' ) allocate(corq2(nlat,isig))
+      allocate(hwllin(nlat,isig))
+      if ( isig>1 ) allocate(vscalesin(nlat,isig))
+      if ( var/='sst' ) then
+         if ( var=='q' .or. var=='Q' .or. (var=='cw' .and. cwoption==2) ) then
+            read(inerr,iostat=ier) corzin,corq2
+            call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (corzin,corq2)')
+         else
+            read(inerr,iostat=ier) corzin
+            call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (corzin)')
+         endif
+         read(inerr,iostat=ier) hwllin
+         call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (hwllin)')
+         if ( isig>1 ) then
+            read(inerr,iostat=ier) vscalesin
+            call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (vscalein)')
+         endif
+      else
+         read(inerr,iostat=ier) corsst
+         call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (corsst)')
+         read(inerr,iostat=ier) hsst
+         call check_iostat(ier,myname_,'read('//trim(berror_stats)//') for (hsst)')
+      endif
+
+
+      deallocate(corzin,hwllin)
+      if (isig>1) deallocate(vscalesin)
+      if ( var=='q' .or. var=='cw' ) deallocate(corq2)
+
+   enddo readloop
+   close(inerr)
+
+   return
+  end subroutine read_hsst_
+
+  subroutine berror_read_hsst_(iope,hsst)
+!$$$  subprogram documentation block
+!                .      .    .                                       .
+! subprogram:    berror_read_hsst_     read hsst from berror_stats file
+!   prgmmr: xuli          org: np23                date: 2020-03-13
+!
+! program history log:
+!
+!   input argument list:
+!     iope        - mpi task handling i/o
+!
+!   output argument list:
+!     hsst      - sst bg error correlation length (km)
+!
+! attributes:
+!   language: f90
+!   machine:  ibm RS/6000 SP
+!
+!$$$
+    use kinds, only: r_kind,i_kind,r_single
+    use gridmod, only: nlat,nlon
+    use mpimod, only: mpi_itype,mpi_rtype4,mpi_comm_world,mype
+    use jfunc, only: cwoption
+    implicit none
+
+!   Declare passed variables
+    integer(i_kind),                       intent(in)  :: iope
+    real(r_single),  dimension(nlat,nlon), intent(out) :: hsst
+
+!   Declare local variables
+    integer(i_kind):: iret,npts
+
+!-----------------------------------------------------------------------------
+!   Read surface history file on processor iope
+    if(mype == iope)then
+       call read_hsst_(hsst,cwoption)
+       write(*,*) 'read hsst from berror_stats'
+    endif
+
+!   Load onto all processors
+    npts=nlat*nlon
+    call mpi_bcast(hsst,npts,mpi_rtype4,iope,mpi_comm_world,iret)
+  end subroutine berror_read_hsst_
+
   subroutine read_nst_ (tref,dt_cool,z_c,dt_warm,z_w,c_0,c_d,w_0,w_d)
 
 !$$$  subprogram documentation block

src/gsi/qcmod.f90

@@ -580,7 +580,7 @@ subroutine setup_tzr_qc(obstype)
       tzchk = 0.50_r_kind
     elseif ( obstype == 'amsua' .or. obstype == 'ssmis' .or. obstype == 'ssmi' ) then
       tzchk = 0.12_r_kind
-    elseif (  obstype == 'avhrr' .or. obstype == 'avhrr_navy' ) then 
+    elseif (  obstype == 'avhrr' .or. obstype == 'avhrr_navy' .or. obstype == 'viirs') then 
       tzchk = 0.85_r_kind
     elseif (  obstype == 'hirs2' .or. obstype == 'hirs3' .or. obstype == 'hirs4' .or. & 
               obstype == 'sndr' .or. obstype == 'sndrd1' .or. obstype == 'sndrd2'.or. &

src/gsi/radiance_mod.f90

@@ -322,6 +322,7 @@ subroutine radiance_obstype_init
 !   2018-04-04  zhu -- move rad_type_info(k)%cclr and rad_type_info(k)%ccld to this subroutine
 !   2018-04-06  derber -- change rad_type_info(k)%cclr default value from zero to a large number
 !   2019-03-21  Wei/Martin - fix capabilities for AOD assimilation
+!   2021-03-05  X.Li - add viirs (viirs-m with sstviirs data set)
 !
 !   input argument list:
 !
@@ -401,6 +402,7 @@ subroutine radiance_obstype_init
        if (index(dtype(i),'sndr') /= 0)   rtype(i)='sndr'
        if (index(dtype(i),'hirs') /= 0)   rtype(i)='hirs'
        if (index(dtype(i),'avhrr') /= 0)  rtype(i)='avhrr'
+       if (index(dtype(i),'viirs-m') /= 0)  rtype(i)='viirs'
        if (index(dtype(i),'modis') /= 0)  rtype(i)='modis'
        if (index(dtype(i),'seviri') /= 0) rtype(i)='seviri'
 
@@ -411,7 +413,7 @@ subroutine radiance_obstype_init
           rtype(i) == 'avhrr'  .or. rtype(i) == 'amsre'    .or.  rtype(i) == 'ssmis'  .or. & 
           rtype(i) == 'ssmi'   .or. rtype(i) == 'atms'     .or.  rtype(i) == 'cris'   .or. & 
           rtype(i) == 'amsr2'  .or. rtype(i) == 'gmi'      .or.  rtype(i) == 'saphir' .or. &
-          rtype(i) == 'cris-fsr' .or. rtype(i) == 'abi' ) then
+          rtype(i) == 'cris-fsr' .or. rtype(i) == 'abi'    .or.  rtype(i) == 'viirs' ) then
           drtype(i)='rads'
        end if
     end do

src/gsi/radinfo.f90

@@ -1690,7 +1690,7 @@ subroutine init_predx
    logical mean_only
    logical ssmi,ssmis,amsre,amsre_low,amsre_mid,amsre_hig,tmi,gmi,amsr2,saphir
    logical ssmis_las,ssmis_uas,ssmis_env,ssmis_img
-   logical avhrr,avhrr_navy,goessndr,goes_img,ahi,seviri,abi
+   logical avhrr,avhrr_navy,goessndr,goes_img,ahi,seviri,abi,viirs
 
    character(len=20):: obstype,platid
    character(len=20):: satsens,satsens_id
@@ -1872,6 +1872,7 @@ subroutine init_predx
       ahi        = obstype == 'ahi'
       abi        = obstype == 'abi'
       avhrr      = obstype == 'avhrr'
+      viirs      = obstype == 'viirs-m'
       avhrr_navy = obstype == 'avhrr_navy'
       ssmi       = obstype == 'ssmi'
       amsre_low  = obstype == 'amsre_low'