pkg/ggl90/ggl90_calc.F

C $Header: /u/gcmpack/MITgcm/pkg/ggl90/ggl90_calc.F,v 1.1 2004/09/16 11:27:18 mlosch Exp $
C $Name:  $

#include "GGL90_OPTIONS.h"

CBOP
C !ROUTINE: GGL90_CALC

C !INTERFACE: ======================================================
      subroutine GGL90_CALC(
     I     bi, bj, myTime, myThid )

C !DESCRIPTION: \bv
C     /==========================================================\
C     | SUBROUTINE GGL90_CALC                                    |
C     | o Compute all GGL90 fields defined in GGL90.h            |
C     |==========================================================|
C     | Equation numbers refer to                                |
C     | Gaspar et al. (1990), JGR 95 (C9), pp 16,179             |
C     | Some parts of the implementation follow Blanke and       |
C     | Delecuse (1993), JPO, and OPA code, in particular the    |
C     | computation of the                                       |
C     | mixing length = max(min(lk,depth),lkmin)                 |
C     \==========================================================/
      IMPLICIT NONE
C
C--------------------------------------------------------------------

C global parameters updated by ggl90_calc
C     GGL90TKE     - sub-grid turbulent kinetic energy           (m^2/s^2)
C     GGL90viscAz  - GGL90 eddy viscosity coefficient              (m^2/s)
C     GGL90diffKzT - GGL90 diffusion coefficient for temperature   (m^2/s)
C
C \ev

C !USES: ============================================================
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GGL90.h"
#include "FFIELDS.h"
#include "GRID.h"

C !INPUT PARAMETERS: ===================================================
c Routine arguments
c     bi, bj - array indices on which to apply calculations
c     myTime - Current time in simulation

      INTEGER bi, bj
      INTEGER myThid
      _RL     myTime

#ifdef ALLOW_GGL90

C !LOCAL VARIABLES: ====================================================
c Local constants
C     iMin, iMax, jMin, jMax, I, J - array computation indices
C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices
C     ab15, ab05       - weights for implicit timestepping
C     uStarSquare      - square of friction velocity
C     verticalShear    - (squared) vertical shear of horizontal velocity
C     Nsquare          - squared buoyancy freqency 
C     RiNumber         - local Richardson number
C     KappaM           - (local) viscosity parameter (eq.10)
C     KappaH           - (local) diffusivity parameter for temperature (eq.11)
C     KappaE           - (local) diffusivity parameter for TKE (eq.15)
C     buoyFreq         - buoyancy freqency
C     TKEdissipation   - dissipation of TKE
C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse
C     totalDepth       - thickness of water column (inverse of recip_Rcol)
C     TKEPrandtlNumber - here, an empirical function of the Richardson number
C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)
C     gTKE             - right hand side of implicit equation
      INTEGER iMin ,iMax ,jMin ,jMax
      INTEGER I, J, K, Kp1, Km1, kSurf, kBottom
      _RL     ab15, ab05
      _RL     uStarSquare
      _RL     verticalShear
      _RL     KappaM, KappaH
      _RL     Nsquare
      _RL     deltaTggl90
      _RL     SQRTTKE
      _RL     RiNumber
      _RL     TKEdissipation
      _RL     tempU, tempV, prTemp
      _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
C     tri-diagonal matrix
      _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
#ifdef ALLOW_GGL90_HORIZDIFF
C     xA, yA   - area of lateral faces
C     dfx, dfy - diffusive flux across lateral faces
      _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
#endif /* ALLOW_GGL90_HORIZDIFF */
CEOP
      iMin = 2-OLx
      iMax = sNx+OLx-1
      jMin = 2-OLy
      jMax = sNy+OLy-1

C     set separate time step (should be deltaTtracer)
      deltaTggl90 = deltaTtracer
C     
      kSurf = 1
C     implicit timestepping weights for dissipation
      ab15 =  1.5 _d 0
      ab05 = -0.5 _d 0
      ab15 =  1. _d 0
      ab05 =  0. _d 0

