!======================================================================
! CROCO is a branch of ROMS developped at IRD, INRIA, 
! Ifremer, CNRS and Univ. Toulouse III  in France
! The two other branches from UCLA (Shchepetkin et al)
! and Rutgers University (Arango et al) are under MIT/X style license.
! CROCO specific routines (nesting) are under CeCILL-C license.
!
! CROCO website : http://www.croco-ocean.org
!======================================================================
!
#include "cppdefs.h"
#ifndef CHILD_SPG
      subroutine uv3dmix (tile)
!====================================================================!
!                                                                    !
!  This subroutine computes biharmonic mixing of momentum, along     !
!  constant  S-surfaces, from the  horizontal divergence  of the     !
!  stress tensor. A transverse isotropy is assumed so the stress     !
!  tensor is splitted into vertical and horizontal subtensors.       !
!  Components of the stress tensor are:                              !
!                                                                    !
!                         du     dv                                  !
!         s_xx = -s_yy = ---- - -----                                !
!                         dx     dy                                  !
!                                                                    !
!                         dv     du                                  !
!         s_xy =  s_yx = ---- + ----                                 !
!                         dx     dy                                  !
!                                                                    !
!  References:                                                       !
!                                                                    !
!      Wajsowicz, R.C, 1993: A consistent formulation of the         !
!         anisotropic stress tensor for use in models of the         !
!         large-scale ocean circulation, JCP, 105, 333-338.          !
!                                                                    !
!      Sadourny, R. and K. Maynard, 1997: Formulations of            !
!         lateral diffusion in geophysical fluid dynamics            !
!         models, In "Numerical Methods of Atmospheric and           !
!         Oceanic Modelling". Lin, Laprise, and Ritchie,             !
!         Eds., NRC Research Press, 547-556.                         !
!                                                                    !
!      Griffies, S.M. and R.W. Hallberg, 2000: Biharmonic            !
!         friction with a Smagorinsky-like viscosity for             !
!         use in large-scale eddy-permitting ocean models,           !
!         Monthly Weather Rev., 128, 8, 2935-2946.                   !
!                                                                    !
!=====================================================================
!
      implicit none
      integer tile, trd, omp_get_thread_num
# include "param.h"
# include "private_scratch.h"
# include "compute_tile_bounds.h"
      trd=omp_get_thread_num()
!
# ifdef MASKING
#  define SWITCH *
# else
#  define SWITCH !
# endif
!
# ifdef AGRIF
      if (AGRIF_Root()) then
      call uv3dmix_tile (Istr,Iend,Jstr,Jend,
     &                      A2d(1, 1,trd),A2d(1, 3,trd),A2d(1, 5,trd),
     &                      A2d(1, 7,trd),A2d(1, 9,trd),A2d(1,11,trd))
      else
      call uv3dmix_child_tile (Istr,Iend,Jstr,Jend,
     &                      A2d(1, 1,trd),A2d(1, 3,trd),A2d(1, 5,trd),
     &                      A2d(1, 7,trd),A2d(1, 9,trd),A2d(1,11,trd))


       endif
      return
      end
# else
      call uv3dmix_tile (Istr,Iend,Jstr,Jend,
     &                      A2d(1, 1,trd),A2d(1, 3,trd),A2d(1, 5,trd),
     &                      A2d(1, 7,trd),A2d(1, 9,trd),A2d(1,11,trd))
      return
      end
# endif   /* AGRIF */
!
!---------------------------------------------------------------------
!*********************************************************************
!---------------------------------------------------------------------
!
!PARENT
!
      subroutine uv3dmix_tile (Istr,Iend,Jstr,Jend,UFx,UFe,LapU,
     &                            VFx,VFe,LapV)
# undef CLIMAT_UV_MIXH_FINE
!
#else
!
! CHILD
!
      subroutine uv3dmix_child_tile (Istr,Iend,Jstr,Jend,UFx,UFe,LapU,
     &                            VFx,VFe,LapV)
!
!Diffusion always applied on U-UCLM in fine grid
# if !defined UV_HADV_RSUP3
#  define CLIMAT_UV_MIXH_FINE
# endif
!
#endif /* CHILD_SPG */
!
!---------------------------------------------------------------------
! ******************************Common Code***************************
!---------------------------------------------------------------------
!
      implicit none
