GitHub Issue NOAA-EMC/GSI#237: merge release/gfsda.v16.x at a723db1 into master

scripts/exgdas_enkf_update.sh

@@ -349,6 +349,7 @@ 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',
    $SATOBS_ENKF
 /
 &ozobs_enkf
@@ -362,6 +363,9 @@ cat > enkf.nml << EOFnml
    sattypes_oz(8) = 'mls30_aura',
    sattypes_oz(9) = 'ompsnp_npp',
    sattypes_oz(10) = 'ompstc8_npp',
+   sattypes_oz(11) = 'ompstc8_n20',
+   sattypes_oz(12) = 'ompsnp_n20',
+   sattypes_oz(13) = 'ompslp_npp',
    $OZOBS_ENKF
 /
 EOFnml

scripts/exglobal_atmos_analysis.sh

@@ -165,10 +165,10 @@ SSMISBF=${SSMISBF:-${COMIN_OBS}/${OPREFIX}ssmisu.tm00.bufr_d${OSUFFIX}}
 SBUVBF=${SBUVBF:-${COMIN_OBS}/${OPREFIX}osbuv8.tm00.bufr_d${OSUFFIX}}
 OMPSNPBF=${OMPSNPBF:-${COMIN_OBS}/${OPREFIX}ompsn8.tm00.bufr_d${OSUFFIX}}
 OMPSTCBF=${OMPSTCBF:-${COMIN_OBS}/${OPREFIX}ompst8.tm00.bufr_d${OSUFFIX}}
+OMPSLPBF=${OMPSLPBF:-${COMIN_OBS}/${OPREFIX}ompslp.tm00.bufr_d${OSUFFIX}}
 GOMEBF=${GOMEBF:-${COMIN_OBS}/${OPREFIX}gome.tm00.bufr_d${OSUFFIX}}
 OMIBF=${OMIBF:-${COMIN_OBS}/${OPREFIX}omi.tm00.bufr_d${OSUFFIX}}
 MLSBF=${MLSBF:-${COMIN_OBS}/${OPREFIX}mls.tm00.bufr_d${OSUFFIX}}
-OMPSLPBF=${OMPSLPBF:-${COMIN_OBS}/${OPREFIX}ompslp.tm00.bufr_d${OSUFFIX}}
 SMIPCP=${SMIPCP:-${COMIN_OBS}/${OPREFIX}spssmi.tm00.bufr_d${OSUFFIX}}
 TMIPCP=${TMIPCP:-${COMIN_OBS}/${OPREFIX}sptrmm.tm00.bufr_d${OSUFFIX}}
 GPSROBF=${GPSROBF:-${COMIN_OBS}/${OPREFIX}gpsro.tm00.bufr_d${OSUFFIX}}
@@ -499,6 +499,7 @@ $NLN $AMUBDB           amsubbufr_db
 #$NLN $MHSDB            mhsbufr_db
 $NLN $SBUVBF           sbuvbufr
 $NLN $OMPSNPBF         ompsnpbufr
+$NLN $OMPSLPBF         ompslpbufr
 $NLN $OMPSTCBF         ompstcbufr
 $NLN $GOMEBF           gomebufr
 $NLN $OMIBF            omibufr
@@ -904,7 +905,10 @@ OBS_INPUT::
    abibufr        abi         g17         abi_g17             0.0     1     0
    rapidscatbufr  uv          null        uv                  0.0     0     0
    ompsnpbufr     ompsnp      npp         ompsnp_npp          0.0     0     0
+   ompslpbufr     ompslp      npp         ompslp_npp          0.0     0     0
    ompstcbufr     ompstc8     npp         ompstc8_npp         0.0     2     0
+   ompsnpbufr     ompsnp      n20         ompsnp_n20          0.0     0     0
+   ompstcbufr     ompstc8     n20         ompstc8_n20         0.0     2     0
    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

scripts/exglobal_diag.sh