C     Initialize local fields
      DO K = 1, Nr
       DO J=1-Oly,sNy+Oly
        DO I=1-Olx,sNx+Olx
         gTKE(I,J,K)              = 0. _d 0
         KappaE(I,J,K)            = 0. _d 0
         TKEPrandtlNumber(I,J,K)  = 0. _d 0
         GGL90mixingLength(I,J,K) = 0. _d 0
        ENDDO
       ENDDO    
      ENDDO
      DO J=1-Oly,sNy+Oly
       DO I=1-Olx,sNx+Olx
        rhoK    (I,J) = 0. _d 0
        rhoKm1  (I,J) = 0. _d 0
        totalDepth(I,J)   = 0. _d 0
        IF ( recip_Rcol(I,J,bi,bj) .NE. 0. ) 
     &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)
       ENDDO
      ENDDO

C     start k-loop
      DO K = 2, Nr
       Km1 = K-1
       Kp1 = MIN(Nr,K+1)
       CALL FIND_RHO(
     I      bi, bj, iMin, iMax, jMin, jMax, Km1, K,
     I      theta, salt,
     O      rhoKm1,
     I      myThid )
       CALL FIND_RHO(
     I      bi, bj, iMin, iMax, jMin, jMax, K, K,
     I      theta, salt,
     O      rhoK,
     I      myThid )
       DO J=jMin,jMax
        DO I=iMin,iMax
         SQRTTKE=SQRT( GGL90TKE(I,J,K,bi,bj) )
C
C     buoyancy frequency
C
         Nsquare = - gravity*recip_rhoConst*recip_drC(K)
     &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)
C     vertical shear term (dU/dz)^2+(dV/dz)^2
         tempu= .5*(  uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)
     &             - (uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )
     &        *recip_drC(K)
         tempv= .5*(  vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)
     &             - (vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )
     &        *recip_drC(K)
         verticalShear = tempU*tempU + tempV*tempV
         RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)
C     compute Prandtl number (always greater than 0)
         prTemp = 1. _d 0
         IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber
         TKEPrandtlNumber(I,J,K) = MIN(10.,prTemp)
C     mixing length
         GGL90mixingLength(I,J,K) = 
     &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) ) 
C     impose upper bound for mixing length (total depth)
         GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K), 
     &        totalDepth(I,J))
C     impose minimum mixing length (to avoid division by zero)
         GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K), 
     &        GGL90mixingLengthMin)
C     viscosity of last timestep
         KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE
         KappaE(I,J,K) = KappaM*GGL90alpha
C     dissipation term
         TKEdissipation = ab05*GGL90ceps
     &        *SQRTTKE/GGL90mixingLength(I,J,K)
     &        *GGL90TKE(I,J,K,bi,bj)      
C     sum up contributions to form the right hand side
         gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj) 
     &        + deltaTggl90*(
     &        + KappaM*verticalShear
     &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)
     &        - TKEdissipation 
     &        )
        ENDDO   
       ENDDO
      ENDDO
C     horizontal diffusion of TKE (requires an exchange in 
C       do_fields_blocking_exchanges)
#ifdef ALLOW_GGL90_HORIZDIFF
      IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
       DO K=2,Nr
C     common factors
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          xA(i,j) = _dyG(i,j,bi,bj)
     &         *drF(k)*_hFacW(i,j,k,bi,bj)
          yA(i,j) = _dxG(i,j,bi,bj)
     &         *drF(k)*_hFacS(i,j,k,bi,bj)
         ENDDO
        ENDDO   
C     Compute diffusive fluxes
C     ... across x-faces
        DO j=1-Oly,sNy+Oly
         dfx(1-Olx,j)=0.
         DO i=1-Olx+1,sNx+Olx
          dfx(i,j) = -GGL90diffTKEh*xA(i,j)
     &      *_recip_dxC(i,j,bi,bj)
     &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
     &      *CosFacU(j,bi,bj)
         ENDDO
        ENDDO
