#include "cppdefs.h"
MODULE zeta_balance_mod
#if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
!git $Id$
!svn $Id: zeta_balance.F 1180 2023-07-13 02:42:10Z arango $
!=================================================== Andrew M. Moore ===
! Copyright (c) 2002-2023 The ROMS/TOMS Group Hernan G. Arango !
! Licensed under a MIT/X style license !
! See License_ROMS.md !
!=======================================================================
! !
! This module contains routines to solve the elliptic equation for !
! free surface (SSH) that is used as part of the balance operator !
! during 4D-Var. !
! !
! The balanced (baroclinic) SSH 4DVar increment is approximated as !
! the sum of two terms: a baroclinic term that depends on density !
! and a barotropic term that depends on the depth-integrated !
! transport. !
! !
! div (h grad(zeta)) = - div (int{int{grad(rho)z'/rho0) dz'} dz}) !
! !
! References: !
! !
! Fukumori, I., R. Raghunath and L. Fu, 1998: Nature of global !
! large-scale sea level variability in relation to atmospheric !
! forcing: a modeling study, J. Geophys. Res., 103, 5493-5512. !
! !
! Weaver, A.T., C. Deltel, E. Machu, S. Ricci, and N. Daget, 2005: !
! A multivariate balance operator for variational data assimila- !
! tion, Q. J. R. Meteorol. Soc., 131, 3605-3625. !
! !
!=======================================================================
!
implicit none
PRIVATE
PUBLIC :: balance_ref
PUBLIC :: biconj
PUBLIC :: ad_biconj_tile
PUBLIC :: tl_biconj_tile
PUBLIC :: r2d_bc
PUBLIC :: ad_r2d_bc
PUBLIC :: u2d_bc
PUBLIC :: v2d_bc
CONTAINS
!
!***********************************************************************
SUBROUTINE balance_ref (ng, tile, Lbck)
!***********************************************************************
!
USE mod_param
USE mod_grid
# ifdef SOLVE3D
USE mod_coupling
USE mod_mixing
# endif
USE mod_fourdvar
USE mod_ocean
USE mod_stepping
# ifdef SOLVE3D
!
USE rho_eos_mod
USE set_depth_mod
# endif
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, Lbck
!
! Local variable declarations.
!
# include "tile.h"
# ifdef SOLVE3D
!
! Compute background state thickness, depth arrays, and density fields.
! Use zero value free-surface.
!
COUPLING(ng) % Zt_avg1 = 0.0_r8
CALL set_depth (ng, tile, iNLM)
nrhs(ng)=Lbck
CALL rho_eos (ng, tile, iNLM)
!
# endif
CALL balance_ref_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& Lbck, &
& GRID(ng) % pm, &
& GRID(ng) % pn, &
& GRID(ng) % pmon_r, &
& GRID(ng) % pnom_r, &
& GRID(ng) % pmon_u, &
& GRID(ng) % pnom_v, &
& GRID(ng) % h, &
# ifdef SOLVE3D
& GRID(ng) % Hz, &
& GRID(ng) % z_r, &
& GRID(ng) % z_w, &
# endif
# ifdef MASKING
& GRID(ng) % rmask, &
& GRID(ng) % umask, &
& GRID(ng) % vmask, &
# endif
# ifdef SOLVE3D
& MIXING(ng) % alpha, &
& MIXING(ng) % beta, &
& OCEAN(ng) % rho, &
# endif
& OCEAN(ng) % zeta, &
& FOURDVAR(ng) % rhs_r2d)
RETURN
END SUBROUTINE balance_ref
!
!***********************************************************************
SUBROUTINE balance_ref_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& Lbck, &
& pm, pn, pmon_r, pnom_r, &
& pmon_u, pnom_v, h, &
# ifdef SOLVE3D
& Hz, z_r, z_w, &
# endif
# ifdef MASKING
& rmask, umask, vmask, &
# endif
# ifdef SOLVE3D
& alpha, beta, rho, &
# endif
& zeta, &
& rhs_r2d)
!***********************************************************************
!
USE mod_param
USE mod_parallel
USE mod_scalars
USE mod_grid
!
USE exchange_2d_mod
USE exchange_3d_mod
# ifdef DISTRIBUTE
USE mp_exchange_mod, ONLY : mp_exchange2d
# endif
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
integer, intent(in) :: LBi, UBi, LBj, UBj
integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
integer, intent(in) :: Lbck
!
# ifdef ASSUMED_SHAPE
real(r8), intent(in) :: pm(LBi:,LBj:)
real(r8), intent(in) :: pn(LBi:,LBj:)
real(r8), intent(in) :: pmon_r(LBi:,LBj:)
real(r8), intent(in) :: pnom_r(LBi:,LBj:)
real(r8), intent(in) :: pmon_u(LBi:,LBj:)
real(r8), intent(in) :: pnom_v(LBi:,LBj:)
real(r8), intent(in) :: h(LBi:,LBj:)
# ifdef SOLVE3D
real(r8), intent(in) :: Hz(LBi:,LBj:,:)
real(r8), intent(in) :: z_r(LBi:,LBj:,:)
real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
# endif
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:,LBj:)
real(r8), intent(in) :: umask(LBi:,LBj:)
real(r8), intent(in) :: vmask(LBi:,LBj:)
# endif
# ifdef SOLVE3D
real(r8), intent(in) :: alpha(LBi:,LBj:)
real(r8), intent(in) :: beta(LBi:,LBj:)
real(r8), intent(in) :: rho(LBi:,LBj:,:)
# endif
real(r8), intent(inout) :: rhs_r2d(LBi:,LBj:)
real(r8), intent(in) :: zeta(LBi:,LBj:,:)
# else
real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pmon_r(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pnom_r(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pmon_u(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pnom_v(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
# ifdef SOLVE3D
real(r8), intent(in) :: Hz(LBi:UBi,LBj:UBj,N(ng))
real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
# endif
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
# endif
# ifdef SOLVE3D
real(r8), intent(in) :: alpha(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: beta(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
# endif
real(r8), intent(inout) :: rhs_r2d(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,:)
# endif
!
! Local variable declarations.
!
integer :: i, j, k
real(r8) :: fac, fac1, fac2, fac3, gamma
real(r8) :: cff, cff1, cff2, cff3, cff4
real(r8), dimension(IminS:ImaxS) :: phie
real(r8), dimension(IminS:ImaxS) :: phix
real(r8), dimension(IminS:ImaxS,1:N(ng)) :: phi
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gradPx
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: gradPy
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
! Compute the RHS side term of the elliptic equation for the
! biconjugate gradient solver.
!-----------------------------------------------------------------------
!
! NOTE: fac2=0 because the balanced component should consist of the
! baroclinic pressure gradient only.
!
fac1=0.5_r8*g/rho0
!! fac2=1000.0_r8*g/rho0
fac2=0.0_r8
fac3=0.25_r8*g/rho0
DO j=JminS,JmaxS
DO i=IminS,ImaxS
gradPx(i,j)=0.0_r8
gradPy(i,j)=0.0_r8
END DO
END DO
!
! Compute balanced, surface U-momentum from baroclinic and barotropic
! surface pressure gradient.
!
DO j=Jstr,Jend+1
DO i=Istr-1,Iend
cff1=z_w(i,j ,N(ng))-z_r(i,j ,N(ng))+ &
& z_w(i,j-1,N(ng))-z_r(i,j-1,N(ng))
phie(i)=fac1*(rho(i,j,N(ng))-rho(i,j-1,N(ng)))*cff1+ &
& fac2*(zeta(i,j,Lbck)-zeta(i,j-1,Lbck))
phi(i,N(ng))=phie(i)
END DO
!
! Compute balance, interior U-momentum from baroclinic pressure
! gradient (differentiate and then vertically integrate).
!
DO k=N(ng)-1,1,-1
DO i=Istr-1,Iend
cff1=1.0_r8/((z_r(i,j ,k+1)-z_r(i,j ,k))* &
& (z_r(i,j-1,k+1)-z_r(i,j-1,k)))
cff2=z_r(i,j ,k )-z_r(i,j-1,k )+ &
& z_r(i,j ,k+1)-z_r(i,j-1,k+1)
cff3=z_r(i,j ,k+1)-z_r(i,j ,k )- &
& z_r(i,j-1,k+1)+z_r(i,j-1,k )
gamma=0.125_r8*cff1*cff2*cff3
!
cff1=(1.0_r8+gamma)*(rho(i,j,k+1)-rho(i,j-1,k+1))+ &
& (1.0_r8-gamma)*(rho(i,j,k )-rho(i,j-1,k ))
cff2=rho(i,j,k+1)+rho(i,j-1,k+1)- &
& rho(i,j,k )-rho(i,j-1,k )
cff3=z_r(i,j,k+1)+z_r(i,j-1,k+1)- &
& z_r(i,j,k )-z_r(i,j-1,k )
cff4=(1.0_r8+gamma)*(z_r(i,j,k+1)-z_r(i,j-1,k+1))+ &
& (1.0_r8-gamma)*(z_r(i,j,k )-z_r(i,j-1,k ))
phie(i)=phie(i)+fac3*(cff1*cff3-cff2*cff4)
phi(i,k)=phie(i)
END DO
END DO
!
! Integrate from bottom to surface.
!
DO k=1,N(ng)
DO i=Istr-1,Iend
cff=0.5_r8*(Hz(i,j-1,k)+Hz(i,j,k))*phi(i,k)
# ifdef MASKING
cff=cff*vmask(i,j)
# endif
gradPy(i,j)=gradPy(i,j)+cff
END DO
END DO
END DO
!
! Compute balanced, surface V-momentum from baroclinic and barotropic
! surface pressure gradient.
!
DO j=Jstr-1,Jend
DO i=Istr,Iend+1
cff1=z_w(i ,j,N(ng))-z_r(i ,j,N(ng))+ &
& z_w(i-1,j,N(ng))-z_r(i-1,j,N(ng))
phix(i)=fac1*(rho(i,j,N(ng))-rho(i-1,j,N(ng)))*cff1+ &
& fac2*(zeta(i,j,Lbck)-zeta(i-1,j,Lbck))
phi(i,N(ng))=phix(i)
END DO
!
! Compute balance, interior V-momentum from baroclinic pressure
! gradient (differentiate and then vertically integrate).
!
DO k=N(ng)-1,1,-1
DO i=Istr,Iend+1
cff1=1.0_r8/((z_r(i ,j,k+1)-z_r(i ,j,k))* &
& (z_r(i-1,j,k+1)-z_r(i-1,j,k)))
cff2=z_r(i ,j,k )-z_r(i-1,j,k )+ &
& z_r(i ,j,k+1)-z_r(i-1,j,k+1)
cff3=z_r(i ,j,k+1)-z_r(i ,j,k )- &
& z_r(i-1,j,k+1)+z_r(i-1,j,k )
gamma=0.125_r8*cff1*cff2*cff3
!
cff1=(1.0_r8+gamma)*(rho(i,j,k+1)-rho(i-1,j,k+1))+ &
& (1.0_r8-gamma)*(rho(i,j,k )-rho(i-1,j,k ))
cff2=rho(i,j,k+1)+rho(i-1,j,k+1)- &
& rho(i,j,k )-rho(i-1,j,k )
cff3=z_r(i,j,k+1)+z_r(i-1,j,k+1)- &
& z_r(i,j,k )-z_r(i-1,j,k )
cff4=(1.0_r8+gamma)*(z_r(i,j,k+1)-z_r(i-1,j,k+1))+ &
& (1.0_r8-gamma)*(z_r(i,j,k )-z_r(i-1,j,k ))
phix(i)=phix(i)+fac3*(cff1*cff3-cff2*cff4)
phi(i,k)=phix(i)
END DO
END DO
!
! Integrate from bottom to surface.
!
DO k=1,N(ng)
DO i=Istr,Iend+1
cff=0.5_r8*(Hz(i-1,j,k)+Hz(i,j,k))*phi(i,k)
# ifdef MASKING
cff=cff*umask(i,j)
# endif
gradPx(i,j)=gradPx(i,j)+cff
END DO
END DO
END DO
!
! Impose zero gradient conditions at open boundaries.
!
CALL u2d_bc (ng, tile, &
& IminS, ImaxS, JminS, JmaxS, &
& gradPx)
CALL v2d_bc (ng, tile, &
& IminS, ImaxS, JminS, JmaxS, &
& gradPy)
!
! Compute the RHS term, -div(phi_bar), for the biconjugate gradient
! solver.
!
DO j=Jstr,Jend
DO i=Istr,Iend
rhs_r2d(i,j)=-pm(i,j)*pn(i,j)* &
& (pmon_u(i+1,j)*gradPx(i+1,j)- &
& pmon_u(i ,j)*gradPx(i ,j)+ &
& pnom_v(i,j+1)*gradPy(i,j+1)- &
& pnom_v(i,j )*gradPy(i,j ))
# ifdef MASKING
rhs_r2d(i,j)=rhs_r2d(i,j)*rmask(i,j)
# endif
END DO
END DO
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& rhs_r2d)
# ifdef DISTRIBUTE
CALL mp_exchange2d (ng, tile, iNLM, 1, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& rhs_r2d)
# endif
RETURN
END SUBROUTINE balance_ref_tile
!
!***********************************************************************
SUBROUTINE biconj (ng, tile, model, Lbck)
!***********************************************************************
!
USE mod_param
USE mod_grid
# ifdef SOLVE3D
USE mod_coupling
USE mod_mixing
# endif
USE mod_fourdvar
USE mod_ocean
USE mod_stepping
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, Lbck, model
!
! Local variable declarations.
!
# include "tile.h"
!
! Call the biconjugate gradient solver to compute "zeta_ref".
!
CALL biconj_tile (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& Lbck, &
& GRID(ng)% h, &
& GRID(ng) % pmon_u, &
& GRID(ng) % pnom_v, &
& GRID(ng) % pm, &
& GRID(ng) % pn, &
# ifdef MASKING
& GRID(ng) % umask, &
& GRID(ng) % vmask, &
& GRID(ng) % rmask, &
# endif
& FOURDVAR(ng) % bc_ak, &
& FOURDVAR(ng) % bc_bk, &
& FOURDVAR(ng) % zdf1, &
& FOURDVAR(ng) % zdf2, &
& FOURDVAR(ng) % zdf3, &
& FOURDVAR(ng) % pc_r2d, &
& FOURDVAR(ng) % r_r2d, &
& FOURDVAR(ng) % br_r2d, &
& FOURDVAR(ng) % p_r2d, &
& FOURDVAR(ng) % bp_r2d, &
& OCEAN(ng) % zeta, &
& FOURDVAR(ng) % zeta_ref, &
& FOURDVAR(ng) % rhs_r2d)
RETURN
END SUBROUTINE biconj
!
!***********************************************************************
SUBROUTINE biconj_tile (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& Lbck, &
& h, pmon_u, pnom_v, pm, pn, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& bc_ak, bc_bk, zdf1, zdf2, zdf3, &
& pc_r2d, r_r2d, br_r2d, p_r2d, bp_r2d, &
& zeta, r2d_ref, rhs_r2d)
!***********************************************************************
!
USE mod_param
USE mod_scalars
USE mod_grid
USE mod_iounits
USE mod_parallel
USE mod_scalars
!
# ifdef DISTRIBUTE
USE mp_exchange_mod, ONLY : ad_mp_exchange2d
USE mp_exchange_mod, ONLY : mp_exchange2d
# endif
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, model
integer, intent(in) :: LBi, UBi, LBj, UBj
integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
integer, intent(in) :: Lbck
!
# ifdef ASSUMED_SHAPE
# ifdef MASKING
real(r8), intent(in) :: umask(LBi:,LBj:)
real(r8), intent(in) :: vmask(LBi:,LBj:)
real(r8), intent(in) :: rmask(LBi:,LBj:)
# endif
real(r8), intent(in) :: h(LBi:,LBj:)
real(r8), intent(in) :: pmon_u(LBi:,LBj:)
real(r8), intent(in) :: pnom_v(LBi:,LBj:)
real(r8), intent(in) :: pm(LBi:,LBj:)
real(r8), intent(in) :: pn(LBi:,LBj:)
real(r8), intent(in) :: zeta(LBi:,LBj:,:)
real(r8), intent(inout) :: bc_ak(:)
real(r8), intent(inout) :: bc_bk(:)
real(r8), intent(inout) :: zdf1(:)
real(r8), intent(inout) :: zdf2(:)
real(r8), intent(inout) :: zdf3(:)
real(r8), intent(inout) :: pc_r2d(LBi:,LBj:)
real(r8), intent(inout) :: r_r2d(LBi:,LBj:,:)
real(r8), intent(inout) :: br_r2d(LBi:,LBj:,:)
real(r8), intent(inout) :: p_r2d(LBi:,LBj:,:)
real(r8), intent(inout) :: bp_r2d(LBi:,LBj:,:)
real(r8), intent(inout) :: rhs_r2d(LBi:,LBj:)
real(r8), intent(inout) :: r2d_ref(LBi:,LBj:)
# else
# ifdef MASKING
real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pmon_u(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pnom_v(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,:)
real(r8), intent(inout) :: bc_ak(Nbico(ng))
real(r8), intent(inout) :: bc_bk(Nbico(ng))
real(r8), intent(inout) :: zdf1(Nbico(ng))
real(r8), intent(inout) :: zdf2(Nbico(ng))
real(r8), intent(inout) :: zdf3(Nbico(ng))
real(r8), intent(inout) :: pc_r2d(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: r_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
real(r8), intent(inout) :: br_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
real(r8), intent(inout) :: p_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
real(r8), intent(inout) :: bp_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
real(r8), intent(inout) :: rhs_r2d(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: r2d_ref(LBi:UBi,LBj:UBj)
# endif
!
! Local variable declarations.
!
logical :: Ltrans
integer :: i, j, it
real(r8) :: error, zdf4, zdf5
real(r8), dimension(LBi:UBi,LBj:UBj) :: bv_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: v_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: z1_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: z2_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: z3_r2d
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
! Biconjugate gradient algorithm.
!-----------------------------------------------------------------------
!
! Initialize local arrays.
!
DO j=LBj,UBj
DO i=LBi,UBi
v_r2d(i,j)=0.0_r8
bv_r2d(i,j)=0.0_r8
z1_r2d(i,j)=0.0_r8
z2_r2d(i,j)=0.0_r8
z3_r2d(i,j)=0.0_r8
END DO
END DO
!
! Choose the starting value of "zeta_ref". Use background free-surface.
!
DO j=Jstr,Jend
DO i=Istr,Iend
r2d_ref(i,j)=zeta(i,j,Lbck)
# ifdef MASKING
r2d_ref(i,j)=r2d_ref(i,j)*rmask(i,j)
# endif
END DO
END DO
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& r2d_ref)
# ifdef DISTRIBUTE
CALL mp_exchange2d (ng, tile, model, 1, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& r2d_ref)
# endif
!
! Compute starting value of divergence operator:
! z1_r2d = div[h grad(r2d_ref)].
!
Ltrans=.FALSE.
CALL r2d_oper (ng, tile, &
& Ltrans, Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, r2d_ref, z1_r2d)
!
! Set the initial values for residual vectors "r" and "br". Then, use
! recurrence relationship to compute direction vectors "p" and "bp".
!
DO j=Jstr,Jend
DO i=Istr,Iend
r_r2d(i,j,1)=rhs_r2d(i,j)-z1_r2d(i,j)
br_r2d(i,j,1)=r_r2d(i,j,1)
END DO
END DO
DO j=Jstr,Jend
DO i=Istr,Iend
p_r2d(i,j,1)=r_r2d(i,j,1)/pc_r2d(i,j)
bp_r2d(i,j,1)=br_r2d(i,j,1)/pc_r2d(i,j)
END DO
END DO
!
!=======================================================================
! Iterate.
!=======================================================================
!
ITERATE : DO it=1,Nbico(ng)-1
DO j=Jstr,Jend
DO i=Istr,Iend
z1_r2d(i,j)=p_r2d(i,j,it)
END DO
END DO
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& z1_r2d)
# ifdef DISTRIBUTE
CALL mp_exchange2d (ng, tile, model, 1, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& z1_r2d)
# endif
!
! Compute divergence operator: v_r2d = div[h grad(z1_r2d)].
!
Ltrans=.FALSE.
CALL r2d_oper (ng, tile, &
& Ltrans, Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, z1_r2d, v_r2d)
!
! Compute dot products and "bc_ak" coefficient.
!
DO j=Jstr,Jend
DO i=Istr,Iend
z1_r2d(i,j)=r_r2d(i,j,it)/pc_r2d(i,j)
z2_r2d(i,j)=br_r2d(i,j,it)
z3_r2d(i,j)=bp_r2d(i,j,it)
END DO
END DO
CALL r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& zdf1(it), &
# ifdef MASKING
& rmask, &
# endif
& z2_r2d, z1_r2d)
CALL r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& zdf2(it), &
# ifdef MASKING
& rmask, &
# endif
& z3_r2d, v_r2d)
bc_ak(it)=zdf1(it)/zdf2(it)
!
! Solve for new iterate of "r2d_ref".
!
DO j=Jstr,Jend
DO i=Istr,Iend
r2d_ref(i,j)=r2d_ref(i,j)+bc_ak(it)*p_r2d(i,j,it)
END DO
END DO
!
!-----------------------------------------------------------------------
! Test for convergence: use "bv_r2d" as temporary storage.
!-----------------------------------------------------------------------
!
IF (it.eq.Nbico(ng)-1) THEN
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& r2d_ref)
# ifdef DISTRIBUTE
CALL mp_exchange2d (ng, tile, model, 1, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& r2d_ref)
# endif
!
! Compute divergence operator: bv_r2d = div[h grad(r2d_ref)].
!
Ltrans=.FALSE.
CALL r2d_oper (ng, tile, &
& Ltrans, Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, r2d_ref, bv_r2d)
!
! Compute dot products and report convergence value.
!
DO j=Jstr,Jend
DO i=Istr,Iend
bv_r2d(i,j)=bv_r2d(i,j)-rhs_r2d(i,j)
END DO
END DO
CALL r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& zdf4, &
# ifdef MASKING
& rmask, &
# endif
& bv_r2d, bv_r2d)
CALL r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& zdf5, &
# ifdef MASKING
& rmask, &
# endif
& rhs_r2d, rhs_r2d)
IF (Master) THEN
error=SQRT(zdf4/zdf5)
WRITE (stdout,10) error
10 FORMAT (/,' BICONJ - balance operator, error in ', &
& 'reference free-surface = ',1p,e14.7)
END IF
END IF
!
!-----------------------------------------------------------------------
! Compute new (it+1) residual and direction vectors.
!-----------------------------------------------------------------------
!
DO j=Jstr,Jend
DO i=Istr,Iend
z1_r2d(i,j)=bp_r2d(i,j,it)
END DO
END DO
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& z1_r2d)
# ifdef DISTRIBUTE
CALL mp_exchange2d (ng, tile, model, 1, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& z1_r2d)
# endif
!
! Compute divergence operator, tranpose: bv_r2d = div[h grad(z1_r2d)].
! Need to call "ad_r2d_bc" here since Ltrans is TRUE.
# ifdef DISTRIBUTE
! Need to call "ad_mp_exchange" here since Ltrans is TRUE.
# endif
!
Ltrans=.TRUE.
CALL r2d_oper (ng, tile, &
& Ltrans, Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, z1_r2d, bv_r2d)
# ifdef DISTRIBUTE
CALL ad_mp_exchange2d (ng, tile, model, 1, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& bv_r2d)
# endif
CALL ad_r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& bv_r2d)
!
! Compute new residual vectors "r" and "br".
!
DO j=Jstr,Jend
DO i=Istr,Iend
r_r2d(i,j,it+1)=r_r2d(i,j,it)-bc_ak(it)*v_r2d(i,j)
br_r2d(i,j,it+1)=br_r2d(i,j,it)-bc_ak(it)*bv_r2d(i,j)
END DO
END DO
!
! Compute dot product and "bc_ak" coefficient.
!
DO j=Jstr,Jend
DO i=Istr,Iend
z1_r2d(i,j)=r_r2d(i,j,it+1)/pc_r2d(i,j)
z2_r2d(i,j)=br_r2d(i,j,it+1)
END DO
END DO
CALL r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& zdf3(it), &
# ifdef MASKING
& rmask, &
# endif
& z1_r2d, z2_r2d)
bc_bk(it)=zdf3(it)/zdf1(it)
!
! Use recurrence relationship to compute new direction vectors
! "p" and "bp".
!
DO j=Jstr,Jend
DO i=Istr,Iend
p_r2d(i,j,it+1)=r_r2d(i,j,it+1)/pc_r2d(i,j)+ &
& bc_bk(it)*p_r2d(i,j,it)
bp_r2d(i,j,it+1)=br_r2d(i,j,it+1)/pc_r2d(i,j)+ &
& bc_bk(it)*bp_r2d(i,j,it)
END DO
END DO
END DO ITERATE
RETURN
END SUBROUTINE biconj_tile
!
!***********************************************************************
SUBROUTINE r2d_oper (ng, tile, &
& Ltrans, Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, &
& r2d_in, r2d_out)
!***********************************************************************
!
USE mod_param
USE mod_scalars
USE mod_ocean
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, Lbck
integer, intent(in) :: LBi, UBi, LBj, UBj
integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
logical, intent(in) :: Ltrans
!
# ifdef ASSUMED_SHAPE
# ifdef MASKING
real(r8), intent(in) :: umask(LBi:,LBj:)
real(r8), intent(in) :: vmask(LBi:,LBj:)
real(r8), intent(in) :: rmask(LBi:,LBj:)
# endif
real(r8), intent(in) :: h(LBi:,LBj:)
real(r8), intent(in) :: pmon_u(LBi:,LBj:)
real(r8), intent(in) :: pnom_v(LBi:,LBj:)
real(r8), intent(in) :: pm(LBi:,LBj:)
real(r8), intent(in) :: pn(LBi:,LBj:)
real(r8), intent(inout) :: pc_r2d(LBi:,LBj:)
real(r8), intent(inout) :: r2d_in(LBi:,LBj:)
real(r8), intent(out) :: r2d_out(LBi:,LBj:)
# else
# ifdef MASKING
real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pmon_u(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pnom_v(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: r2d_in(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: r2d_out(LBi:UBi,LBj:UBj)
# endif
!
! Local variable declarations.
!
integer :: i, j
real(r8) :: cff, fac
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FE
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: FX
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
! Compute divergence XI- and ETA-components: transposed operator.
! (Ltrans=.TRUE.)
!-----------------------------------------------------------------------
!
cff=0.5_r8*g
IF (Ltrans) THEN
DO j=Jstr,Jend+1
DO i=Istr,Iend+1
FX(i,j)=0.0_r8
FE(i,j)=0.0_r8
r2d_out(i,j)=0.0_r8
END DO
END DO
DO j=Jstr,Jend
DO i=Istr,Iend
!^ r2d_out(i,j)=pm(i,j)*pn(i,j)* &
!^ & (FX(i+1,j)-FX(i,j)+ &
!^ & FE(i,j+1)-FE(i,j))
!^
fac=pm(i,j)*pn(i,j)*r2d_in(i,j)
# ifdef MASKING
fac=fac*rmask(i,j)
# endif
FX(i ,j )=FX(i ,j )-fac
FX(i+1,j )=FX(i+1,j )+fac
FE(i ,j )=FE(i ,j )-fac
FE(i ,j+1)=FE(i ,j+1)+fac
END DO
END DO
DO j=Jstr,Jend+1
DO i=Istr,Iend
!^ FE(i,j)=cff*pnom_v(i,j)*(h(i,j)+h(i,j-1))* &
!^ & (r2d_in(i,j)-r2d_in(i,j-1))
!^
fac=cff*pnom_v(i,j)*(h(i,j)+h(i,j-1))*FE(i,j)
# ifdef MASKING
fac=fac*vmask(i,j)
# endif
r2d_out(i,j-1)=r2d_out(i,j-1)-fac
r2d_out(i,j )=r2d_out(i,j )+fac
END DO
END DO
DO j=Jstr,Jend
DO i=Istr,Iend+1
!^ FX(i,j)=cff*pmon_u(i,j)*(h(i,j)+h(i-1,j))* &
!^ & (r2d_in(i,j)-r2d_in(i-1,j))
!^
fac=cff*pmon_u(i,j)*(h(i,j)+h(i-1,j))*FX(i,j)
# ifdef MASKING
fac=fac*umask(i,j)
# endif
r2d_out(i-1,j)=r2d_out(i-1,j)-fac
r2d_out(i ,j)=r2d_out(i ,j)+fac
END DO
END DO
!
!-----------------------------------------------------------------------
! Compute divergence XI- and ETA-components: regular operator
! (Ltrans=.FALSE.).
!-----------------------------------------------------------------------
!
ELSE
DO j=Jstr,Jend
DO i=Istr,Iend+1
FX(i,j)=cff*pmon_u(i,j)*(h(i,j)+h(i-1,j))* &
& (r2d_in(i,j)-r2d_in(i-1,j))
# ifdef MASKING
FX(i,j)=FX(i,j)*umask(i,j)
# endif
END DO
END DO
DO j=Jstr,Jend+1
DO i=Istr,Iend
FE(i,j)=cff*pnom_v(i,j)*(h(i,j)+h(i,j-1))* &
& (r2d_in(i,j)-r2d_in(i,j-1))
# ifdef MASKING
FE(i,j)=FE(i,j)*vmask(i,j)
# endif
END DO
END DO
DO j=Jstr,Jend
DO i=Istr,Iend
r2d_out(i,j)=pm(i,j)*pn(i,j)* &
& (FX(i+1,j)-FX(i,j)+ &
& FE(i,j+1)-FE(i,j))
# ifdef MASKING
r2d_out(i,j)=r2d_out(i,j)*rmask(i,j)
# endif
END DO
END DO
END IF
!
!-----------------------------------------------------------------------
! Compute scale coefficient.
!-----------------------------------------------------------------------
!
DO j=Jstr,Jend
DO i=Istr,Iend
pc_r2d(i,j)=-pm(i,j)*pn(i,j)* &
& cff*(pnom_v(i ,j+1)*(h(i ,j+1)+h(i ,j ))+ &
& pnom_v(i ,j )*(h(i ,j )+h(i ,j-1))+ &
& pmon_u(i+1,j )*(h(i+1,j )+h(i ,j ))+ &
& pmon_u(i ,j )*(h(i ,j )+h(i-1,j )))
END DO
END DO
RETURN
END SUBROUTINE r2d_oper
!
!***********************************************************************
SUBROUTINE r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& A)
!***********************************************************************
!
USE mod_param
USE mod_scalars
!
USE exchange_2d_mod, ONLY : exchange_r2d_tile
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
integer, intent(in) :: LBi, UBi, LBj, UBj
# ifdef ASSUMED_SHAPE
real(r8), intent(inout) :: A(LBi:,LBj:)
# else
real(r8), intent(inout) :: A(LBi:UBi,LBj:UBj)
# endif
!
! Local variable declarations.
!
integer :: i, j
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
! East-West gradient boundary conditions.
!-----------------------------------------------------------------------
!
IF (.not.EWperiodic(ng)) THEN
IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
DO j=Jstr,Jend
# ifdef R2D_GRADIENT
A(Iend+1,j)=A(Iend,j)
# else
A(Iend+1,j)=0.0_r8
# endif
END DO
END IF
IF (DOMAIN(ng)%Western_Edge(tile)) THEN
DO j=Jstr,Jend
# ifdef R2D_GRADIENT
A(Istr-1,j)=A(Istr,j)
# else
A(Istr-1,j)=0.0_r8
# endif
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
! North-South gradient boundary conditions.
!-----------------------------------------------------------------------
!
IF (.not.NSperiodic(ng)) THEN
IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
DO i=Istr,Iend
# ifdef R2D_GRADIENT
A(i,Jend+1)=A(i,Jend)
# else
A(i,Jend+1)=0.0_r8
# endif
END DO
END IF
IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
DO i=Istr,Iend
# ifdef R2D_GRADIENT
A(i,Jstr-1)=A(i,Jstr)
# else
A(i,Jstr-1)=0.0_r8
# endif
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
! Boundary corners.
!-----------------------------------------------------------------------
!
IF (.not.(EWperiodic(ng).or.NSperiodic(ng))) THEN
IF (DOMAIN(ng)%SouthWest_Corner(tile)) THEN
A(Istr-1,Jstr-1)=0.5_r8*(A(Istr ,Jstr-1)+ &
& A(Istr-1,Jstr ))
END IF
IF (DOMAIN(ng)%SouthEast_Corner(tile)) THEN
A(Iend+1,Jstr-1)=0.5_r8*(A(Iend ,Jstr-1)+ &
& A(Iend+1,Jstr ))
END IF
IF (DOMAIN(ng)%NorthWest_Corner(tile)) THEN
A(Istr-1,Jend+1)=0.5_r8*(A(Istr-1,Jend )+ &
& A(Istr ,Jend+1))
END IF
IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN
A(Iend+1,Jend+1)=0.5_r8*(A(Iend+1,Jend )+ &
& A(Iend ,Jend+1))
END IF
END IF
!
!-----------------------------------------------------------------------
! Exchange boundary data.
!-----------------------------------------------------------------------
!
IF (EWperiodic(ng).or.NSperiodic(ng)) THEN
CALL exchange_r2d_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& A)
END IF
RETURN
END SUBROUTINE r2d_bc
!
!***********************************************************************
SUBROUTINE u2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& A)
!***********************************************************************
!
USE mod_param
USE mod_grid
USE mod_scalars
!
USE exchange_2d_mod, ONLY : exchange_u2d_tile
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
integer, intent(in) :: LBi, UBi, LBj, UBj
#ifdef ASSUMED_SHAPE
real(r8), intent(inout) :: A(LBi:,LBj:)
#else
real(r8), intent(inout) :: A(LBi:UBi,LBj:UBj)
#endif
!
! Local variable declarations.
!
integer :: i, j
#include "set_bounds.h"
!
!-----------------------------------------------------------------------
! East-West boundary conditions: Closed or gradient
!-----------------------------------------------------------------------
!
IF (.not.EWperiodic(ng)) THEN
IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
DO j=Jstr,Jend
A(Iend+1,j)=0.0_r8
END DO
END IF
IF (DOMAIN(ng)%Western_Edge(tile)) THEN
DO j=Jstr,Jend
A(Istr,j)=0.0_r8
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
! North-South boundary conditions: Closed (free-slip/no-slip) or
! gradient.
!-----------------------------------------------------------------------
!
IF (.not.EWperiodic(ng)) THEN
IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
DO i=IstrU,Iend
A(i,Jend+1)=0.0_r8
END DO
END IF
IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
DO i=IstrU,Iend
A(i,Jstr-1)=0.0_r8
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
! Boundary corners.
!-----------------------------------------------------------------------
!
IF (.not.(EWperiodic(ng).or.NSperiodic(ng))) THEN
IF (DOMAIN(ng)%SouthWest_Corner(tile)) THEN
A(Istr ,Jstr-1)=0.5_r8*(A(Istr+1,Jstr-1)+ &
& A(Istr ,Jstr ))
END IF
IF (DOMAIN(ng)%SouthEast_Corner(tile)) THEN
A(Iend+1,Jstr-1)=0.5_r8*(A(Iend ,Jstr-1)+ &
& A(Iend+1,Jstr ))
END IF
IF (DOMAIN(ng)%NorthWest_Corner(tile)) THEN
A(Istr ,Jend+1)=0.5_r8*(A(Istr ,Jend )+ &
& A(Istr+1,Jend+1))
END IF
IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN
A(Iend+1,Jend+1)=0.5_r8*(A(Iend+1,Jend )+ &
& A(Iend ,Jend+1))
END IF
END IF
!
!-----------------------------------------------------------------------
! Exchange boundary data.
!-----------------------------------------------------------------------
!
IF (EWperiodic(ng).or.NSperiodic(ng)) THEN
CALL exchange_u2d_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& A)
END IF
RETURN
END SUBROUTINE u2d_bc
!
!***********************************************************************
SUBROUTINE v2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& A)
!***********************************************************************
!
USE mod_param
USE mod_grid
USE mod_scalars
!
USE exchange_2d_mod, ONLY : exchange_v2d_tile
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
integer, intent(in) :: LBi, UBi, LBj, UBj
#ifdef ASSUMED_SHAPE
real(r8), intent(inout) :: A(LBi:,LBj:)
#else
real(r8), intent(inout) :: A(LBi:UBi,LBj:UBj)
#endif
!
! Local variable declarations.
!
integer :: i, j
#include "set_bounds.h"
!
!-----------------------------------------------------------------------
! East-West boundary conditions: Closed (free-slip/no-slip) or
! gradient.
!-----------------------------------------------------------------------
!
IF (.not.EWperiodic(ng)) THEN
IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
DO j=JstrV,Jend
A(Iend+1,j)=0.0_r8
END DO
END IF
IF (DOMAIN(ng)%Western_Edge(tile)) THEN
DO j=JstrV,Jend
A(Istr-1,j)=0.0_r8
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
! North-South boundary conditions: Closed or Gradient.
!-----------------------------------------------------------------------
!
IF (.not.NSperiodic(ng)) THEN
IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
DO i=Istr,Iend
A(i,Jend+1)=0.0_r8
END DO
END IF
IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
DO i=Istr,Iend
A(i,Jstr)=0.0_r8
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
! Boundary corners.
!-----------------------------------------------------------------------
!
IF (.not.(EWperiodic(ng).or.NSperiodic(ng))) THEN
IF (DOMAIN(ng)%SouthWest_Corner(tile)) THEN
A(Istr-1,Jstr )=0.5_r8*(A(Istr ,Jstr )+ &
& A(Istr-1,Jstr+1))
END IF
IF (DOMAIN(ng)%SouthEast_Corner(tile)) THEN
A(Iend+1,Jstr )=0.5_r8*(A(Iend ,Jstr )+ &
& A(Iend+1,Jstr+1))
END IF
IF (DOMAIN(ng)%NorthWest_Corner(tile)) THEN
A(Istr-1,Jend+1)=0.5_r8*(A(Istr-1,Jend )+ &
& A(Istr ,Jend+1))
END IF
IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN
A(Iend+1,Jend+1)=0.5_r8*(A(Iend+1,Jend )+ &
& A(Iend ,Jend+1))
END IF
END IF
!
!-----------------------------------------------------------------------
! Exchange boundary data.
!-----------------------------------------------------------------------
!
IF (EWperiodic(ng).or.NSperiodic(ng)) THEN
CALL exchange_v2d_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& A)
END IF
RETURN
END SUBROUTINE v2d_bc
!
!***********************************************************************
SUBROUTINE r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& DotProd, &
# ifdef MASKING
& rmask, &
# endif
& s1_zeta, s2_zeta)
!***********************************************************************
!
USE mod_param
USE mod_parallel
USE mod_ncparam
# ifdef DISTRIBUTE
!
USE distribute_mod, ONLY : mp_reduce
# endif
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, model
integer, intent(in) :: LBi, UBi, LBj, UBj
!
# ifdef ASSUMED_SHAPE
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:,LBj:)
# endif
real(r8), intent(in) :: s1_zeta(LBi:,LBj:)
real(r8), intent(in) :: s2_zeta(LBi:,LBj:)
# else
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: s1_zeta(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: s2_zeta(LBi:UBi,LBj:UBj)
# endif
!
real(r8), intent(out) :: DotProd
!
! Local variable declarations.
!
integer :: NSUB, i, j
real(r8) :: cff
real(r8) :: my_DotProd
# ifdef DISTRIBUTE
character (len=3) :: op_handle
# endif
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
! Compute dot product between "s1_zeta" and "s2_zeta".
!-----------------------------------------------------------------------
!
my_DotProd=0.0_r8
DO j=Jstr,Jend
DO i=Istr,Iend
cff=s1_zeta(i,j)*s2_zeta(i,j)
# ifdef MASKING
cff=cff*rmask(i,j)
# endif
my_DotProd=my_DotProd+cff
END DO
END DO
!
!-----------------------------------------------------------------------
! Perform parallel global reduction operations.
!-----------------------------------------------------------------------
!
# ifdef DISTRIBUTE
NSUB=1 ! distributed-memory
# else
IF (DOMAIN(ng)%SouthWest_Corner(tile).and. &
& DOMAIN(ng)%NorthEast_Corner(tile)) THEN
NSUB=1 ! non-tiled application
ELSE
NSUB=NtileX(ng)*NtileE(ng) ! tiled application
END IF
# endif
!$OMP CRITICAL (DOT_PROD)
IF (tile_count.eq.0) THEN
DotProd=0.0_r8
END IF
DotProd=DotProd+my_DotProd
tile_count=tile_count+1
IF (tile_count.eq.NSUB) THEN
tile_count=0
# ifdef DISTRIBUTE
op_handle='SUM'
CALL mp_reduce (ng, model, 1, DotProd, op_handle)
# endif
END IF
!$OMP END CRITICAL (DOT_PROD)
RETURN
END SUBROUTINE r2d_dotp
!
!***********************************************************************
SUBROUTINE tl_biconj_tile (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& Lbck, &
& h, pmon_u, pnom_v, pm, pn, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& bc_ak, bc_bk, zdf1, zdf2, zdf3, &
& pc_r2d, r_r2d, br_r2d, p_r2d, bp_r2d, &
& tl_r2d_ref, tl_rhs_r2d)
!***********************************************************************
!
USE mod_param
USE mod_scalars
USE mod_grid
USE mod_parallel
USE mod_scalars
!
# ifdef DISTRIBUTE
USE mp_exchange_mod, ONLY : ad_mp_exchange2d
USE mp_exchange_mod, ONLY : mp_exchange2d
# endif
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, model
integer, intent(in) :: LBi, UBi, LBj, UBj
integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
integer, intent(in) :: Lbck
!
# ifdef ASSUMED_SHAPE
# ifdef MASKING
real(r8), intent(in) :: umask(LBi:,LBj:)
real(r8), intent(in) :: vmask(LBi:,LBj:)
real(r8), intent(in) :: rmask(LBi:,LBj:)
# endif
real(r8), intent(in) :: h(LBi:,LBj:)
real(r8), intent(in) :: pmon_u(LBi:,LBj:)
real(r8), intent(in) :: pnom_v(LBi:,LBj:)
real(r8), intent(in) :: pm(LBi:,LBj:)
real(r8), intent(in) :: pn(LBi:,LBj:)
real(r8), intent(in) :: bc_ak(:)
real(r8), intent(in) :: bc_bk(:)
real(r8), intent(in) :: zdf1(:)
real(r8), intent(in) :: zdf2(:)
real(r8), intent(in) :: zdf3(:)
real(r8), intent(inout) :: pc_r2d(LBi:,LBj:)
real(r8), intent(inout) :: r_r2d(LBi:,LBj:,:)
real(r8), intent(inout) :: br_r2d(LBi:,LBj:,:)
real(r8), intent(inout) :: p_r2d(LBi:,LBj:,:)
real(r8), intent(inout) :: bp_r2d(LBi:,LBj:,:)
real(r8), intent(inout) :: tl_rhs_r2d(LBi:,LBj:)
real(r8), intent(inout) :: tl_r2d_ref(LBi:,LBj:)
# else
# ifdef MASKING
real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pmon_u(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pnom_v(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: bc_ak(Nbico(ng))
real(r8), intent(inout) :: bc_bk(Nbico(ng))
real(r8), intent(inout) :: zdf1(Nbico(ng))
real(r8), intent(inout) :: zdf2(Nbico(ng))
real(r8), intent(inout) :: zdf3(Nbico(ng))
real(r8), intent(inout) :: pc_r2d(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: r_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
real(r8), intent(inout) :: br_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
real(r8), intent(inout) :: p_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
real(r8), intent(inout) :: tl_bp_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
real(r8), intent(inout) :: tl_rhs_r2d(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: tl_r2d_ref(LBi:UBi,LBj:UBj)
# endif
!
! Local variable declarations.
!
logical :: Ltrans
integer :: i, j, it
real(r8) :: tl_zdf1, tl_zdf2, tl_zdf3, tl_bc_ak, tl_bc_bk
real(r8), dimension(LBi:UBi,LBj:UBj) :: bv_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: v_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: z1_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: z2_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: z3_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: tl_bp_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: tl_br_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: tl_bv_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: tl_p_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: tl_r_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: tl_v_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: tl_z1_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: tl_z2_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: tl_z3_r2d
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
! Tangent linear of biconjugate gradient algorithm.
!-----------------------------------------------------------------------
!
! Initialize local arrays.
!
DO j=LBj,UBj
DO i=LBi,UBi
bv_r2d(i,j)=0.0_r8
v_r2d(i,j)=0.0_r8
z1_r2d(i,j)=0.0_r8
z2_r2d(i,j)=0.0_r8
z3_r2d(i,j)=0.0_r8
tl_bp_r2d(i,j)=0.0_r8
tl_br_r2d(i,j)=0.0_r8
tl_bv_r2d(i,j)=0.0_r8
tl_p_r2d(i,j)=0.0_r8
tl_r_r2d(i,j)=0.0_r8
tl_v_r2d(i,j)=0.0_r8
tl_z1_r2d(i,j)=0.0_r8
tl_z2_r2d(i,j)=0.0_r8
tl_z3_r2d(i,j)=0.0_r8
END DO
END DO
!
! Compute starting value of divergence TLM operator:
! tl_z1_r2d = div[h grad(tl_r2d_ref)].
!
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& tl_r2d_ref)
# ifdef DISTRIBUTE
CALL mp_exchange2d (ng, tile, model, 1, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& tl_r2d_ref)
# endif
CALL tl_r2d_oper (ng, tile, &
& Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, tl_r2d_ref, tl_z1_r2d)
!
! Set the initial values for residual vectors "tl_r" and "tl_br". Then,
! use recurrence relationship to compute direction vectors "tl_p" and
! "tl_bp".
!
DO j=Jstr,Jend
DO i=Istr,Iend
tl_r_r2d(i,j)=tl_rhs_r2d(i,j)-tl_z1_r2d(i,j)
tl_br_r2d(i,j)=tl_r_r2d(i,j)
END DO
END DO
DO j=Jstr,Jend
DO i=Istr,Iend
tl_p_r2d(i,j)=tl_r_r2d(i,j)/pc_r2d(i,j)
tl_bp_r2d(i,j)=tl_br_r2d(i,j)/pc_r2d(i,j)
END DO
END DO
!
!=======================================================================
! Iterate.
!=======================================================================
!
ITERATE : DO it=1,Nbico(ng)-1
DO j=Jstr,Jend
DO i=Istr,Iend
z1_r2d(i,j)=p_r2d(i,j,it)
tl_z1_r2d(i,j)=tl_p_r2d(i,j)
END DO
END DO
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& z1_r2d)
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& tl_z1_r2d)
# ifdef DISTRIBUTE
CALL mp_exchange2d (ng, tile, model, 2, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& z1_r2d, &
& tl_z1_r2d)
# endif
!
! Compute divergence NLM operator: v_r2d = div[h grad(z1_r2d)].
!
Ltrans=.FALSE.
CALL r2d_oper (ng, tile, &
& Ltrans, Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, z1_r2d, v_r2d)
!
! Compute divergence TLM operator: tl_v_r2d = div[h grad(tl_z1_r2d)].
!
CALL tl_r2d_oper (ng, tile, &
& Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, tl_z1_r2d, tl_v_r2d)
!
! Compute dot products and "tl_bc_ak" coefficient.
!
DO j=Jstr,Jend
DO i=Istr,Iend
z1_r2d(i,j)=r_r2d(i,j,it)/pc_r2d(i,j)
z2_r2d(i,j)=br_r2d(i,j,it)
z3_r2d(i,j)=bp_r2d(i,j,it)
tl_z1_r2d(i,j)=tl_r_r2d(i,j)/pc_r2d(i,j)
tl_z2_r2d(i,j)=tl_br_r2d(i,j)
tl_z3_r2d(i,j)=tl_bp_r2d(i,j)
END DO
END DO
CALL tl_r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& tl_zdf1, &
# ifdef MASKING
& rmask, &
# endif
& z2_r2d, z1_r2d, tl_z2_r2d, tl_z1_r2d)
CALL tl_r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& tl_zdf2, &
# ifdef MASKING
& rmask, &
# endif
& z3_r2d, v_r2d, tl_z3_r2d, tl_v_r2d)
tl_bc_ak=tl_zdf1/zdf2(it)-tl_zdf2*bc_ak(it)/zdf2(it)
!
! Solve for new iterate of "tl_r2d_ref".
!
DO j=Jstr,Jend
DO i=Istr,Iend
tl_r2d_ref(i,j)=tl_r2d_ref(i,j)+ &
& tl_bc_ak*p_r2d(i,j,it)+ &
& bc_ak(it)*tl_p_r2d(i,j)
END DO
END DO
!
!-----------------------------------------------------------------------
! Compute new (it+1) residual and direction vectors.
!-----------------------------------------------------------------------
!
DO j=Jstr,Jend
DO i=Istr,Iend
z1_r2d(i,j)=bp_r2d(i,j,it)
tl_z1_r2d(i,j)=tl_bp_r2d(i,j)
tl_bv_r2d(i,j)=0.0_r8
END DO
END DO
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& z1_r2d)
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& tl_z1_r2d)
# ifdef DISTRIBUTE
CALL mp_exchange2d (ng, tile, model, 1, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& z1_r2d)
# endif
!
! Compute divergence operator, tranpose: bv_r2d = div[h grad(z1_r2d)]
! for basic state and TLM (need to use "ad_r2d_oper").
! Need to call "ad_r2d_bc" here since Ltrans is TRUE.
# ifdef DISTRIBUTE
! Need to call "ad_mp_exchange" here since Ltrans is TRUE.
# endif
!
Ltrans=.TRUE.
CALL r2d_oper (ng, tile, &
& Ltrans, Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, z1_r2d, bv_r2d)
CALL ad_r2d_oper (ng, tile, &
& Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, tl_bv_r2d, tl_z1_r2d)
# ifdef DISTRIBUTE
CALL ad_mp_exchange2d (ng, tile, model, 2, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& bv_r2d, &
& tl_bv_r2d)
# endif
CALL ad_r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& bv_r2d)
CALL ad_r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& tl_bv_r2d)
!
! Compute new residual vectors "tl_r" and "tl_br".
!
DO j=Jstr,Jend
DO i=Istr,Iend
tl_r_r2d(i,j)=tl_r_r2d(i,j)- &
& tl_bc_ak*v_r2d(i,j)-bc_ak(it)*tl_v_r2d(i,j)
tl_br_r2d(i,j)=tl_br_r2d(i,j)- &
& tl_bc_ak*bv_r2d(i,j)-bc_ak(it)*tl_bv_r2d(i,j)
END DO
END DO
!
! Compute dot product and "tl_bc_ak" coefficient.
!
DO j=Jstr,Jend
DO i=Istr,Iend
z1_r2d(i,j)=r_r2d(i,j,it+1)/pc_r2d(i,j)
z2_r2d(i,j)=br_r2d(i,j,it+1)
tl_z1_r2d(i,j)=tl_r_r2d(i,j)/pc_r2d(i,j)
tl_z2_r2d(i,j)=tl_br_r2d(i,j)
END DO
END DO
CALL tl_r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& tl_zdf3, &
# ifdef MASKING
& rmask, &
# endif
& z1_r2d, z2_r2d, tl_z1_r2d, tl_z2_r2d)
tl_bc_bk=tl_zdf3/zdf1(it)-tl_zdf1*bc_bk(it)/zdf1(it)
!
! Use recurrence relationship to compute new direction vectors
! "tl_p" and "tl_bp".
!
DO j=Jstr,Jend
DO i=Istr,Iend
tl_p_r2d(i,j)=tl_r_r2d(i,j)/pc_r2d(i,j)+ &
& tl_bc_bk*p_r2d(i,j,it)+ &
& bc_bk(it)*tl_p_r2d(i,j)
tl_bp_r2d(i,j)=tl_br_r2d(i,j)/pc_r2d(i,j)+ &
& tl_bc_bk*bp_r2d(i,j,it)+ &
& bc_bk(it)*tl_bp_r2d(i,j)
END DO
END DO
END DO ITERATE
RETURN
END SUBROUTINE tl_biconj_tile
!
!***********************************************************************
SUBROUTINE tl_r2d_oper (ng, tile, &
& Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, &
& tl_r2d_in, tl_r2d_out)
!***********************************************************************
!
USE mod_param
USE mod_scalars
USE mod_ocean
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, Lbck
integer, intent(in) :: LBi, UBi, LBj, UBj
integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
!
# ifdef ASSUMED_SHAPE
# ifdef MASKING
real(r8), intent(in) :: umask(LBi:,LBj:)
real(r8), intent(in) :: vmask(LBi:,LBj:)
real(r8), intent(in) :: rmask(LBi:,LBj:)
# endif
real(r8), intent(in) :: h(LBi:,LBj:)
real(r8), intent(in) :: pmon_u(LBi:,LBj:)
real(r8), intent(in) :: pnom_v(LBi:,LBj:)
real(r8), intent(in) :: pm(LBi:,LBj:)
real(r8), intent(in) :: pn(LBi:,LBj:)
real(r8), intent(inout) :: pc_r2d(LBi:,LBj:)
real(r8), intent(inout) :: tl_r2d_in(LBi:,LBj:)
real(r8), intent(out) :: tl_r2d_out(LBi:,LBj:)
# else
# ifdef MASKING
real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pmon_u(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pnom_v(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: tl_r2d_in(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: tl_r2d_out(LBi:UBi,LBj:UBj)
# endif
!
! Local variable declarations.
!
integer :: i, j
real(r8) :: cff, tl_fac
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FE
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tl_FX
# include "set_bounds.h"
!
!----------------------------------------------------------------------
! Compute tl_r2d_out = div ( h grad(tl_r2d_in) ).
!----------------------------------------------------------------------
!
cff=0.5_r8*g
DO j=Jstr,Jend
DO i=Istr,Iend+1
tl_FX(i,j)=cff*pmon_u(i,j)*(h(i,j)+h(i-1,j))* &
& (tl_r2d_in(i,j)-tl_r2d_in(i-1,j))
# ifdef MASKING
tl_FX(i,j)=tl_FX(i,j)*umask(i,j)
# endif
END DO
END DO
DO j=Jstr,Jend+1
DO i=Istr,Iend
tl_FE(i,j)=cff*pnom_v(i,j)*(h(i,j)+h(i,j-1))* &
& (tl_r2d_in(i,j)-tl_r2d_in(i,j-1))
# ifdef MASKING
tl_FE(i,j)=tl_FE(i,j)*vmask(i,j)
# endif
END DO
END DO
DO j=Jstr,Jend
DO i=Istr,Iend
tl_r2d_out(i,j)=pm(i,j)*pn(i,j)* &
& (tl_FX(i+1,j)-tl_FX(i,j)+ &
& tl_FE(i,j+1)-tl_FE(i,j))
# ifdef MASKING
tl_r2d_out(i,j)=tl_r2d_out(i,j)*rmask(i,j)
# endif
END DO
END DO
RETURN
END SUBROUTINE tl_r2d_oper
!
!***********************************************************************
SUBROUTINE tl_r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& tl_DotProd, &
# ifdef MASKING
& rmask, &
# endif
& s1_zeta, s2_zeta, &
& tl_s1_zeta, tl_s2_zeta)
!***********************************************************************
!
USE mod_param
USE mod_parallel
USE mod_ncparam
# ifdef DISTRIBUTE
!
USE distribute_mod, ONLY : mp_reduce
# endif
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, model
integer, intent(in) :: LBi, UBi, LBj, UBj
!
# ifdef ASSUMED_SHAPE
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:,LBj:)
# endif
real(r8), intent(in) :: s1_zeta(LBi:,LBj:)
real(r8), intent(in) :: s2_zeta(LBi:,LBj:)
real(r8), intent(in) :: tl_s1_zeta(LBi:,LBj:)
real(r8), intent(in) :: tl_s2_zeta(LBi:,LBj:)
# else
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: s1_zeta(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: s2_zeta(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: tl_s1_zeta(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: tl_s2_zeta(LBi:UBi,LBj:UBj)
# endif
!
real(r8), intent(out) :: tl_DotProd
!
! Local variable declarations.
!
integer :: NSUB, i, j
real(r8) :: tl_cff
real(r8) :: tl_my_DotProd
# ifdef DISTRIBUTE
character (len=3) :: op_handle
# endif
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
! Compute dot product between "tl_s1_zeta" and "tl_s2_zeta".
!-----------------------------------------------------------------------
!
tl_my_DotProd=0.0_r8
DO j=Jstr,Jend
DO i=Istr,Iend
tl_cff=tl_s1_zeta(i,j)*s2_zeta(i,j)+ &
& s1_zeta(i,j)*tl_s2_zeta(i,j)
# ifdef MASKING
tl_cff=tl_cff*rmask(i,j)
# endif
tl_my_DotProd=tl_my_DotProd+tl_cff
END DO
END DO
!
!-----------------------------------------------------------------------
! Perform parallel global reduction operations.
!-----------------------------------------------------------------------
!
# ifdef DISTRIBUTE
NSUB=1 ! distributed-memory
# else
IF (DOMAIN(ng)%SouthWest_Corner(tile).and. &
& DOMAIN(ng)%NorthEast_Corner(tile)) THEN
NSUB=1 ! non-tiled application
ELSE
NSUB=NtileX(ng)*NtileE(ng) ! tiled application
END IF
# endif
!$OMP CRITICAL (TL_DOT_PROD)
IF (tile_count.eq.0) THEN
tl_DotProd=0.0_r8
END IF
tl_DotProd=tl_DotProd+tl_my_DotProd
tile_count=tile_count+1
IF (tile_count.eq.NSUB) THEN
tile_count=0
# ifdef DISTRIBUTE
op_handle='SUM'
CALL mp_reduce (ng, model, 1, tl_DotProd, op_handle)
# endif
END IF
!$OMP END CRITICAL (TL_DOT_PROD)
RETURN
END SUBROUTINE tl_r2d_dotp
!
!***********************************************************************
SUBROUTINE ad_biconj_tile (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& Lbck, &
& h, pmon_u, pnom_v, pm, pn, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& bc_ak, bc_bk, zdf1, zdf2, zdf3, &
& pc_r2d, r_r2d, br_r2d, p_r2d, bp_r2d, &
& ad_r2d_ref, ad_rhs_r2d)
!***********************************************************************
!
USE mod_param
USE mod_scalars
USE mod_grid
USE mod_parallel
USE mod_scalars
!
# ifdef DISTRIBUTE
USE distribute_mod, ONLY : mp_reduce
USE mp_exchange_mod, ONLY : ad_mp_exchange2d
USE mp_exchange_mod, ONLY : mp_exchange2d
# endif
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, model
integer, intent(in) :: LBi, UBi, LBj, UBj
integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
integer, intent(in) :: Lbck
!
# ifdef ASSUMED_SHAPE
# ifdef MASKING
real(r8), intent(in) :: umask(LBi:,LBj:)
real(r8), intent(in) :: vmask(LBi:,LBj:)
real(r8), intent(in) :: rmask(LBi:,LBj:)
# endif
real(r8), intent(in) :: h(LBi:,LBj:)
real(r8), intent(in) :: pmon_u(LBi:,LBj:)
real(r8), intent(in) :: pnom_v(LBi:,LBj:)
real(r8), intent(in) :: pm(LBi:,LBj:)
real(r8), intent(in) :: pn(LBi:,LBj:)
real(r8), intent(in) :: bc_ak(:)
real(r8), intent(in) :: bc_bk(:)
real(r8), intent(in) :: zdf1(:)
real(r8), intent(in) :: zdf2(:)
real(r8), intent(in) :: zdf3(:)
real(r8), intent(inout) :: pc_r2d(LBi:,LBj:)
real(r8), intent(inout) :: r_r2d(LBi:,LBj:,:)
real(r8), intent(inout) :: br_r2d(LBi:,LBj:,:)
real(r8), intent(inout) :: p_r2d(LBi:,LBj:,:)
real(r8), intent(inout) :: bp_r2d(LBi:,LBj:,:)
real(r8), intent(inout) :: ad_rhs_r2d(LBi:,LBj:)
real(r8), intent(inout) :: ad_r2d_ref(LBi:,LBj:)
# else
# ifdef MASKING
real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pmon_u(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pnom_v(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: bc_ak(Nbico(ng))
real(r8), intent(in) :: bc_bk(Nbico(ng))
real(r8), intent(in) :: zdf1(Nbico(ng))
real(r8), intent(in) :: zdf2(Nbico(ng))
real(r8), intent(in) :: zdf3(Nbico(ng))
real(r8), intent(inout) :: pc_r2d(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: r_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
real(r8), intent(inout) :: br_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
real(r8), intent(inout) :: p_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
real(r8), intent(inout) :: ad_bp_r2d(LBi:UBi,LBj:UBj,Nbico(ng))
real(r8), intent(inout) :: ad_rhs_r2d(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: ad_r2d_ref(LBi:UBi,LBj:UBj)
# endif
!
! Local variable declarations.
!
# ifdef DISTRIBUTE
character (len=3) :: op_handle
# endif
logical :: Ltrans
integer :: i, j, it, NSUB
real(r8) :: ad_zdf1, ad_zdf2, ad_zdf3, ad_bc_ak, ad_bc_bk
real(r8) :: my_ad_bc_ak, my_ad_bc_bk
real(r8), dimension(LBi:UBi,LBj:UBj) :: bv_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: v_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: z1_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: z2_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: z3_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: ad_bp_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: ad_br_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: ad_bv_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: ad_p_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: ad_r_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: ad_v_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: ad_z1_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: ad_z2_r2d
real(r8), dimension(LBi:UBi,LBj:UBj) :: ad_z3_r2d
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
! Adjoint of biconjugate gradient algorithm.
!-----------------------------------------------------------------------
!
! Initialize local variables.
!
DO j=LBj,UBj
DO i=LBi,UBi
bv_r2d(i,j)=0.0_r8
v_r2d(i,j)=0.0_r8
z1_r2d(i,j)=0.0_r8
z2_r2d(i,j)=0.0_r8
z3_r2d(i,j)=0.0_r8
ad_bp_r2d(i,j)=0.0_r8
ad_br_r2d(i,j)=0.0_r8
ad_bv_r2d(i,j)=0.0_r8
ad_p_r2d(i,j)=0.0_r8
ad_r_r2d(i,j)=0.0_r8
ad_v_r2d(i,j)=0.0_r8
ad_z1_r2d(i,j)=0.0_r8
ad_z2_r2d(i,j)=0.0_r8
ad_z3_r2d(i,j)=0.0_r8
END DO
END DO
ad_bc_ak=0.0_r8
ad_bc_bk=0.0_r8
ad_zdf1=0.0_r8
ad_zdf2=0.0_r8
ad_zdf3=0.0_r8
my_ad_bc_ak=0.0_r8
my_ad_bc_bk=0.0_r8
!
!=======================================================================
! Iterate in reverse order.
!=======================================================================
!
ITERATE : DO it=Nbico(ng)-1,1,-1
!
! Compute "v_r2d" and "bv_r2d"
!
DO j=Jstr,Jend
DO i=Istr,Iend
z1_r2d(i,j)=p_r2d(i,j,it)
z2_r2d(i,j)=bp_r2d(i,j,it)
END DO
END DO
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& z1_r2d)
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& z2_r2d)
# ifdef DISTRIBUTE
CALL mp_exchange2d (ng, tile, model, 2, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& z1_r2d, &
& z2_r2d)
# endif
!
! Compute divergence operator: v_r2d = div[h grad(z1_r2d)].
!
Ltrans=.FALSE.
CALL r2d_oper (ng, tile, &
& Ltrans, Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, z1_r2d, v_r2d)
!
! Compute divergence operator, tranpose: bv_r2d = div[h grad(z2_r2d)].
!
Ltrans=.TRUE.
CALL r2d_oper (ng, tile, &
& Ltrans, Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, z2_r2d, bv_r2d)
# ifdef DISTRIBUTE
CALL ad_mp_exchange2d (ng, tile, model, 1, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& bv_r2d)
# endif
CALL ad_r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& bv_r2d)
!
! Adjoint of use recurrence relationship to compute new direction
! vectors "tl_p" and "tl_bp".
!
DO j=Jstr,Jend
DO i=Istr,Iend
!^ tl_bp_r2d(i,j)=tl_br_r2d(i,j)/pc_r2d(i,j)+ &
!^ & tl_bc_bk*bp_r2d(i,j,it)+ &
!^ & bc_bk(it)*tl_bp_r2d(i,j)
!^
my_ad_bc_bk=my_ad_bc_bk+bp_r2d(i,j,it)*ad_bp_r2d(i,j)
ad_br_r2d(i,j)=ad_br_r2d(i,j)+ad_bp_r2d(i,j)/pc_r2d(i,j)
ad_bp_r2d(i,j)=bc_bk(it)*ad_bp_r2d(i,j)
!^ tl_p_r2d(i,j)=tl_r_r2d(i,j)/pc_r2d(i,j)+ &
!^ & tl_bc_bk*p_r2d(i,j,it)+ &
!^ & bc_bk(it)*tl_p_r2d(i,j)
!^
my_ad_bc_bk=my_ad_bc_bk+p_r2d(i,j,it)*ad_p_r2d(i,j)
ad_r_r2d(i,j)=ad_r_r2d(i,j)+ad_p_r2d(i,j)/pc_r2d(i,j)
ad_p_r2d(i,j)=bc_bk(it)*ad_p_r2d(i,j)
END DO
END DO
!
! Perform parallel global reduction operation on "my_ad_bc_bk".
!
# ifdef DISTRIBUTE
NSUB=1 ! distributed-memory
# else
IF (DOMAIN(ng)%SouthWest_Corner(tile).and. &
& DOMAIN(ng)%NorthEast_Corner(tile)) THEN
NSUB=1 ! non-tiled application
ELSE
NSUB=NtileX(ng)*NtileE(ng) ! tiled application
END IF
# endif
!$OMP CRITICAL (AD_DOT_PROD1)
IF (tile_count.eq.0) THEN
ad_bc_bk=0.0_r8
END IF
ad_bc_bk=ad_bc_bk+my_ad_bc_bk
tile_count=tile_count+1
IF (tile_count.eq.NSUB) THEN
tile_count=0
# ifdef DISTRIBUTE
op_handle='SUM'
CALL mp_reduce (ng, model, 1, ad_bc_bk, op_handle)
# endif
END IF
my_ad_bc_bk=0.0_r8
!$OMP END CRITICAL (AD_DOT_PROD1)
!
! Adjoint of compute dot product and "tl_bc_ak" coefficient.
!
!^ tl_bc_bk=tl_zdf3/zdf1(it)- &
!^ & tl_zdf1*bc_bk(it)/zdf1(it)
!^
ad_zdf1=ad_zdf1-ad_bc_bk*bc_bk(it)/zdf1(it)
ad_zdf3=ad_zdf3+ad_bc_bk/zdf1(it)
ad_bc_bk=0.0
DO j=Jstr,Jend
DO i=Istr,Iend
z1_r2d(i,j)=r_r2d(i,j,it+1)/pc_r2d(i,j)
z2_r2d(i,j)=br_r2d(i,j,it+1)
END DO
END DO
CALL ad_r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& ad_zdf3, &
# ifdef MASKING
& rmask, &
# endif
& z1_r2d, z2_r2d, ad_z1_r2d, ad_z2_r2d)
DO j=Jstr,Jend
DO i=Istr,Iend
!^ tl_z2_r2d(i,j)=tl_br_r2d(i,j)
!^
ad_br_r2d(i,j)=ad_br_r2d(i,j)+ad_z2_r2d(i,j)
ad_z2_r2d(i,j)=0.0_r8
!^ tl_z1_r2d(i,j)=tl_r_r2d(i,j)/pc_r2d(i,j)
!^
ad_r_r2d(i,j)=ad_r_r2d(i,j)+ad_z1_r2d(i,j)/pc_r2d(i,j)
ad_z1_r2d(i,j)=0.0_r8
END DO
END DO
!
! Adjoint of compute new residual vectors "tl_r" and "tl_br".
!
DO j=Jstr,Jend
DO i=Istr,Iend
!^ tl_br_r2d(i,j)=tl_br_r2d(i,j)-tl_bc_ak*bv_r2d(i,j)- &
!^ & bc_ak(it)*tl_bv_r2d(i,j)
!^
ad_bv_r2d(i,j)=ad_bv_r2d(i,j)-bc_ak(it)*ad_br_r2d(i,j)
my_ad_bc_ak=my_ad_bc_ak-bv_r2d(i,j)*ad_br_r2d(i,j)
!^ tl_r_r2d(i,j)=tl_r_r2d(i,j)-tl_bc_ak*v_r2d(i,j)- &
!^ & bc_ak(it)*tl_v_r2d(i,j)
!^
ad_v_r2d(i,j)=ad_v_r2d(i,j)-bc_ak(it)*ad_r_r2d(i,j)
my_ad_bc_ak=my_ad_bc_ak-v_r2d(i,j)*ad_r_r2d(i,j)
END DO
END DO
!
! Adjoint of compute divergence operator, tranpose:
! tl_bv_r2d = div[h grad(tl_z1_r2d)]. Need to use "tl_r2d_oper" here.
!
CALL r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& ad_bv_r2d)
# ifdef DISTRIBUTE
CALL mp_exchange2d (ng, tile, model, 1, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& ad_bv_r2d)
# endif
CALL tl_r2d_oper (ng, tile, &
& Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, ad_bv_r2d, ad_z1_r2d)
!
! Adjoint of loading "tl_z1_r2d".
!
DO j=Jstr,Jend
DO i=Istr,Iend
!^ tl_z1_r2d(i,j)=tl_bp_r2d(i,j)
!^
ad_bp_r2d(i,j)=ad_bp_r2d(i,j)+ad_z1_r2d(i,j)
ad_z1_r2d(i,j)=0.0_r8
ad_bv_r2d(i,j)=0.0_r8 ! cleared because it was used in
! "tl_r2d_oper"
END DO
END DO
!
! Adjoint of solve for new iterate "tl_r2d_ref".
!
DO j=Jstr,Jend
DO i=Istr,Iend
!^ tl_r2d_ref(i,j)=tl_r2d_ref(i,j)+tl_bc_ak*p_r2d(i,j,it)+ &
!^ & bc_ak(it)*tl_p_r2d(i,j)
!^
ad_p_r2d(i,j)=ad_p_r2d(i,j)+bc_ak(it)*ad_r2d_ref(i,j)
my_ad_bc_ak=my_ad_bc_ak+p_r2d(i,j,it)*ad_r2d_ref(i,j)
END DO
END DO
!
! Perform parallel global reduction operation on "my_ad_bc_ak".
!
# ifdef DISTRIBUTE
NSUB=1 ! distributed-memory
# else
IF (DOMAIN(ng)%SouthWest_Corner(tile).and. &
& DOMAIN(ng)%NorthEast_Corner(tile)) THEN
NSUB=1 ! non-tiled application
ELSE
NSUB=NtileX(ng)*NtileE(ng) ! tiled application
END IF
# endif
!$OMP CRITICAL (AD_DOT_PROD2)
IF (tile_count.eq.0) THEN
ad_bc_ak=0.0_r8
END IF
ad_bc_ak=ad_bc_ak+my_ad_bc_ak
tile_count=tile_count+1
IF (tile_count.eq.NSUB) THEN
tile_count=0
# ifdef DISTRIBUTE
op_handle='SUM'
CALL mp_reduce (ng, model, 1, ad_bc_ak, op_handle)
# endif
END IF
my_ad_bc_ak=0.0_r8
!$OMP END CRITICAL (AD_DOT_PROD2)
!
! Adjoint of compute dot products and "tl_bc_ak" coefficient.
!
!^ tl_bc_ak=tl_zdf1/zdf2(it)- &
!^ & tl_zdf2*bc_ak(it)/zdf2(it)
!^
ad_zdf2=ad_zdf2-ad_bc_ak*bc_ak(it)/zdf2(it)
ad_zdf1=ad_zdf1+ad_bc_ak/zdf2(it)
ad_bc_ak=0.0
DO j=Jstr,Jend
DO i=Istr,Iend
z1_r2d(i,j)=r_r2d(i,j,it)/pc_r2d(i,j)
z2_r2d(i,j)=br_r2d(i,j,it)
z3_r2d(i,j)=bp_r2d(i,j,it)
END DO
END DO
CALL ad_r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& ad_zdf2, &
# ifdef MASKING
& rmask, &
# endif
& z3_r2d, v_r2d, ad_z3_r2d, ad_v_r2d)
CALL ad_r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& ad_zdf1, &
# ifdef MASKING
& rmask, &
# endif
& z2_r2d, z1_r2d, ad_z2_r2d, ad_z1_r2d)
DO j=Jstr,Jend
DO i=Istr,Iend
!^ tl_z1_r2d(i,j)=tl_r_r2d(i,j)/pc_r2d(i,j)
!^
ad_r_r2d(i,j)=ad_r_r2d(i,j)+ad_z1_r2d(i,j)/pc_r2d(i,j)
ad_z1_r2d(i,j)=0.0_r8
!^ tl_z2_r2d(i,j)=tl_br_r2d(i,j)
!^
ad_br_r2d(i,j)=ad_br_r2d(i,j)+ad_z2_r2d(i,j)
ad_z2_r2d(i,j)=0.0_r8
!^ tl_z3_r2d(i,j)=tl_bp_r2d(i,j)
!^
ad_bp_r2d(i,j)=ad_bp_r2d(i,j)+ad_z3_r2d(i,j)
ad_z3_r2d(i,j)=0.0_r8
END DO
END DO
!
! Adjoint of compute divergence operator:
! tl_v_r2d = div[h grad(tl_z1_r2d)].
!
CALL ad_r2d_oper (ng, tile, &
& Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, ad_z1_r2d, ad_v_r2d)
# ifdef DISTRIBUTE
CALL ad_mp_exchange2d (ng, tile, model, 1, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& ad_z1_r2d)
# endif
CALL ad_r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& ad_z1_r2d)
DO j=Jstr,Jend
DO i=Istr,Iend
!^ tl_z1_r2d(i,j)=tl_p_r2d(i,j)
!^
ad_p_r2d(i,j)=ad_p_r2d(i,j)+ad_z1_r2d(i,j)
ad_z1_r2d(i,j)=0.0_r8
END DO
END DO
END DO ITERATE
!
! Adjoint of set the initial values for residual vectors "tl_r" and
! "tl_br". Then, use recurrence relationship to compute direction
! vectors "tl_p" and "tl_bp".
!
DO j=Jstr,Jend
DO i=Istr,Iend
!^ tl_p_r2d(i,j)=tl_r_r2d(i,j)/pc_r2d(i,j)
!^
ad_r_r2d(i,j)=ad_r_r2d(i,j)+ad_p_r2d(i,j)/pc_r2d(i,j)
ad_p_r2d(i,j)=0.0
!^ tl_bp_r2d(i,j)=tl_br_r2d(i,j)/pc_r2d(i,j)
!^
ad_br_r2d(i,j)=ad_br_r2d(i,j)+ad_bp_r2d(i,j)/pc_r2d(i,j)
ad_bp_r2d(i,j)=0.0
END DO
END DO
DO j=Jstr,Jend
DO i=Istr,Iend
!^ tl_br_r2d(i,j)=tl_r_r2d(i,j)
!^
ad_r_r2d(i,j)=ad_r_r2d(i,j)+ad_br_r2d(i,j)
ad_br_r2d(i,j)=0.0
!^ tl_r_r2d(i,j)=tl_rhs_r2d(i,j)-tl_z1_r2d(i,j)
!^
ad_rhs_r2d(i,j)=ad_rhs_r2d(i,j)+ad_r_r2d(i,j)
ad_z1_r2d(i,j)=ad_z1_r2d(i,j)-ad_r_r2d(i,j)
ad_r_r2d(i,j)=0.0_r8
END DO
END DO
!
! Adjoint of compute starting value of divergence TLM operator:
! tl_z1_r2d = div[h grad(tl_r2d_ref)].
!
CALL ad_r2d_oper (ng, tile, &
& Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, ad_r2d_ref, ad_z1_r2d)
# ifdef DISTRIBUTE
CALL ad_mp_exchange2d (ng, tile, model, 1, &
& LBi, UBi, LBj, UBj, &
& NghostPoints, &
& EWperiodic(ng), NSperiodic(ng), &
& ad_r2d_ref)
# endif
CALL ad_r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& ad_r2d_ref)
RETURN
END SUBROUTINE ad_biconj_tile
!
!***********************************************************************
SUBROUTINE ad_r2d_oper (ng, tile, &
& Lbck, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
# ifdef MASKING
& umask, vmask, rmask, &
# endif
& h, &
& pmon_u, pnom_v, pm, pn, &
& pc_r2d, &
& ad_r2d_in, ad_r2d_out)
!***********************************************************************
!
USE mod_param
USE mod_scalars
USE mod_ocean
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
integer, intent(in) :: LBi, UBi, LBj, UBj
integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
integer, intent(in) :: Lbck
!
# ifdef ASSUMED_SHAPE
# ifdef MASKING
real(r8), intent(in) :: umask(LBi:,LBj:)
real(r8), intent(in) :: vmask(LBi:,LBj:)
real(r8), intent(in) :: rmask(LBi:,LBj:)
# endif
real(r8), intent(in) :: h(LBi:,LBj:)
real(r8), intent(in) :: pmon_u(LBi:,LBj:)
real(r8), intent(in) :: pnom_v(LBi:,LBj:)
real(r8), intent(in) :: pm(LBi:,LBj:)
real(r8), intent(in) :: pn(LBi:,LBj:)
real(r8), intent(inout) :: pc_r2d(LBi:,LBj:)
real(r8), intent(inout) :: ad_r2d_in(LBi:,LBj:)
real(r8), intent(inout) :: ad_r2d_out(LBi:,LBj:)
# else
# ifdef MASKING
real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pmon_u(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pnom_v(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: ad_r2d_in(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: ad_r2d_out(LBi:UBi,LBj:UBj)
# endif
!
! Local variable declarations.
!
integer :: i, j
real(r8) :: adfac, cff
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_FE
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ad_FX
# include "set_bounds.h"
!
!----------------------------------------------------------------------
! Adjoint of compute tl_r2d_out = div( h grad(ad_r2d_in) ).
!----------------------------------------------------------------------
!
! Initialize local arrays.
!
DO j=Jstr,Jend+1
DO i=Istr,Iend+1
ad_FX(i,j)=0.0_r8
ad_FE(i,j)=0.0_r8
END DO
END DO
!
! Compute adjoint of divergence XI- and ETA-components.
!
cff=0.5_r8*g
DO j=Jstr,Jend
DO i=Istr,Iend
# ifdef MASKING
!^ tl_r2d_out(i,j)=tl_r2d_out(i,j)*rmask(i,j)
!^
ad_r2d_out(i,j)=ad_r2d_out(i,j)*rmask(i,j)
# endif
!^ tl_r2d_out(i,j)=pm(i,j)*pn(i,j)* &
!^ & (tl_FX(i+1,j)-tl_FX(i,j)+ &
!^ & tl_FE(i,j+1)-tl_FE(i,j))
!^
adfac=pm(i,j)*pn(i,j)*ad_r2d_out(i,j)
ad_FX(i ,j)=ad_FX(i ,j)-adfac
ad_FX(i+1,j)=ad_FX(i+1,j)+adfac
ad_FE(i,j+1)=ad_FE(i,j+1)+adfac
ad_FE(i,j )=ad_FE(i,j )-adfac
ad_r2d_out(i,j)=0.0_r8
END DO
END DO
DO j=Jstr,Jend+1
DO i=Istr,Iend
# ifdef MASKING
!^ tl_FE(i,j)=tl_FE(i,j)*vmask(i,j)
!^
ad_FE(i,j)=ad_FE(i,j)*vmask(i,j)
# endif
!^ tl_FE(i,j)=cff*pnom_v(i,j)*(h(i,j)+h(i,j-1))* &
!^ & (tl_r2d_in(i,j)-tl_r2d_in(i,j-1))
!^
adfac=cff*pnom_v(i,j)*(h(i,j)+h(i,j-1))*ad_FE(i,j)
ad_r2d_in(i,j-1)=ad_r2d_in(i,j-1)-adfac
ad_r2d_in(i,j )=ad_r2d_in(i,j )+adfac
ad_FE(i,j)=0.0_r8
END DO
END DO
DO j=Jstr,Jend
DO i=Istr,Iend+1
# ifdef MASKING
!^ tl_FX(i,j)=tl_FX(i,j)*umask(i,j)
!^
ad_FX(i,j)=ad_FX(i,j)*umask(i,j)
# endif
!^ tl_FX(i,j)=cff*pmon_u(i,j)*(h(i,j)+h(i-1,j))* &
!^ & (tl_r2d_in(i,j)-tl_r2d_in(i-1,j))
!^
adfac=cff*pmon_u(i,j)*(h(i,j)+h(i-1,j))*ad_FX(i,j)
ad_r2d_in(i-1,j)=ad_r2d_in(i-1,j)-adfac
ad_r2d_in(i ,j)=ad_r2d_in(i ,j)+adfac
ad_FX(i,j)=0.0_r8
END DO
END DO
RETURN
END SUBROUTINE ad_r2d_oper
!
!***********************************************************************
SUBROUTINE ad_r2d_bc (ng, tile, &
& LBi, UBi, LBj, UBj, &
& ad_A)
!***********************************************************************
!
USE mod_param
USE mod_scalars
!
USE ad_exchange_2d_mod, ONLY : ad_exchange_r2d_tile
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
integer, intent(in) :: LBi, UBi, LBj, UBj
!
# ifdef ASSUMED_SHAPE
real(r8), intent(inout) :: ad_A(LBi:,LBj:)
# else
real(r8), intent(inout) :: ad_A(LBi:UBi,LBj:UBj)
# endif
!
! Local variable declarations.
!
integer :: i, j
real(r8) :: adfac
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
! Set adjoint periodic boundary conditons.
!-----------------------------------------------------------------------
!
IF (EWperiodic(ng).or.NSperiodic(ng)) THEN
CALL ad_exchange_r2d_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& ad_A)
END IF
!
!-----------------------------------------------------------------------
! Boundary corners.
!-----------------------------------------------------------------------
!
IF (.not.(EWperiodic(ng).or.NSperiodic(ng))) THEN
IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN
!^ tl_A(Iend+1,Jend+1)=0.5_r8*(tl_A(Iend+1,Jend )+ &
!^ & tl_A(Iend ,Jend+1))
!^
adfac=0.5_r8*ad_A(Iend+1,Jend+1)
ad_A(Iend+1,Jend )=ad_A(Iend+1,Jend )+adfac
ad_A(Iend ,Jend+1)=ad_A(Iend ,Jend+1)+adfac
ad_A(Iend+1,Jend+1)=0.0_r8
END IF
IF (DOMAIN(ng)%NorthWest_Corner(tile)) THEN
!^ tl_A(Istr-1,Jend+1)=0.5_r8*(tl_A(Istr-1,Jend )+ &
!^ & tl_A(Istr ,Jend+1))
!^
adfac=0.5_r8*ad_A(Istr-1,Jend+1)
ad_A(Istr-1,Jend )=ad_A(Istr-1,Jend )+adfac
ad_A(Istr ,Jend+1)=ad_A(Istr ,Jend+1)+adfac
ad_A(Istr-1,Jend+1)=0.0_r8
END IF
IF (DOMAIN(ng)%SouthEast_Corner(tile)) THEN
!^ tl_A(Iend+1,Jstr-1)=0.5_r8*(tl_A(Iend ,Jstr-1)+ &
!^ & tl_A(Iend+1,Jstr ))
!^
adfac=0.5_r8*ad_A(Iend+1,Jstr-1)
ad_A(Iend ,Jstr-1)=ad_A(Iend ,Jstr-1)+adfac
ad_A(Iend+1,Jstr )=ad_A(Iend+1,Jstr )+adfac
ad_A(Iend+1,Jstr-1)=0.0_r8
END IF
IF (DOMAIN(ng)%SouthWest_Corner(tile)) THEN
!^ tl_A(Istr-1,Jstr-1)=0.5_r8*(tl_A(Istr ,Jstr-1)+ &
!^ & tl_A(Istr-1,Jstr ))
!^
adfac=0.5_r8*ad_A(Istr-1,Jstr-1)
ad_A(Istr ,Jstr-1)=ad_A(Istr ,Jstr-1)+adfac
ad_A(Istr-1,Jstr )=ad_A(Istr-1,Jstr )+adfac
ad_A(Istr-1,Jstr-1)=0.0_r8
END IF
END IF
!
!-----------------------------------------------------------------------
! Adjoint North-South gradient boundary conditions.
!-----------------------------------------------------------------------
!
IF (.not.NSperiodic(ng)) THEN
IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
DO i=Istr,Iend
# ifdef R2D_GRADIENT
!^ tl_A(i,Jstr-1)=tl_A(i,Jstr)
!^
ad_A(i,Jstr )=ad_A(i,Jstr)+ad_A(i,Jstr-1)
ad_A(i,Jstr-1)=0.0_r8
# else
!^ tl_A(i,Jstr-1)=0.0_r8
!^
ad_A(i,Jstr-1)=0.0_r8
# endif
END DO
END IF
IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
DO i=Istr,Iend
# ifdef R2D_GRADIENT
!^ tl_A(i,Jend+1)=tl_A(i,Jend)
!^
ad_A(i,Jend )=ad_A(i,Jend)+ad_A(i,Jend+1)
ad_A(i,Jend+1)=0.0_r8
# else
!^ tl_A(i,Jend+1)=0.0_r8
!^
ad_A(i,Jend+1)=0.0_r8
# endif
END DO
END IF
END IF
!
!-----------------------------------------------------------------------
! Adjoint East-West gradient boundary conditions.
!-----------------------------------------------------------------------
!
IF (.not.EWperiodic(ng)) THEN
IF (DOMAIN(ng)%Western_Edge(tile)) THEN
DO j=Jstr,Jend
# ifdef R2D_GRADIENT
!^ tl_A(Istr-1,j)=tl_A(Istr,j)
!^
ad_A(Istr ,j)=ad_A(Istr,j)+ad_A(Istr-1,j)
ad_A(Istr-1,j)=0.0_r8
# else
!^ tl_A(Istr-1,j)=0.0_r8
!^
ad_A(Istr-1,j)=0.0_r8
# endif
END DO
END IF
IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
DO j=Jstr,Jend
# ifdef R2D_GRADIENT
!^ tl_A(Iend+1,j)=tl_A(Iend,j)
!^
ad_A(Iend ,j)=ad_A(Iend,j)+ad_A(Iend+1,j)
ad_A(Iend+1,j)=0.0_r8
# else
!^ tl_A(Iend+1,j)=0.0_r8
!^
ad_A(Iend+1,j)=0.0_r8
# endif
END DO
END IF
END IF
RETURN
END SUBROUTINE ad_r2d_bc
!
!***********************************************************************
SUBROUTINE ad_r2d_dotp (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& ad_DotProd, &
# ifdef MASKING
& rmask, &
# endif
& s1_zeta, s2_zeta, &
& ad_s1_zeta, ad_s2_zeta)
!***********************************************************************
!
USE mod_param
USE mod_parallel
USE mod_ncparam
# ifdef DISTRIBUTE
!
USE distribute_mod, ONLY : mp_reduce
# endif
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, model
integer, intent(in) :: LBi, UBi, LBj, UBj
!
# ifdef ASSUMED_SHAPE
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:,LBj:)
# endif
real(r8), intent(in) :: s1_zeta(LBi:,LBj:)
real(r8), intent(in) :: s2_zeta(LBi:,LBj:)
real(r8), intent(inout) :: ad_s1_zeta(LBi:,LBj:)
real(r8), intent(inout) :: ad_s2_zeta(LBi:,LBj:)
# else
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: s1_zeta(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: s2_zeta(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: ad_s1_zeta(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: ad_s2_zeta(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(inout) :: ad_DotProd
!
! Local variable declarations.
!
integer :: NSUB, i, j
real(r8) :: ad_cff
real(r8) :: ad_my_DotProd
# ifdef DISTRIBUTE
character (len=3) :: op_handle
# endif
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
! Compute adjoint dot product between ad_s1_zeta and ad_s2_zeta.
!-----------------------------------------------------------------------
!
ad_my_DotProd=ad_DotProd
ad_DotProd=0.0_r8
ad_cff=0.0_r8
!
DO j=Jstr,Jend
DO i=Istr,Iend
!^ tl_my_DotProd=tl_my_DotProd+tl_cff
!^
ad_cff=ad_cff+ad_my_DotProd
# ifdef MASKING
!^ tl_cff=tl_cff*rmask(i,j)
!^
ad_cff=ad_cff*rmask(i,j)
# endif
!^ tl_cff=tl_s1_zeta(i,j)*s2_zeta(i,j)+ &
!^ & s1_zeta(i,j)*tl_s2_zeta(i,j)
!^
ad_s1_zeta(i,j)=ad_s1_zeta(i,j)+ad_cff*s2_zeta(i,j)
ad_s2_zeta(i,j)=ad_s2_zeta(i,j)+ad_cff*s1_zeta(i,j)
ad_cff=0.0_r8
END DO
END DO
!^ tl_my_DotProd=0.0_r8
!^
ad_my_DotProd=0.0_r8
RETURN
END SUBROUTINE ad_r2d_dotp
#endif
END MODULE zeta_balance_mod