@@ -115,7 +115,7 @@ if [ $GENDIAG = "YES" ] ; then
 
    diagtype[0]="conv conv_gps conv_ps conv_pw conv_q conv_sst conv_t conv_tcp conv_uv conv_spd"
    diagtype[1]="pcp_ssmi_dmsp pcp_tmi_trmm"
-   diagtype[2]="sbuv2_n16 sbuv2_n17 sbuv2_n18 sbuv2_n19 gome_metop-a gome_metop-b omi_aura mls30_aura ompsnp_npp ompstc8_npp gome_metop-c"
+   diagtype[2]="sbuv2_n16 sbuv2_n17 sbuv2_n18 sbuv2_n19 gome_metop-a gome_metop-b omi_aura mls30_aura ompsnp_npp ompstc8_npp  ompstc8_n20 ompsnp_n20 ompslp_npp gome_metop-c"
    diagtype[3]="hirs2_n14 msu_n14 sndr_g08 sndr_g11 sndr_g12 sndr_g13 sndr_g08_prep sndr_g11_prep sndr_g12_prep sndr_g13_prep sndrd1_g11 sndrd2_g11 sndrd3_g11 sndrd4_g11 sndrd1_g12 sndrd2_g12 sndrd3_g12 sndrd4_g12 sndrd1_g13 sndrd2_g13 sndrd3_g13 sndrd4_g13 sndrd1_g14 sndrd2_g14 sndrd3_g14 sndrd4_g14 sndrd1_g15 sndrd2_g15 sndrd3_g15 sndrd4_g15 hirs3_n15 hirs3_n16 hirs3_n17 amsua_n15 amsua_n16 amsua_n17 amsub_n15 amsub_n16 amsub_n17 hsb_aqua airs_aqua amsua_aqua imgr_g08 imgr_g11 imgr_g12 imgr_g14 imgr_g15 ssmi_f13 ssmi_f15 hirs4_n18 hirs4_metop-a amsua_n18 amsua_metop-a mhs_n18 mhs_metop-a amsre_low_aqua amsre_mid_aqua amsre_hig_aqua ssmis_f16 ssmis_f17 ssmis_f18 ssmis_f19 ssmis_f20 iasi_metop-a hirs4_n19 amsua_n19 mhs_n19 seviri_m08 seviri_m09 seviri_m10 seviri_m11 cris_npp cris-fsr_npp cris-fsr_n20 atms_npp atms_n20 hirs4_metop-b amsua_metop-b mhs_metop-b iasi_metop-b avhrr_metop-b avhrr_n18 avhrr_n19 avhrr_metop-a amsr2_gcom-w1 gmi_gpm saphir_meghat ahi_himawari8 abi_g16 abi_g17 amsua_metop-c mhs_metop-c iasi_metop-c avhrr_metop-c"
 
    diaglist[0]=listcnv

src/enkf/innovstats.f90

@@ -219,7 +219,7 @@ subroutine print_innovstats(obfit,obsprd)
    call printstats('  all  rw',sumrw_nh,biasq_nh,sumrw_spread_nh,sumrw_oberr_nh,nobsrw_nh,&
         sumrw_sh,biasrw_sh,sumrw_spread_sh,sumrw_oberr_sh,nobsrw_sh,&
         sumrw_tr,biasrw_tr,sumrw_spread_tr,sumrw_oberr_tr,nobsrw_tr)
