/[MITgcm]/MITgcm/pkg/ggl90/ggl90_calc.F
ViewVC logotype

Diff of /MITgcm/pkg/ggl90/ggl90_calc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.15 by dfer, Sat Nov 14 14:27:56 2009 UTC revision 1.33 by heimbach, Thu Jul 2 04:42:43 2015 UTC
# Line 7  CBOP Line 7  CBOP
7  C !ROUTINE: GGL90_CALC  C !ROUTINE: GGL90_CALC
8    
9  C !INTERFACE: ======================================================  C !INTERFACE: ======================================================
10        subroutine GGL90_CALC(        SUBROUTINE GGL90_CALC(
11       I     bi, bj, myTime, myThid )       I                 bi, bj, sigmaR, myTime, myIter, myThid )
12    
13  C !DESCRIPTION: \bv  C !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
# Line 35  C !USES: =============================== Line 35  C !USES: ===============================
35  #include "EEPARAMS.h"  #include "EEPARAMS.h"
36  #include "PARAMS.h"  #include "PARAMS.h"
37  #include "DYNVARS.h"  #include "DYNVARS.h"
 #include "GGL90.h"  
38  #include "FFIELDS.h"  #include "FFIELDS.h"
39  #include "GRID.h"  #include "GRID.h"
40    #include "GGL90.h"
41    
42  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
43  C Routine arguments  C Routine arguments
44  C     bi, bj :: array indices on which to apply calculations  C     bi, bj :: Current tile indices
45    C     sigmaR :: Vertical gradient of iso-neutral density
46  C     myTime :: Current time in simulation  C     myTime :: Current time in simulation
47    C     myIter :: Current time-step number
48  C     myThid :: My Thread Id number  C     myThid :: My Thread Id number
49        INTEGER bi, bj        INTEGER bi, bj
50          _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51        _RL     myTime        _RL     myTime
52          INTEGER myIter
53        INTEGER myThid        INTEGER myThid
 CEOP  