C     ... across y-faces
        DO i=1-Olx,sNx+Olx
         dfy(i,1-Oly)=0.
        ENDDO
        DO j=1-Oly+1,sNy+Oly
         DO i=1-Olx,sNx+Olx
          dfy(i,j) = -GGL90diffTKEh*yA(i,j)
     &      *_recip_dyC(i,j,bi,bj)
     &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
#ifdef ISOTROPIC_COS_SCALING
     &      *CosFacV(j,bi,bj)
#endif /* ISOTROPIC_COS_SCALING */
         ENDDO
        ENDDO   
C     Compute divergence of fluxes
        DO j=1-Oly,sNy+Oly-1
         DO i=1-Olx,sNx+Olx-1
          gTKE(i,j,k)=gTKE(i,j,k)
     &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
     &         *( (dfx(i+1,j)-dfx(i,j))
     &           +(dfy(i,j+1)-dfy(i,j)) 
     &           )
         ENDDO  
        ENDDO
C       end of k-loop
       ENDDO
C     end if GGL90diffTKEh .eq. 0.
      ENDIF
#endif /* ALLOW_GGL90_HORIZDIFF */

C     ============================================
C     Implicit time step to update TKE for k=1,Nr;
C     TKE(Nr+1)=0 by default
C     ============================================
C     set up matrix
C--   Lower diagonal
      DO j=jMin,jMax
       DO i=iMin,iMax
         a(i,j,1) = 0. _d 0 
       ENDDO
      ENDDO
      DO k=2,Nr
       km1=MAX(1,k-1)
       DO j=jMin,jMax
        DO i=iMin,iMax
         a(i,j,k) = -deltaTggl90
     &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)
     &        *.5*(KappaE(i,j, k )+KappaE(i,j,km1))
     &        *recip_drC(k)
          IF (recip_hFacI(i,j,km1,bi,bj).EQ.0.) a(i,j,k)=0.
        ENDDO
       ENDDO
      ENDDO
C--   Upper diagonal
      DO j=jMin,jMax
       DO i=iMin,iMax
         c(i,j,1)  = 0. _d 0
         c(i,j,Nr) = 0. _d 0
       ENDDO
      ENDDO
CML      DO k=1,Nr-1
      DO k=2,Nr-1
       kp1=min(Nr,k+1)
       DO j=jMin,jMax
        DO i=iMin,iMax
          c(i,j,k) = -deltaTggl90
     &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)
     &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))
     &        *recip_drC(k)
          IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0.) c(i,j,k)=0.
        ENDDO
       ENDDO
      ENDDO
C--   Center diagonal
      DO k=1,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax
          b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
     &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))
     &        /GGL90mixingLength(I,J,K)
         ENDDO
       ENDDO
      ENDDO
C     end set up matrix

C
C     Apply boundary condition
C     
      DO J=jMin,jMax
       DO I=iMin,iMax
C     estimate friction velocity uStar from surface forcing
        uStarSquare = SQRT( 
     &         ( .5*( surfaceForcingU(I,  J,  bi,bj)
     &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2
     &       + ( .5*( surfaceForcingV(I,  J,  bi,bj) 
     &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2
     &                     )
C     Dirichlet surface boundary condition for TKE
        gTKE(I,J,kSurf) = MAX(GGL90TKEmin,GGL90m2*uStarSquare)
     &                     *maskC(I,J,kSurf,bi,bj)
C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
        kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)
        gTKE(I,J,kBottom) = gTKE(I,J,kBottom) 
     &       - GGL90TKEbottom*c(I,J,kBottom)
       ENDDO
      ENDDO     
C
C     solve tri-diagonal system, and store solution on gTKE (previously rhs)
C
      CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,
     I     a, b, c,
     U     gTKE,
     I     myThid )
C
C     now update TKE
C     
      DO K=1,Nr
       DO J=jMin,jMax
        DO I=iMin,iMax
C     impose minimum TKE to avoid numerical undershoots below zero
         GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin ) 
     &        * maskC(I,J,K,bi,bj)
        ENDDO
       ENDDO
      ENDDO     
C
C     end of time step
C     ===============================
C     compute viscosity coefficients
C     
      DO K=2,Nr
       DO J=jMin,jMax
        DO I=iMin,iMax
C     Eq. (11), (18) and (21)
         KappaM = GGL90ck*GGL90mixingLength(I,J,K)*
     &                  SQRT( GGL90TKE(I,J,K,bi,bj) )
         KappaH = KappaM/TKEPrandtlNumber(I,J,K)
C     Set a minium (= background) value
         KappaM = MAX(KappaM,viscAr)
         KappaH = MAX(KappaH,diffKrT)
C     Set a maximum and mask land point
         GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)
     &        * maskC(I,J,K,bi,bj)
         GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)
     &        * maskC(I,J,K,bi,bj)
        ENDDO
       ENDDO 
