Merge pull request #354 from ADCollard/RO_Updates

GitHub Issue #353. RO Updates

src/gsi/gsimod.F90

@@ -89,7 +89,7 @@ module gsimod
   use turblmod, only: use_pbl,init_turbl
   use qcmod, only: dfact,dfact1,create_qcvars,destroy_qcvars,&
       erradar_inflate,tdrerr_inflate,use_poq7,qc_satwnds,&
-      init_qcvars,vadfile,noiqc,c_varqc,qc_noirjaco3,qc_noirjaco3_pole,&
+      init_qcvars,vadfile,noiqc,c_varqc,gps_jacqc,qc_noirjaco3,qc_noirjaco3_pole,&
       buddycheck_t,buddydiag_save,njqc,vqc,nvqc,hub_norm,vadwnd_l2rw_qc, &
       pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres,cao_check
   use qcmod, only: troflg,lat_c,nrand
@@ -935,6 +935,7 @@ module gsimod
 !     tcp_width  - parameter for tcps oberr inflation (width, mb)
 !     tcp_ermin  - parameter for tcps oberr inflation (minimum oberr, mb)
 !     tcp_ermax  - parameter for tcps oberr inflation (maximum oberr, mb)
+!     gps_jacqc  - logical to turn on GNSSRO Jacobian QC (default is off)
 !     qc_noirjaco3 - controls whether to use O3 Jac from IR instruments
 !     qc_noirjaco3_pole - controls wheter to use O3 Jac from IR instruments near poles
 !     qc_satwnds - allow bypass sat-winds qc normally removing lots of mid-tropo obs
@@ -1009,7 +1010,7 @@ module gsimod
 
   namelist/obsqc/dfact,dfact1,erradar_inflate,tdrerr_inflate,oberrflg,&
        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,&
+       tcp_ermin,tcp_ermax,gps_jacqc,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,ompslp_mult_fact,  &
        pvis,pcldch,scale_cv,estvisoe,estcldchoe,vis_thres,cldch_thres,cld_det_dec2bin, &

src/gsi/lagmod.f90

@@ -12,7 +12,8 @@ module lagmod
 !   2005-12-01  cucurull - implement tangent linear and adjoint code 
 !                        - adapt the code to GSI standards
 !   2008-05-09  safford  - complete documentation block
-! 
+!   2021-11-04  cucurull - bug fix
+!
 ! subroutines included:
 !   setq
 !   setq_TL
@@ -717,6 +718,13 @@ SUBROUTINE slagdw_AD(x,x_AD,xt,aq,aq_AD,bq,bq_AD,w_AD,dw,dw_AD,n)
    dwb_AD=zero
 ENDIF
 
+!Criterion for target lying outside the central interval:                                                                                                
+if(dwb==zero)then
+   wb_AD =zero
+   dwb_AD=zero
+endif
+
+
 xa_AD        =xa_AD-wb_AD*dwb
 dwb_AD       =dwb_AD+(xt-xa)*wb_AD
 xb_AD        =xb_AD-(dwb_AD/(xb-xa)**2)

src/gsi/pcgsoi.f90

@@ -377,7 +377,7 @@ subroutine pcgsoi()
      if (.not. restart .or. iter > 0) then
         if (iter > 1 .or. .not. read_success)then
            if (gsave>1.e-16_r_kind) b=gnorm(2)/gsave
-           if (b<zero .or. b>10.0_r_kind) then
+           if (b<zero .or. b>30.0_r_kind) then
               if (mype==0) then
                  if (iout_6) write(6,105) gnorm(2),gsave,b
                  write(iout_iter,105) gnorm(2),gsave,b

src/gsi/qcmod.f90

@@ -184,7 +184,7 @@ module qcmod
   public :: qc_saphir
 ! set passed variables to public
   public :: npres_print,nlnqc_iter,varqc_iter,pbot,ptop,c_varqc,njqc,vqc,nvqc,hub_norm
-  public :: use_poq7,noiqc,vadfile,dfact1,dfact,erradar_inflate
+  public :: use_poq7,noiqc,vadfile,dfact1,dfact,erradar_inflate,gps_jacqc
   public :: pboto3,ptopo3,pbotq,ptopq,newvad,tdrerr_inflate
   public :: igood_qc,ifail_crtm_qc,ifail_satinfo_qc,ifail_interchan_qc,&
             ifail_gross_qc,ifail_cloud_qc,ifail_outside_range,&
@@ -205,6 +205,7 @@ module qcmod
   logical use_poq7
   logical qc_noirjaco3
   logical qc_noirjaco3_pole
