module intozmod
!$$$ module documentation block
! . . . .
! module: intozmod module for intoz and its tangent linear intoz_tl
! prgmmr:
!
! abstract: module for intoz and its tangent linear intoz_tl
!
! program history log:
! 2005-05-13 Yanqiu zhu - wrap intoz and its tangent linear intoz_tl into one module
! 2005-11-16 Derber - remove interfaces
! 2008-11-26 Todling - remove intoz_tl; add interface back
! 2009-08-13 lueken - update documentation
! 2012-09-14 Syed RH Rizvi, NCAR/NESL/MMM/DAS - implemented obs adjoint test
! 2016-05-18 guo - replaced ob_type with polymorphic obsNode through type casting
! 2016-08-29 J. Guo - added individual interfaces, intozlay() and intozlev()
! 2018-07-13 J. Guo - splitted original module intozmod into this intozmod with
! subroutine intozlay() only, and module into3lmod below.
!
! subroutines included:
! sub intozlay_
!
! variable definitions:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use m_obsNode, only: obsNode
implicit none
PRIVATE
public:: intozlay
interface intozlay; module procedure intozlay_; end interface
contains
subroutine intozlay_(ozhead,rval,sval)
!$$$ subprogram documentation block
! . . . .
! subprogram: intoz apply nonlin qc obs operator for ozone
! prgmmr: derber org: np23 date: 1995-07-11
!
! abstract: This routine applies the observation operator (forward
! model) and adjoint of this operator for ozone observations
! with the addition of nonlinear qc.
!
! program history log:
! 1995-07-11 derber
! 1999-03-01 wu - port cray90 code to ibm-sp (mpi version)
! 2004-06-16 treadon - update documentation
! 2004-08-02 treadon - add only to module use, add intent in/out
! 2004-10-08 parrish - add nonlinear qc
! 2005-03-01 parrish - nonlinear qc change to account for inflated obs error
! 2005-04-11 treadon - merge intoz and intoz_qc into single routine
! 2005-06-14 wu - add OMI total ozone
! 2005-09-28 derber - consolidate location and weight arrays
! 2006-07-28 derber - modify to use new inner loop obs data structure
! - unify NL qc
! 2007-02-16 sienkiewicz - add call to routine for level ozone contrib.
! 2007-03-19 tremolet - binning of observations
! 2007-05-30 h.liu - move interpolation weights w1-w4 inside k loop
! 2007-06-05 tremolet - use observation diagnostics structure
! 2007-07-09 tremolet - observation sensitivity
! 2008-01-04 tremolet - Don't apply H^T if l_do_adjoint is false
! 2008-??-?? ?????? - remove nonlinear qc gradient; folded OMI within layer O3
! 2008-11-28 todling - turn FOTO optional; changed ptr%time handle
! 2009-01-18 todling - treat val in quad precision (revisit later)
! 2010-05-13 todling - update to use gsi_bundle; update interface
! 2012-09-08 wargan - add OMI with efficiency factors
! 2012-09-14 Syed RH Rizvi, NCAR/NESL/MMM/DAS - introduced ladtest_obs
! 2014-12-03 derber - modify so that use of obsdiags can be turned off
!
! input argument list:
! ozhead - layer ozone obs type pointer to obs structure
! soz - ozone increment in grid space
! roz
!
! output argument list:
! roz - ozone results from observation operator
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
!--------
use kinds, only: r_kind,i_kind,r_quad
use obsmod, only: lsaveobsens,l_do_adjoint,nloz_omi,luse_obsdiag
use gridmod, only: lat2,lon2,nsig
use jfunc, only: jiter
use constants, only: one,zero,r3600,zero_quad
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_4dvar, only: ladtest_obs
use m_ozNode , only: ozNode
use m_ozNode , only: ozNode_typecast
use m_ozNode , only: ozNode_nextcast
use m_obsdiagNode, only: obsdiagNode_get
use m_obsdiagNode, only: obsdiagNode_set
implicit none
! Declare passed variables
class(obsNode) ,pointer, intent(in ) :: ozhead
type(gsi_bundle), intent(in ) :: sval
type(gsi_bundle), intent(inout) :: rval
! Declare local variables
integer(i_kind) i,j,ij,ier,istatus
integer(i_kind) k,j1,j2,j3,j4,kk,iz1,iz2,kl
real(r_kind) dz1,pob,delz
real(r_quad) val1,valx
!-- real(r_kind) valx_
real(r_kind) w1,w2,w3,w4
real(r_kind),pointer,dimension(:,:,:) :: sozp
real(r_kind),pointer,dimension(:,:,:) :: rozp
real(r_kind),allocatable,dimension(:,:) :: soz
real(r_kind),allocatable,dimension(:,:) :: roz
type(ozNode), pointer :: ozptr
real(r_kind),dimension(nloz_omi):: val_lay
! If no data, return
if(.not. associated(ozhead))return
! Retrieve pointers
! Simply return if any pointer not found
ier=0
call gsi_bundlegetpointer(sval,'oz',sozp,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'oz',rozp,istatus);ier=istatus+ier
if(ier/=0)return
! Can't do rank-2 pointer into rank-2, therefore, allocate work space
allocate(soz(lat2*lon2,nsig),roz(lat2*lon2,nsig))
do k=1,nsig
ij=0
do j=1,lon2
do i=1,lat2
ij=ij+1
soz(ij,k) = sozp(i,j,k)
roz(ij,k) = rozp(i,j,k)
enddo
enddo
enddo
!
! SBUV OZONE: LAYER O3 and TOTAL O3
!
! Loop over ozone observations.
ozptr => ozNode_typecast(ozhead)
do while (associated(ozptr))
! Set location
j1=ozptr%ij(1)
j2=ozptr%ij(2)
j3=ozptr%ij(3)
j4=ozptr%ij(4)
! Accumulate contribution from layer observations
dz1=nsig+1
if ( ozptr%nloz >= 1 ) then
do k=1,ozptr%nloz
val1= zero_quad
pob = ozptr%prs(k)
iz1=dz1
if (iz1 > nsig) iz1=nsig
iz2=pob
do kk=iz1,iz2,-1
delz=one
if (kk==iz1) delz=dz1-iz1
if (kk==iz2) delz=delz-pob+iz2
w1=ozptr%wij(1,kk)
w2=ozptr%wij(2,kk)
w3=ozptr%wij(3,kk)
w4=ozptr%wij(4,kk)
val1=val1 + ( &
w1* soz(j1,kk)+ &
w2* soz(j2,kk)+ &
w3* soz(j3,kk)+ &
w4* soz(j4,kk))*delz
enddo
if(luse_obsdiag)then
if (lsaveobsens) then
valx=val1*ozptr%err2(k)*ozptr%raterr2(k)
!-- ozptr%diags(k)%ptr%obssen(jiter)=valx
call obsdiagNode_set(ozptr%diags(k)%ptr,jiter=jiter,obssen=real(valx,r_kind))
else
!-- if (ozptr%luse) ozptr%diags(k)%ptr%tldepart(jiter)=val1
if (ozptr%luse) call obsdiagNode_set(ozptr%diags(k)%ptr,jiter=jiter,tldepart=real(val1,r_kind))
endif
endif
if (l_do_adjoint) then
if (.not. lsaveobsens) then
if(ladtest_obs) then
valx=val1
else
val1=val1-ozptr%res(k)
valx = val1*ozptr%err2(k)
valx = valx*ozptr%raterr2(k)
endif
endif
do kk=iz1,iz2,-1
delz=one
if(kk==iz1)delz=dz1-iz1
if(kk==iz2)delz=delz-pob+iz2
w1=ozptr%wij(1,kk)
w2=ozptr%wij(2,kk)
w3=ozptr%wij(3,kk)
w4=ozptr%wij(4,kk)
roz(j1,kk) = roz(j1,kk) + valx*w1*delz
roz(j2,kk) = roz(j2,kk) + valx*w2*delz
roz(j3,kk) = roz(j3,kk) + valx*w3*delz
roz(j4,kk) = roz(j4,kk) + valx*w4*delz
enddo
dz1=pob
endif
end do
end if ! (ozptr%nloz >= 1)
! Add contribution from total column observation
k=ozptr%nloz+1
if (ozptr%apriori(1) .lt. zero) then ! non-OMIEFF ozone
val1= zero
do kk=nsig,1,-1
w1=ozptr%wij(1,kk)
w2=ozptr%wij(2,kk)
w3=ozptr%wij(3,kk)
w4=ozptr%wij(4,kk)
val1=val1 + &
w1* soz(j1,kk)+ &
w2* soz(j2,kk)+ &
w3* soz(j3,kk)+ &
w4* soz(j4,kk)
enddo
else ! OMI ozone with efficiency factor
! Integrate ozone within each layer
dz1=nsig+1
do kl=1,nloz_omi
val_lay(kl)= zero_quad
pob = ozptr%prs(kl)
iz1=dz1
if (iz1 > nsig) iz1=nsig
iz2=pob
do kk=iz1,iz2,-1
delz=one
if (kk==iz1) delz=dz1-iz1
if (kk==iz2) delz=delz-pob+iz2
w1=ozptr%wij(1,kk)
w2=ozptr%wij(2,kk)
w3=ozptr%wij(3,kk)
w4=ozptr%wij(4,kk)
val_lay(kl)=val_lay(kl) + ( &
w1* soz(j1,kk)+ &
w2* soz(j2,kk)+ &
w3* soz(j3,kk)+ &
w4* soz(j4,kk))*delz
enddo
dz1=pob
enddo
! Apply the efficiency factor
val1=zero_quad
do j=1,nloz_omi
val1=val1+ozptr%efficiency(j)*val_lay(j)
enddo
endif ! OMI ozone with efficiency factor
if(luse_obsdiag)then
if (lsaveobsens) then
valx=val1*ozptr%err2(k)*ozptr%raterr2(k)
!-- ozptr%diags(k)%ptr%obssen(jiter)=valx
call obsdiagNode_set(ozptr%diags(k)%ptr,jiter=jiter,obssen=real(valx,r_kind))
else
!-- if (ozptr%luse) ozptr%diags(k)%ptr%tldepart(jiter)=val1
if (ozptr%luse) call obsdiagNode_set(ozptr%diags(k)%ptr,jiter=jiter,tldepart=real(val1,r_kind))
endif
endif
if (l_do_adjoint) then
if (ozptr%apriori(1) .lt. zero) then ! non-OMI ozone
if (.not. lsaveobsens) then
if(ladtest_obs) then
valx=val1
else
val1=val1-ozptr%res(k)
valx = val1*ozptr%err2(k)
valx = valx*ozptr%raterr2(k)
endif
endif
do kk=nsig,1,-1
w1=ozptr%wij(1,kk)
w2=ozptr%wij(2,kk)
w3=ozptr%wij(3,kk)
w4=ozptr%wij(4,kk)
roz(j1,kk) = roz(j1,kk) + valx*w1
roz(j2,kk) = roz(j2,kk) + valx*w2
roz(j3,kk) = roz(j3,kk) + valx*w3
roz(j4,kk) = roz(j4,kk) + valx*w4
enddo
else ! OMI ozone with efficiency factor
if (lsaveobsens) then
! Precondition: luse_obsdiag .or. .not.lsaveobsens
! -------------------------------------------------
! lsaveobsens implies luse_obsdiag in a valid configuration. So
! there is no need to get valx back from %diags(k)%ptr%obssen(jiter).
! Also, this operation to get the value back from %obssen(jiter) will
! result an accuracy lost due to the kind difference between valx and
! %obssen(:).
!-- valx = ozptr%diags(k)%ptr%obssen(jiter) ! or ...
!-- call obsdiagNode_get(ozptr%diags(k)%ptr,jiter=jiter,obssen=valx_)
!-- valx=valx_ ! see vlax_ declaration on the top)
else
if(ladtest_obs) then
valx=val1
else
val1=val1-ozptr%res(k)
valx = val1*ozptr%err2(k)
valx = valx*ozptr%raterr2(k)
endif
endif
! Apply the transpose of the efficiency factor
do j=1,nloz_omi
val_lay(j)=ozptr%efficiency(j)*valx
enddo
! spread the info over GSI levels
do kl=1,nloz_omi
pob = ozptr%prs(kl)
iz1=dz1
if (iz1 > nsig) iz1=nsig
iz2=pob
do kk=iz1,iz2,-1
delz=one
if(kk==iz1)delz=dz1-iz1
if(kk==iz2)delz=delz-pob+iz2
w1=ozptr%wij(1,kk)
w2=ozptr%wij(2,kk)
w3=ozptr%wij(3,kk)
w4=ozptr%wij(4,kk)
roz(j1,kk) = roz(j1,kk) + val_lay(kl)*w1*delz
roz(j2,kk) = roz(j2,kk) + val_lay(kl)*w2*delz
roz(j3,kk) = roz(j3,kk) + val_lay(kl)*w3*delz
roz(j4,kk) = roz(j4,kk) + val_lay(kl)*w4*delz
enddo
dz1=pob
end do
endif ! OMI
endif ! do adjoint
ozptr => ozNode_nextcast(ozptr)
! End loop over observations
enddo
! Copy output and clean up
do k=1,nsig
ij=0
do j=1,lon2
do i=1,lat2
ij=ij+1
rozp(i,j,k) = roz(ij,k)
enddo
enddo
enddo
deallocate(soz,roz)
! End of routine
return
end subroutine intozlay_
end module intozmod
module into3lmod
!$$$ module documentation block
! . . . .
! module: intozlevmod module for intozlev()
! prgmmr:
!
! abstract: module for intozlev()
!
! program history log:
! 2018-07-13 J. Guo - splitted from original module intozmod into this into3lmod
! with subroutine intozlev(). See intozmod for more
! about earlier history logs.
!
! subroutines included:
! sub intozlev_
!
! variable definitions:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use m_obsNode, only: obsNode
implicit none
PRIVATE
public:: intozlev
interface intozlev; module procedure intozlev_; end interface
contains
subroutine intozlev_(o3lhead,rval,sval)
!$$$ subprogram documentation block
! . . . .
! subprogram: into3l apply nonlin qc obs operator for o3 level
! prgmmr: sienkiewicz org: GMAO date: 2006-09-14
!
! abstract: This routine applies the observation operator (forward
! model) and adjoint of this operator for ozone level
! observations with the addition of nonlinear qc.
!
! to do: add time derivatives correctly (Todling)
!
! program history log:
! 2006-09-14 sienkiewicz - add level ozone
! 2007-01-02 sienkiewicz - separate from intoz (again)
! 2007-02-16 sienkiewicz - changes for new inner loop obs data structure
! 2009-01-08 todling - remove nonlinear qc
! 2009-01-22 sienkiewicz - add time derivative
! 2010-05-13 todling - update to use gsi_bundle; update interface
! 2012-09-14 Syed RH Rizvi, NCAR/NESL/MMM/DAS - introduced ladtest_obs
! 2014-12-03 derber - modify so that use of obsdiags can be turned off
! 2015-12-01 todling - modify so that use of obsdiags can be turned off
!
! input argument list:
! o3lhead - level ozone obs type pointer to obs structure
! soz1d - ozone increment in grid space
! roz1d
!
! output argument list:
! roz1d - results from observation operator (0 for no data)
!
! attributes:
! language: f90
! machine: -
!
!$$$
!--------
use kinds, only: r_kind,i_kind
use obsmod, only: lsaveobsens, l_do_adjoint,luse_obsdiag
use constants, only: r3600
use jfunc, only: jiter
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_4dvar, only: ladtest_obs
use m_o3lNode, only: o3lNode
use m_o3lNode, only: o3lNode_typecast
use m_o3lNode, only: o3lNode_nextcast
use m_obsdiagNode, only: obsdiagNode_set
implicit none
! Declare passed variables
class(obsNode) ,pointer, intent(in ) :: o3lhead
type(gsi_bundle), intent(in ) :: sval
type(gsi_bundle), intent(inout) :: rval
! Declare local variables
integer(i_kind) ier,istatus
integer(i_kind) j1,j2,j3,j4,j5,j6,j7,j8
real(r_kind) val,grad
real(r_kind) w1,w2,w3,w4,w5,w6,w7,w8
real(r_kind),pointer,dimension(:) :: soz1d
real(r_kind),pointer,dimension(:) :: roz1d
type(o3lNode), pointer :: o3lptr
! If no data, return
if(.not. associated(o3lhead))return
! Retrieve pointers
! Simply return if any pointer not found
ier=0
call gsi_bundlegetpointer(sval,'oz',soz1d,istatus);ier=istatus+ier
call gsi_bundlegetpointer(rval,'oz',roz1d,istatus);ier=istatus+ier
if(ier/=0)return
! LEVEL-OZONE OBSERVATIONS
! Loop over ozone observations.
o3lptr => o3lNode_typecast(o3lhead)
do while (associated(o3lptr))
j1=o3lptr%ij(1)
j2=o3lptr%ij(2)
j3=o3lptr%ij(3)
j4=o3lptr%ij(4)
j5=o3lptr%ij(5)
j6=o3lptr%ij(6)
j7=o3lptr%ij(7)
j8=o3lptr%ij(8)
w1=o3lptr%wij(1)
w2=o3lptr%wij(2)
w3=o3lptr%wij(3)
w4=o3lptr%wij(4)
w5=o3lptr%wij(5)
w6=o3lptr%wij(6)
w7=o3lptr%wij(7)
w8=o3lptr%wij(8)
! Forward model
val=w1*soz1d(j1)+w2*soz1d(j2)+w3*soz1d(j3)+w4*soz1d(j4)+ &
w5*soz1d(j5)+w6*soz1d(j6)+w7*soz1d(j7)+w8*soz1d(j8)
if (luse_obsdiag ) then ! need to save either obssen=grad or tldepart=val
if (lsaveobsens) then
grad = val*o3lptr%raterr2*o3lptr%err2
!-- o3lptr%diags%obssen(jiter) = grad
call obsdiagNode_set(o3lptr%diags,jiter=jiter,obssen=grad)
else
!-- if (o3lptr%luse) o3lptr%diags%tldepart(jiter)=val
if (o3lptr%luse) call obsdiagNode_set(o3lptr%diags,jiter=jiter,tldepart=val)
endif
endif
if (l_do_adjoint) then
if (.not. lsaveobsens) then
if(ladtest_obs) then
grad = val
else
val=val-o3lptr%res
grad = val*o3lptr%raterr2*o3lptr%err2
endif
endif
! Adjoint
roz1d(j1)=roz1d(j1)+w1*grad
roz1d(j2)=roz1d(j2)+w2*grad
roz1d(j3)=roz1d(j3)+w3*grad
roz1d(j4)=roz1d(j4)+w4*grad
roz1d(j5)=roz1d(j5)+w5*grad
roz1d(j6)=roz1d(j6)+w6*grad
roz1d(j7)=roz1d(j7)+w7*grad
roz1d(j8)=roz1d(j8)+w8*grad
endif
o3lptr => o3lNode_nextcast(o3lptr)
end do
! End of routine
return
end subroutine intozlev_
end module into3lmod