module stpradmod !$$$ module documentation block ! . . . . ! module: stpradmod module for stprad and its tangent linear stprad_tl ! prgmmr: ! ! abstract: module for stprad and its tangent linear stprad_tl ! ! program history log: ! 2005-05-20 Yanqiu zhu - wrap stprad and its tangent linear stprad_tl into one module ! 2005-11-16 Derber - remove interfaces ! 2008-12-02 Todling - remove stprad_tl ! 2009-08-12 lueken - update documentation ! 2011-05-17 todling - add internal routine set_ ! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting ! ! subroutines included: ! sub stprad ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use kinds, only: i_kind implicit none PRIVATE PUBLIC stprad contains subroutine stprad(radhead,dval,xval,rpred,spred,out,sges,nstep) !$$$ subprogram documentation block ! . . . . ! subprogram: stprad compute contribution to penalty and stepsize ! from rad, using nonlinear qc. ! prgmmr: parrish org: np22 date: 1990-10-11 ! ! abstract: compute contribution to penalty and stepsize from radiances. ! ! program history log: ! 1990-10-11 parrish ! 1992-07-21 ! 1995-07-17 derber ! 1997-03-10 wu ! 1998-02-02 weiyu yang ! 1999-08-24 derber, j., treadon, r., yang, w., first frozen mpp version ! 2004-07-30 treadon - add only to module use, add intent in/out ! 2004-10-07 parrish - add nonlinear qc option ! 2005-01-20 okamoto - add wind components ! 2005-04-11 treadon - merge stprad and stprad_qc into single routine ! 2005-09-28 derber - modify var qc and change location and weight arrays ! 2007-03-19 tremolet - binning of observations ! 2007-07-28 derber - modify to use new inner loop obs data structure ! - unify NL qc ! 2007-06-04 derber - use quad precision to get reproducability over number of processors ! 2008-04-09 safford - rm unused vars and uses ! 2008-12-03 todling - changed handling of ptr%time ! 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 sst,oz cases; add pointer_state ! 2010-05-13 todling - update to use gsi_bundle ! - on-the-spot handling of non-essential vars ! 2010-07-10 todling - remove omp directives (per merge w/ r8741; Derber?) ! 2011-05-04 todling - merge in Min-Jeong Kim's cloud clear assimilation (connect to Metguess) ! 2011-05-16 todling - generalize entries in radiance jacobian ! 2011-05-17 augline/todling - add hydrometeors ! 2016-07-19 kbathmann- adjustment to bias correction when using correlated obs ! ! input argument list: ! radhead ! rt - search direction for temperature ! rq - search direction for moisture ! roz - search direction for ozone ! ru - search direction for zonal wind ! rv - search direction for meridional wind ! rst - search direction for skin temperature ! st - input temperature correction field ! sq - input q correction field ! soz - input ozone correction field ! su - input u correction field ! sv - input v correction field ! sst - input skin temp. vector ! rpred - search direction for predictors ! spred - input predictor values ! sges - step size estimates(nstep) ! nstep - number of stepsizes (==0 means use outer iteration value) ! ! output argument list: ! out(1:nstep) - penalty for radiance data sges(1:nstep) ! ! attributes: ! language: f90 ! machine: ibm RS/6000 SP ! !$$$ use kinds, only: r_kind,i_kind,r_quad use radinfo, only: npred,jpch_rad,b_rad,pg_rad use radinfo, only: nsigradjac use qcmod, only: nlnqc_iter,varqc_iter use constants, only: zero,half,one,two,tiny_r_kind,cg_term,r3600,zero_quad,one_quad use gridmod, only: nsig,latlon11 use gsi_bundlemod, only: gsi_bundle use gsi_bundlemod, only: gsi_bundlegetpointer use gsi_metguess_mod, only: gsi_metguess_get use mpeu_util, only: getindex use intradmod, only: luseu,lusev,luset,luseq,lusecw,luseoz,luseqg,luseqh,luseqi,luseql, & luseqr,luseqs use intradmod, only: itv,iqv,ioz,icw,ius,ivs,isst,iqg,iqh,iqi,iql,iqr,iqs,lgoback use m_obsNode, only: obsNode use m_radNode, only: radNode use m_radNode, only: radNode_typecast use m_radNode, only: radNode_nextcast implicit none ! Declare passed variables class(obsNode), pointer ,intent(in ) :: radhead integer(i_kind) ,intent(in ) :: nstep real(r_quad),dimension(max(1,nstep)) ,intent(inout) :: out real(r_kind),dimension(npred,jpch_rad) ,intent(in ) :: rpred,spred real(r_kind),dimension(max(1,nstep)) ,intent(in ) :: sges type(gsi_bundle),intent(in) :: dval type(gsi_bundle),intent(in) :: xval ! Declare local variables integer(i_kind) istatus integer(i_kind) nn,n,ic,k,nx,j1,j2,j3,j4,kk, mm, ic1 real(r_kind) val2,val,w1,w2,w3,w4 real(r_kind),dimension(nsigradjac):: tdir,rdir real(r_kind) cg_rad,wgross,wnotgross integer(i_kind),dimension(nsig) :: j1n,j2n,j3n,j4n real(r_kind),dimension(max(1,nstep)) :: term,rad type(radNode), pointer :: radptr real(r_kind),pointer,dimension(:) :: rt,rq,rcw,roz,ru,rv,rqg,rqh,rqi,rql,rqr,rqs real(r_kind),pointer,dimension(:) :: st,sq,scw,soz,su,sv,sqg,sqh,sqi,sql,sqr,sqs real(r_kind),pointer,dimension(:) :: rst,sst out=zero_quad ! If no rad data return if(.not. associated(radhead))return if(lgoback)return ! Retrieve pointers call gsi_bundlegetpointer(xval,'u', su, istatus) call gsi_bundlegetpointer(xval,'v', sv, istatus) call gsi_bundlegetpointer(xval,'tv' ,st, istatus) call gsi_bundlegetpointer(xval,'q', sq, istatus) call gsi_bundlegetpointer(xval,'cw' ,scw,istatus) call gsi_bundlegetpointer(xval,'oz' ,soz,istatus) call gsi_bundlegetpointer(xval,'sst',sst,istatus) call gsi_bundlegetpointer(xval,'qg' ,sqg,istatus) call gsi_bundlegetpointer(xval,'qh' ,sqh,istatus) call gsi_bundlegetpointer(xval,'qi' ,sqi,istatus) call gsi_bundlegetpointer(xval,'ql' ,sql,istatus) call gsi_bundlegetpointer(xval,'qr' ,sqr,istatus) call gsi_bundlegetpointer(xval,'qs' ,sqs,istatus) call gsi_bundlegetpointer(dval,'u', ru, istatus) call gsi_bundlegetpointer(dval,'v', rv, istatus) call gsi_bundlegetpointer(dval,'tv' ,rt, istatus) call gsi_bundlegetpointer(dval,'q', rq, istatus) call gsi_bundlegetpointer(dval,'cw' ,rcw,istatus) call gsi_bundlegetpointer(dval,'oz' ,roz,istatus) call gsi_bundlegetpointer(dval,'sst',rst,istatus) call gsi_bundlegetpointer(dval,'qg' ,rqg,istatus) call gsi_bundlegetpointer(dval,'qh' ,rqh,istatus) call gsi_bundlegetpointer(dval,'qi' ,rqi,istatus) call gsi_bundlegetpointer(dval,'ql' ,rql,istatus) call gsi_bundlegetpointer(dval,'qr' ,rqr,istatus) call gsi_bundlegetpointer(dval,'qs' ,rqs,istatus) tdir=zero rdir=zero radptr=> radNode_typecast(radhead) do while(associated(radptr)) if(radptr%luse)then if(nstep > 0)then j1=radptr%ij(1) j2=radptr%ij(2) j3=radptr%ij(3) j4=radptr%ij(4) w1=radptr%wij(1) w2=radptr%wij(2) w3=radptr%wij(3) w4=radptr%wij(4) if(luseu)then tdir(ius+1)=w1* su(j1) + w2* su(j2) + w3* su(j3) + w4* su(j4) rdir(ius+1)=w1* ru(j1) + w2* ru(j2) + w3* ru(j3) + w4* ru(j4) endif if(lusev)then tdir(ivs+1)=w1* sv(j1) + w2* sv(j2) + w3* sv(j3) + w4* sv(j4) rdir(ivs+1)=w1* rv(j1) + w2* rv(j2) + w3* rv(j3) + w4* rv(j4) endif if (isst>=0) then tdir(isst+1)=w1*sst(j1) + w2*sst(j2) + w3*sst(j3) + w4*sst(j4) rdir(isst+1)=w1*rst(j1) + w2*rst(j2) + w3*rst(j3) + w4*rst(j4) end if j1n(1) = j1 j2n(1) = j2 j3n(1) = j3 j4n(1) = j4 do n=2,nsig j1n(n) = j1n(n-1)+latlon11 j2n(n) = j2n(n-1)+latlon11 j3n(n) = j3n(n-1)+latlon11 j4n(n) = j4n(n-1)+latlon11 enddo do n=1,nsig j1 = j1n(n) j2 = j2n(n) j3 = j3n(n) j4 = j4n(n) ! Input state vector ! Input search direction vector if(luset)then tdir(itv+n)=w1* st(j1) +w2* st(j2) + w3* st(j3) +w4* st(j4) rdir(itv+n)=w1* rt(j1) +w2* rt(j2) + w3* rt(j3) +w4* rt(j4) endif if(luseq)then tdir(iqv+n)=w1* sq(j1) +w2* sq(j2) + w3* sq(j3) +w4* sq(j4) rdir(iqv+n)=w1* rq(j1) +w2* rq(j2) + w3* rq(j3) +w4* rq(j4) endif if (luseoz) then tdir(ioz+n)=w1*soz(j1)+w2*soz(j2)+ w3*soz(j3)+w4*soz(j4) rdir(ioz+n)=w1*roz(j1)+w2*roz(j2)+ w3*roz(j3)+w4*roz(j4) end if if (lusecw) then tdir(icw+n)=w1*scw(j1)+w2*scw(j2)+ w3*scw(j3)+w4*scw(j4) rdir(icw+n)=w1*rcw(j1)+w2*rcw(j2)+ w3*rcw(j3)+w4*rcw(j4) end if if (luseqg) then tdir(iqg+n)=w1*sqg(j1)+w2*sqg(j2)+ w3*sqg(j3)+w4*sqg(j4) rdir(iqg+n)=w1*rqg(j1)+w2*rqg(j2)+ w3*rqg(j3)+w4*rqg(j4) end if if (luseqh) then tdir(iqh+n)=w1*sqh(j1)+w2*sqh(j2)+ w3*sqh(j3)+w4*sqh(j4) rdir(iqh+n)=w1*rqh(j1)+w2*rqh(j2)+ w3*rqh(j3)+w4*rqh(j4) end if if (luseqi) then tdir(iqi+n)=w1*sqi(j1)+w2*sqi(j2)+ w3*sqi(j3)+w4*sqi(j4) rdir(iqi+n)=w1*rqi(j1)+w2*rqi(j2)+ w3*rqi(j3)+w4*rqi(j4) end if if (luseql) then tdir(iql+n)=w1*sql(j1)+w2*sql(j2)+ w3*sql(j3)+w4*sql(j4) rdir(iql+n)=w1*rql(j1)+w2*rql(j2)+ w3*rql(j3)+w4*rql(j4) end if if (luseqr) then tdir(iqr+n)=w1*sqr(j1)+w2*sqr(j2)+ w3*sqr(j3)+w4*sqr(j4) rdir(iqr+n)=w1*rqr(j1)+w2*rqr(j2)+ w3*rqr(j3)+w4*rqr(j4) end if if (luseqs) then tdir(iqs+n)=w1*sqs(j1)+w2*sqs(j2)+ w3*sqs(j3)+w4*sqs(j4) rdir(iqs+n)=w1*rqs(j1)+w2*rqs(j2)+ w3*rqs(j3)+w4*rqs(j4) end if end do end if do nn=1,radptr%nchan val2=-radptr%res(nn) if(nstep > 0)then val = zero ! contribution from bias corection ic=radptr%icx(nn) if (radptr%use_corr_obs) then do mm=1,radptr%nchan do nx=1,npred ic1=radptr%icx(mm) val2=val2+spred(nx,ic1)*radptr%rsqrtinv(mm,nn)*radptr%pred(nx,mm) val=val+rpred(nx,ic1)*radptr%rsqrtinv(mm,nn)*radptr%pred(nx,mm) end do end do else do nx=1,npred val2=val2+spred(nx,ic)*radptr%pred(nx,nn) val =val +rpred(nx,ic)*radptr%pred(nx,nn) end do end if ! end do ! contribution from atmosphere do k=1,nsigradjac val2=val2+tdir(k)*radptr%dtb_dvar(k,nn) val =val +rdir(k)*radptr%dtb_dvar(k,nn) end do ! calculate radiances for each stepsize do kk=1,nstep rad(kk)=val2+sges(kk)*val end do else rad(kk)= val2 end if ! calculate contribution to J do kk=1,max(1,nstep) term(kk) = radptr%err2(nn)*rad(kk)*rad(kk) end do ! Modify penalty term if nonlinear QC if(nlnqc_iter .and. pg_rad(ic) > tiny_r_kind .and. & b_rad(ic) > tiny_r_kind)then cg_rad=cg_term/b_rad(ic) wnotgross= one-pg_rad(ic)*varqc_iter wgross = varqc_iter*pg_rad(ic)*cg_rad/wnotgross do kk=1,max(1,nstep) term(kk) = -two*log((exp(-half*term(kk) ) + wgross)/(one+wgross)) end do endif out(1) = out(1) + term(1)*radptr%raterr2(nn) do kk=2,nstep out(kk) = out(kk) + (term(kk)-term(1))*radptr%raterr2(nn) end do end do end if radptr => radNode_nextcast(radptr) end do return end subroutine stprad end module stpradmod