#include "param.h"
#include "grid.h"
#include "mixing.h"
#include "coupling.h"
#include "ocean3d.h"
#include "scalars.h"
#ifdef CLIMAT_UV_MIXH
# include "climat.h"
#endif
#ifdef DIAGNOSTICS_UV
# include "diagnostics.h"
#else
# if defined DIAGNOSTICS_VRT
#  include "diags_vrt.h"
# endif
# if defined DIAGNOSTICS_EK
#  include "diags_ek.h"
# endif
# if defined DIAGNOSTICS_PV
#  include "diags_pv.h"
# endif
#endif
!
      integer  Iend, Istr, Jend, Jstr, i, j, k, indx,
     &         imin,imax,jmin,jmax, IminU,JminV
      real     cff,cff1,
     &         LapU(PRIVATE_2D_SCRATCH_ARRAY),
     &         LapV(PRIVATE_2D_SCRATCH_ARRAY),
     &          UFe(PRIVATE_2D_SCRATCH_ARRAY),
     &          UFx(PRIVATE_2D_SCRATCH_ARRAY),
     &          VFe(PRIVATE_2D_SCRATCH_ARRAY),
     &          VFx(PRIVATE_2D_SCRATCH_ARRAY)
!
#ifdef AGRIF
# include "zoom.h"
#endif
!
#include "compute_auxiliary_bounds.h"
!
#ifdef CHILD_SPG
# define UCLM usponge
# define VCLM vsponge
#else
# define UCLM uclm
# define VCLM vclm
#endif
!
#ifndef EW_PERIODIC
      if (WESTERN_EDGE) then
        imin=istr
        iminU=istr+1
      else
        imin=istr-1
        iminU=imin
      endif
      if (EASTERN_EDGE) then
        imax=iend
      else
        imax=iend+1
      endif
#else
      imin=istr-1
      imax=iend+1
      iminU=imin
#endif
#ifndef NS_PERIODIC
      if (SOUTHERN_EDGE) then
        jmin=jstr
        jminV=jstr+1
      else
        jmin=jstr-1
        jminV=jmin
      endif
      if (NORTHERN_EDGE) then
        jmax=jend
      else
        jmax=jend+1
      endif
#else
      jmin=jstr-1
      jmax=jend+1
      jminV=jmin
#endif

      indx=3-nstp    !--> time index for target arrays;

      do k=1,N
!
!---------------------------------------------------------------------
!  Compute horizontal biharmonic viscosity along constant S-surfaces.
!  The biharmonic operator is computed by applying the harmonic
!  operator twice.
!---------------------------------------------------------------------
!
!  Compute flux-components of the horizontal divergence of the stress
!  tensor (m4 s^-3/2) in XI- and ETA-directions. For momentum balance
!  purposes, the thickness "Hz" appears only when computing the second
!  harmonic operator.
!                         du     dv
!         s_xx = -s_yy = ---- - -----
!                         dx     dy
!
        do j=jminV-1,jmax
          do i=iminU-1,imax
            cff=pmon_r(i,j)*
#if defined CLIMAT_UV_MIXH || defined CLIMAT_UV_MIXH_FINE
     &           (pn_u(i+1,j)*(u(i+1,j,k,nstp)-UCLM(i+1,j,k))-
     &            pn_u(i  ,j)*(u(i  ,j,k,nstp)-UCLM(i  ,j,k)) ) -
#else
     &           (pn_u(i+1,j)*u(i+1,j,k,nstp)-
     &            pn_u(i  ,j)*u(i  ,j,k,nstp) ) -
#endif
     &          pnom_r(i,j)*
#if defined CLIMAT_UV_MIXH || defined CLIMAT_UV_MIXH_FINE
     &           (pm_v(i,j+1)*(v(i,j+1,k,nstp)-VCLM(i,j+1,k))-
     &            pm_v(i,j  )*(v(i,j  ,k,nstp)-VCLM(i,j  ,k)) )
#else
     &           (pm_v(i,j+1)*v(i,j+1,k,nstp)-
     &            pm_v(i,j  )*v(i,j  ,k,nstp))
#endif

            UFx(i,j)=on_r(i,j)*on_r(i,j)*cff
            VFe(i,j)=om_r(i,j)*om_r(i,j)*cff
          enddo
        enddo
!                         dv     du
!         s_xy =  s_yx = ---- + ----
!                         dx     dy
        do j=jmin,jmax+1
          do i=imin,imax+1
            cff=pmon_p(i,j)*
#if defined CLIMAT_UV_MIXH || defined CLIMAT_UV_MIXH_FINE
     &           (pn_v(i  ,j)*(v(i  ,j,k,nstp)-VCLM(i  ,j,k)) -
     &            pn_v(i-1,j)*(v(i-1,j,k,nstp)-VCLM(i-1,j,k)) )+