54    
55  #ifdef ALLOW_GGL90  #ifdef ALLOW_GGL90
56    
57  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
58  C Local constants  C Local constants
59  C     iMin, iMax, jMin, jMax, I, J - array computation indices  C     iMin,iMax,jMin,jMax :: index boundaries of computation domain
60  C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices  C     i, j, k, kp1,km1 :: array computation indices
61  C     ab15, ab05       - weights for implicit timestepping  C     kSurf, kBottom   :: vertical indices of domain boundaries
62  C     uStarSquare      - square of friction velocity  C     hFac/hFacI       :: fractional thickness of W-cell
63  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     explDissFac      :: explicit Dissipation Factor (in [0-1])
64  C     Nsquare          - squared buoyancy freqency  C     implDissFac      :: implicit Dissipation Factor (in [0-1])
65  C     RiNumber         - local Richardson number  C     uStarSquare      :: square of friction velocity
66  C     KappaM           - (local) viscosity parameter (eq.10)  C     verticalShear    :: (squared) vertical shear of horizontal velocity
67  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     Nsquare          :: squared buoyancy freqency
68  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     RiNumber         :: local Richardson number
69  C     TKEdissipation   - dissipation of TKE  C     KappaM           :: (local) viscosity parameter (eq.10)
70  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
71  C         rMixingLength- inverse of mixing length  C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
72  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C     TKEdissipation   :: dissipation of TKE
73  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C     GGL90mixingLength:: mixing length of scheme following Banke+Delecuse
74  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C         rMixingLength:: inverse of mixing length
75  C     gTKE             - right hand side of implicit equation  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
76    C     TKEPrandtlNumber :: here, an empirical function of the Richardson number
77        INTEGER iMin ,iMax ,jMin ,jMax        INTEGER iMin ,iMax ,jMin ,jMax
78        INTEGER I, J, K, Kp1, Km1, kSurf, kBottom        INTEGER i, j, k, kp1, km1, kSurf, kBottom
79        _RL     ab15, ab05        _RL     explDissFac, implDissFac
80        _RL     uStarSquare        _RL     uStarSquare
81        _RL     verticalShear        _RL     verticalShear(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82        _RL     KappaM, KappaH        _RL     KappaM(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83          _RL     KappaH
84  c     _RL     Nsquare  c     _RL     Nsquare
85        _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
86        _RL     deltaTggl90        _RL     deltaTggl90
87  c     _RL     SQRTTKE  c     _RL     SQRTTKE
88        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89        _RL     RiNumber        _RL     RiNumber
90    #ifdef ALLOW_GGL90_IDEMIX
91          _RL     IDEMIX_RiNumber
92    #endif
93        _RL     TKEdissipation        _RL     TKEdissipation
94        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
95        _RL     MaxLength, tmpmlx        _RL     MaxLength, tmpmlx, tmpVisc
96        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
98        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
99        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100        _RL     KappaE           (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)  
101        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103    #ifdef ALLOW_DIAGNOSTICS
104          _RL     surf_flx_tke     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105    #endif /* ALLOW_DIAGNOSTICS */
106    C     hFac(I)  :: fractional thickness of W-cell
107          _RL       hFac
108    #ifdef ALLOW_GGL90_IDEMIX
109          _RL       hFacI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
110    #endif /* ALLOW_GGL90_IDEMIX */
111          _RL recip_hFacI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
112  C-    tri-diagonal matrix  C-    tri-diagonal matrix
113        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
114        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
115        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
116        INTEGER errCode        INTEGER errCode
117  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
118  C-    xA, yA   - area of lateral faces  C     xA, yA   :: area of lateral faces
119  C-    dfx, dfy - diffusive flux across lateral faces  C     dfx, dfy :: diffusive flux across lateral faces
120    C     gTKE     :: right hand side of diffusion equation
121        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125          _RL    gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
127  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
128        _RL p4, p8, p16, tmpdiffKrS        _RL p4, p8, p16
129    CEOP
130        p4=0.25 _d 0        p4=0.25 _d 0
131        p8=0.125 _d 0        p8=0.125 _d 0
132        p16=0.0625 _d 0        p16=0.0625 _d 0
# Line 122  C     set separate time step (should be Line 140  C     set separate time step (should be
140        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
141    
142        kSurf = 1        kSurf = 1
143  C     implicit timestepping weights for dissipation  C     explicit/implicit timestepping weights for dissipation
144        ab15 =  1.5 _d 0        explDissFac = 0. _d 0
145        ab05 = -0.5 _d 0        implDissFac = 1. _d 0 - explDissFac
146        ab15 =  1. _d 0  
147        ab05 =  0. _d 0  C     For nonlinear free surface and especially with r*-coordinates, the
148    C     hFacs change every timestep, so we need to update them here in the
149    C     case of using IDEMIX.
150           DO K=1,Nr
151            Km1 = MAX(K-1,1)
152            DO j=1-OLy,sNy+OLy
153             DO i=1-OLx,sNx+OLx
154              hFac =
155         &         MIN(.5 _d 0,_hFacC(i,j,km1,bi,bj) ) +
156         &         MIN(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
157              recip_hFacI(I,J,K)=0. _d 0
158              IF ( hFac .NE. 0. _d 0 )
159         &         recip_hFacI(I,J,K)=1. _d 0/hFac
160    #ifdef ALLOW_GGL90_IDEMIX
161              hFacI(i,j,k) = hFac
162    #endif /* ALLOW_GGL90_IDEMIX */
163             ENDDO
164            ENDDO
165           ENDDO
166    
167  C     Initialize local fields  C     Initialize local fields
168        DO K = 1, Nr        DO k = 1, Nr
169         DO J=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
170          DO I=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
171           gTKE(I,J,K)              = 0. _d 0           rMixingLength(i,j,k)     = 0. _d 0
172           KappaE(I,J,K)            = 0. _d 0           mxLength_Dn(i,j,k)       = 0. _d 0
173           TKEPrandtlNumber(I,J,K)  = 1. _d 0           GGL90visctmp(i,j,k)      = 0. _d 0
174           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin           KappaE(i,j,k)            = 0. _d 0
175             TKEPrandtlNumber(i,j,k)  = 1. _d 0
176             GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
177             GGL90visctmp(i,j,k)      = 0. _d 0
178    #ifndef SOLVE_DIAGONAL_LOWMEMORY
179             a3d(i,j,k) = 0. _d 0
180             b3d(i,j,k) = 1. _d 0
181             c3d(i,j,k) = 0. _d 0
182    #endif
183             Nsquare(i,j,k) = 0. _d 0
184             SQRTTKE(i,j,k) = 0. _d 0
185          ENDDO          ENDDO
186         ENDDO         ENDDO
187        ENDDO        ENDDO
188        DO J=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
189         DO I=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
190          rhoK(I,J)          = 0. _d 0          KappaM(i,j)        = 0. _d 0
191          rhoKm1(I,J)        = 0. _d 0          verticalShear(i,j) = 0. _d 0
192          totalDepth(I,J)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)          totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
193          rMixingLength(i,j,1) = 0. _d 0          rMixingLength(i,j,1) = 0. _d 0
194          mxLength_Dn(I,J,1) = GGL90mixingLengthMin          mxLength_Dn(i,j,1) = GGL90mixingLengthMin
195          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )          SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
196    #ifdef ALLOW_GGL90_HORIZDIFF
197            xA(i,j)  = 0. _d 0
198            yA(i,j)  = 0. _d 0
199            dfx(i,j) = 0. _d 0
200            dfy(i,j) = 0. _d 0
201            gTKE(i,j) = 0. _d 0
202    #endif /* ALLOW_GGL90_HORIZDIFF */
203         ENDDO         ENDDO
204        ENDDO        ENDDO
205    
206  C     start k-loop  #ifdef ALLOW_GGL90_IDEMIX
207        DO K = 2, Nr        IF ( useIDEMIX) CALL GGL90_IDEMIX(
208         Km1 = K-1       &   bi, bj, hFacI, recip_hFacI, sigmaR, myTime, myIter, myThid )
209  c      Kp1 = MIN(Nr,K+1)  #endif /* ALLOW_GGL90_IDEMIX */
210         CALL FIND_RHO_2D(  
211       I      iMin, iMax, jMin, jMax, K,        DO k = 2, Nr
212       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),         DO j=jMin,jMax
213       O      rhoKm1,          DO i=iMin,iMax
214       I      Km1, bi, bj, myThid )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
   
        CALL FIND_RHO_2D(  
      I      iMin, iMax, jMin, jMax, K,  
      I      theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),  
      O      rhoK,  
      I      K, bi, bj, myThid )  
        DO J=jMin,jMax  
         DO I=iMin,iMax  
          SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )  
