module intlagmod
!$$$ module documentation block
! . . . .
! module: intlagmod module for intlag
! prgmmr: meunier
!
! abstract: module for intlag
!
! program history log:
! 2008-03-23 lmeunier - initial code
! 2009-08-13 lueken - update documentation
!
! subroutines included:
! sub intlag
!
! variable definitions:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
PRIVATE
PUBLIC intlag
contains
subroutine intlag(laghead,ru,rv,su,sv,obsbin)
!$$$ subprogram documentation block
! . . . .
! subprogram: intlag apply nonlin qc obs operator for conv. lag
! prgmmr: meunier org: np23 date: 2009-03-23
!
! abstract: apply observation operator and adjoint for conventional lag
! observations with nonlinear qc operator
!
! program history log:
! 2009-03-23 meunier
!
! input argument list:
! laghead - obs type pointer to appropriate obs structure
! su - zonal wind increment in grid space
! sv - meridian wind increment in grid space
! obsbin
!
! output argument list:
! su - result for zonal wind increment in grid space
! sv - result fo meridian wind increment in grid space
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: half,one,zero,izero,ione,tiny_r_kind,cg_term,rad2deg
use obsmod, only: lag_ob_type, lsaveobsens, l_do_adjoint
use qcmod, only: nlnqc_iter
use gridmod, only: latlon1n,iglobal
use jfunc, only: jiter
use lag_fields, only: lag_gather_stateuv,lag_ADscatter_stateuv
use lag_fields, only: lag_u_full,lag_v_full,lag_uv_alloc,lag_uv_fill
use lag_fields, only: lag_tl_vec,lag_ad_vec
use lag_fields, only: orig_lag_num,ntotal_orig_lag,lag_kcount
use lag_traj, only: lag_rk2iter_tl,lag_rk2iter_ad
! use lag_traj, only: lag_rk4iter_tl,lag_rk4iter_ad
use mpimod, only: mype
implicit none
! Declare passed variables
type(lag_ob_type),pointer ,intent(in ) :: laghead
real(r_kind),dimension(latlon1n),intent(in ) :: su,sv
real(r_kind),dimension(latlon1n),intent(inout) :: ru,rv
integer(i_kind) ,intent(in ) :: obsbin
! Print level
integer(i_kind),parameter::iv_debug = izero
! real(r_kind) penalty
real(r_kind) cg_lag,p0,wnotgross,wgross
real(r_kind) lon_tl,lat_tl,p_tl
real(r_kind) grad_lon,grad_lat,grad_p
type(lag_ob_type), pointer :: lagptr
integer(i_kind) i,j
real(r_kind),dimension(:,:),allocatable:: adu_tmp,adv_tmp
! If no balloons at all exit (save time)
if (ntotal_orig_lag==izero) return
! Gather the wind fields if not already done
if (.not.lag_uv_alloc) then
call lag_gather_stateuv(su,sv,obsbin)
else
if (.not.lag_uv_fill(obsbin)) then
call lag_gather_stateuv(su,sv,obsbin)
end if
end if
! Allocate AD fields
if (l_do_adjoint) then
allocate(adu_tmp(iglobal,lag_kcount))
allocate(adv_tmp(iglobal,lag_kcount))
adu_tmp=zero; adv_tmp=zero
end if
lagptr => laghead
do while (associated(lagptr))
! Forward model
lon_tl=lag_tl_vec(orig_lag_num(lagptr%intnum,3),obsbin,1)
lat_tl=lag_tl_vec(orig_lag_num(lagptr%intnum,3),obsbin,2)
p_tl =lag_tl_vec(orig_lag_num(lagptr%intnum,3),obsbin,3)
if (iv_debug>=2_i_kind) then
print *,'TL Orig:',lon_tl*rad2deg,lat_tl*rad2deg
print *,'SPECI',lagptr%speci
print *,'SPECR',lagptr%specr
end if
if (iv_debug>=ione) then
print *,'MAX INC U',maxval(abs(lag_u_full(:,:,obsbin)))
print *,'MAX INC V',maxval(abs(lag_v_full(:,:,obsbin)))
end if
call lag_rk2iter_tl(lagptr%speci,lagptr%specr,&
&lon_tl,lat_tl,p_tl,&
&lag_u_full(:,:,obsbin),lag_v_full(:,:,obsbin))
if (iv_debug>=ione) print *,'TL correction:',lon_tl*rad2deg,lat_tl*rad2deg
if (lsaveobsens) then
lagptr%diag_lon%obssen(jiter) = lon_tl*lagptr%raterr2*lagptr%err2_lon
lagptr%diag_lat%obssen(jiter) = lat_tl*lagptr%raterr2*lagptr%err2_lat
else
if (lagptr%luse) lagptr%diag_lon%tldepart(jiter)=lon_tl
if (lagptr%luse) lagptr%diag_lat%tldepart(jiter)=lat_tl
endif
if (l_do_adjoint) then
if (lsaveobsens) then
grad_lon = lagptr%diag_lon%obssen(jiter)
grad_lat = lagptr%diag_lat%obssen(jiter)
grad_p = zero
else
lon_tl=lon_tl-lagptr%res_lon
lat_tl=lat_tl-lagptr%res_lat
! gradient of nonlinear operator
if (nlnqc_iter .and. lagptr%pg > tiny_r_kind .and. &
lagptr%b > tiny_r_kind) then
cg_lag=cg_term/lagptr%b
wnotgross= one-lagptr%pg
wgross = lagptr%pg*cg_lag/wnotgross
p0 = wgross/(wgross+&
&exp(-half*(lagptr%err2_lon*lon_tl**2+lagptr%err2_lat*lat_tl**2)))
lon_tl = lon_tl*(one-p0)
lat_tl = lat_tl*(one-p0)
if (iv_debug>=ione) print *,'Do nlnqc_iter'
end if
grad_lon = lon_tl*lagptr%raterr2*lagptr%err2_lon
grad_lat = lat_tl*lagptr%raterr2*lagptr%err2_lat
grad_p = zero
endif
if (iv_debug>=2_i_kind) then
print *,'Residual a. cor:',lon_tl*rad2deg,lat_tl*rad2deg
print *,'R application :',grad_lon*rad2deg,grad_lat*rad2deg
end if
! Adjoint
call lag_rk2iter_ad(lagptr%speci,lagptr%specr,&
&grad_lon,grad_lat,grad_p,adu_tmp,adv_tmp)
lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,1)=&
&lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,1)+grad_lon
lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,2)=&
&lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,2)+grad_lat
lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,3)=&
&lag_ad_vec(orig_lag_num(lagptr%intnum,3),obsbin,3)+grad_p
if (iv_debug>=2) then
do i=1,iglobal
do j=1,lag_kcount
if (adu_tmp(i,j)/=zero) print *,'IntLag Mype',mype,'GradU',adu_tmp(i,j)
if (adv_tmp(i,j)/=zero) print *,'IntLag Mype',mype,'GradV',adv_tmp(i,j)
end do
end do
end if
endif
lagptr => lagptr%llpoint
end do
if (l_do_adjoint) then
lag_u_full(:,:,obsbin)=adu_tmp
lag_v_full(:,:,obsbin)=adv_tmp
! Scater back the UV adjoint values back to state vector
call lag_ADscatter_stateuv(ru,rv,obsbin)
lag_u_full(:,:,obsbin)=zero
lag_v_full(:,:,obsbin)=zero
deallocate(adu_tmp)
deallocate(adv_tmp)
end if
! Force the reload of the new increment for the next minimization step
lag_uv_fill(obsbin)=.false.
return
end subroutine intlag
end module intlagmod