#else
     &           (pn_v(i  ,j)*v(i  ,j,k,nstp)-
     &            pn_v(i-1,j)*v(i-1,j,k,nstp))+
#endif
     &          pnom_p(i,j)*
#if defined CLIMAT_UV_MIXH || defined CLIMAT_UV_MIXH_FINE
     &           (pm_u(i,j  )*( u(i,j  ,k,nstp)-UCLM(i,j  ,k))-
     &            pm_u(i,j-1)*( u(i,j-1,k,nstp)-UCLM(i,j-1,k)) )
#else
     &           (pm_u(i,j  )*u(i,j  ,k,nstp)-
     &            pm_u(i,j-1)*u(i,j-1,k,nstp))
#endif
#ifdef MASKING
     &                               * pmask(i,j)
#endif

            UFe(i,j)=om_p(i,j)*om_p(i,j)*cff
            VFx(i,j)=on_p(i,j)*on_p(i,j)*cff
          enddo
        enddo
!
!  Compute first harmonic operator (m s^-3/2).
!
        do j=jmin,jmax
          do i=iminU,imax
            LapU(i,j)=0.125*
     &                (pm(i-1,j)+pm(i,j))*(pn(i-1,j)+pn(i,j))*
     &                ((pn(i-1,j)+pn(i,j))*(UFx(i,j  )-UFx(i-1,j))+
     &                 (pm(i-1,j)+pm(i,j))*(UFe(i,j+1)-UFe(i  ,j)))
          enddo
        enddo
        do j=jminV,jmax
          do i=imin,imax
            LapV(i,j)=0.125*
     &                (pm(i,j)+pm(i,j-1))*(pn(i,j)+pn(i,j-1))*
     &                ((pn(i,j-1)+pn(i,j))*(VFx(i+1,j)-VFx(i,j  ))-
     &                 (pm(i,j-1)+pm(i,j))*(VFe(i  ,j)-VFe(i,j-1)))
          enddo
        enddo
!
!  Apply boundary conditions (other than periodic) to the first
!  harmonic operator. These are gradient or closed (free slip or
!  no slip) boundary conditions.
!
#ifndef EW_PERIODIC
        if (WESTERN_EDGE) then
          do j=jmin,jmax
# ifndef OBC_WEST
            LapU(Istr,j)=0.0
# else
            LapU(Istr,j)=LapU(Istr+1,j)
# endif
          enddo
          do j=jminV,jmax
# ifndef OBC_WEST
            LapV(Istr-1,j)=gamma2*LapV(Istr,j)
# else
            LapV(Istr-1,j)=0.0
# endif
          enddo
        endif
        if (EASTERN_EDGE) then
          do j=jmin,jmax
# ifndef OBC_EAST
            LapU(Iend+1,j)=0.0
# else
            LapU(Iend+1,j)=LapU(Iend,j)
# endif
          enddo
          do j=jminV,jmax
# ifndef OBC_EAST
            LapV(Iend+1,j)=gamma2*LapV(Iend,j)
# else
            LapV(Iend+1,j)=0.0
# endif
          enddo
        endif
#endif /* !EW_PERIODIC */
#ifndef NS_PERIODIC
        if (SOUTHERN_EDGE) then
          do i=iminU,imax
# ifndef OBC_SOUTH
            LapU(i,Jstr-1)=gamma2*LapU(i,Jstr)
# else
            LapU(i,Jstr-1)=0.0
# endif
          enddo
          do i=imin,imax
# ifndef OBC_SOUTH
            LapV(i,Jstr)=0.0
# else
            LapV(i,Jstr)=LapV(i,Jstr+1)
# endif
          enddo
        endif
        if (NORTHERN_EDGE) then
          do i=iminU,imax
# ifndef OBC_NORTH
            LapU(i,Jend+1)=gamma2*LapU(i,Jend)
# else
            LapU(i,Jend+1)=0.0
# endif
          enddo
          do i=imin,imax
# ifndef OBC_NORTH
            LapV(i,Jend+1)=0.0
# else
            LapV(i,Jend+1)=LapV(i,Jend)
# endif
          enddo
        endif