215    
216  C     buoyancy frequency  C     buoyancy frequency
217           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
218       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &                  * sigmaR(i,j,k)
219  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  C     vertical shear term (dU/dz)^2+(dV/dz)^2 is computed later
220  c         tempU= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)  C     to save some memory
 c     &                 -( uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )  
 c     &        *recip_drC(K)  
 c         tempV= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)  
 c     &                 -( vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )  
 c     &        *recip_drC(K)  
 c         verticalShear = tempU*tempU + tempV*tempV  
 c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  
 cC     compute Prandtl number (always greater than 0)  
 c         prTemp = 1. _d 0  
 c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber  
 c         TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  
221  C     mixing length  C     mixing length
222           GGL90mixingLength(I,J,K) = SQRTTWO *           GGL90mixingLength(i,j,k) = SQRTTWO *
223       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )       &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
224          ENDDO          ENDDO
225         ENDDO         ENDDO
226        ENDDO        ENDDO
227    
228  C- Impose upper and lower bound for mixing length  C- ensure mixing between first and second level
229          IF (mxlSurfFlag) THEN
230           DO j=jMin,jMax
231            DO i=iMin,iMax
232             GGL90mixingLength(i,j,2)=drF(1)
233            ENDDO
234           ENDDO
235          ENDIF
236    
237    C--   Impose upper and lower bound for mixing length
238    C--   Impose minimum mixing length to avoid division by zero
239        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
240  C-  
241         DO k=2,Nr         DO k=2,Nr
242          DO J=jMin,jMax          DO j=jMin,jMax
243           DO I=iMin,iMax           DO i=iMin,iMax
244            MaxLength=totalDepth(I,J)            MaxLength=totalDepth(i,j)
245            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
246       &                                   MaxLength)       &                                   MaxLength)
247           ENDDO           ENDDO
248          ENDDO          ENDDO
249         ENDDO         ENDDO
250    
251         DO k=2,Nr         DO k=2,Nr
252          DO J=jMin,jMax          DO j=jMin,jMax
253           DO I=iMin,iMax           DO i=iMin,iMax
254            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
255       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
256            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
257           ENDDO           ENDDO
258          ENDDO          ENDDO
259         ENDDO         ENDDO
260    
261        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
262  C-  
263         DO k=2,Nr         DO k=2,Nr
264          DO J=jMin,jMax          DO j=jMin,jMax
265           DO I=iMin,iMax           DO i=iMin,iMax
266            MaxLength=MIN(Ro_surf(I,J,bi,bj)-rF(k),rF(k)-R_low(I,J,bi,bj))            MaxLength=MIN(Ro_surf(i,j,bi,bj)-rF(k),rF(k)-R_low(i,j,bi,bj))
267  c         MaxLength=MAX(MaxLength,20. _d 0)  c         MaxLength=MAX(MaxLength,20. _d 0)
268            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
269       &                                   MaxLength)       &                                   MaxLength)
270           ENDDO           ENDDO
271          ENDDO          ENDDO
272         ENDDO         ENDDO
273    
274         DO k=2,Nr         DO k=2,Nr
275          DO J=jMin,jMax          DO j=jMin,jMax
276           DO I=iMin,iMax           DO i=iMin,iMax
277            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
278       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
279            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
280           ENDDO           ENDDO
281          ENDDO          ENDDO
282         ENDDO         ENDDO
283    
284        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
285  C-  
286         DO k=2,Nr         DO k=2,Nr
287          DO J=jMin,jMax          DO j=jMin,jMax
288           DO I=iMin,iMax           DO i=iMin,iMax
289            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
290       &        GGL90mixingLength(I,J,K-1)+drF(k-1))       &        GGL90mixingLength(i,j,k-1)+drF(k-1))
291           ENDDO           ENDDO
292          ENDDO          ENDDO
293         ENDDO         ENDDO
294         DO J=jMin,jMax         DO j=jMin,jMax
295          DO I=iMin,iMax          DO i=iMin,iMax
296            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
297       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
298          ENDDO          ENDDO
299         ENDDO         ENDDO
300         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
301          DO J=jMin,jMax          DO j=jMin,jMax
302           DO I=iMin,iMax           DO i=iMin,iMax
303            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
304       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
305           ENDDO           ENDDO
306          ENDDO          ENDDO
307         ENDDO         ENDDO
308    
309         DO k=2,Nr         DO k=2,Nr
310          DO J=jMin,jMax          DO j=jMin,jMax
311           DO I=iMin,iMax           DO i=iMin,iMax
312            GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
313       &                                   GGL90mixingLengthMin)       &                                   GGL90mixingLengthMin)
314            rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)            rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
315           ENDDO           ENDDO
316          ENDDO          ENDDO
317         ENDDO         ENDDO
318    
319        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
320  C-  
321         DO k=2,Nr         DO k=2,Nr
322          DO J=jMin,jMax          DO j=jMin,jMax
323           DO I=iMin,iMax           DO i=iMin,iMax
324            mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),            mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
325       &        mxLength_Dn(I,J,K-1)+drF(k-1))       &        mxLength_Dn(i,j,k-1)+drF(k-1))
326           ENDDO           ENDDO
327          ENDDO          ENDDO
328         ENDDO         ENDDO
329         DO J=jMin,jMax         DO j=jMin,jMax
330          DO I=iMin,iMax          DO i=iMin,iMax
331            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),            GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
332       &       GGL90mixingLengthMin+drF(Nr))       &       GGL90mixingLengthMin+drF(Nr))
333          ENDDO          ENDDO
334         ENDDO         ENDDO
335         DO k=Nr-1,2,-1         DO k=Nr-1,2,-1
336          DO J=jMin,jMax          DO j=jMin,jMax
337           DO I=iMin,iMax           DO i=iMin,iMax
338            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
339       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(i,j,k+1)+drF(k))
340           ENDDO           ENDDO
341          ENDDO          ENDDO
342         ENDDO         ENDDO
343    
344         DO k=2,Nr         DO k=2,Nr
345          DO J=jMin,jMax          DO j=jMin,jMax
346           DO I=iMin,iMax           DO i=iMin,iMax
347            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
348       &                                  mxLength_Dn(I,J,K))       &                                  mxLength_Dn(i,j,k))
349            tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )            tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
350            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)            tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
351            rMixingLength(I,J,K) = 1. _d 0 / tmpmlx            rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
352           ENDDO           ENDDO
353          ENDDO          ENDDO
354         ENDDO         ENDDO
# Line 314  C- Line 357  C-
357         STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'         STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
358        ENDIF        ENDIF
359    
360  C- Impose minimum mixing length (to avoid division by zero)  C     start "proper" k-loop (the code above was moved out and up to
361  c      DO k=2,Nr  C     implemement various mixing parameters efficiently)
 c      DO J=jMin,jMax  
 c       DO I=iMin,iMax  
 c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  
 c    &        GGL90mixingLengthMin)  
 c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)  
 c       ENDDO  
 c      ENDDO  
 c     ENDDO  
   
