From 16b28c0664142fc0fcd4bfb7f252c72fa11a8007 Mon Sep 17 00:00:00 2001 From: Kristen Bathmann Date: Thu, 19 Aug 2021 15:31:20 +0000 Subject: [PATCH] updates to setupbend to change commercial ro obs errors/qc and add Lidia's fixes --- src/gsi/setupbend.f90 | 124 ++++++++++++++++++++++-------------------- 1 file changed, 65 insertions(+), 59 deletions(-) diff --git a/src/gsi/setupbend.f90 b/src/gsi/setupbend.f90 index 885bcfd29..28e74fb31 100644 --- a/src/gsi/setupbend.f90 +++ b/src/gsi/setupbend.f90 @@ -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 @@ -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 @@ -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. @@ -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 @@ -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 @@ -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 @@ -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)) @@ -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 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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 @@ -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