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

updates to setupbend #201

Merged
merged 1 commit into from
Aug 22, 2021
Merged

updates to setupbend #201

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
124 changes: 65 additions & 59 deletions src/gsi/setupbend.f90
View file
Open in desktop
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