-   call printstats(' sbuv2 oz',sumoz_nh,biasoz_nh,sumoz_spread_nh,sumoz_oberr_nh,nobsoz_nh,&
+   call printstats('  all  oz',sumoz_nh,biasoz_nh,sumoz_spread_nh,sumoz_oberr_nh,nobsoz_nh,&
         sumoz_sh,biasoz_sh,sumoz_spread_sh,sumoz_oberr_sh,nobsoz_sh,&
         sumoz_tr,biasoz_tr,sumoz_spread_tr,sumoz_oberr_tr,nobsoz_tr)
 end if ! nobs_conv+nobs_oz > 0

src/gsi/correlated_obsmod.F90

@@ -55,11 +55,11 @@
 \begin{verbatim}
 correlated_observations::
 ! isis       method   kreq   kmut  type    cov_file
-  airs281_aqua  1      60.   1.0   ice     airs_rcov.bin
-  airs281_aqua  1      60.   1.0   land    airs_rcov.bin
-  airs281_aqua  1      60.   1.0   sea     airs_rcov.bin
-  airs281_aqua  1      60.   1.0   snow    airs_rcov.bin
-  airs281_aqua  1      60.   1.0   mixed   airs_rcov.bin
+  airs_aqua  1      60.   1.0   ice     airs_rcov.bin
+  airs_aqua  1      60.   1.0   land    airs_rcov.bin
+  airs_aqua  1      60.   1.0   sea     airs_rcov.bin
+  airs_aqua  1      60.   1.0   snow    airs_rcov.bin
+  airs_aqua  1      60.   1.0   mixed   airs_rcov.bin
 # cris_npp      1     -99.   1.0   snow    cris_rcov.bin
 # cris_npp      1     -99.   1.0   land    cris_rcov.bin
 # cris_npp      1     -99.   1.0   sea     cris_rcov.bin
@@ -201,13 +201,13 @@ module correlated_obsmod
 type ObsErrorCov
      character(len=40) :: name                        ! R covariance name
      character(len=20) :: instrument                  ! instrument
-     integer(i_kind)   :: nch_active=-1               ! active channels
-     integer(i_kind)   :: nctot=-1                    ! total number of channels (active+passive)
+     integer(i_kind)   :: nch_active=-1               ! number of channels actively assimilated, according to the satinfo
+     integer(i_kind)   :: nctot=-1                    ! total number of channels (active+passive), according to the covariance file
      integer(i_kind)   :: method    =-1               ! define method of computation
      real(r_kind)      :: kreq      =-99._r_kind      ! Weston et al-like spectrum adjustment factor
      real(r_kind)      :: kmut      =-99._r_kind      ! multiplicative inflation factor
      character(len=20) :: mask      ='global'         ! Apply covariance for profiles over all globe
-     integer(i_kind),pointer :: indxR(:)   =>NULL()   ! indexes of active channels in between 1 and nchanl
+     integer(i_kind),pointer :: indxR(:)   =>NULL()   ! indexes of active channels in between 1 and nchanl, according to the satinfo
      real(r_kind),   pointer :: R(:,:)     =>NULL()   ! nch_active x nch_active
      real(r_kind),   pointer :: Revals(:)  =>NULL()   ! eigenvalues of R
 end type
@@ -333,7 +333,8 @@ end subroutine ini_
 ! !INTERFACE:
 !
 subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov)
-use radinfo, only: nusis,iuse_rad,jpch_rad
+use radinfo, only: nusis,iuse_rad,jpch_rad,varch,nuchan
+use constants, only: zero
 implicit none
 
 ! !INPUT PARAMETERS:
@@ -368,10 +369,13 @@ subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov)
 !BOC
 
 character(len=*),parameter :: myname_=myname//'*set'
-integer(i_kind) nch_active,lu,ii,ioflag,iprec,nctot,coun
-
-real(r_single),allocatable, dimension(:,:) :: readR4  ! nch_active x nch_active x ninstruments
-real(r_double),allocatable, dimension(:,:) :: readR8  ! nch_active x nch_active x ninstruments
+integer(i_kind) nch_active   !number of channels accounted for in the covariance file
+integer(i_kind) nctot        !the total number of channels (active+passive), according to the covariance file
+integer(i_kind) lu,ii,jj,ioflag,iprec,coun,couns,istart,indR,nctotf
+integer(i_kind),dimension(:),allocatable:: indxR,indxRf !channel indices read in from the covariance file
+real(r_kind),dimension(:,:),allocatable:: Rcov   !covariance matrix read in from the covariance file
+real(r_single),allocatable, dimension(:,:) :: readR4  ! nch_active x nch_active 
+real(r_double),allocatable, dimension(:,:) :: readR8  ! nch_active x nch_active 
 real(r_kind),allocatable, dimension(:) :: diag
 logical :: corr_obs
 
