C $Header: /u/gcmpack/MITgcm/pkg/mom_vecinv/mom_vi_del2uv.F,v 1.17 2015/09/10 18:08:51 jmc Exp $
C $Name:  $

#include "MOM_VECINV_OPTIONS.h"

      SUBROUTINE MOM_VI_DEL2UV(
     I        bi,bj,k,
     I        hDiv,vort3,hFacZ,
     O        del2u,del2v,
     I        myThid)
      IMPLICIT NONE
C
C     Calculate del^2 of (u,v) in terms of hDiv and vort3
C

C     == Global variables ==
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"

C     == Routine arguments ==
      INTEGER bi,bj,k
      _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL del2u(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL del2v(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER myThid

C     == Local variables ==
      INTEGER i,j
c     _RL Zip,Zij,Zpj,Dim,Dij,Dmj,uDij

C     - bi-harmonic viscosity :

c     DO j=2-OLy,sNy+OLy-1
c      DO i=2-OLx,sNx+OLx-1

c       Dim=dyF( i ,j-1,bi,bj)*hFacC( i ,j-1,k,bi,bj)*hDiv( i ,j-1)
c       Dij=dyF( i , j ,bi,bj)*hFacC( i , j ,k,bi,bj)*hDiv( i , j )
c       Dmj=dyF(i-1, j ,bi,bj)*hFacC(i-1, j ,k,bi,bj)*hDiv(i-1, j )
c       Dim=dyF( i ,j-1,bi,bj)*                       hDiv( i ,j-1)
c       Dij=dyF( i , j ,bi,bj)*                       hDiv( i , j )
c       Dmj=dyF(i-1, j ,bi,bj)*                       hDiv(i-1, j )
c       Dim=                                          hDiv( i ,j-1)
c       Dij=                                          hDiv( i , j )
c       Dmj=                                          hDiv(i-1, j )

c       Zip=dxV( i ,j+1,bi,bj)*hFacZ( i ,j+1)*vort3( i ,j+1)
c       Zij=dxV( i , j ,bi,bj)*hFacZ( i , j )*vort3( i , j )
c       Zpj=dxV(i+1, j ,bi,bj)*hFacZ(i+1, j )*vort3(i+1, j )
c       Zip=                   hFacZ( i ,j+1)*vort3( i ,j+1)
c       Zij=                   hFacZ( i , j )*vort3( i , j )
c       Zpj=                   hFacZ(i+1, j )*vort3(i+1, j )

c       del2u(i,j) = recip_rAw(i,j,bi,bj)*(
c    &   +recip_hFacW(i,j,k,bi,bj)*( Dij-Dmj )
c    &   -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij ) )
c       del2u(i,j) = recip_rAw(i,j,bi,bj)*(
c    &   +                         ( Dij-Dmj )
c    &   -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij ) )
c       del2u(i,j) =
c    &   +                         ( Dij-Dmj )*recip_DXC(i,j,bi,bj)
c    &   -recip_hFacW(i,j,k,bi,bj)*( Zip-Zij )*recip_DYG(i,j,bi,bj)

c       del2v(i,j) = recip_rAs(i,j,bi,bj)*(
c    &    recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )
c    &   +recip_hFacS(i,j,k,bi,bj)*( Dij-Dim ) )
c       del2v(i,j) = recip_rAs(i,j,bi,bj)*(
c    &    recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )
c    &   +                         ( Dij-Dim ) )
c       del2v(i,j) =
c    &    recip_hFacS(i,j,k,bi,bj)*( Zpj-Zij )*recip_DXG(i,j,bi,bj)
c    &   +                         ( Dij-Dim )*recip_DYC(i,j,bi,bj)

c      ENDDO
c     ENDDO

C     - bi-harmonic viscosity : del^2(U_component)

cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
      IF (useCubedSphereExchange) THEN
C     to compute d/dx(hDiv), fill corners with appropriate values:
        CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
     &                             hDiv, bi,bj, myThid )
      ENDIF
cph-exch2#endif
c     DO j=1,sNy
c      DO i=1,sNx+1
      DO j=2-OLy,sNy+OLy-1
       DO i=2-OLx,sNx+OLx-1
        del2u(i,j) =
     &      (  ( hDiv(i,j) - hDiv(i-1,j) )*recip_dxC(i,j,bi,bj)
     &        -_recip_hFacW(i,j,k,bi,bj)*
     &         ( hFacZ(i,j+1)*vort3(i,j+1) - hFacZ(i,j)*vort3(i,j) )
     &                                    *recip_dyG(i,j,bi,bj)
     &      )*maskW(i,j,k,bi,bj)*recip_deepFacC(k)
#ifdef ALLOW_OBCS
     &       *maskInW(i,j,bi,bj)
#endif
       ENDDO
      ENDDO

C     - bi-harmonic viscosity : del^2(V_component)

cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
      IF (useCubedSphereExchange) THEN
C     to compute d/dy(hDiv), fill corners with appropriate values:
        CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
     &                             hDiv, bi,bj, myThid )
      ENDIF
cph-exch2#endif
c     DO j=1,sNy+1
c      DO i=1,sNx
      DO j=2-OLy,sNy+OLy-1
       DO i=2-OLx,sNx+OLx-1
        del2v(i,j) =
     &      (  ( hDiv(i,j) - hDiv(i,j-1) )*recip_dyC(i,j,bi,bj)
     &        +_recip_hFacS(i,j,k,bi,bj)*
     &         ( hFacZ(i+1,j)*vort3(i+1,j) - hFacZ(i,j)*vort3(i,j) )
     &                                    *recip_dxG(i,j,bi,bj)
     &      )*maskS(i,j,k,bi,bj)*recip_deepFacC(k)
#ifdef ALLOW_OBCS
     &       *maskInS(i,j,bi,bj)
#endif
       ENDDO
      ENDDO

      RETURN
      END