#include "GMREDI_OPTIONS.h"

CBOP
C     !ROUTINE: SUBMESO_CALC_PSI
C     !INTERFACE:
      SUBROUTINE SUBMESO_CALC_PSI(
     I             bi, bj, iMin, iMax, jMin, jMax,
     I             sigmaX, sigmaY, sigmaR,
     I             locMixLayer,
     I             myIter, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SUBMESO_CALC_PSI
C     | o Calculate stream-functions for Sub-Meso bolus velocity
C     *==========================================================*
C     | Ref: B. Fox-Kemper etal, Oce.Model., 39:61-78, 2011
C     |      B. Fox-Kemper etal, JPO, 38(6):1145-1165, 2008
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

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

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
      _RL sigmaX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL sigmaY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL locMixLayer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER myIter
      INTEGER myThid
CEOP

#ifndef GM_EXCLUDE_SUBMESO

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER i,j,k
      _RL mixLayerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL mixLayerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL dBuoyX_Hu   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL dBuoyY_Hv   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL NHmixLay    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL MsquareH    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL lengthScaleF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL fcorLoc     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL PsiLoc      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL dzLoc
#ifdef GM_BOLUS_ADVEC
      _RL z2H, mu_z
#endif
      _RL five_ov21
      PARAMETER( five_ov21 = 5. _d 0 / 21. _d 0 )

C--   parameter to move to GMREDI.h
c     _RL subMeso_invTau, subMeso_LfMin, subMeso_Ceff
c     _RS subMeso_Lmax

c     subMeso_invTau = 1.6 _d -6  ! ~ 1/(7.2 days)
c     subMeso_LfMin  = 1000. _d 0
c     subMeso_Ceff   = 0.07 _d 0
c     subMeso_Lmax   = 111. _d 3

C-    Initialization : <= done in S/R gmredi_init

c     IF ( GM_useSubMeso ) THEN
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx+1,sNx+OLx
         mixLayerU(i,j) = op5*( locMixLayer(i-1,j)+locMixLayer(i,j) )
         mixLayerU(i,j) = MIN( mixLayerU(i,j), -rLowW(i,j,bi,bj) )
        ENDDO
       ENDDO
       DO j=1-OLy+1,sNy+OLy
        DO i=1-OLx,sNx+OLx
         mixLayerV(i,j)=op5*( locMixLayer(i,j-1)+locMixLayer(i,j) )
         mixLayerV(i,j) = MIN( mixLayerV(i,j), -rLowS(i,j,bi,bj) )
        ENDDO
       ENDDO

C--    Integrate buoyancy gradient over the Mixed-Layer
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
          dBuoyX_Hu(i,j)= 0.
          dBuoyY_Hv(i,j)= 0.
          NHmixLay(i,j) = 0.
          fcorLoc(i,j) = SQRT( fCori(i,j,bi,bj)*fCori(i,j,bi,bj)
     &                       + subMeso_invTau*subMeso_invTau )
        ENDDO
       ENDDO
       DO k=1,Nr
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx+1,sNx+OLx
          dzLoc = MAX( 0. _d 0, MIN( drF(k), mixLayerU(i,j)+rF(k) ) )
          dBuoyX_Hu(i,j) = dBuoyX_Hu(i,j) + sigmaX(i,j,k)*dzLoc
         ENDDO
        ENDDO
        DO j=1-OLy+1,sNy+OLy
         DO i=1-OLx,sNx+OLx
          dzLoc = MAX( 0. _d 0, MIN( drF(k), mixLayerV(i,j)+rF(k) ) )
          dBuoyY_Hv(i,j) = dBuoyY_Hv(i,j) + sigmaY(i,j,k)*dzLoc
         ENDDO
        ENDDO
       ENDDO
       DO k=2,Nr
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          dzLoc = 0.
          IF ( locMixLayer(i,j)+rC(k-1).GE.0. ) dzLoc = drC(k)
          NHmixLay(i,j) = NHmixLay(i,j)
     &                  + dzLoc*MAX( -sigmaR(i,j,k), 0. _d 0 )
         ENDDO
        ENDDO
       ENDDO
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
          dBuoyX_Hu(i,j)= -dBuoyX_Hu(i,j)*gravity*recip_rhoConst
          dBuoyY_Hv(i,j)= -dBuoyY_Hv(i,j)*gravity*recip_rhoConst
          NHmixLay(i,j) = SQRT( NHmixLay(i,j)*gravity*recip_rhoConst
     &                         *locMixLayer(i,j) )
        ENDDO
       ENDDO
       DO j=2-OLy,sNy+OLy-1
        DO i=2-OLx,sNx+OLx-1
          MsquareH(i,j)= SQRT( op25*(
     &            (dBuoyX_Hu(i,j) + dBuoyX_Hu(i+1,j))**2
     &          + (dBuoyY_Hv(i,j) + dBuoyY_Hv(i,j+1))**2
     &                     ) )
        ENDDO
       ENDDO
C-     Compute Lf at grid-cell center
       DO j=2-OLy,sNy+OLy-1
        DO i=2-OLx,sNx+OLx-1
          lengthScaleF(i,j)= MAX(
     &        MsquareH(i,j)/(fcorLoc(i,j)*fcorLoc(i,j)) ,
     &        NHmixLay(i,j)/fcorLoc(i,j) ,
     &        subMeso_LfMin )
        ENDDO
       ENDDO

C      Mix-Layer Eddies contribution to Bolus Transport in X dir.
       DO j=2-OLy,sNy+OLy-1
        DO i=3-OLx,sNx+OLx-1
         PsiLoc(i,j) = -subMeso_Ceff*dBuoyX_Hu(i,j)
     &                 *mixLayerU(i,j)
     &                 *MIN( dxC(i,j,bi,bj), subMeso_Lmax )
     &                 *2. _d 0/(lengthScaleF(i-1,j)+lengthScaleF(i,j))
     &                 *2. _d 0/(fcorLoc(i-1,j)+fcorLoc(i,j))
        ENDDO
       ENDDO
#ifdef GM_BOLUS_ADVEC
       DO k=2,Nr
        DO j=2-OLy,sNy+OLy-1
         DO i=3-OLx,sNx+OLx-1
          IF ( mixLayerU(i,j).GT.0. _d 0 ) THEN
            z2H = 2. _d 0*rF(k)/mixLayerU(i,j)
          ELSE
            z2H = 0. _d 0
          ENDIF
          mu_z = ( z2H + 1. _d 0 )*( z2H + 1. _d 0 )
          mu_z = ( 1. _d 0 - mu_z )*(1. _d 0 + mu_z*five_ov21 )
          mu_z = MAX( 0. _d 0, mu_z )
          GM_PsiX(i,j,k,bi,bj) = GM_PsiX(i,j,k,bi,bj)
     &                         + mu_z*PsiLoc(i,j)
         ENDDO
        ENDDO
       ENDDO
#endif /* GM_BOLUS_ADVEC */
#ifdef ALLOW_DIAGNOSTICS
       IF ( useDiagnostics ) THEN
         CALL DIAGNOSTICS_FILL( lengthScaleF, 'SubMesLf',
     &                           0, 1, 2, bi, bj, myThid )
         CALL DIAGNOSTICS_FILL( PsiLoc, 'SubMpsiX',
     &                          0, 1, 2, bi, bj, myThid )
       ENDIF
#endif
       IF ( debugLevel.GE.debLevD ) THEN
         CALL WRITE_LOCAL_RL( 'subMeso_Lf','I10',1,lengthScaleF,
     &                         bi,bj,1,myIter,myThid )
         CALL WRITE_LOCAL_RL( 'subMeso_psiX','I10',1,PsiLoc,
     &                         bi,bj,1,myIter,myThid )
       ENDIF

C      Mix-Layer Eddies contribution to Bolus Transport in Y dir.
       DO j=3-OLy,sNy+OLy-1
        DO i=2-OLx,sNx+OLx-1
         PsiLoc(i,j) = -subMeso_Ceff*dBuoyY_Hv(i,j)
     &                 *mixLayerV(i,j)
     &                 *MIN( dyC(i,j,bi,bj), subMeso_Lmax )
     &                 *2. _d 0/(lengthScaleF(i,j-1)+lengthScaleF(i,j))
     &                 *2. _d 0/(fcorLoc(i,j-1)+fcorLoc(i,j))
        ENDDO
       ENDDO
#ifdef GM_BOLUS_ADVEC
       DO k=2,Nr
        DO j=3-OLy,sNy+OLy-1
         DO i=2-OLx,sNx+OLx-1
          IF ( mixLayerV(i,j).GT.0. _d 0 ) THEN
            z2H = 2. _d 0*rF(k)/mixLayerV(i,j)
          ELSE
            z2H = 0. _d 0
          ENDIF
          mu_z = ( z2H + 1. _d 0 )*( z2H + 1. _d 0 )
          mu_z = ( 1. _d 0 - mu_z )*(1. _d 0 + mu_z*five_ov21 )
          mu_z = MAX( 0. _d 0, mu_z )
          GM_PsiY(i,j,k,bi,bj) = GM_PsiY(i,j,k,bi,bj)
     &                         + mu_z*PsiLoc(i,j)
         ENDDO
        ENDDO
       ENDDO
#endif /* GM_BOLUS_ADVEC */
#ifdef ALLOW_DIAGNOSTICS
       IF ( useDiagnostics ) THEN
         CALL DIAGNOSTICS_FILL( PsiLoc, 'SubMpsiY',
     &                          0, 1, 2, bi, bj, myThid )
       ENDIF
#endif
       IF ( debugLevel.GE.debLevD ) THEN
         CALL WRITE_LOCAL_RL( 'subMeso_psiY','I10',1,PsiLoc,
     &                         bi,bj,1,myIter,myThid )
       ENDIF

c     ENDIF
#endif /* ndef GM_EXCLUDE_SUBMESO */

      RETURN
      END