+  logical gps_jacqc
   logical newvad
   logical tdrerr_inflate
   logical qc_satwnds
@@ -429,6 +430,8 @@ subroutine init_qcvars
     use_poq7 = .false.
     cao_check = .false.
 
+    gps_jacqc = .false.     ! Jacobian QC for GNSS RO is off by default
+
     qc_noirjaco3 = .false.  ! when .f., use O3 Jac from IR instruments
     qc_noirjaco3_pole = .false. ! true=do not use O3 Jac from IR instruments near poles
 

src/gsi/setupbend.f90

@@ -104,6 +104,10 @@ subroutine setupbend(obsLL,odiagLL, &
 !   2020-08-26  Shao/Bathmann - add Jacobian QC
 !   2021-07-29  cucurull - revert gross error check to default values
 !   2021-07-29  cucurull - fix forward operator issues identified with L127
+!   2021-11-04  cucurull - turn off Jacbian QC
+!   2021-11-05  cucurull - update QCs and optimize/improve forward operator; bug fixes
+!   2022-01-28  cucurull - add Sentinel-6, PAZ
+!   2022-04-06  collard  - reintroduce Jacbian QC as an option (default off)
 !
 !   input argument list:
 !     lunin    - unit from which to read observations
@@ -125,6 +129,7 @@ subroutine setupbend(obsLL,odiagLL, &
   use obsmod , only: nprof_gps,lobsdiag_allocated,&
       lobsdiagsave,nobskeep,&
       time_offset,lobsdiag_forenkf
+  use qcmod, only: gps_jacqc
   use m_obsNode, only: obsNode
   use m_gpsNode , only: gpsNode
   use m_gpsNode , only: gpsNode_appendto
@@ -198,6 +203,7 @@ subroutine setupbend(obsLL,odiagLL, &
   real(r_kind),parameter:: r12=12.0_r_kind
   real(r_kind),parameter:: r20=20.0_r_kind
   real(r_kind),parameter:: r40=40.0_r_kind
+  real(r_kind),parameter:: r80=80.0_r_kind
   real(r_kind),parameter:: r1em3 = 1.0e-3_r_kind
   real(r_kind),parameter:: r1em6 = 1.0e-6_r_kind
   character(len=*),parameter :: myname='setupbend'
@@ -294,6 +300,7 @@ subroutine setupbend(obsLL,odiagLL, &
 !267 => PlanetiQ GNOMES-A
 !268 => PlanetiQ GNOMES-B
 !269 => Spire Lemur 3U CubeSat
+!66 => Sentinel-6 
 
 ! Check to see if required guess fields are available
   call check_vars_(proceed)
@@ -332,7 +339,7 @@ subroutine setupbend(obsLL,odiagLL, &
   mm1=mype+1
   ns=nsig/two
   nsigstart=nint(ns)
-  ns=(r61/r63)*nsig+r18
+  ns=r80
   grids_dim=nint(ns)  ! grid points for integration of GPS bend
   ds=r10000
   allocate(ddnj(grids_dim),grid_s(grids_dim),ref_rad_s(grids_dim)) 
@@ -468,7 +475,6 @@ subroutine setupbend(obsLL,odiagLL, &
      top_layer_SR=0
      bot_layer_SR=0
 
-     alt=(tpdpres(i)-rocprof)*r1em3
 !$omp parallel do  schedule(dynamic,1) private(k,qmean,tmean,fact,pw,pressure,nrefges1,nrefges2,nrefges3)
      do k=1,nsig 
         zges(k) = (termr*hges(k)) / (termrg-hges(k))  ! eq (23) at interface (topo corrected)
@@ -511,21 +517,18 @@ subroutine setupbend(obsLL,odiagLL, &
            (k4/(tmean**2*pw))*fact*qmean*pressure
 
      end do
-     alt=(tpdpres(i)-rocprof)*r1em3
-     if (alt<= six) then 
+     alt=(tpdpres(i)-rocprof-zsges-unprof)*r1em3
+     if (alt<= five) then 
         do k=nsigstart,1,-1 
 !       check for model SR layer at obs location
            grad_mod=1000.0_r_kind*(nrefges(k+1,i)-nrefges(k,i))/(rges(k+1,i)-rges(k,i))
            if (abs(grad_mod)>= half*crit_grad) then  ! SR - likely, to be used in obs SR qc
               qc_layer_SR=.true.   !SR-likely layer detected
            endif
-!           if (((ref_rad(k+1)-ref_rad(k))/(rges(k+1,i)-rges(k,i))) < zero) then
-           if (abs(grad_mod)>= 0.75_r_kind*crit_grad) then  !relax to close-to-SR conditions
-              count_SR=count_SR+1 ! layers of SR
+           if (abs(grad_mod) >= 0.75_r_kind*crit_grad) then  !relax to close-to-SR conditions
+              count_SR=count_SR+1 ! layers of SR                                           
               if (count_SR > 1 ) then
-!                 if(abs(bot_layer_SR-k) > 1) write(6,*) 'WARNING GPSRO: non-consecutive SR layers'
                  bot_layer_SR=k
-
               else
                  top_layer_SR=k
                  bot_layer_SR=top_layer_SR
@@ -543,7 +546,7 @@ subroutine setupbend(obsLL,odiagLL, &
      error(i)=tiny_r_kind
      error_adjst(i)=tiny_r_kind
 
-     if (hob<one .or. hob>rsig) then
+     if (hob<three .or. hob>rsig) then
         data(ier,i) = zero
         ratio_errors(i) = zero
         muse(i)=.false.
@@ -582,23 +585,22 @@ subroutine setupbend(obsLL,odiagLL, &
 
      if(ratio_errors(i) > tiny_r_kind)  then ! obs inside model grid
 
-       if (alt <= six) then
-          if (top_layer_SR >= 1) then ! SR exists for at least one layer. Check if obs is inside 
-             if (tpdpres(i)==ref_rad(top_layer_SR+1)) then !obs inside model SR layer
-                qcfail(i)=.true.
-             elseif (tpdpres(i) < ref_rad(top_layer_SR+1)) then !obs below model close-to-SR layer
-                qcfail(i)=.true.
-             elseif (tpdpres(i) >= ref_rad(top_layer_SR+1) .and. tpdpres(i) <= ref_rad(top_layer_SR+5)) then !source too close
-                qcfail(i)=.true.
-             else !above
-                qcfail(i)=.false.
-                if(hob < top_layer_SR+1) then !correct if obs location is below non-monotonic section
-                   hob = tpdpres(i)
-                   call grdcrd1(hob,ref_rad(top_layer_SR+1),nsig-top_layer_SR-1,1)
-                   data(ihgt,i) = hob+top_layer_SR
-                   hob = hob+top_layer_SR
-                   rdiagbuf(19,i) = hob
-                endif
+       if (alt <= five) then
+          if (top_layer_SR >= 1) then ! SR exists for at least one layer. Check if obs is inside                                                                       
+             if ((tpdpres(i)<ref_rad(top_layer_SR+1)) .and. &
+                    (tpdpres(i)<ref_rad(bot_layer_SR))) then !obs below model SR/close-to-SR layer
+                    qcfail(i)=.true.
+             elseif (tpdpres(i)>ref_rad(top_layer_SR+5)) then ! obs above SR/close-to-SR layer
+                    qcfail(i)=.false.
+                    if(hob < top_layer_SR+1) then !location might be aliased to the lower section of the non-monotonicity 
+                      hob = tpdpres(i)
+                      call grdcrd1(hob,ref_rad(top_layer_SR+1),nsig-top_layer_SR-1,1) ! only non-monotonic section above SR layer 
+                      data(ihgt,i) = hob+top_layer_SR
+                      hob = hob+top_layer_SR
+                      rdiagbuf(19,i) = hob
+                    endif
+             else ! obs inside model SR/shadow or close-to-SR layer                                                                                                     
+                    qcfail(i)=.true.
              endif
           endif
 
@@ -609,6 +611,7 @@ subroutine setupbend(obsLL,odiagLL, &
           endif
        endif 
 
+       alt=(tpdpres(i)-rocprof)*r1em3
 !      get pressure (in mb), temperature and moisture at obs location
        call tintrp31(ges_lnprsi(:,:,1:nsig,:),dpressure,dlat,dlon,hob,&
               dtime,hrdifsig,mype,nfldsig)
@@ -636,7 +639,7 @@ subroutine setupbend(obsLL,odiagLL, &
             (data(isatid,i)==42) .or.(data(isatid,i)==3)  .or. &
             (data(isatid,i)==821).or.(data(isatid,i)==421).or. &
             (data(isatid,i)==440).or.(data(isatid,i)==43) .or. &
-            (data(isatid,i)==5)) then
+            (data(isatid,i)==5).or.(data(isatid,i)==66)) then
 
            if((data(ilate,i)> r40).or.(data(ilate,i)< -r40)) then
               if(alt>r12) then
@@ -720,10 +723,12 @@ subroutine setupbend(obsLL,odiagLL, &
            xj(j,i)=ref_rad_s(j)
            hob_s=ref_rad_s(j)
            call grdcrd1(hob_s,ref_rad(1),nsig_up,1)
-           if(hob_s < top_layer_SR+1) then !correct if wrong location
-              hob_s = ref_rad_s(j)
-              call grdcrd1(hob_s,ref_rad(top_layer_SR+1),nsig_up-top_layer_SR-1,1)
-              hob_s = hob_s+top_layer_SR
+           if(top_layer_SR >=1) then
+             if(hob_s < top_layer_SR+1) then !correct if wrong location
+                hob_s = ref_rad_s(j)
+                call grdcrd1(hob_s,ref_rad(top_layer_SR+1),nsig_up-top_layer_SR-1,1)
+                hob_s = hob_s+top_layer_SR
+             endif
            endif
            dbend_loc(j,i)=hob_s  !location of x_j with respect to extended x_i
 
@@ -746,8 +751,14 @@ subroutine setupbend(obsLL,odiagLL, &
                  dw4(1)=dw4(1)+dw4(4); dw4(2:4)=dw4(1:3);dw4(1)=zero
                  ihob=ihob-1
               endif
-              ddnj(j)=dot_product(dw4,nrefges(ihob-1:ihob+2,i))!derivative (dN/dx)_j
-              ddnj(j)=max(zero,abs(ddnj(j)))
+              ddnj(j)=dot_product(dw4,nrefges(ihob-1:ihob+2,i))!derivative (dN/dx)_j                                                                      
+              if(ddnj(j)>zero) then
+                 qcfail(i)=.true.
+                 data(ier,i) = zero
+                 ratio_errors(i) = zero
+                 muse(i)=.false.
+                 cycle loopoverobs1 ! reject observation                                                                                                  
+              endif
            else
               obs_check=.true.
               if (obs_check) then ! only once per obs here
@@ -780,12 +791,21 @@ subroutine setupbend(obsLL,odiagLL, &
             ddbend=ds*ddnj(j)/ref_rad_s(j)
             dbend=dbend+two*ddbend
          end do
-         dbend=r1em6*tpdpres(i)*dbend  
+         dbend=-r1em6*tpdpres(i)*dbend  
 
 !        Accumulate diagnostic information
          rdiagbuf( 5,i)  = (data(igps,i)-dbend)/data(igps,i) ! incremental bending angle (x100 %)
 
          data(igps,i)=data(igps,i)-dbend !innovation vector
+
+!        Remove very large simulated values
+         if(dbend > 0.05_r_kind) then
+           data(ier,i) = zero
+           ratio_errors(i) = zero
+           qcfail(i)=.true.
+           muse(i)=.false.
+         endif
+
          if (alt <= gpstop) then ! go into qc checks
             if ((alt <= commgpstop) .or. (.not.commdat)) then 
                cgrossuse=cgross(ikx)
@@ -1069,7 +1089,7 @@ subroutine setupbend(obsLL,odiagLL, &
 
            allocate(my_head)
            call gpsNode_appendto(my_head,gpshead(ibin))
-           
+
            my_head%idv = is
            my_head%iob = ioid(i)
            my_head%elat= data(ilate,i)
@@ -1194,14 +1214,13 @@ subroutine setupbend(obsLL,odiagLL, &
                                                    dbenddn(j) * dndp(j,k)
               end do
 
-              if ((abs(my_head%jac_t(k)) > 0.0016_r_kind).or.(abs(my_head%jac_q(k)) > 7.5_r_kind).or. &
-                  (abs(my_head%jac_p(k)) > 0.004_r_kind)) then
-                 qcfail_jac(i) = one
-              end if
+              if (gps_jacqc .and. ((abs(my_head%jac_t(k)) > 0.0016_r_kind).or.(abs(my_head%jac_q(k)) > 7.5_r_kind).or. &
+                  (abs(my_head%jac_p(k)) > 0.004_r_kind))) qcfail_jac(i) = one
+
            end do 
 
            my_head%jac_p(nsig+1) = zero
-   
+
            if (qcfail_jac(i) == one) then
               do k=1,nsig
                  my_head%jac_t(k) = zero
@@ -1212,8 +1231,7 @@ subroutine setupbend(obsLL,odiagLL, &
               data(ier,i) = zero
               rdiagbuf(12,i) = -one
               rdiagbuf(10,i) = six
-           end if 
-
+           end if
 
            if (save_jacobian) then
               ! fill in the jacobian