@@ -393,56 +397,120 @@ subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov)
          read(lu,IOSTAT=ioflag) nch_active, nctot, iprec
       endif
       if(ioflag/=0) call die(myname_,' failed to read nch from '//trim(fname))
-!     if no data available, turn off Correlated Error
       coun=0
+      couns=0
+      istart=0 
+      nctotf=0
       do ii=1,jpch_rad
         if (nusis(ii)==ErrorCov%instrument) then
-           if (iuse_rad(ii)>0) coun=coun+1
+           if (couns==0) then 
+              istart=ii-1
+              couns=1
+           endif
+           if (iuse_rad(ii)>0) then
+              coun=coun+1
+           endif
+           nctotf=nctotf+1
         endif
       enddo
-      if (coun<nch_active) then
+!     if no data available, turn off Correlated Error
+      if (coun==0) then
          if (iamroot_) write(6,*) 'WARNING: ',trim(ErrorCov%instrument), &
                        ' is not initiallized. Turning off Correlated Error'
          return
       endif
-      ErrorCov%nch_active = nch_active
+      ErrorCov%nch_active = coun
       if (.not.GMAO_ObsErrorCov) ErrorCov%nctot = nctot
-      call create_(nch_active,ErrorCov)
+      call create_(coun,ErrorCov)
+      allocate(indxRf(nch_active),indxR(nch_active),Rcov(nctot,nctot))
 
 !     Read GSI-like channel numbers used in estimating R for this instrument
-      read(lu,IOSTAT=ioflag) ErrorCov%indxR
+      read(lu,IOSTAT=ioflag) indxR
       if(ioflag/=0) call die(myname_,' failed to read indx from '//trim(fname))
 
 !     Read estimate of observation error covariance
+      Rcov=0.0_r_kind
       if(iprec==4) then
         allocate(readR4(nch_active,nch_active))
         read(lu,IOSTAT=ioflag) readR4
         if(ioflag/=0) call die(myname_,' failed to read R from '//trim(fname))
-        ErrorCov%R = readR4
+        Rcov(1:nch_active,1:nch_active)=readR4
         deallocate(readR4)
       endif
       if(iprec==8) then
         allocate(readR8(nch_active,nch_active))
         read(lu,IOSTAT=ioflag) readR8
         if(ioflag/=0) call die(myname_,' failed to read R from '//trim(fname))
-        ErrorCov%R = readR8
+        Rcov(1:nch_active,1:nch_active)=readR8
         deallocate(readR8)
       endif
-
+      if (GMAO_ObsErrorCov) then
+         ErrorCov%indxR(1:nch_active)=indxR(1:nch_active)
+         ErrorCov%nch_active=nch_active
+      else
+         coun=0
+         ErrorCov%R=zero
+         do ii=1,nctotf 
+            if (iuse_rad(ii+istart)>0) then
+               coun=coun+1
+               ErrorCov%indxR(coun)=ii
+               indxRf(coun)=nuchan(ii+istart)
+            endif
+         enddo
+!Add rows and columns for active channels in the satinfo that are not in the covariance file
+        couns=1
+        do ii=1,ErrorCov%nch_active
+           indR=0
+           do jj=couns,nch_active 
+              if (indxR(jj)==indxRf(ii)) then
+                 indR=jj
+                 couns=jj+1
+                 exit
+              endif
+           enddo
+           if (indR==0) then
+              do jj=nctot-1,ii,-1
+                 Rcov(jj+1,:)=Rcov(jj,:)
+                 Rcov(:,jj+1)=Rcov(:,jj)
+              enddo
+              Rcov(ii,:)=zero
+              Rcov(:,ii)=zero
+              Rcov(ii,ii)=varch(istart+ErrorCov%indxR(ii))*varch(istart+ErrorCov%indxR(ii))
+           endif
+        enddo
+!Remove rows and columns that are in the covariance file, but not in the satinfo
+         couns=1
+         do ii=1,nch_active
+           indR=0
+           do jj=couns,ErrorCov%nch_active
+              if (indxRf(jj)==indxR(ii)) then
+                 indR=jj
+                 couns=jj+1
+                 exit
+              endif
+           enddo
+           if (indR==0) then
+              do jj=ii,nctot-1
+                 Rcov(jj,:)=Rcov(jj+1,:)
+                 Rcov(:,jj)=Rcov(:,jj+1) 
+              enddo
+            endif
+         enddo
+      endif
+      ErrorCov%R(1:ErrorCov%nch_active,1:ErrorCov%nch_active)=Rcov(1:ErrorCov%nch_active,1:ErrorCov%nch_active)
 !     Done reading file
       close(lu)
    else
       if (iamroot_) write(6,*) 'No Rcov files found.  Turning off Correlated Error'
       return
    end if
-
 !  If method<0 there is really nothing else to do
 !  ----------------------------------------------
    if (method<0) then
       initialized_=.true.
       return
    endif
-
+   nch_active=ErrorCov%nch_active
    if (VERBOSE_) then
        allocate(diag(nch_active))
        do ii=1,nch_active
@@ -471,7 +539,7 @@ subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov)
        endif
        deallocate(diag)
    endif
-
+   deallocate(indxR,Rcov,indxRf)
    initialized_=.true.
 end subroutine set_
 !EOC
@@ -906,9 +974,6 @@ subroutine upd_varch_
                do jj=1,nch_active
                   nn=GSI_BundleErrorCov(itbl)%indxR(jj)
                   mm=ich1(nn)
-                  if( iuse_rad(mm)<1 ) then
-                    call die(myname_,' active channels used in R do not match those used in GSI, aborting')
-                  endif
                   if(isurf==1) then 
                     if(iamroot_)write(6,'(1x,a6,a20,2i6,2f20.15)')'>>>',idnames(itbl),jj,nn,varch(mm),sqrt(GSI_BundleErrorCov(itbl)%R(jj,jj))
                     varch_sea(mm)=sqrt(GSI_BundleErrorCov(itbl)%R(jj,jj))
@@ -945,9 +1010,6 @@ subroutine upd_varch_
                   endif 
                enddo
                ncp=count(ircv>0) ! number of active channels in profile
-               if(ncp/=nch_active) then
-                  call die(myname_,'serious inconsistency in handling correlated obs')
-               endif
                allocate(IRsubset(ncp)) ! these indexes apply to the matrices/vec in ErrorCov
                allocate(IJsubset(ncp)) ! these indexes in 1 to nchanl
                iii=0;jjj=0

src/gsi/gsimod.F90

@@ -34,7 +34,7 @@ module gsimod
      obsmod_init_instr_table,obsmod_final_instr_table
   use obsmod, only: luse_obsdiag
   use obsmod, only: netcdf_diag, binary_diag
-  use obsmod, only: l_wcp_cwm
+  use obsmod, only: l_wcp_cwm,ompslp_mult_fact
   use obsmod, only: aircraft_recon, &
 
        ! The following variables are the coefficients that describe
@@ -936,6 +936,7 @@ module gsimod
 !     vis_thres  - threshold value for both vis observation and input first guess
 !     cldch_thres  - threshold value for both cldch observation and input first guess
 !     cld_det_dec2bin - re-interprets cld_det in satinfo as binary entries
+!     ompslp_mult_fact - multiplication factor for OMPS LP obserror read in
 
 ! The following variables are the coefficients that describe the
 ! linear regression fits that are used to define the dynamic
@@ -986,7 +987,7 @@ module gsimod
        vadfile,noiqc,c_varqc,blacklst,use_poq7,hilbert_curve,tcp_refps,tcp_width,&
        tcp_ermin,tcp_ermax,qc_noirjaco3,qc_noirjaco3_pole,qc_satwnds,njqc,vqc,nvqc,hub_norm,troflg,lat_c,nrand,&
        aircraft_t_bc_pof,aircraft_t_bc,aircraft_t_bc_ext,biaspredt,upd_aircraft,cleanup_tail,&
-       hdist_aircraft,buddycheck_t,buddydiag_save,vadwnd_l2rw_qc,  &
+       hdist_aircraft,buddycheck_t,buddydiag_save,vadwnd_l2rw_qc,ompslp_mult_fact,  &
        pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres,cld_det_dec2bin, &
        q_doe_a_136,q_doe_a_137,q_doe_b_136,q_doe_b_137, &
        t_doe_a_136,t_doe_a_137,t_doe_b_136,t_doe_b_137, &

src/gsi/obsmod.F90

@@ -530,6 +530,10 @@ module obsmod
 
   public :: vad_near_analtime
 
+  ! The following correspond to OMPS LP observations:
+
+  public :: ompslp_mult_fact
+
   interface obsmod_init_instr_table
           module procedure init_instr_table_
   end interface
@@ -679,6 +683,10 @@ module obsmod
 
   logical vad_near_analtime 
 
+ ! The following correspond to OMPS LP observations:
+
+  real(r_kind) :: ompslp_mult_fact
+
 contains
 
   subroutine init_obsmod_dflts
@@ -951,6 +959,10 @@ subroutine init_obsmod_dflts
     ! see 'read_prepbufr.f90'
     vad_near_analtime = .false.
 
+    ! The following correspond to OMPS LP  observations:
+
+    ompslp_mult_fact = 2.0_r_kind
+
     return
   end subroutine init_obsmod_dflts
 

src/gsi/read_nsstbufr.f90

@@ -266,6 +266,8 @@ subroutine read_nsstbufr(nread,ndata,nodata,gstime,infile,obstype,lunout, &
         elseif ( trim(subset) == 'NC001007' ) then            ! CSTGD
            msst = 16.0_r_kind                                 ! for CSTGD, assign to be 16
            call ufbint(lunin,sst,1,1,iret,'SST1')             ! read SST
+        else
+           cycle read_loop
         endif
 
           call ufbint(lunin,loc,4,1,iret,'CLAT CLATH CLON CLONH')

src/gsi/read_ozone.f90

@@ -103,7 +103,7 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, &
   use satthin, only: radthin_time_info,tdiff2crit
   use gridmod, only: nlat,nlon,regional,tll2xy,rlats,rlons
   use constants, only: deg2rad,zero,one_tenth,r60inv,two
-  use obsmod, only: nloz_v6,nloz_v8
+  use obsmod, only: nloz_v6,nloz_v8, ompslp_mult_fact
   use obsmod, only: time_window_max
   use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen
   use radinfo, only: dec2bin
@@ -1219,7 +1219,7 @@ subroutine read_ozone(nread,ndata,nodata,jsatid,infile,gstime,lunout, &
 
          ozout(8,ndata)=usage1(k)          ! 
          ozout(9,ndata)=log(press(k))      ! ompslp pressure in log(cb)
-         ozout(10,ndata)=omrstd(k)         ! ozone mixing ratio precision in ppmv
+         ozout(10,ndata)=omrstd(k)*ompslp_mult_fact   ! ozone mixing ratio precision in ppmv
          ozout(11,ndata)=float(ipos(k))    ! pointer of obs level index in 
                                            ! ozinfo.txt
          ozout(12,ndata)=j !nloz              ! # of ompslp vertical levels