C $Header: /u/gcmpack/MITgcm/pkg/ggl90/ggl90_idemix.F,v 1.7 2016/10/26 00:49:05 jmc Exp $
C $Name: $
#include "GGL90_OPTIONS.h"
CBOP
C !ROUTINE: GGL90_IDEMIX
C !INTERFACE: ======================================================
SUBROUTINE GGL90_IDEMIX(
I bi, bj, hFacI, recip_hFacI, sigmaR, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | S/R GGL90_IDEMIX
C |
C | IDEMIX1 model as described in
C | - Olbers, D. and Eden, C. (2013), JPO, doi:10.1175/JPO-D-12-0207.1
C | in a nutshell:
C | computes contribution of internal wave field to vertical mixing
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GGL90.h"
#include "FFIELDS.h"
#include "GRID.h"
#ifdef ALLOW_GMREDI
#include "GMREDI_OPTIONS.h"
#include "GMREDI.h"
#endif
C !INPUT/OUTPUT PARAMETERS:
C === Routine arguments ===
C bi, bj :: Current tile indices
C hFacI :: thickness factors for w-cells (interface)
C with reciprocal of hFacI = recip_hFacI
C sigmaR :: Vertical gradient of iso-neutral density
C myTime :: Current time in simulation
C myIter :: Current time-step number
C myThid :: My Thread Id number
INTEGER bi, bj
_RL hFacI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL recip_hFacI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL myTime
INTEGER myIter
INTEGER myThid
#ifdef ALLOW_GGL90_IDEMIX
C !LOCAL VARIABLES :
C === Local variables ===
INTEGER iMin ,iMax ,jMin ,jMax
INTEGER i, j, k, kp1, km1, kBottom
INTEGER errCode
_RL Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL deltaTggl90
_RL fxa,fxb,fxc,cstar
_RL dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL IDEMIX_gofx2,IDEMIX_hofx1
_RL delta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL bN0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL osborn_diff(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL c0(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL forc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL gm_forc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP
iMin = 2-OLx
iMax = sNx+OLx-1
jMin = 2-OLy
jMax = sNy+OLy-1
C set separate time step (should be deltaTtracer)
deltaTggl90 = dTtracerLev(1)
C Initialize local fields
DO k = 1, Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
Nsquare(i,j,k) = 0. _d 0
delta(i,j,k) = 0. _d 0
a3d(i,j,k) = 0. _d 0
b3d(i,j,k) = 1. _d 0
c3d(i,j,k) = 0. _d 0
osborn_diff(i,j,k) = 0. _d 0
c0(i,j,k) = 0. _d 0
forc(i,j,k) = 0. _d 0
ENDDO
ENDDO
ENDDO
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
dfx(i,j) = 0. _d 0
dfy(i,j) = 0. _d 0
bN0(i,j) = 0. _d 0
gm_forc(i,j) = 0. _d 0
ENDDO
ENDDO
c-----------------------------------------------------------------------
c allow for IW everywhere by limiting buoyancy freq.
c-----------------------------------------------------------------------
DO k=2,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
& * sigmaR(i,j,k)
fxb = max( 1. _d -6, abs( fCori(i,j,bi,bj) ))
Nsquare(i,j,k)= max( 100.*fxb*fxb, Nsquare(i,j,k) )
& *maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
ENDDO
ENDDO
ENDDO
c-----------------------------------------------------------------------
c vertically integrated N
c-----------------------------------------------------------------------
DO k=2,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
bN0(i,j)=bN0(i,j)
& +SQRT(Nsquare(i,j,k))*drC(k)*hFacI(i,j,k)
ENDDO
ENDDO
ENDDO
c-----------------------------------------------------------------------
c vertical and horizontal group velocities
c and constant for dissipation
c-----------------------------------------------------------------------
DO k=2,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
fxb = max( 1. _d -6, abs( fCori(i,j,bi,bj) ))
fxa = SQRT(Nsquare(i,j,k))/fxb
cstar = bN0(i,j)/(pi*IDEMIX_jstar)
c0(i,j,k)=max(0. _d 0,
& cstar*IDEMIX_gamma*IDEMIX_gofx2(fxa))
IDEMIX_V0(i,j,k,bi,bj)=max(0. _d 0,
& cstar*IDEMIX_gamma*IDEMIX_hofx1(fxa))
fxc = max( 1. _d 0 , fxa )
fxc = log( fxc + sqrt( fxc*fxc -1.))
IDEMIX_tau_d(i,j,k,bi,bj) = IDEMIX_mu0*fxb*fxc*
& (IDEMIX_jstar*pi/(GGL90eps+bN0(i,j)) )**2
ENDDO
ENDDO
ENDDO
IF ( IDEMIX_tau_h .GT. 0. _d 0 ) THEN
C horizontal diffusion of IW energy can become unstable for long
C time steps, so limit horizontal group velocity to satisfy simple
C CFL-like criterion:
C tau_h V0**2 *dt/dx**2 < 0.25 <=> V0 < sqrt( 0.25 * dx**2/(dt*tau_h) )
fxa = sqrt( 1. _d 0/( deltaTggl90 * IDEMIX_tau_h ) )
DO k=2,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
fxb = 0.5*min( _dxF(i,j,bi,bj), _dyF(i,j,bi,bj) )*fxa
IDEMIX_V0(i,j,k,bi,bj) = min( IDEMIX_V0(i,j,k,bi,bj), fxb )
ENDDO
ENDDO
ENDDO
ENDIF
c-----------------------------------------------------------------------
c forcing by mesoscale GM
c-----------------------------------------------------------------------
c vertically integrated forcing
#ifdef ALLOW_GMREDI
if (useGmredi) then
#ifdef GM_EG_PROGNOSTIC
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
gm_forc(i,j) = gm_forc(i,j)
& +GM_EG_diss(i,j,k,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
#else
DO k=2,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
gm_forc(i,j) = gm_forc(i,j)
& +max( 0. _d 0,Kwz(i,j,k,bi,bj)*Nsquare(i,j,k) )
& *drC(k)*hFacI(i,j,k)
ENDDO
ENDDO
ENDDO
#endif
endif
if (IDEMIX_include_GM .and. useGmredi) then
c inject locally
#ifdef GM_EG_PROGNOSTIC
DO k=2,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
forc(i,j,k) = forc(i,j,k)
& +.5 _d 0*(GM_EG_diss(i,j,k,bi,bj)+
& GM_EG_diss(i,j,k-1,bi,bj))
ENDDO
ENDDO
ENDDO
#else
DO k=2,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
forc(i,j,k) = forc(i,j,k)
& +max( 0. _d 0,Kwz(i,j,k,bi,bj)*Nsquare(i,j,k) )
ENDDO
ENDDO
ENDDO
#endif
endif
if (IDEMIX_include_GM_bottom .and. useGmredi) then
c inject at bottom box only
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
kBottom = MAX(kLowC(i,j,bi,bj),1)
forc(i,j,kbottom)=forc(i,j,kbottom)
& + gm_forc(i,j)*recip_drC(kbottom)
& *recip_hFacI(i,j,kbottom)
ENDDO
ENDDO
endif
#endif
c-----------------------------------------------------------------------
c horizontal diffusion of IW energy
c-----------------------------------------------------------------------
DO k=2,Nr
DO j=1-OLy,sNy+OLy
dfx(1-OLx,j)=0. _d 0
DO i=1-OLx+1,sNx+OLx
fxa = IDEMIX_tau_h*0.5 _d 0*(
& IDEMIX_V0(i-1,j,k,bi,bj)*maskC(i-1,j,k,bi,bj)
& +IDEMIX_V0(i ,j,k,bi,bj)*maskC(i ,j,k,bi,bj))
dfx(i,j) = -fxa*_dyG(i,j,bi,bj)*drC(k)
& *(min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
& min(.5 _d 0,_hFacW(i,j,k ,bi,bj) ) )
& *_recip_dxC(i,j,bi,bj)
& *(IDEMIX_V0(i ,j,k,bi,bj)*IDEMIX_E(i ,j,k,bi,bj)
& -IDEMIX_V0(i-1,j,k,bi,bj)*IDEMIX_E(i-1,j,k,bi,bj))
& *maskW(i,j,k,bi,bj) ! paranoia setting
ENDDO
ENDDO
DO i=1-OLx,sNx+OLx
dfy(i,1-OLy)=0. _d 0
ENDDO
DO j=1-OLy+1,sNy+OLy
DO i=1-OLx,sNx+OLx
fxa = IDEMIX_tau_h*0.5 _d 0*(
& IDEMIX_V0(i,j ,k,bi,bj)*maskC(i,j ,k,bi,bj)
& +IDEMIX_V0(i,j-1,k,bi,bj)*maskC(i,j-1,k,bi,bj) )
dfy(i,j) = -fxa*_dxG(i,j,bi,bj)*drC(k)
& *(min(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
& min(.5 _d 0,_hFacS(i,j,k ,bi,bj) ) )
& *_recip_dyC(i,j,bi,bj)
& *(IDEMIX_V0(i,j ,k,bi,bj)*IDEMIX_E(i,j ,k,bi,bj)
& -IDEMIX_V0(i,j-1,k,bi,bj)*IDEMIX_E(i,j-1,k,bi,bj))
& *maskS(i,j,k,bi,bj) ! paranoia setting
ENDDO
ENDDO
c-----------------------------------------------------------------------
C Compute divergence of fluxes, add time tendency
c-----------------------------------------------------------------------
DO j=jMin,jMax
DO i=iMin,iMax
IDEMIX_E(i,j,k,bi,bj) = IDEMIX_E(i,j,k,bi,bj)
& + deltaTggl90*(-recip_drC(k)*recip_rA(i,j,bi,bj)
& *recip_hFacI(i,j,k)
& *((dfx(i+1,j)-dfx(i,j))+(dfy(i,j+1)-dfy(i,j)) ) )
& *maskC(i,j,k,bi,bj) ! paranoia setting
ENDDO
ENDDO
ENDDO ! k loop
c-----------------------------------------------------------------------
c add interior forcing e.g. by mesoscale GM
c-----------------------------------------------------------------------
DO k=2,Nr
DO j=jMin,jMax
DO i=iMin,iMax
IDEMIX_E(i,j,k,bi,bj) = IDEMIX_E(i,j,k,bi,bj)
& + forc(i,j,k)*deltaTggl90
ENDDO
ENDDO
ENDDO
c-----------------------------------------------------------------------
c solve vertical diffusion implicitly
c-----------------------------------------------------------------------
C delta_k = dt tau_v /drF_k (c_k+c_k+1)/2
DO k=2,Nr-1
DO j=jMin,jMax
DO i=iMin,iMax
delta(i,j,k) = deltaTggl90*IDEMIX_tau_v
& *recip_drF(k)*recip_hFacC(i,j,k,bi,bj)
& *.5 _d 0*(c0(i,j,k)+c0(i,j,k+1))
ENDDO
ENDDO
ENDDO
DO j=jMin,jMax
DO i=iMin,iMax
delta(i,j,1) = 0. _d 0
delta(i,j,Nr) = 0. _d 0
kBottom = MAX(kLowC(i,j,bi,bj),1)
delta(i,j,kBottom) = 0. _d 0
ENDDO
ENDDO
C-- Lower diagonal for E_(k-1) : -delta_k-1 c_k-1/drC_k
DO j=jMin,jMax
DO i=iMin,iMax
a3d(i,j,1) = 0. _d 0
a3d(i,j,2) = 0. _d 0
ENDDO
ENDDO
DO k=3,Nr
km1=MAX(2,k-1)
DO j=jMin,jMax
DO i=iMin,iMax
C- No need for maskC(k-1) with recip_hFacC(k-1) in delta(k-1)
a3d(i,j,k) = -delta(i,j,k-1)*c0(i,j,km1)
& *recip_drC(k)*recip_hFacI(i,j,k)
& *maskC(i,j,k,bi,bj)!*maskC(i,j,km1,bi,bj)
ENDDO
ENDDO
ENDDO
C-- Upper diagonal for E_(k+1): delta_k c_k+1/drC_k
DO k=2,Nr
DO j=jMin,jMax
DO i=iMin,iMax
C- No need for maskC(k) with recip_hFacC(k) in delta(k)
kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
c3d(i,j,k) = -delta(i,j,k)*c0(i,j,kp1)
& *recip_drC(k)*recip_hFacI(i,j,k)
& *maskC(i,j,k-1,bi,bj)
! & *maskC(i,j,k,bi,bj)*maskC(i,j,kp1,bi,bj)
ENDDO
ENDDO
ENDDO
DO j=jMin,jMax ! c3d at bottom is zero
DO i=iMin,iMax
c3d(i,j,1) = 0. _d 0
kBottom = MAX(kLowC(i,j,bi,bj),1)
c3d(i,j,kBottom) = 0. _d 0
ENDDO
ENDDO
C-- Center diagonal
DO j=jMin,jMax
DO i=iMin,iMax
b3d(i,j,1) = 1. _d 0
ENDDO
ENDDO
DO k=2,Nr
km1 = MAX(k-1,2)
DO j=jMin,jMax
DO i=iMin,iMax
b3d(i,j,k) = 1. _d 0 + deltaTggl90*IDEMIX_tau_d(i,j,k,bi,bj)
& *IDEMIX_E(i,j,k,bi,bj)
& *maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
C- No need for maskC(k) with recip_hFacC(k) in delta(k)
b3d(i,j,k) = b3d(i,j,k) + delta(i,j,k)*c0(i,j,k)
& *recip_drC(k)*recip_hFacI(i,j,k)
& *maskC(i,j,km1,bi,bj)
C- No need for maskC(k-1) with recip_hFacC(k-1) in delta(k-1)
b3d(i,j,k) = b3d(i,j,k) + delta(i,j,km1)*c0(i,j,k)
& *recip_drC(k)*recip_hFacI(i,j,k)
& *maskC(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
c at surface and bottom
DO j=jMin,jMax
DO i=iMin,iMax
k = MAX(kLowC(i,j,bi,bj),1)
km1 = MAX(k-1,2)
b3d(i,j,k) = 1. _d 0 + deltaTggl90*IDEMIX_tau_d(i,j,k,bi,bj)
& *IDEMIX_E(i,j,k,bi,bj)
& *maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
C- No need for maskC(k-1) with recip_hFacC(k-1) in delta(k-1)
& + delta(i,j,km1 )*c0(i,j,k)
& *recip_drC(k)*recip_hFacI(i,j,k)
& *maskC(i,j,k,bi,bj)
k=2
b3d(i,j,k) = 1. _d 0 + deltaTggl90*IDEMIX_tau_d(i,j,k,bi,bj)
& *IDEMIX_E(i,j,k,bi,bj)
& *maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
C- No need for maskC(k) with recip_hFacC(k) in delta(k)
& + delta(i,j,k)*c0(i,j,k)
& *recip_drC(k)*recip_hFacI(i,j,k)
& *maskC(i,j,km1,bi,bj)
ENDDO
ENDDO
C Apply flux boundary condition
DO j=jMin,jMax
DO i=iMin,iMax
k=2
IDEMIX_E(i,j,k,bi,bj) = IDEMIX_E(i,j,k,bi,bj)
& +deltaTggl90*IDEMIX_F_s(i,j,bi,bj)
& *recip_drC(k)*recip_hFacI(i,j,k)
& *maskC(i,j,k,bi,bj)
k = MAX(kLowC(i,j,bi,bj),1)
IDEMIX_E(i,j,k,bi,bj) = IDEMIX_E(i,j,k,bi,bj)
& -deltaTggl90*IDEMIX_F_b(i,j,bi,bj)
& *recip_drC(k)*recip_hFacI(i,j,k)
& *maskC(i,j,k,bi,bj)
ENDDO
ENDDO
C solve tri-diagonal system
errCode = -1
CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
I a3d, b3d, c3d,
U IDEMIX_E(1-OLx,1-OLy,1,bi,bj),
O errCode,
I bi, bj, myThid )
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
c-----------------------------------------------------------------------
c compute diffusivity due to internal wave breaking
c assuming local Osborn-Cox balance model
c kept for diagnostics only
c-----------------------------------------------------------------------
DO k=2,Nr
DO j=jMin,jMax
DO i=iMin,iMax
osborn_diff(i,j,k) = IDEMIX_mixing_efficiency
& *IDEMIX_tau_d(i,j,k,bi,bj)
& *IDEMIX_E(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)
& /max(1. _d -12,Nsquare(i,j,k))*maskC(i,j,k,bi,bj)
osborn_diff(i,j,k) = min(IDEMIX_diff_max,osborn_diff(i,j,k))
ENDDO
ENDDO
ENDDO
CALL DIAGNOSTICS_FILL( IDEMIX_E ,'IDEMIX_E',
& 0,Nr, 1, bi, bj, myThid )
CALL DIAGNOSTICS_FILL( IDEMIX_V0 ,'IDEMIX_v',
& 0,Nr, 1, bi, bj, myThid )
CALL DIAGNOSTICS_FILL( IDEMIX_tau_d ,'IDEMIX_t',
& 0,Nr, 1, bi, bj, myThid )
CALL DIAGNOSTICS_FILL( c0 ,'IDEMIX_c',
& 0,Nr, 2, bi, bj, myThid )
CALL DIAGNOSTICS_FILL( osborn_diff ,'IDEMIX_K',
& 0,Nr, 2, bi, bj, myThid )
CALL DIAGNOSTICS_FILL( forc ,'IDEMIX_F',
& 0,Nr, 2, bi, bj, myThid )
CALL DIAGNOSTICS_FILL(IDEMIX_F_b,'IDEM_F_b',0,1,1,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(IDEMIX_F_s,'IDEM_F_s',0,1,1,bi,bj,myThid)
CALL DIAGNOSTICS_FILL(gm_forc,'IDEM_F_g',
& 0,1,2,bi,bj,myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
#endif /* ALLOW_GGL90_IDEMIX */
RETURN
END
#ifdef ALLOW_GGL90_IDEMIX
C helper functions
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
_RL FUNCTION IDEMIX_gofx2(xx)
IMPLICIT NONE
_RL x,c,xx
_RL pi
PARAMETER( pi = 3.14159265358979323846264338327950588d0 )
x=MAX(3.d0,xx)
c= 1.d0-(2.d0/pi)*ASIN(1.d0/x)
IDEMIX_gofx2 = 2.d0/pi/c*0.9d0*x**(-2.d0/3.d0)*(1.-EXP(-x/4.3d0))
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
_RL FUNCTION IDEMIX_hofx1(x)
IMPLICIT NONE
_RL x
_RL pi
PARAMETER( pi = 3.14159265358979323846264338327950588d0 )
IDEMIX_hofx1 = (2.d0/pi)/(1.d0-(2.d0/pi)*
& ASIN(1.d0/MAX(1.01d0,x)))*(x-1.d0)/(x+1.d0)
END
#endif /* ALLOW_GGL90_IDEMIX */