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 */