362        DO k=2,Nr        DO k=2,Nr
363         Km1 = K-1         km1 = k-1
        DO J=jMin,jMax  
         DO I=iMin,iMax  
 C     vertical shear term (dU/dz)^2+(dV/dz)^2  
          tempU= .5 _d 0*( 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 _d 0*( 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(i,j,k),0. _d 0)/(verticalShear+GGL90eps)  
 C     compute Prandtl number (always greater than 0)  
          prTemp = 1. _d 0  
          IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber  
          TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)  
 c         TKEPrandtlNumber(I,J,K) = 1. _d 0  
   
 C     viscosity and diffusivity  
          KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)  
          KappaH = KappaM/TKEPrandtlNumber(I,J,K)  
364    
 C     Set a minium (= background) and maximum value  
          KappaM = MAX(KappaM,viscArNr(k))  
          KappaH = MAX(KappaH,diffKrNrT(k))  
          KappaM = MIN(KappaM,GGL90viscMax)  
          KappaH = MIN(KappaH,GGL90diffMax)  
   
 C     Mask land points and storage  
          GGL90viscAr(I,J,K,bi,bj) = KappaM * maskC(I,J,K,bi,bj)  
          GGL90diffKr(I,J,K,bi,bj) = KappaH * maskC(I,J,K,bi,bj)  
          KappaE(I,J,K) = GGL90alpha * GGL90viscAr(I,J,K,bi,bj)  
   
 C     dissipation term  
          TKEdissipation = ab05*GGL90ceps  
      &        *SQRTTKE(i,j,k)*rMixingLength(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  
      &        - KappaH*Nsquare(i,j,k)  
      &        - TKEdissipation  
      &        )  
         ENDDO  
        ENDDO  
       ENDDO  
   
 C     horizontal diffusion of TKE (requires an exchange in  
 C      do_fields_blocking_exchanges)  