#endif /* !NS_PERIODIC */
!
! Corner points
!
#if !defined EW_PERIODIC && !defined NS_PERIODIC
        IF ((SOUTHERN_EDGE).and.(WESTERN_EDGE)) THEN
          LapU(Istr  ,Jstr-1)=0.5*(LapU(Istr+1,Jstr-1)+
     &                             LapU(Istr  ,Jstr  ))
          LapV(Istr-1,Jstr  )=0.5*(LapV(Istr-1,Jstr+1)+
     &                             LapV(Istr  ,Jstr  ))
        END IF
        IF ((SOUTHERN_EDGE).and.(EASTERN_EDGE)) THEN
          LapU(Iend+1,Jstr-1)=0.5*(LapU(Iend  ,Jstr-1)+
     &                             LapU(Iend+1,Jstr  ))
          LapV(Iend+1,Jstr  )=0.5*(LapV(Iend  ,Jstr  )+
     &                             LapV(Iend+1,Jstr+1))
        END IF
        IF ((NORTHERN_EDGE).and.(WESTERN_EDGE)) THEN
          LapU(Istr  ,Jend+1)=0.5*(LapU(Istr+1,Jend+1)+
     &                             LapU(Istr  ,Jend  ))
          LapV(Istr-1,Jend+1)=0.5*(LapV(Istr  ,Jend+1)+
     &                             LapV(Istr-1,Jend  ))
        END IF
        IF ((NORTHERN_EDGE).and.(EASTERN_EDGE)) THEN
          LapU(Iend+1,Jend+1)=0.5*(LapU(Iend  ,Jend+1)+
     &                             LapU(Iend+1,Jend  ))
          LapV(Iend+1,Jend+1)=0.5*(LapV(Iend  ,Jend+1)+
     &                             LapV(Iend+1,Jend  ))
        END IF
#endif
!
!---------------------------------------------------------------------
!  Compute horizontal gradients associated with the
!  second harmonic operator.
!---------------------------------------------------------------------
!
!  Compute flux-components of the horizontal divergence of the
!  harmonic stress tensor (m4/s2) in XI- and ETA-directions.
!
        do j=JstrV-1,Jend
          do i=IstrU-1,Iend
            cff=Hz(i,j,k)*
     &          (pmon_r(i,j)*
     &           (pn_u(i+1,j)*LapU(i+1,j)-
     &            pn_u(i  ,j)*LapU(i  ,j))-
     &           pnom_r(i,j)*
     &           (pm_v(i,j+1)*LapV(i,j+1)-
     &            pm_v(i,j  )*LapV(i,j  )))
            UFx(i,j)=on_r(i,j)*on_r(i,j)*cff
            VFe(i,j)=om_r(i,j)*om_r(i,j)*cff
#ifdef VIS_COEF_3D
# ifdef UV_HADV_RSUP3
            UFx(i,j)=viscU_r(i,j,k)*UFx(i,j)
            VFe(i,j)=viscV_r(i,j,k)*VFe(i,j)
# else
            UFx(i,j)=visc3d_r(i,j,k)*UFx(i,j)
            VFe(i,j)=visc3d_r(i,j,k)*VFe(i,j)
# endif
#else
            UFx(i,j)=visc4_r(i,j)*UFx(i,j)
            VFe(i,j)=visc4_r(i,j)*VFe(i,j)
#endif
          enddo
        enddo
        do j=Jstr,Jend+1
          do i=Istr,Iend+1
            cff=0.25*(Hz(i-1,j  ,k)+Hz(i,j  ,k)+
     &                Hz(i-1,j-1,k)+Hz(i,j-1,k))*
     &          (pmon_p(i,j)*
     &           (pn_v(i  ,j)*LapV(i  ,j)-
     &            pn_v(i-1,j)*LapV(i-1,j))+
     &           pnom_p(i,j)*
     &           (pm_u(i,j  )*LapU(i,j  )-
     &            pm_u(i,j-1)*LapU(i,j-1)))
#ifdef MASKING
     &                             * pmask(i,j)
#endif

            UFe(i,j)=om_p(i,j)*om_p(i,j)*cff
            VFx(i,j)=on_p(i,j)*on_p(i,j)*cff
#ifdef VIS_COEF_3D
# ifdef UV_HADV_RSUP3
            UFe(i,j)=viscU_p(i,j,k)*UFe(i,j)
            VFx(i,j)=viscV_p(i,j,k)*VFx(i,j)
# else
            UFe(i,j)=visc3d_p(i,j,k)*UFe(i,j)
            VFx(i,j)=visc3d_p(i,j,k)*VFx(i,j)
# endif
#else
            UFe(i,j)=visc4_p(i,j)*UFe(i,j)
            VFx(i,j)=visc4_p(i,j)*VFx(i,j)
