module stppcpmod !$$$ module documentation block ! . . . . ! module: stppcpmod module for stppcp and its tangent linear stppcp_tl ! prgmmr: ! ! abstract: module for stppcp and its tangent linear stppcp_tl ! ! program history log: ! 2005-05-17 Yanqiu zhu - wrap stppcp and its tangent linear stppcp_tl into one module ! 2005-11-16 Derber - remove interfaces ! 2008-12-02 Todling - remove stppcp_tl ! 2009-08-12 lueken - update documentation ! 2010-05-13 todling - uniform interface across stp routines ! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting ! ! subroutines included: ! sub stppcp ! ! attributes: ! language: f90 ! machine: ! !$$$ end documnetation block implicit none PRIVATE PUBLIC stppcp contains subroutine stppcp(pcphead,dval,xval,out,sges,nstep) !$$$ subprogram documentation block ! . . . . ! subprogram: stppcp compute contribution to penalty and ! stepsize from pcp, with nonlinear qc ! prgmmr: treadon org: np23 date: 2003-09-13 ! ! abstract: compute contribution to penalty and stepsize from precipitation ! observations ! ! program history log: ! 2003-12-18 treadon - initial routine ! 2004-06-15 treadon - update documentation ! 2004-07-29 treadon - add only to module use, add intent in/out ! 2004-10-07 parrish - add nonlinear qc option ! 2005-04-11 treadon - merge stppcp and stppcp_qc into single routine ! 2005-09-28 derber - modify var qc and change location and weight arrays ! 2006-07-28 derber - modify to use new inner loop obs data structure ! - unify NL qc ! 2006-09-18 derber - modify output b1 and b3 ! 2007-01-19 derber - limit pcp_ges* > zero ! 2007-03-19 tremolet - binning of observations ! 2007-05-10 tremolet - add opt to run as linear procedure ! 2007-06-04 derber - use quad precision to get reproducability over number of processors ! 2008-06-02 safford - rm unused var and uses ! 2008-12-03 todling - changed handling of ptr%time ! 2009-01-26 todling - re-implement Tremolet's linearization for q1fy10 ! 2010-01-04 zhang,b - bug fix: accumulate penalty for multiple obs bins ! 2010-03-25 zhu - use state_vector in the interface; ! - add handlings of cw case; add pointer_state ! 2010-05-13 todling - update to use gsi_bundle ! - on-the-spot handling of non-essential vars ! 2010-09-25 todling - fix linearization ! 2011-11-01 eliu - add handling for ql and qi increments ! ! input argument list: ! pcphead ! rt - search direction for temperature ! rq - search direction for moisture ! ru - search direction for zonal wind ! rv - search direction for meridional wind ! rql - search direction for cloud liquid water mixing ratio ! rqi - search direction for cloud ice water mixing ratio ! st - input temperature correction field ! sq - input q correction field ! su - input u correction field ! sv - input v correction field ! sql - input ql correction field ! sqi - input qi correction field ! sges - step size estimates (nstep) ! nstep - number of stepsizes (==0 means use outer iteration values) ! ! output argument list: ! out(1:nstep) - contribution to penalty from precipitation rate - sges(1:nstep) ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind,r_quad use pcpinfo, only: b_pcp,pg_pcp,tinym1_obs use constants, only: zero,one,half,two,tiny_r_kind,cg_term,zero_quad,r3600 use qcmod, only: nlnqc_iter,varqc_iter use gridmod, only: latlon11,nsig use gsi_4dvar, only: ltlint use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer use m_obsNode, only: obsNode use m_pcpNode, only: pcpNode use m_pcpNode, only: pcpNode_typecast use m_pcpNode, only: pcpNode_nextcast implicit none ! Declare passed variables class(obsNode), pointer ,intent(in ) :: pcphead integer(i_kind) ,intent(in ) :: nstep real(r_kind),dimension(max(1,nstep)),intent(in ) :: sges real(r_quad),dimension(max(1,nstep)),intent(inout) :: out type(gsi_bundle),intent(in) :: dval type(gsi_bundle),intent(in) :: xval ! Declare local variables logical:: lcld integer(i_kind) n,ncwm,nq,nt,nu,nv,kx,ier,istatus,icw,iql,iqi integer(i_kind) j1,j2,j3,j4,kk real(r_kind) dt,dt0,w1,w2,w3,w4 real(r_kind) dq,dq0 real(r_kind) du,du0 real(r_kind) dv,dv0 real(r_kind) dcwm,dcwm0 real(r_kind) pcp_gest,pcp_ges0,pcp_ges,obsges,termges,termgtl,obsgtl real(r_kind),dimension(max(1,nstep)):: pen real(r_kind) cg_pcp,wgross,wnotgross,pentl,pencur type(pcpNode), pointer :: pcpptr real(r_kind),pointer,dimension(:):: rt,st,rq,sq,ru,su,rv,sv,rcwm,scwm real(r_kind),pointer,dimension(:):: rql,rqi,sql,sqi ! Initialize penalty, b1, and b3 to zero out=zero_quad ! If no pcp data return if(.not. associated(pcphead))return ! Retrieve pointers ier=0; icw=0; iql=0; iqi=0 call gsi_bundlegetpointer(xval,'u', su,istatus);ier=istatus+ier call gsi_bundlegetpointer(xval,'v', sv,istatus);ier=istatus+ier call gsi_bundlegetpointer(xval,'tsen' ,st,istatus);ier=istatus+ier call gsi_bundlegetpointer(xval,'q', sq,istatus);ier=istatus+ier call gsi_bundlegetpointer(xval,'cw', scwm,istatus);icw=istatus+icw call gsi_bundlegetpointer(xval,'ql', sql,istatus);iql=istatus+iql call gsi_bundlegetpointer(xval,'qi', sqi,istatus);iqi=istatus+iqi if(ier/=0)return call gsi_bundlegetpointer(dval,'u', ru,istatus);ier=istatus+ier call gsi_bundlegetpointer(dval,'v', rv,istatus);ier=istatus+ier call gsi_bundlegetpointer(dval,'tsen' ,rt,istatus);ier=istatus+ier call gsi_bundlegetpointer(dval,'q', rq,istatus);ier=istatus+ier call gsi_bundlegetpointer(dval,'cw', rcwm,istatus);icw=istatus+icw call gsi_bundlegetpointer(dval,'ql', rql,istatus);iql=istatus+iql call gsi_bundlegetpointer(dval,'qi', rqi,istatus);iqi=istatus+iqi if(ier/=0)return lcld = (icw==0 .or. (iql+iqi)==0) ! Loop over number of observations. pcpptr => pcpNode_typecast(pcphead) do while(associated(pcpptr)) if(pcpptr%luse)then pcp_ges0 = pcpptr%ges if(nstep > 0)then j1=pcpptr%ij(1) j2=pcpptr%ij(2) j3=pcpptr%ij(3) j4=pcpptr%ij(4) w1=pcpptr%wij(1) w2=pcpptr%wij(2) w3=pcpptr%wij(3) w4=pcpptr%wij(4) pcp_gest = zero ! Compute updates to simulated precipitation do n=1,nsig dt =w1* st(j1)+w2* st(j2)+ w3* st(j3)+w4* st(j4) dq =w1* sq(j1)+w2* sq(j2)+ w3* sq(j3)+w4* sq(j4) du =w1* su(j1)+w2* su(j2)+ w3* su(j3)+w4* su(j4) dv =w1* sv(j1)+w2* sv(j2)+ w3* sv(j3)+w4* sv(j4) if (lcld) then if (icw==0) then dcwm=w1* scwm(j1)+w2* scwm(j2)+ w3* scwm(j3)+w4* scwm(j4) else dcwm=w1* (sql(j1)+sqi(j1))+ & w2* (sql(j2)+sqi(j2))+ & w3* (sql(j3)+sqi(j3))+ & w4* (sql(j4)+sqi(j4)) end if else dcwm=zero endif dt0 =w1* rt(j1)+w2* rt(j2)+ w3* rt(j3)+w4* rt(j4) dq0 =w1* rq(j1)+w2* rq(j2)+ w3* rq(j3)+w4* rq(j4) du0 =w1* ru(j1)+w2* ru(j2)+ w3* ru(j3)+w4* ru(j4) dv0 =w1* rv(j1)+w2* rv(j2)+ w3* rv(j3)+w4* rv(j4) if (lcld) then if (icw==0) then dcwm0=w1* rcwm(j1)+w2* rcwm(j2)+ w3* rcwm(j3)+w4* rcwm(j4) else dcwm0=w1* (rql(j1)+rqi(j1))+ & w2* (rql(j2)+rqi(j2))+ & w3* (rql(j3)+rqi(j3))+ & w4* (rql(j4)+rqi(j4)) end if else dcwm0=zero endif nt=n; nq=nt+nsig; nu=nq+nsig; nv=nu+nsig; ncwm=nv+nsig pcp_ges0 = pcp_ges0 + pcpptr%dpcp_dvar(nt) *dt + & pcpptr%dpcp_dvar(nq) *dq + & pcpptr%dpcp_dvar(nu) *du + & pcpptr%dpcp_dvar(nv) *dv + & pcpptr%dpcp_dvar(ncwm)*dcwm pcp_gest = pcp_gest + pcpptr%dpcp_dvar(nt) *dt0+ & pcpptr%dpcp_dvar(nq) *dq0+ & pcpptr%dpcp_dvar(nu) *du0+ & pcpptr%dpcp_dvar(nv) *dv0+ & pcpptr%dpcp_dvar(ncwm)*dcwm0 j1=j1+latlon11 j2=j2+latlon11 j3=j3+latlon11 j4=j4+latlon11 end do if (ltlint) then pcp_ges = pcp_ges0 ! Logrithmic formulation. Ensure pcp_ges > zero pcp_ges = max(pcp_ges,zero) termges = log(one+pcp_ges) obsges= pcpptr%obs - termges pencur = pcpptr%err2*obsges*obsges do kk=1,nstep if (pcp_ges>tinym1_obs) then termgtl = pcp_gest/(one+pcp_ges) else termgtl = zero endif obsgtl= - termgtl pentl = two*pcpptr%err2*obsges*obsgtl pen(kk) = pencur+sges(kk)*pentl enddo else do kk=1,nstep pcp_ges = pcp_ges0 + sges(kk)*pcp_gest ! Logrithmic formulation. Ensure pcp_ges > zero pcp_ges = max(pcp_ges,zero) termges = log(one+pcp_ges) obsges= pcpptr%obs - termges pen(kk) = pcpptr%err2*obsges*obsges enddo end if else pen(1)=pcpptr%err2*pcp_ges0*pcp_ges0 end if kx=pcpptr%icxp ! Modify penalty term if nonlinear QC if (nlnqc_iter .and. pg_pcp(kx) > tiny_r_kind .and. & b_pcp(kx) > tiny_r_kind) then cg_pcp=cg_term/b_pcp(kx) wnotgross= one-pg_pcp(kx)*varqc_iter wgross = varqc_iter*pg_pcp(kx)*cg_pcp/wnotgross do kk=1,max(1,nstep) pen(kk)= -two*log((exp(-half*pen(kk)) + wgross)/(one+wgross)) end do endif ! Accumulate stepsize terms out(1) = out(1)+ pen(1) * pcpptr%raterr2 do kk=2,nstep out(kk) = out(kk)+ (pen(kk)-pen(1)) * pcpptr%raterr2 end do end if ! pcpptr => pcpNode_nextcast(pcpptr) end do return end subroutine stppcp end module stppcpmod