365  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
366        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
367         DO K=2,Nr  C     horizontal diffusion of TKE (requires an exchange in
368    C      do_fields_blocking_exchanges)
369  C     common factors  C     common factors
370          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
371           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
372            xA(i,j) = _dyG(i,j,bi,bj)            xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
373       &         *drF(k)*_hFacW(i,j,k,bi,bj)       &                 (min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
374            yA(i,j) = _dxG(i,j,bi,bj)       &                  min(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
375       &         *drF(k)*_hFacS(i,j,k,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)*drC(k)*
376         &                 (min(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
377         &                  min(.5 _d 0,_hFacS(i,j,k  ,bi,bj) ) )
378           ENDDO           ENDDO
379          ENDDO          ENDDO
380  C     Compute diffusive fluxes  C     Compute diffusive fluxes
381  C     ... across x-faces  C     ... across x-faces
382          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
383           dfx(1-Olx,j)=0. _d 0           dfx(1-OLx,j)=0. _d 0
384           DO i=1-Olx+1,sNx+Olx           DO i=1-OLx+1,sNx+OLx
385            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
386       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
387       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
388    #ifdef ISOTROPIC_COS_SCALING
389       &      *CosFacU(j,bi,bj)       &      *CosFacU(j,bi,bj)
390    #endif /* ISOTROPIC_COS_SCALING */
391           ENDDO           ENDDO
392          ENDDO          ENDDO
393  C     ... across y-faces  C     ... across y-faces
394          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
395           dfy(i,1-Oly)=0. _d 0           dfy(i,1-OLy)=0. _d 0
396          ENDDO          ENDDO
397          DO j=1-Oly+1,sNy+Oly          DO j=1-OLy+1,sNy+OLy
398           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
399            dfy(i,j) = -GGL90diffTKEh*yA(i,j)            dfy(i,j) = -GGL90diffTKEh*yA(i,j)
400       &      *_recip_dyC(i,j,bi,bj)       &      *_recip_dyC(i,j,bi,bj)
401       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))       &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
# Line 414  C     ... across y-faces Line 405  C     ... across y-faces
405           ENDDO           ENDDO
406          ENDDO          ENDDO
407  C     Compute divergence of fluxes  C     Compute divergence of fluxes
408          DO j=1-Oly,sNy+Oly-1          DO j=1-OLy,sNy+OLy-1
409           DO i=1-Olx,sNx+Olx-1           DO i=1-OLx,sNx+OLx-1
410            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)
411       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &         *recip_hFacI(i,j,k)
412       &         *( (dfx(i+1,j)-dfx(i,j))       &         *((dfx(i+1,j)-dfx(i,j))
413       &           +(dfy(i,j+1)-dfy(i,j))       &         + (dfy(i,j+1)-dfy(i,j)) )
      &           )*deltaTggl90  
414           ENDDO           ENDDO
415          ENDDO          ENDDO
 C       end of k-loop  
        ENDDO  
