Skip to content

Commit

Permalink
Merge pull request #201 from KristenBathmann-NOAA/release/gfsda.v16.1.3
Browse files Browse the repository at this point in the history 
Issue #200:  updates to setupbend
  • Loading branch information
@RussTreadon-NOAA
RussTreadon-NOAA authored Aug 22, 2021
2 parents 3616ecb + 16b28c0 commit 75feef4
Showing 1 changed file with 65 additions and 59 deletions.
124 changes: 65 additions & 59 deletions src/gsi/setupbend.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,8 @@ subroutine setupbend(obsLL,odiagLL, &
! 2020-04-13 Shao - update the statistis QC for COSMIC-2
! 2020-05-21 Shao - add comments to include commercial data ID information
! 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
!
! input argument list:
! lunin - unit from which to read observations
Expand Down Expand Up @@ -586,14 +588,22 @@ subroutine setupbend(obsLL,odiagLL, &
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+2)) then !source too close
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
endif
endif

! check for SR in obs, will be updated in genstats.
if ( data(igps,i) >= 0.03 .and. qc_layer_SR) then
if ( data(igps,i) >= 0.03_r_kind .and. qc_layer_SR) then
kprof = data(iprof,i)
toss_gps_sub(kprof) = max (toss_gps_sub(kprof),data(igps,i))
endif
Expand All @@ -613,7 +623,8 @@ subroutine setupbend(obsLL,odiagLL, &
rdiagbuf( 6,i) = ten*exp(dpressure) ! pressure at obs location (hPa) if monotone grid
rdiagbuf(18,i) = trefges ! temperature at obs location (Kelvin) if monotone grid
rdiagbuf(21,i) = qrefges ! specific humidity at obs location (kg/kg) if monotone grid

commdat=.false.
if (data(isatid,i)>=265 .and. data(isatid,i)<=269) commdat=.true.
if (.not. qcfail(i)) then ! not SR

! Modify error to account for representativeness error.
Expand Down Expand Up @@ -642,19 +653,19 @@ subroutine setupbend(obsLL,odiagLL, &
endif
else
! CDAAC-type processing
if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
if (((data(isatid,i) > 749).and.(data(isatid,i) < 756)).or.commdat) then
if ((data(ilate,i)> r40).or.(data(ilate,i)< -r40)) then
if (alt.le.8.0_r_kind) then
if (alt <= 8.0_r_kind) then
repe_gps=-1.0304261_r_kind+0.3203316_r_kind*alt+0.0141337_r_kind*alt**2
elseif (alt.gt.8.0_r_kind.and.alt.le.r12) then
elseif (alt > 8.0_r_kind.and.alt <= r12) then
repe_gps=2.1750271_r_kind+0.0431177_r_kind*alt-0.0008567_r_kind*alt**2
else
repe_gps=-0.3447429_r_kind+0.2829981_r_kind*alt-0.0028545_r_kind*alt**2
endif
else
if (alt.le.4.0_r_kind) then
if (alt <= 4.0_r_kind) then
repe_gps=0.7285212_r_kind-1.1138755_r_kind*alt+0.2311123_r_kind*alt**2
elseif (alt.le.r18.and.alt.gt.4.0_r_kind) then
elseif (alt <= r18.and.alt > 4.0_r_kind) then
repe_gps=-3.3878629_r_kind+0.8691249_r_kind*alt-0.0297196_r_kind*alt**2
else
repe_gps=-2.3875749_r_kind+0.3667211_r_kind*alt-0.0037542_r_kind*alt**2
Expand All @@ -680,11 +691,7 @@ subroutine setupbend(obsLL,odiagLL, &

repe_gps=exp(repe_gps) ! one/modified error in (rad-1*1E3)
repe_gps= r1em3*(one/abs(repe_gps)) ! modified error in rad
commdat=.false.
if (data(isatid,i)>=265 .and. data(isatid,i)<=269) then
commdat=.true.
repe_gps=commgpserrinf*repe_gps ! Inflate error for commercial data
endif
if (commdat) repe_gps=commgpserrinf*repe_gps ! Inflate error for commercial data
ratio_errors(i) = data(ier,i)/abs(repe_gps)

error(i)=one/data(ier,i) ! one/original error
Expand Down Expand Up @@ -713,6 +720,11 @@ 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
endif
dbend_loc(j,i)=hob_s !location of x_j with respect to extended x_i


Expand Down Expand Up @@ -779,11 +791,6 @@ subroutine setupbend(obsLL,odiagLL, &
cgrossuse=cgross(ikx)
cermaxuse=cermax(ikx)
cerminuse=cermin(ikx)
if (alt > five) then
cgrossuse=cgrossuse*r400
cermaxuse=cermaxuse*r400
cerminuse=cerminuse*r100
endif
! Gross error check
obserror = one/max(ratio_errors(i)*data(ier,i),tiny_r_kind)
obserrlm = max(cerminuse,min(cermaxuse,obserror))
Expand All @@ -801,7 +808,7 @@ subroutine setupbend(obsLL,odiagLL, &
else
! Statistics QC check if obs passed gross error check
cutoff=zero
if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
if (((data(isatid,i) > 749).and.(data(isatid,i) < 756)).or.commdat) then
cutoff1=(-4.725_r_kind+0.045_r_kind*alt+0.005_r_kind*alt**2)*one/two
else
cutoff1=(-4.725_r_kind+0.045_r_kind*alt+0.005_r_kind*alt**2)*two/three
Expand All @@ -812,12 +819,12 @@ subroutine setupbend(obsLL,odiagLL, &
else
cutoff3=0.005_r_kind*trefges**2-2.3_r_kind*trefges+266_r_kind
endif
if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
if (((data(isatid,i) > 749).and.(data(isatid,i) < 756)).or.commdat) then
cutoff3=cutoff3*one/two
else
cutoff3=cutoff3*two/three
end if
if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
if (((data(isatid,i) > 749).and.(data(isatid,i) < 756)).or.commdat) then
cutoff4=(four+eight*cos(data(ilate,i)*deg2rad))*one/two
else
cutoff4=(four+eight*cos(data(ilate,i)*deg2rad))*two/three
Expand All @@ -836,7 +843,7 @@ subroutine setupbend(obsLL,odiagLL, &
if((alt<=six).and.(alt>four)) cutoff=cutoff34
if(alt<=four) cutoff=cutoff4

if ((data(isatid,i) > 749).and.(data(isatid,i) < 756)) then
if (((data(isatid,i) > 749).and.(data(isatid,i) < 756)).or.commdat) then
cutoff=two*cutoff*r0_01
else
cutoff=three*cutoff*r0_01
Expand Down Expand Up @@ -1017,43 +1024,43 @@ subroutine setupbend(obsLL,odiagLL, &

! Fill obs diagnostics structure
if (luse_obsdiag) then
call obsdiagNode_set(my_diag,wgtjo=(data(ier,i)*ratio_errors(i))**2, &
call obsdiagNode_set(my_diag,wgtjo=(data(ier,i)*ratio_errors(i))**2, &
jiter=jiter,muse=muse(i),nldepart=data(igps,i) )
endif

! Load additional obs diagnostic structure
ioff = mreal
if (lobsdiagsave) then
associate(odiag => my_diag )
do jj=1,miter
ioff=ioff+1
if (odiag%muse(jj)) then
rdiagbuf(ioff,i) = one
else
rdiagbuf(ioff,i) = -one
endif
enddo
do jj=1,miter+1
ioff=ioff+1
rdiagbuf(ioff,i) = odiag%nldepart(jj)
enddo
do jj=1,miter
ioff=ioff+1
rdiagbuf(ioff,i) = odiag%tldepart(jj)
enddo
do jj=1,miter
ioff=ioff+1
rdiagbuf(ioff,i) = odiag%obssen(jj)
enddo
end associate ! odiag
endif

! if obs is not "acceptable" and jacobian is not computed, fill jacobian
! with zeros
if (save_jacobian) then
dhx_dx%val = 0._r_kind
call writearray(dhx_dx, rdiagbuf(ioff+1:nreal,i))
endif
associate(odiag => my_diag )
do jj=1,miter
ioff=ioff+1
if (odiag%muse(jj)) then
rdiagbuf(ioff,i) = one
else
rdiagbuf(ioff,i) = -one
endif
enddo
do jj=1,miter+1
ioff=ioff+1
rdiagbuf(ioff,i) = odiag%nldepart(jj)
enddo
do jj=1,miter
ioff=ioff+1
rdiagbuf(ioff,i) = odiag%tldepart(jj)
enddo
do jj=1,miter
ioff=ioff+1
rdiagbuf(ioff,i) = odiag%obssen(jj)
enddo
end associate ! odiag
endif

! if obs is not "acceptable" and jacobian is not computed, fill jacobian
! with zeros
if (save_jacobian) then
dhx_dx%val = 0._r_kind
call writearray(dhx_dx, rdiagbuf(ioff+1:nreal,i))
endif

! If obs is "acceptable", load array with obs info for use
! in inner loop minimization (int* and stp* routines)
Expand Down Expand Up @@ -1197,9 +1204,9 @@ subroutine setupbend(obsLL,odiagLL, &

if (qcfail_jac(i) == one) then
do k=1,nsig
my_head%jac_t(k) = zero
my_head%jac_q(k) = zero
my_head%jac_p(k) = zero
my_head%jac_t(k) = zero
my_head%jac_q(k) = zero
my_head%jac_p(k) = zero
end do
ratio_errors(i) = zero
data(ier,i) = zero
Expand All @@ -1220,7 +1227,6 @@ subroutine setupbend(obsLL,odiagLL, &
ioff = ioff + size(dhx_dx)
endif

my_head%jac_p(nsig+1) = zero
my_head%raterr2= ratio_errors(i)**2
my_head%res = data(igps,i)
my_head%err2 = data(ier,i)**2
Expand All @@ -1243,7 +1249,7 @@ subroutine setupbend(obsLL,odiagLL, &
gps_alltail(ibin)%head%obserr = data(ier,i)
gps_alltail(ibin)%head%dataerr = data(ier,i)*data(igps,i)
gps_alltail(ibin)%head%muse = muse(i) ! logical
endif ! (last_pass)
endif ! (last_pass)
end do ! i=1,nobs
deallocate(ddnj,grid_s,ref_rad_s)
! Release memory of local guess arrays
Expand Down

0 comments on commit 75feef4

Please sign in to comment.