#include "cppdefs.h"
      MODULE zeta_balance_mod

#if defined BALANCE_OPERATOR  && defined ZETA_ELLIPTIC
!svn $Id: zeta_balance.F 995 2020-01-10 04:01:28Z arango $
!=================================================== Andrew M. Moore ===
!  Copyright (c) 2002-2020 The ROMS/TOMS Group      Hernan G. Arango   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.txt                                              !
!=======================================================================
!                                                                      !
!  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,3)
# 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,3)
      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