416  C     end if GGL90diffTKEh .eq. 0.  C     end if GGL90diffTKEh .eq. 0.
417        ENDIF         ENDIF
418  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
419    
420    C     viscosity and diffusivity
421           DO j=jMin,jMax
422            DO i=iMin,iMax
423             KappaM(i,j) = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
424             GGL90visctmp(i,j,k) = MAX(KappaM(i,j),diffKrNrS(k))
425         &                            * maskC(i,j,k,bi,bj)
426    C        note: storing GGL90visctmp like this, and using it later to compute
427    C              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
428             KappaM(i,j) = MAX(KappaM(i,j),viscArNr(k)) * maskC(i,j,k,bi,bj)
429            ENDDO
430           ENDDO
431    
432    C     compute Prandtl number (always greater than 0)
433    #ifdef ALLOW_GGL90_IDEMIX
434           IF ( useIDEMIX ) THEN
435           DO j=jMin,jMax
436            DO i=iMin,iMax
437    C     vertical shear term (dU/dz)^2+(dV/dz)^2
438             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
439         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
440         &        *recip_drC(k)*recip_hFacI(i,j,k)
441             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
442         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
443         &        *recip_drC(k)*recip_hFacI(i,j,k)
444             verticalShear(i,j) = tempU*tempU + tempV*tempV
445             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
446         &        /(verticalShear(i,j)+GGL90eps)
447    CML         IDEMIX_RiNumber = 1./GGL90eps
448             IDEMIX_RiNumber = MAX( KappaM(i,j)*Nsquare(i,j,k), 0. _d 0)/
449         &    (GGL90eps+IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2)
450             prTemp         = MIN(5.*RiNumber, 6.6 _d 0*IDEMIX_RiNumber)
451             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
452             TKEPrandtlNumber(i,j,k) = MAX( 1. _d 0,TKEPrandtlNumber(i,j,k))
453            ENDDO
454           ENDDO
455           ELSE
456    #else /* ndef ALLOW_GGL90_IDEMIX */
457           IF (.TRUE.) THEN
458    #endif /* ALLOW_GGL90_IDEMIX */
459           DO j=jMin,jMax
460            DO i=iMin,iMax
461             tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
462         &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
463         &        *recip_drC(k)
464             tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
465         &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
466         &        *recip_drC(k)
467             verticalShear(i,j) = tempU*tempU + tempV*tempV
468             RiNumber = MAX(Nsquare(i,j,k),0. _d 0)
469         &        /(verticalShear(i,j)+GGL90eps)
470             prTemp = 1. _d 0
471             IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
472             TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
473            ENDDO
474           ENDDO
475           ENDIF
476    
477           DO j=jMin,jMax
478            DO i=iMin,iMax
479    C        diffusivity
480             KappaH = KappaM(i,j)/TKEPrandtlNumber(i,j,k)
481             KappaE(i,j,k) = GGL90alpha * KappaM(i,j) * maskC(i,j,k,bi,bj)
482    
483    C     dissipation term
484             TKEdissipation = explDissFac*GGL90ceps
485         &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
486         &        *GGL90TKE(i,j,k,bi,bj)
487    C     partial update with sum of explicit contributions
488             GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
489         &        + deltaTggl90*(
490         &        + KappaM(i,j)*verticalShear(i,j)
491         &        - KappaH*Nsquare(i,j,k)
492         &        - TKEdissipation
493         &        )
494            ENDDO
495           ENDDO
496    
497    #ifdef ALLOW_GGL90_IDEMIX
498           IF ( useIDEMIX ) THEN
499    C     add IDEMIX contribution to the turbulent kinetic energy
500            DO j=jMin,jMax
501             DO i=iMin,iMax
502              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
503         &         + deltaTggl90*(
504         &         + IDEMIX_tau_d(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)**2
505         &         )
506             ENDDO
507            ENDDO
508           ENDIF
509    #endif /* ALLOW_GGL90_IDEMIX */
510    
511    #ifdef ALLOW_GGL90_HORIZDIFF
512           IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
513    C--    Add horiz. diffusion tendency
514            DO j=jMin,jMax
515             DO i=iMin,iMax
516              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
517         &                          + gTKE(i,j)*deltaTggl90
518             ENDDO
519            ENDDO
520           ENDIF
521    #endif /* ALLOW_GGL90_HORIZDIFF */
522    
523    C--   end of k loop
524          ENDDO
525    
526  C     ============================================  C     ============================================
527  C     Implicit time step to update TKE for k=1,Nr;  C     Implicit time step to update TKE for k=1,Nr;
528  C     TKE(Nr+1)=0 by default  C     TKE(Nr+1)=0 by default
# Line 437  C     set up matrix Line 531  C     set up matrix
531  C--   Lower diagonal  C--   Lower diagonal
532        DO j=jMin,jMax        DO j=jMin,jMax
533         DO i=iMin,iMax         DO i=iMin,iMax
534           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
535         ENDDO         ENDDO
536        ENDDO        ENDDO
537        DO k=2,Nr        DO k=2,Nr
# Line 447  C--   Lower diagonal Line 541  C--   Lower diagonal
541  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
542  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
543  C-    No need for maskC(k-1) with recip_hFacC(k-1)  C-    No need for maskC(k-1) with recip_hFacC(k-1)
544           a(i,j,k) = -deltaTggl90           a3d(i,j,k) = -deltaTggl90
545       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
546       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
547       &        *recip_drC(k)*maskC(i,j,k,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
# Line 457  C-    No need for maskC(k-1) with recip_ Line 551  C-    No need for maskC(k-1) with recip_
551  C--   Upper diagonal  C--   Upper diagonal
552        DO j=jMin,jMax        DO j=jMin,jMax
553         DO i=iMin,iMax         DO i=iMin,iMax
554           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
555         ENDDO         ENDDO
556        ENDDO        ENDDO
557        DO k=2,Nr        DO k=2,Nr
558         DO j=jMin,jMax         DO j=jMin,jMax
559          DO i=iMin,iMax          DO i=iMin,iMax
560            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))           kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
561  C-    We keep recip_hFacC in the diffusive flux calculation,  C-    We keep recip_hFacC in the diffusive flux calculation,
562  C-    but no hFacC in TKE volume control  C-    but no hFacC in TKE volume control
563  C-    No need for maskC(k) with recip_hFacC(k)  C-    No need for maskC(k) with recip_hFacC(k)
564            c(i,j,k) = -deltaTggl90           c3d(i,j,k) = -deltaTggl90
565       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
566       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
567       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
568          ENDDO          ENDDO
569         ENDDO         ENDDO
570        ENDDO        ENDDO
571    
572    #ifdef ALLOW_GGL90_IDEMIX
573          IF ( useIDEMIX ) THEN
574           DO k=2,Nr
575            DO j=jMin,jMax
576             DO i=iMin,iMax
577              a3d(i,j,k) = a3d(i,j,k)*recip_hFacI(i,j,k)
578              c3d(i,j,k) = c3d(i,j,k)*recip_hFacI(i,j,k)
579             ENDDO
580            ENDDO
581           ENDDO
582          ENDIF
583    #endif /* ALLOW_GGL90_IDEMIX */
584    
585          IF (.NOT.GGL90_dirichlet) THEN
586    C      Neumann bottom boundary condition for TKE: no flux from bottom
587           DO j=jMin,jMax
588            DO i=iMin,iMax
589             kBottom   = MAX(kLowC(i,j,bi,bj),1)
590             c3d(i,j,kBottom) = 0. _d 0
591            ENDDO
592           ENDDO
593          ENDIF
594    
595  C--   Center diagonal  C--   Center diagonal
596        DO k=1,Nr        DO k=1,Nr
597         km1 = MAX(k-1,1)         km1 = MAX(k-1,1)
598         DO j=jMin,jMax         DO j=jMin,jMax
599          DO i=iMin,iMax          DO i=iMin,iMax
600            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)           b3d(i,j,k) = 1. _d 0 - c3d(i,j,k) - a3d(i,j,k)
601       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
602       &        * rMixingLength(I,J,K)       &        * rMixingLength(i,j,k)
603       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)       &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
604           ENDDO          ENDDO
605         ENDDO         ENDDO
606        ENDDO        ENDDO
607  C     end set up matrix  C     end set up matrix
608    
609  C     Apply boundary condition  C     Apply boundary condition
610        kp1 = MIN(Nr,kSurf+1)        kp1 = MIN(Nr,kSurf+1)
611        DO J=jMin,jMax        DO j=jMin,jMax
612         DO I=iMin,iMax         DO i=iMin,iMax
613  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
614          uStarSquare = SQRT(          uStarSquare = SQRT(
615       &    ( .5 _d 0*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
616       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
617       &  + ( .5 _d 0*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
618       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
619       &                     )       &                     )
620  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
621          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
622       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
623          gTKE(i,j,kp1) = gTKE(i,j,kp1)          GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
624       &                - a(i,j,kp1)*gTKE(i,j,kSurf)       &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
625          a(i,j,kp1) = 0. _d 0          a3d(i,j,kp1) = 0. _d 0
 C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  
         kBottom   = MAX(kLowC(I,J,bi,bj),1)  
         gTKE(I,J,kBottom) = gTKE(I,J,kBottom)  
      &                    - GGL90TKEbottom*c(I,J,kBottom)  
         c(I,J,kBottom) = 0. _d 0  
626         ENDDO         ENDDO
627        ENDDO        ENDDO
628    
629  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)        IF (GGL90_dirichlet) THEN
630    C      Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
631           DO j=jMin,jMax
632            DO i=iMin,iMax
633             kBottom   = MAX(kLowC(i,j,bi,bj),1)
634             GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
635         &                              - GGL90TKEbottom*c3d(i,j,kBottom)
636             c3d(i,j,kBottom) = 0. _d 0
637            ENDDO
638           ENDDO
639          ENDIF
640    
641    C     solve tri-diagonal system
642        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
643       I                        a, b, c,       I                        a3d, b3d, c3d,
644       U                        gTKE,       U                        GGL90TKE(1-OLx,1-OLy,1,bi,bj),
645       O                        errCode,       O                        errCode,
646       I                        bi, bj, myThid )       I                        bi, bj, myThid )
647  c     CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,  
648  c    I     a, b, c,        DO k=1,Nr
649  c    U     gTKE,         DO j=jMin,jMax
650  c    I     myThid )          DO i=iMin,iMax
   
 C     now update TKE  
       DO K=1,Nr  
        DO J=jMin,jMax  
         DO I=iMin,iMax  