C     end third k-loop
      ENDDO     

#endif /* ALLOW_GGL90 */

      RETURN
      END

1	mlosch	1.2	C $Header: /u/gcmpack/MITgcm/pkg/ggl90/ggl90_calc.F,v 1.1 2004/09/16 11:27:18 mlosch Exp $
2	mlosch	1.1	C $Name: $
3
4			#include "GGL90_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: GGL90_CALC
8
9			C !INTERFACE: ======================================================
10			subroutine GGL90_CALC(
11			I bi, bj, myTime, myThid )
12
13			C !DESCRIPTION: \bv
14			C /==========================================================\
15			C \| SUBROUTINE GGL90_CALC \|
16			C \| o Compute all GGL90 fields defined in GGL90.h \|
17			C \|==========================================================\|
18			C \| Equation numbers refer to \|
19			C \| Gaspar et al. (1990), JGR 95 (C9), pp 16,179 \|
20			C \| Some parts of the implementation follow Blanke and \|
21			C \| Delecuse (1993), JPO, and OPA code, in particular the \|
22			C \| computation of the \|
23			C \| mixing length = max(min(lk,depth),lkmin) \|
24			C \==========================================================/
25			IMPLICIT NONE
26			C
27			C--------------------------------------------------------------------
28
29			C global parameters updated by ggl90_calc
30			C GGL90TKE - sub-grid turbulent kinetic energy (m^2/s^2)
31			C GGL90viscAz - GGL90 eddy viscosity coefficient (m^2/s)
32			C GGL90diffKzT - GGL90 diffusion coefficient for temperature (m^2/s)
33			C
34			C \ev
35
36			C !USES: ============================================================
37			#include "SIZE.h"
38			#include "EEPARAMS.h"
39			#include "PARAMS.h"
40			#include "DYNVARS.h"
41			#include "GGL90.h"
42			#include "FFIELDS.h"
43			#include "GRID.h"
44
45			C !INPUT PARAMETERS: ===================================================
46			c Routine arguments
47			c bi, bj - array indices on which to apply calculations
48			c myTime - Current time in simulation
49
50			INTEGER bi, bj
51			INTEGER myThid
52			_RL myTime
53
54			#ifdef ALLOW_GGL90
55
56			C !LOCAL VARIABLES: ====================================================
57			c Local constants
58			C iMin, iMax, jMin, jMax, I, J - array computation indices
59			C K, Kp1, km1, kSurf, kBottom - vertical loop indices
60			C ab15, ab05 - weights for implicit timestepping
61			C uStarSquare - square of friction velocity
62			C verticalShear - (squared) vertical shear of horizontal velocity
63			C Nsquare - squared buoyancy freqency
64			C RiNumber - local Richardson number
65			C KappaM - (local) viscosity parameter (eq.10)
66			C KappaH - (local) diffusivity parameter for temperature (eq.11)
67			C KappaE - (local) diffusivity parameter for TKE (eq.15)
68			C buoyFreq - buoyancy freqency
69			C TKEdissipation - dissipation of TKE
70			C GGL90mixingLength- mixing length of scheme following Banke+Delecuse
71			C totalDepth - thickness of water column (inverse of recip_Rcol)
72			C TKEPrandtlNumber - here, an empirical function of the Richardson number
73			C rhoK, rhoKm1 - density at layer K and Km1 (relative to K)
74			C gTKE - right hand side of implicit equation
75			INTEGER iMin ,iMax ,jMin ,jMax
76			INTEGER I, J, K, Kp1, Km1, kSurf, kBottom
77			_RL ab15, ab05
78			_RL uStarSquare
79			_RL verticalShear
80			_RL KappaM, KappaH
81			_RL Nsquare
82			_RL deltaTggl90
83			_RL SQRTTKE
84			_RL RiNumber
85			_RL TKEdissipation
86			_RL tempU, tempV, prTemp
87			_RL TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
88			_RL GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89			_RL KappaE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
90			_RL rhoK (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91			_RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92			_RL totalDepth (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93			_RL gTKE (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
94			C tri-diagonal matrix
95			_RL a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96			_RL b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97			_RL c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
98	mlosch	1.2	#ifdef ALLOW_GGL90_HORIZDIFF
99			C xA, yA - area of lateral faces
100			C dfx, dfy - diffusive flux across lateral faces
101			_RL xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102			_RL yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103			_RL dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104			_RL dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105			#endif /* ALLOW_GGL90_HORIZDIFF */
106	mlosch	1.1	CEOP
107			iMin = 2-OLx
108			iMax = sNx+OLx-1
109			jMin = 2-OLy
110			jMax = sNy+OLy-1
111
112			C set separate time step (should be deltaTtracer)
113			deltaTggl90 = deltaTtracer
114			C
115			kSurf = 1
116			C implicit timestepping weights for dissipation
117			ab15 = 1.5 _d 0
118			ab05 = -0.5 _d 0
119			ab15 = 1. _d 0
120			ab05 = 0. _d 0
121
122			C Initialize local fields
123			DO K = 1, Nr
124			DO J=1-Oly,sNy+Oly
125			DO I=1-Olx,sNx+Olx
126			gTKE(I,J,K) = 0. _d 0
127			KappaE(I,J,K) = 0. _d 0
128			TKEPrandtlNumber(I,J,K) = 0. _d 0
129			GGL90mixingLength(I,J,K) = 0. _d 0
130			ENDDO
131			ENDDO
132			ENDDO
133			DO J=1-Oly,sNy+Oly
134			DO I=1-Olx,sNx+Olx
135			rhoK (I,J) = 0. _d 0
136			rhoKm1 (I,J) = 0. _d 0
137			totalDepth(I,J) = 0. _d 0
138			IF ( recip_Rcol(I,J,bi,bj) .NE. 0. )
139			& totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)
140			ENDDO
141			ENDDO
142
143			C start k-loop
144			DO K = 2, Nr
145			Km1 = K-1
146			Kp1 = MIN(Nr,K+1)
147			CALL FIND_RHO(
148			I bi, bj, iMin, iMax, jMin, jMax, Km1, K,
149			I theta, salt,
150			O rhoKm1,
151			I myThid )
152			CALL FIND_RHO(
153			I bi, bj, iMin, iMax, jMin, jMax, K, K,
154			I theta, salt,
155			O rhoK,
156			I myThid )
157			DO J=jMin,jMax
158			DO I=iMin,iMax
159			SQRTTKE=SQRT( GGL90TKE(I,J,K,bi,bj) )
160			C
161			C buoyancy frequency
162			C
163			Nsquare = - gravityrecip_rhoConstrecip_drC(K)
164			& * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)
165			C vertical shear term (dU/dz)^2+(dV/dz)^2
166			tempu= .5*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)
167			& - (uVel(I,J,K ,bi,bj)+uVel(I+1,J,K ,bi,bj)) )
168			& *recip_drC(K)
169			tempv= .5*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)
170			& - (vVel(I,J,K ,bi,bj)+vVel(I,J+1,K ,bi,bj)) )
171			& *recip_drC(K)
172			verticalShear = tempUtempU + tempVtempV
173			RiNumber = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)
174			C compute Prandtl number (always greater than 0)
175			prTemp = 1. _d 0
176			IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber
177			TKEPrandtlNumber(I,J,K) = MIN(10.,prTemp)
178			C mixing length
179			GGL90mixingLength(I,J,K) =
180			& SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )
181			C impose upper bound for mixing length (total depth)
182			GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
183			& totalDepth(I,J))
184			C impose minimum mixing length (to avoid division by zero)
185			GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
186			& GGL90mixingLengthMin)
187			C viscosity of last timestep
188			KappaM = GGL90ckGGL90mixingLength(I,J,K)SQRTTKE
189			KappaE(I,J,K) = KappaM*GGL90alpha
190			C dissipation term
191			TKEdissipation = ab05*GGL90ceps
192			& *SQRTTKE/GGL90mixingLength(I,J,K)
193			& *GGL90TKE(I,J,K,bi,bj)
194			C sum up contributions to form the right hand side
195			gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)
196			& + deltaTggl90*(
197			& + KappaM*verticalShear
198			& - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)
199			& - TKEdissipation
200			& )
201			ENDDO
202			ENDDO
203	mlosch	1.2	ENDDO
204			C horizontal diffusion of TKE (requires an exchange in
205			C do_fields_blocking_exchanges)
206			#ifdef ALLOW_GGL90_HORIZDIFF
207			IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
208			DO K=2,Nr
209			C common factors
210			DO j=1-Oly,sNy+Oly
211			DO i=1-Olx,sNx+Olx
212			xA(i,j) = _dyG(i,j,bi,bj)
213			& drF(k)_hFacW(i,j,k,bi,bj)
214			yA(i,j) = _dxG(i,j,bi,bj)
215			& drF(k)_hFacS(i,j,k,bi,bj)
216			ENDDO
217			ENDDO
218			C Compute diffusive fluxes
219			C ... across x-faces
220			DO j=1-Oly,sNy+Oly
221			dfx(1-Olx,j)=0.
222			DO i=1-Olx+1,sNx+Olx
223			dfx(i,j) = -GGL90diffTKEh*xA(i,j)
224			& *_recip_dxC(i,j,bi,bj)
225			& *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
226			& *CosFacU(j,bi,bj)
227			ENDDO
228			ENDDO
229			C ... across y-faces
230			DO i=1-Olx,sNx+Olx
231			dfy(i,1-Oly)=0.
232			ENDDO
233			DO j=1-Oly+1,sNy+Oly
234			DO i=1-Olx,sNx+Olx
235			dfy(i,j) = -GGL90diffTKEh*yA(i,j)
236			& *_recip_dyC(i,j,bi,bj)
237			& *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
238			#ifdef ISOTROPIC_COS_SCALING
239			& *CosFacV(j,bi,bj)
240			#endif /* ISOTROPIC_COS_SCALING */
241			ENDDO
242			ENDDO
243			C Compute divergence of fluxes
244			DO j=1-Oly,sNy+Oly-1
245			DO i=1-Olx,sNx+Olx-1
246			gTKE(i,j,k)=gTKE(i,j,k)
247			& -_recip_hFacC(i,j,k,bi,bj)recip_drF(k)recip_rA(i,j,bi,bj)
248			& *( (dfx(i+1,j)-dfx(i,j))
249			& +(dfy(i,j+1)-dfy(i,j))
250			& )
251			ENDDO
252			ENDDO
253			C end of k-loop
254			ENDDO
255			C end if GGL90diffTKEh .eq. 0.
256			ENDIF
257			#endif /* ALLOW_GGL90_HORIZDIFF */
258
259			C ============================================
260			C Implicit time step to update TKE for k=1,Nr;
261			C TKE(Nr+1)=0 by default
262			C ============================================
263	mlosch	1.1	C set up matrix
264	mlosch	1.2	C-- Lower diagonal
265	mlosch	1.1	DO j=jMin,jMax
266			DO i=iMin,iMax
267			a(i,j,1) = 0. _d 0
268			ENDDO
269			ENDDO
270			DO k=2,Nr
271			km1=MAX(1,k-1)
272			DO j=jMin,jMax
273			DO i=iMin,iMax
274			a(i,j,k) = -deltaTggl90
275			& recip_drF(km1)recip_hFacI(i,j,k,bi,bj)
276			& .5(KappaE(i,j, k )+KappaE(i,j,km1))
277			& *recip_drC(k)
278			IF (recip_hFacI(i,j,km1,bi,bj).EQ.0.) a(i,j,k)=0.
279			ENDDO
280			ENDDO
281			ENDDO
282	mlosch	1.2	C-- Upper diagonal
283	mlosch	1.1	DO j=jMin,jMax
284			DO i=iMin,iMax
285			c(i,j,1) = 0. _d 0
286			c(i,j,Nr) = 0. _d 0
287			ENDDO
288			ENDDO
289			CML DO k=1,Nr-1
290	mlosch	1.2	DO k=2,Nr-1
291	mlosch	1.1	kp1=min(Nr,k+1)
292			DO j=jMin,jMax
293			DO i=iMin,iMax
294			c(i,j,k) = -deltaTggl90
295			& recip_drF( k )recip_hFacI(i,j,k,bi,bj)
296			& .5(KappaE(i,j,k)+KappaE(i,j,kp1))
297			& *recip_drC(k)
298	mlosch	1.2	IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0.) c(i,j,k)=0.
299	mlosch	1.1	ENDDO
300			ENDDO
301			ENDDO
302	mlosch	1.2	C-- Center diagonal
303	mlosch	1.1	DO k=1,Nr
304			DO j=jMin,jMax
305			DO i=iMin,iMax
306			b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
307			& + ab15deltaTggl90GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))
308			& /GGL90mixingLength(I,J,K)
309			ENDDO
310			ENDDO
311			ENDDO
312			C end set up matrix
313
314			C
315			C Apply boundary condition
316			C
317			DO J=jMin,jMax
318			DO I=iMin,iMax
319			C estimate friction velocity uStar from surface forcing
320			uStarSquare = SQRT(
321			& ( .5*( surfaceForcingU(I, J, bi,bj)
322			& + surfaceForcingU(I+1,J, bi,bj) ) )**2
323			& + ( .5*( surfaceForcingV(I, J, bi,bj)
324			& + surfaceForcingV(I, J+1,bi,bj) ) )**2
325			& )
326			C Dirichlet surface boundary condition for TKE
327			gTKE(I,J,kSurf) = MAX(GGL90TKEmin,GGL90m2*uStarSquare)
328			& *maskC(I,J,kSurf,bi,bj)
329			C Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
330			kBottom = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)
331			gTKE(I,J,kBottom) = gTKE(I,J,kBottom)
332			& - GGL90TKEbottom*c(I,J,kBottom)
333			ENDDO
334			ENDDO
335			C
336			C solve tri-diagonal system, and store solution on gTKE (previously rhs)
337			C
338			CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,
339			I a, b, c,
340			U gTKE,
341			I myThid )
342			C
343			C now update TKE
344			C
345			DO K=1,Nr
346			DO J=jMin,jMax
347			DO I=iMin,iMax
348			C impose minimum TKE to avoid numerical undershoots below zero
349			GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )
350			& * maskC(I,J,K,bi,bj)
351			ENDDO
352			ENDDO
353			ENDDO
354			C
355	mlosch	1.2	C end of time step
356			C ===============================
357	mlosch	1.1	C compute viscosity coefficients
358			C
359			DO K=2,Nr
360			DO J=jMin,jMax
361			DO I=iMin,iMax
362			C Eq. (11), (18) and (21)
363			KappaM = GGL90ckGGL90mixingLength(I,J,K)
364			& SQRT( GGL90TKE(I,J,K,bi,bj) )
365			KappaH = KappaM/TKEPrandtlNumber(I,J,K)
366			C Set a minium (= background) value
367			KappaM = MAX(KappaM,viscAr)
368			KappaH = MAX(KappaH,diffKrT)
369			C Set a maximum and mask land point
370			GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)
371			& * maskC(I,J,K,bi,bj)
372			GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)
373			& * maskC(I,J,K,bi,bj)
374			ENDDO
375			ENDDO
376			C end third k-loop
377			ENDDO
378
379			#endif /* ALLOW_GGL90 */
380
381			RETURN
382			END
383