#endif
          enddo
        enddo
!
! Time-step biharmonic, S-surfaces viscosity term.  Notice that
! momentum at this stage is HzU and HzV and has units m2/s.
! Also compute vertical integral of viscous terms and add it
! into coupling terms for the barotropic mode
!
        do j=Jstr,Jend
          do i=IstrU,Iend
            cff=pn_u(i,j)*(UFx(i,j)-UFx(i-1,j))
     &         +pm_u(i,j)*(UFe(i,j+1)-UFe(i,j))
            cff1=pm_u(i,j)*pn_u(i,j)*cff
            rufrc(i,j)=rufrc(i,j) - cff
            u(i,j,k,indx)=u(i,j,k,indx) - dt*cff1
#if defined DIAGNOSTICS_UV || defined DIAGNOSTICS_EK_FULL || defined DIAGNOSTICS_PV
!            if (nnew.ne.3) then
            MHmix(i,j,k,1,indx) = -cff1 * om_u(i,j)*on_u(i,j)
     &                       SWITCH umask(i,j)
!            endif
#elif defined DIAGNOSTICS_EK
!            if (nnew.ne.3) then
              if (k.eq.1) then
                ekwrkHmix(i,j,1,indx) = -cff1 * om_u(i,j)*on_u(i,j)
     &                           * u(i,j,1,nrhs) SWITCH umask(i,j)
              else
                ekwrkHmix(i,j,1,indx) = ekwrkHmix(i,j,1,indx)
     &                           - cff1 * om_u(i,j)*on_u(i,j)
     &                           * u(i,j,1,nrhs) SWITCH umask(i,j)
              endif
!            endif
#endif
#if defined DIAGNOSTICS_VRT && ! defined DIAGNOSTICS_UV
!            if (nnew.ne.3) then
              if (k.eq.1) then
                wrkHmix(i,j,1,indx) = -cff1 * om_u(i,j)*on_u(i,j)
     &                                    SWITCH umask(i,j)
              else
                wrkHmix(i,j,1,indx) = wrkHmix(i,j,1,indx)
     &                          -cff1 * om_u(i,j)*on_u(i,j)
     &                                    SWITCH umask(i,j)
              endif
!            endif
#endif
          enddo
        enddo
        do j=JstrV,Jend
          do i=Istr,Iend
            cff=pn_v(i,j)*(VFx(i+1,j)-VFx(i,j))
     &         -pm_v(i,j)*(VFe(i,j)-VFe(i,j-1))
            cff1=pm_v(i,j)*pn_v(i,j)*cff
            rvfrc(i,j)=rvfrc(i,j) - cff
            v(i,j,k,indx)=v(i,j,k,indx) - dt*cff1
#if defined DIAGNOSTICS_UV || defined DIAGNOSTICS_EK_FULL || defined DIAGNOSTICS_PV
!            if (nnew.ne.3) then
            MHmix(i,j,k,2,indx) =  -cff1 * om_v(i,j)*on_v(i,j)
     &                       SWITCH vmask(i,j)
!            endif
#elif defined DIAGNOSTICS_EK
!            if (nnew.ne.3) then
              if (k.eq.1) then
                ekwrkHmix(i,j,2,indx) = -cff1 * om_v(i,j)*on_v(i,j)
     &                          * v(i,j,1,nrhs) SWITCH vmask(i,j)
              else
                ekwrkHmix(i,j,2,indx) = ekwrkHmix(i,j,2,indx)
     &                           -cff1 * om_v(i,j)*on_v(i,j)
     &                          * v(i,j,1,nrhs) SWITCH vmask(i,j)
              endif
!            endif
#endif
#if defined DIAGNOSTICS_VRT && ! defined DIAGNOSTICS_UV
!            if (nnew.ne.3) then
              if (k.eq.1) then
                wrkHmix(i,j,2,indx) = -cff1 * om_v(i,j)*on_v(i,j)
     &                                    SWITCH vmask(i,j)
              else
                wrkHmix(i,j,2,indx) = wrkHmix(i,j,2,indx)
     &                           - cff1 * om_v(i,j)*on_v(i,j)
     &                                    SWITCH vmask(i,j)
              endif
!            endif
#endif
          enddo
        enddo
      enddo
!
      return
      end
!
#ifndef CHILD_SPG
# undef UCLM VCLM
# define CHILD_SPG
# ifdef AGRIF
#  include "uv3dmix4_S.F"
# endif
# undef CHILD_SPG
#endif  /* !CHILD_SPG */