651  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
652           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
653       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
654          ENDDO          ENDDO
655         ENDDO         ENDDO
656        ENDDO        ENDDO
# Line 538  C     impose minimum TKE to avoid numeri Line 658  C     impose minimum TKE to avoid numeri
658  C     end of time step  C     end of time step
659  C     ===============================  C     ===============================
660    
661          DO k=2,Nr
662           DO j=1,sNy
663            DO i=1,sNx
664  #ifdef ALLOW_GGL90_SMOOTH  #ifdef ALLOW_GGL90_SMOOTH
665        DO K=1,Nr           tmpVisc=
        DO J=jMin,jMax  
         DO I=iMin,iMax  
          tmpdiffKrS=  
666       &  (       &  (
667       &   p4 *  GGL90viscAr(i  ,j  ,k,bi,bj) * mskCor(i  ,j  ,bi,bj)       &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
668       &  +p8 *( GGL90viscAr(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
669       &       + GGL90viscAr(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)       &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
670       &       + GGL90viscAr(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)       &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
671       &       + GGL90viscAr(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))       &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
672       &  +p16*( GGL90viscAr(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)       &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
673       &       + GGL90viscAr(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)       &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
674       &       + GGL90viscAr(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)       &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
675       &       + GGL90viscAr(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))       &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
676       &  )       &  )
677       & /(p4       & /(p4
678       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)       &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
# Line 564  C     =============================== Line 684  C     ===============================
684       &       +       maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)       &       +       maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
685       &       +       maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))       &       +       maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))
686       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)       &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
687       &   /TKEPrandtlNumber(i,j,k)  #else
688           GGL90diffKrS(I,J,K,bi,bj)= MAX( tmpdiffKrS , diffKrNrT(k) )           tmpVisc = GGL90visctmp(i,j,k)
689    #endif
690             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
691             GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrS(k) )
692          ENDDO          ENDDO
693         ENDDO         ENDDO
694        ENDDO        ENDDO
695    
696          DO k=2,Nr
697           DO j=1,sNy
698            DO i=1,sNx+1
699    #ifdef ALLOW_GGL90_SMOOTH
700            tmpVisc =
701         & (
702         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
703         &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
704         &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
705         &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
706         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
707         &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
708         &  )
709         & /(p4 * 2. _d 0
710         &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
711         &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
712         &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
713         &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
714         &  )
715         &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
716         &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
717    #else
718            tmpVisc = _maskW(i,j,k,bi,bj) *
719         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
720         &                            +GGL90visctmp(i-1,j,k))
721         &                   )
722  #endif  #endif
723            tmpVisc = MIN( tmpVisc , GGL90viscMax )
724            GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
725            ENDDO
726           ENDDO
727          ENDDO
728    
729          DO k=2,Nr
730           DO j=1,sNy+1
731            DO i=1,sNx
732    #ifdef ALLOW_GGL90_SMOOTH
733            tmpVisc =
734         & (
735         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
736         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))
737         &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
738         &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
739         &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
740         &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))
741         &  )
742         & /(p4 * 2. _d 0
743         &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
744         &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
745         &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
746         &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))
747         &  )
748         &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
749         &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
750    #else
751            tmpVisc = _maskS(i,j,k,bi,bj) *
752         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
753         &                            +GGL90visctmp(i,j-1,k))
754         &                   )
755    
756    #endif
757            tmpVisc = MIN( tmpVisc , GGL90viscMax )
758            GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
759            ENDDO
760           ENDDO
761          ENDDO
762    
763  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
764        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
765           CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',          CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
766       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
767           CALL DIAGNOSTICS_FILL( GGL90viscAr,'GGL90Ar ',          CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
768       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
769           CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',          CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
770       &                          0,Nr, 1, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
771           CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',          CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
772       &                          0,Nr, 2, bi, bj, myThid )       &                         0,Nr, 1, bi, bj, myThid )
773           CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',          CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
774       &                          0,Nr, 2, bi, bj, myThid )       &                         0,Nr, 2, bi, bj, myThid )
775            CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
776         &                         0,Nr, 2, bi, bj, myThid )
777    
778            kp1 = MIN(Nr,kSurf+1)
779            DO j=jMin,jMax
780             DO i=iMin,iMax
781    C     diagnose surface flux of TKE
782              surf_flx_tke(i,j) =(GGL90TKE(i,j,kSurf,bi,bj)-
783         &                        GGL90TKE(i,j,kp1,bi,bj))
784         &        *recip_drF(kSurf)*recip_hFacC(i,j,kSurf,bi,bj)
785         &        *KappaE(i,j,kp1)
786             ENDDO
787            ENDDO
788            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90flx',
789         &                         0, 1, 2, bi, bj, myThid )
790    
791            k=kSurf
792            DO j=jMin,jMax
793             DO i=iMin,iMax
794    C     diagnose work done by the wind
795              surf_flx_tke(i,j) =
796         &      halfRL*( surfaceForcingU(i,  j,bi,bj)*uVel(i  ,j,k,bi,bj)
797         &              +surfaceForcingU(i+1,j,bi,bj)*uVel(i+1,j,k,bi,bj))
798         &    + halfRL*( surfaceForcingV(i,j,  bi,bj)*vVel(i,j  ,k,bi,bj)
799         &              +surfaceForcingV(i,j+1,bi,bj)*vVel(i,j+1,k,bi,bj))
800             ENDDO
801            ENDDO
802            CALL DIAGNOSTICS_FILL( surf_flx_tke,'GGL90tau',
803         &                         0, 1, 2, bi, bj, myThid )
804    
805        ENDIF        ENDIF
806  #endif  #endif /* ALLOW_DIAGNOSTICS */
807    
808  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
809    
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.22