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