module intgpsmod
!$$$ module documentation block
! . . . .
! module: intrefmod module for intref and its tangent linear intref_tl
! prgmmr:
!
! abstract: module for intref and its tangent linear intref_tl
!
! program history log:
! 2005-05-13 Yanqiu zhu - wrap intref and its tangent linear intref_tl into one module
! 2005-11-16 Derber - remove interfaces
! 2008-11-28 Todling - add interface back
! 2009-08-13 lueken - update documentation
! 2013-10-28 todling - rename p3d to prse
! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting
!
! subroutines included:
! sub intgps_
!
! variable definitions:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use m_obsNode, only: obsNode
use m_gpsNode, only: gpsNode
use m_gpsNode, only: gpsNode_typecast
use m_gpsNode, only: gpsNode_nextcast
implicit none
PRIVATE
PUBLIC intgps
interface intgps; module procedure &
intgps_
end interface
contains
subroutine intgps_(gpshead,rval,sval)
!$$$ subprogram documentation block
! . . . .
! subprogram: intref apply nonlinqc obs operator refractivity
! prgmmr: cucurull, l. org: JCSDA/NCEP date: 2004-04-29
!
! abstract: apply gps local refractivity operator and adjoint with
! addition of nonlinear qc.
!
! program history log:
! 2004-04-29 cucurull- original code
! 2004-06-21 treadon - update documentation
! 2004-08-02 treadon - add only to module use, add intent in/out
! 2004-10-08 parrish - add nonlinear qc option
! 2004-11-19 cucurull- add increments for surface pressure and temperature at levels
! below observation. Install non-linear forward operator.
! 2005-01-26 cucurull- Implement local GPS RO operator
! 2005-03-01 parrish - nonlinear qc change as above; correct bug in zeroing of tl_AD
! 2005-03-23 cucurull- correct bounds for obs below the second level; place
! bounds for k1 and k2
! 2005-04-11 treadon - merge intref and intref_qc into single routine
! 2005-08-02 derber - modify for variational qc parameters for each ob
! 2005-09-28 derber - consolidate location and weight arrays
! 2005-12-02 cucurull - fix bug for dimensions of sp and rp
! 2006-01-03 treadon - include r_kind type in w1,w2,...,w12 declaration
! 2006-07-28 derber - modify to use new inner loop obs data structure
! - unify NL qc
! 2006-09-06 cucurull - generalize code to hybrid vertical coordinate and modify to use
! surface pressure
! 2007-01-13 derber - clean up code and use coding standards
! 2007-03-28 derber - turn intref into generalized intgps
! 2007-07-26 cucurull - in/out 3d pressure to update code to generalized vertical coordinate
! 2008-06-02 safford - rm unused vars
! 2008-11-26 todling - add 4dvar and GSI adjoint capability (obs binning, obsens, etc)
! 2008-11-26 todling - turned FOTO optional; changed ptr%time handle
! 2010-05-13 todling - update to use gsi_bundle; update interface
! 2014-12-03 derber - modify so that use of obsdiags can be turned off
!
! input argument list:
! gpshead - obs type pointer to obs structure
! st - input temperature correction field
! sq - input q correction field
! sp - input (3D) p correction field
!
! output argument list:
! gpshead - obs type pointer to obs structure
! rt - output t vector after inclusion of gps local refractivity
! rq - output q vector after inclusion of gps local refractivity
! rp - output p vector after inclusion of gps local refractivity
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use obsmod, only: lsaveobsens,l_do_adjoint,luse_obsdiag
use qcmod, only: nlnqc_iter,varqc_iter
use gridmod, only: nsig
use constants, only: zero,one,half,tiny_r_kind,cg_term,r3600
use jfunc, only: jiter
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_4dvar, only: ladtest_obs
implicit none
! Declare passed variables
class(obsNode), pointer, intent(in ) :: gpshead
type(gsi_bundle), intent(in ) :: sval
type(gsi_bundle), intent(inout) :: rval
! Declare local variables
integer(i_kind) j,ier,istatus
integer(i_kind),dimension(nsig):: i1,i2,i3,i4
real(r_kind) :: w1,w2,w3,w4
real(r_kind) :: p_TL,p_AD,t_TL,t_AD,q_TL,q_AD
real(r_kind) :: val,pg_gps
real(r_kind) ::cg_gps,grad,p0,wnotgross,wgross
real(r_kind),pointer,dimension(:) :: st,sq
real(r_kind),pointer,dimension(:) :: rt,rq
real(r_kind),pointer,dimension(:) :: sp
real(r_kind),pointer,dimension(:) :: rp
type(gpsNode), pointer :: gpsptr
! If no gps obs return
if(.not. associated(gpshead))return
! Retrieve pointers
! Simply return if any pointer not found
ier=0
call gsi_bundlegetpointer(sval,'tv' ,st,istatus);ier=istatus+ier
call gsi_bundlegetpointer(sval,'q' ,sq,istatus);ier=istatus+ier
call gsi_bundlegetpointer(sval,'prse',sp,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'tv' ,rt,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'q' ,rq,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'prse',rp,istatus);ier=istatus+ier
if(ier/=0)return
!gpsptr => gpshead
gpsptr => gpsNode_typecast(gpshead)
do while (associated(gpsptr))
! Load location information into local variables
do j=1,nsig
i1(j)= gpsptr%ij(1,j)
i2(j)= gpsptr%ij(2,j)
i3(j)= gpsptr%ij(3,j)
i4(j)= gpsptr%ij(4,j)
enddo
w1=gpsptr%wij(1)
w2=gpsptr%wij(2)
w3=gpsptr%wij(3)
w4=gpsptr%wij(4)
val=zero
! local refractivity (linear operator)
do j=1,nsig
t_TL=w1* st(i1(j))+w2* st(i2(j))+w3* st(i3(j))+w4* st(i4(j))
q_TL=w1* sq(i1(j))+w2* sq(i2(j))+w3* sq(i3(j))+w4* sq(i4(j))
p_TL=w1* sp(i1(j))+w2* sp(i2(j))+w3* sp(i3(j))+w4* sp(i4(j))
val = val + p_TL*gpsptr%jac_p(j) + t_TL*gpsptr%jac_t(j)+q_TL*gpsptr%jac_q(j)
end do
if (luse_obsdiag)then
if (lsaveobsens) then
grad = val*gpsptr%raterr2*gpsptr%err2
gpsptr%diags%obssen(jiter) = grad
else
if (gpsptr%luse) gpsptr%diags%tldepart(jiter)=val
endif
endif
! Do adjoint
if (l_do_adjoint) then
if (.not. lsaveobsens) then
if( .not. ladtest_obs) val=val-gpsptr%res
! needed for gradient of nonlinear qc operator
if (nlnqc_iter .and. gpsptr%pg > tiny_r_kind .and. &
gpsptr%b > tiny_r_kind) then
pg_gps=gpsptr%pg*varqc_iter
cg_gps=cg_term/gpsptr%b
wnotgross= one-pg_gps
wgross = pg_gps*cg_gps/wnotgross
p0 = wgross/(wgross+exp(-half*gpsptr%err2*val**2))
val = val*(one-p0)
endif
if( ladtest_obs) then
grad = val
else
grad = val*gpsptr%raterr2*gpsptr%err2
end if
endif
! adjoint
do j=1,nsig
t_AD = grad*gpsptr%jac_t(j)
rt(i1(j))=rt(i1(j))+w1*t_AD
rt(i2(j))=rt(i2(j))+w2*t_AD
rt(i3(j))=rt(i3(j))+w3*t_AD
rt(i4(j))=rt(i4(j))+w4*t_AD
q_AD = grad*gpsptr%jac_q(j)
rq(i1(j))=rq(i1(j))+w1*q_AD
rq(i2(j))=rq(i2(j))+w2*q_AD
rq(i3(j))=rq(i3(j))+w3*q_AD
rq(i4(j))=rq(i4(j))+w4*q_AD
p_AD = grad*gpsptr%jac_p(j)
rp(i1(j))=rp(i1(j))+w1*p_AD
rp(i2(j))=rp(i2(j))+w2*p_AD
rp(i3(j))=rp(i3(j))+w3*p_AD
rp(i4(j))=rp(i4(j))+w4*p_AD
enddo
endif
!gpsptr => gpsptr%llpoint
gpsptr => gpsNode_nextcast(gpsptr)
end do
return
end subroutine intgps_
end module intgpsmod