diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 7177719fde..ca82893398 100644 --- a/src/gsi/gsimod.F90 +++ b/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, & diff --git a/src/gsi/lagmod.f90 b/src/gsi/lagmod.f90 index 965a499fdc..8cd8e2d8e3 100644 --- a/src/gsi/lagmod.f90 +++ b/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) diff --git a/src/gsi/pcgsoi.f90 b/src/gsi/pcgsoi.f90 index 0f7f3d32d8..a4ae2431b1 100644 --- a/src/gsi/pcgsoi.f90 +++ b/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 (b10.0_r_kind) then + if (b30.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 diff --git a/src/gsi/qcmod.f90 b/src/gsi/qcmod.f90 index 525da81979..db6b6072a4 100644 --- a/src/gsi/qcmod.f90 +++ 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 diff --git a/src/gsi/setupbend.f90 b/src/gsi/setupbend.f90 index 9a3cf0a4c9..9bc856d67a 100644 --- a/src/gsi/setupbend.f90 +++ b/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 (hobrsig) then + if (hobrsig) 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+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