Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GitHub Issue NOAA-EMC/GSI#353. RO Updates #354

Merged
merged 1 commit into from
Apr 20, 2022
Merged

GitHub Issue NOAA-EMC/GSI#353. RO Updates #354

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions src/gsi/gsimod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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, &
Expand Down
10 changes: 9 additions & 1 deletion src/gsi/lagmod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion src/gsi/pcgsoi.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
5 changes: 4 additions & 1 deletion src/gsi/qcmod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,&
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand Down
104 changes: 61 additions & 43 deletions src/gsi/setupbend.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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'
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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.
Expand Down Expand Up @@ -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

Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down