/[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.1 by mlosch, Thu Sep 16 11:27:18 2004 UTC revision 1.26 by jmc, Thu Aug 14 16:42:35 2014 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     *==========================================================*
15  C     | SUBROUTINE GGL90_CALC                                    |  C     | SUBROUTINE GGL90_CALC                                    |
16  C     | o Compute all GGL90 fields defined in GGL90.h            |  C     | o Compute all GGL90 fields defined in GGL90.h            |
17  C     |==========================================================|  C     *==========================================================*
18  C     | Equation numbers refer to                                |  C     | Equation numbers refer to                                |
19  C     | Gaspar et al. (1990), JGR 95 (C9), pp 16,179             |  C     | Gaspar et al. (1990), JGR 95 (C9), pp 16,179             |
20  C     | Some parts of the implementation follow Blanke and       |  C     | Some parts of the implementation follow Blanke and       |
21  C     | Delecuse (1993), JPO, and OPA code, in particular the    |  C     | Delecuse (1993), JPO, and OPA code, in particular the    |
22  C     | computation of the                                       |  C     | computation of the                                       |
23  C     | mixing length = max(min(lk,depth),lkmin)                 |  C     | mixing length = max(min(lk,depth),lkmin)                 |
24  C     \==========================================================/  C     *==========================================================*
       IMPLICIT NONE  
 C  
 C--------------------------------------------------------------------  
25    
26  C global parameters updated by ggl90_calc  C global parameters updated by ggl90_calc
27  C     GGL90TKE     - sub-grid turbulent kinetic energy           (m^2/s^2)  C     GGL90TKE     :: sub-grid turbulent kinetic energy          (m^2/s^2)
28  C     GGL90viscAz  - GGL90 eddy viscosity coefficient              (m^2/s)  C     GGL90viscAz  :: GGL90 eddy viscosity coefficient             (m^2/s)
29  C     GGL90diffKzT - GGL90 diffusion coefficient for temperature   (m^2/s)  C     GGL90diffKzT :: GGL90 diffusion coefficient for temperature  (m^2/s)
 C  
30  C \ev  C \ev
31    
32  C !USES: ============================================================  C !USES: ============================================================
33          IMPLICIT NONE
34  #include "SIZE.h"  #include "SIZE.h"
35  #include "EEPARAMS.h"  #include "EEPARAMS.h"
36  #include "PARAMS.h"  #include "PARAMS.h"
# Line 43  C !USES: =============================== Line 40  C !USES: ===============================
40  #include "GRID.h"  #include "GRID.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     myTime - Current time in simulation  C     sigmaR :: Vertical gradient of iso-neutral density
46    C     myTime :: Current time in simulation
47    C     myIter :: Current time-step number
48    C     myThid :: My Thread Id number
49        INTEGER bi, bj        INTEGER bi, bj
50        INTEGER myThid        _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
51        _RL     myTime        _RL     myTime
52          INTEGER myIter
53          INTEGER myThid
54    CEOP
55    
56  #ifdef ALLOW_GGL90  #ifdef ALLOW_GGL90
57    
58  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
59  c Local constants  C Local constants
60  C     iMin, iMax, jMin, jMax, I, J - array computation indices  C     iMin,iMax,jMin,jMax :: index boundaries of computation domain
61  C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices  C     i, j, k, kp1,km1 :: array computation indices
62  C     ab15, ab05       - weights for implicit timestepping  C     kSurf, kBottom   :: vertical indices of domain boundaries
63  C     uStarSquare      - square of friction velocity  C     explDissFac      :: explicit Dissipation Factor (in [0-1])
64  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     implDissFac      :: implicit Dissipation Factor (in [0-1])
65  C     Nsquare          - squared buoyancy freqency  C     uStarSquare      :: square of friction velocity
66  C     RiNumber         - local Richardson number  C     verticalShear    :: (squared) vertical shear of horizontal velocity
67  C     KappaM           - (local) viscosity parameter (eq.10)  C     Nsquare          :: squared buoyancy freqency
68  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     RiNumber         :: local Richardson number
69  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     KappaM           :: (local) viscosity parameter (eq.10)
70  C     buoyFreq         - buoyancy freqency  C     KappaH           :: (local) diffusivity parameter for temperature (eq.11)
71  C     TKEdissipation   - dissipation of TKE  C     KappaE           :: (local) diffusivity parameter for TKE (eq.15)
72  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     TKEdissipation   :: dissipation of TKE
73  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C     GGL90mixingLength:: mixing length of scheme following Banke+Delecuse
74  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C         rMixingLength:: inverse of mixing length
75  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C     totalDepth       :: thickness of water column (inverse of recip_Rcol)
76  C     gTKE             - right hand side of implicit equation  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
82        _RL     KappaM, KappaH        _RL     KappaM, KappaH
83        _RL     Nsquare  c     _RL     Nsquare
84          _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
85        _RL     deltaTggl90        _RL     deltaTggl90
86        _RL     SQRTTKE  c     _RL     SQRTTKE
87          _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
88        _RL     RiNumber        _RL     RiNumber
89        _RL     TKEdissipation        _RL     TKEdissipation
90        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
91          _RL     MaxLength, tmpmlx, tmpVisc
92        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
93        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
94          _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
95          _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
96        _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)  
97        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90visctmp     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
99  C     tri-diagonal matrix  C-    tri-diagonal matrix
100        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c3d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103  CEOP        INTEGER errCode
104    #ifdef ALLOW_GGL90_HORIZDIFF
105    C     hFac     :: fractional thickness of W-cell
106    C     xA, yA   :: area of lateral faces
107    C     dfx, dfy :: diffusive flux across lateral faces
108    C     gTKE     :: right hand side of diffusion equation
109          _RL     hFac
110          _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
111          _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112          _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113          _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114          _RL    gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115    #endif /* ALLOW_GGL90_HORIZDIFF */
116    #ifdef ALLOW_GGL90_SMOOTH
117          _RL p4, p8, p16
118          p4=0.25 _d 0
119          p8=0.125 _d 0
120          p16=0.0625 _d 0
121    #endif
122        iMin = 2-OLx        iMin = 2-OLx
123        iMax = sNx+OLx-1        iMax = sNx+OLx-1
124        jMin = 2-OLy        jMin = 2-OLy
125        jMax = sNy+OLy-1        jMax = sNy+OLy-1
126    
127  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
128        deltaTggl90 = deltaTtracer        deltaTggl90 = dTtracerLev(1)
129  C      
130        kSurf = 1        kSurf = 1
131  C     implicit timestepping weights for dissipation  C     explicit/implicit timestepping weights for dissipation
132        ab15 =  1.5 _d 0        explDissFac = 0. _d 0
133        ab05 = -0.5 _d 0        implDissFac = 1. _d 0 - explDissFac
       ab15 =  1. _d 0  
       ab05 =  0. _d 0  
134    
135  C     Initialize local fields  C     Initialize local fields
136        DO K = 1, Nr        DO k = 1, Nr
137         DO J=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
138          DO I=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
139           gTKE(I,J,K)              = 0. _d 0           KappaE(i,j,k)            = 0. _d 0
140           KappaE(I,J,K)            = 0. _d 0           TKEPrandtlNumber(i,j,k)  = 1. _d 0
141           TKEPrandtlNumber(I,J,K)  = 0. _d 0           GGL90mixingLength(i,j,k) = GGL90mixingLengthMin
142           GGL90mixingLength(I,J,K) = 0. _d 0           GGL90visctmp(i,j,k)      = 0. _d 0
143          ENDDO  #ifndef SOLVE_DIAGONAL_LOWMEMORY
144         ENDDO               a3d(i,j,k) = 0. _d 0
145        ENDDO           b3d(i,j,k) = 1. _d 0
146        DO J=1-Oly,sNy+Oly           c3d(i,j,k) = 0. _d 0
147         DO I=1-Olx,sNx+Olx  #endif
148          rhoK    (I,J) = 0. _d 0          ENDDO
149          rhoKm1  (I,J) = 0. _d 0         ENDDO
150          totalDepth(I,J)   = 0. _d 0        ENDDO
151          IF ( recip_Rcol(I,J,bi,bj) .NE. 0. )        DO j=1-OLy,sNy+OLy
152       &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)         DO i=1-OLx,sNx+OLx
153            totalDepth(i,j)    = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
154            rMixingLength(i,j,1) = 0. _d 0
155            mxLength_Dn(i,j,1) = GGL90mixingLengthMin
156            SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
157         ENDDO         ENDDO
158        ENDDO        ENDDO
159    
160  C     start k-loop  C     start k-loop
161        DO K = 2, Nr        DO k = 2, Nr
162         Km1 = K-1  c      km1 = k-1
163         Kp1 = MIN(Nr,K+1)  c      kp1 = MIN(Nr,k+1)
164         CALL FIND_RHO(         DO j=jMin,jMax
165       I      bi, bj, iMin, iMax, jMin, jMax, Km1, K,          DO i=iMin,iMax
166       I      theta, salt,           SQRTTKE(i,j,k)=SQRT( GGL90TKE(i,j,k,bi,bj) )
167       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  
168  C     buoyancy frequency  C     buoyancy frequency
169  C           Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
170           Nsquare = - gravity*recip_rhoConst*recip_drC(K)       &                  * sigmaR(i,j,k)
171       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
172    c         tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
173    c     &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
174    c     &        *recip_drC(k)
175    c         tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
176    c     &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
177    c     &        *recip_drC(k)
178    c         verticalShear = tempU*tempU + tempV*tempV
179    c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
180    cC     compute Prandtl number (always greater than 0)
181    c         prTemp = 1. _d 0
182    c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
183    c         TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
184    C     mixing length
185             GGL90mixingLength(i,j,k) = SQRTTWO *
186         &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
187            ENDDO
188           ENDDO
189          ENDDO
190    
191    C- ensure mixing between first and second level
192          IF (mxlSurfFlag) THEN
193           DO j=jMin,jMax
194            DO i=iMin,iMax
195             GGL90mixingLength(i,j,2)=drF(1)
196            ENDDO
197           ENDDO
198          ENDIF
199    
200    C- Impose upper and lower bound for mixing length
201          IF ( mxlMaxFlag .EQ. 0 ) THEN
202    
203           DO k=2,Nr
204            DO j=jMin,jMax
205             DO i=iMin,iMax
206              MaxLength=totalDepth(i,j)
207              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
208         &                                   MaxLength)
209             ENDDO
210            ENDDO
211           ENDDO
212    
213           DO k=2,Nr
214            DO j=jMin,jMax
215             DO i=iMin,iMax
216              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
217         &                                   GGL90mixingLengthMin)
218              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
219             ENDDO
220            ENDDO
221           ENDDO
222    
223          ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
224    
225           DO k=2,Nr
226            DO j=jMin,jMax
227             DO i=iMin,iMax
228              MaxLength=MIN(Ro_surf(i,j,bi,bj)-rF(k),rF(k)-R_low(i,j,bi,bj))
229    c         MaxLength=MAX(MaxLength,20. _d 0)
230              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
231         &                                   MaxLength)
232             ENDDO
233            ENDDO
234           ENDDO
235    
236           DO k=2,Nr
237            DO j=jMin,jMax
238             DO i=iMin,iMax
239              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
240         &                                   GGL90mixingLengthMin)
241              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
242             ENDDO
243            ENDDO
244           ENDDO
245    
246          ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
247    
248           DO k=2,Nr
249            DO j=jMin,jMax
250             DO i=iMin,iMax
251              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
252         &        GGL90mixingLength(i,j,k-1)+drF(k-1))
253             ENDDO
254            ENDDO
255           ENDDO
256           DO j=jMin,jMax
257            DO i=iMin,iMax
258              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
259         &       GGL90mixingLengthMin+drF(Nr))
260            ENDDO
261           ENDDO
262           DO k=Nr-1,2,-1
263            DO j=jMin,jMax
264             DO i=iMin,iMax
265              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
266         &        GGL90mixingLength(i,j,k+1)+drF(k))
267             ENDDO
268            ENDDO
269           ENDDO
270    
271           DO k=2,Nr
272            DO j=jMin,jMax
273             DO i=iMin,iMax
274              GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
275         &                                   GGL90mixingLengthMin)
276              rMixingLength(i,j,k) = 1. _d 0 / GGL90mixingLength(i,j,k)
277             ENDDO
278            ENDDO
279           ENDDO
280    
281          ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
282    
283           DO k=2,Nr
284            DO j=jMin,jMax
285             DO i=iMin,iMax
286              mxLength_Dn(i,j,k) = MIN(GGL90mixingLength(i,j,k),
287         &        mxLength_Dn(i,j,k-1)+drF(k-1))
288             ENDDO
289            ENDDO
290           ENDDO
291           DO j=jMin,jMax
292            DO i=iMin,iMax
293              GGL90mixingLength(i,j,Nr) = MIN(GGL90mixingLength(i,j,Nr),
294         &       GGL90mixingLengthMin+drF(Nr))
295            ENDDO
296           ENDDO
297           DO k=Nr-1,2,-1
298            DO j=jMin,jMax
299             DO i=iMin,iMax
300              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
301         &        GGL90mixingLength(i,j,k+1)+drF(k))
302             ENDDO
303            ENDDO
304           ENDDO
305    
306           DO k=2,Nr
307            DO j=jMin,jMax
308             DO i=iMin,iMax
309              GGL90mixingLength(i,j,k) = MIN(GGL90mixingLength(i,j,k),
310         &                                  mxLength_Dn(i,j,k))
311              tmpmlx = SQRT( GGL90mixingLength(i,j,k)*mxLength_Dn(i,j,k) )
312              tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
313              rMixingLength(i,j,k) = 1. _d 0 / tmpmlx
314             ENDDO
315            ENDDO
316           ENDDO
317    
318          ELSE
319           STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
320          ENDIF
321    
322    C- Impose minimum mixing length (to avoid division by zero)
323    c      DO k=2,Nr
324    c      DO j=jMin,jMax
325    c       DO i=iMin,iMax
326    c        GGL90mixingLength(i,j,k) = MAX(GGL90mixingLength(i,j,k),
327    c    &        GGL90mixingLengthMin)
328    c        rMixingLength(i,j,k) = 1. _d 0 /GGL90mixingLength(i,j,k)
329    c       ENDDO
330    c      ENDDO
331    c     ENDDO
332    
333          DO k=2,Nr
334           km1 = k-1
335    
336    #ifdef ALLOW_GGL90_HORIZDIFF
337          IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
338    C     horizontal diffusion of TKE (requires an exchange in
339    C      do_fields_blocking_exchanges)
340    C     common factors
341            DO j=1-OLy,sNy+OLy
342             DO i=1-OLx,sNx+OLx
343              xA(i,j) = _dyG(i,j,bi,bj)*drC(k)*
344         &                 (min(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
345         &                  min(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
346              yA(i,j) = _dxG(i,j,bi,bj)*drC(k)*
347         &                 (min(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
348         &                  min(.5 _d 0,_hFacS(i,j,k  ,bi,bj) ) )
349             ENDDO
350            ENDDO
351    C     Compute diffusive fluxes
352    C     ... across x-faces
353            DO j=1-OLy,sNy+OLy
354             dfx(1-OLx,j)=0. _d 0
355             DO i=1-OLx+1,sNx+OLx
356              dfx(i,j) = -GGL90diffTKEh*xA(i,j)
357         &      *_recip_dxC(i,j,bi,bj)
358         &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i-1,j,k,bi,bj))
359    #ifdef ISOTROPIC_COS_SCALING
360         &      *CosFacU(j,bi,bj)
361    #endif /* ISOTROPIC_COS_SCALING */
362             ENDDO
363            ENDDO
364    C     ... across y-faces
365            DO i=1-OLx,sNx+OLx
366             dfy(i,1-OLy)=0. _d 0
367            ENDDO
368            DO j=1-OLy+1,sNy+OLy
369             DO i=1-OLx,sNx+OLx
370              dfy(i,j) = -GGL90diffTKEh*yA(i,j)
371         &      *_recip_dyC(i,j,bi,bj)
372         &      *(GGL90TKE(i,j,k,bi,bj)-GGL90TKE(i,j-1,k,bi,bj))
373    #ifdef ISOTROPIC_COS_SCALING
374         &      *CosFacV(j,bi,bj)
375    #endif /* ISOTROPIC_COS_SCALING */
376             ENDDO
377            ENDDO
378    C     Compute divergence of fluxes
379            DO j=1-OLy,sNy+OLy-1
380             DO i=1-OLx,sNx+OLx-1
381              hFac = min(.5 _d 0,_hFacC(i,j,k-1,bi,bj) ) +
382         &          min(.5 _d 0,_hFacC(i,j,k  ,bi,bj) )
383              gTKE(i,j) = 0.0
384              if ( hFac .ne. 0.0 )
385         &      gTKE(i,j) = -recip_drC(k)*recip_rA(i,j,bi,bj)/hFac
386         &         *((dfx(i+1,j)-dfx(i,j))
387         &          +(dfy(i,j+1)-dfy(i,j)) )
388             ENDDO
389            ENDDO
390    C      end if GGL90diffTKEh .eq. 0.
391           ENDIF
392    #endif /* ALLOW_GGL90_HORIZDIFF */
393    
394           DO j=jMin,jMax
395            DO i=iMin,iMax
396  C     vertical shear term (dU/dz)^2+(dV/dz)^2  C     vertical shear term (dU/dz)^2+(dV/dz)^2
397           tempu= .5*(  uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)           tempU= .5 _d 0*( uVel(i,j,km1,bi,bj)+uVel(i+1,j,km1,bi,bj)
398       &             - (uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )       &                 -( uVel(i,j,k  ,bi,bj)+uVel(i+1,j,k  ,bi,bj)) )
399       &        *recip_drC(K)       &        *recip_drC(k)
400           tempv= .5*(  vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)           tempV= .5 _d 0*( vVel(i,j,km1,bi,bj)+vVel(i,j+1,km1,bi,bj)
401       &             - (vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )       &                 -( vVel(i,j,k  ,bi,bj)+vVel(i,j+1,k  ,bi,bj)) )
402       &        *recip_drC(K)       &        *recip_drC(k)
403           verticalShear = tempU*tempU + tempV*tempV           verticalShear = tempU*tempU + tempV*tempV
404           RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)           RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
405  C     compute Prandtl number (always greater than 0)  C     compute Prandtl number (always greater than 0)
406           prTemp = 1. _d 0           prTemp = 1. _d 0
407           IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
408           TKEPrandtlNumber(I,J,K) = MIN(10.,prTemp)           TKEPrandtlNumber(i,j,k) = MIN(10. _d 0,prTemp)
409  C     mixing length  c         TKEPrandtlNumber(i,j,k) = 1. _d 0
410           GGL90mixingLength(I,J,K) =  
411       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )  C     viscosity and diffusivity
412  C     impose upper bound for mixing length (total depth)           KappaM = GGL90ck*GGL90mixingLength(i,j,k)*SQRTTKE(i,j,k)
413           GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),           GGL90visctmp(i,j,k) = MAX(KappaM,diffKrNrT(k))
414       &        totalDepth(I,J))       &                            * maskC(i,j,k,bi,bj)
415  C     impose minimum mixing length (to avoid division by zero)  c        note: storing GGL90visctmp like this, and using it later to compute
416           GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),  c              GGL9rdiffKr etc. is robust in case of smoothing (e.g. see OPA)
417       &        GGL90mixingLengthMin)           KappaM = MAX(KappaM,viscArNr(k)) * maskC(i,j,k,bi,bj)
418  C     viscosity of last timestep           KappaH = KappaM/TKEPrandtlNumber(i,j,k)
419           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE           KappaE(i,j,k) = GGL90alpha * KappaM * maskC(i,j,k,bi,bj)
420           KappaE(I,J,K) = KappaM*GGL90alpha  
421  C     dissipation term  C     dissipation term
422           TKEdissipation = ab05*GGL90ceps           TKEdissipation = explDissFac*GGL90ceps
423       &        *SQRTTKE/GGL90mixingLength(I,J,K)       &        *SQRTTKE(i,j,k)*rMixingLength(i,j,k)
424       &        *GGL90TKE(I,J,K,bi,bj)             &        *GGL90TKE(i,j,k,bi,bj)
425  C     sum up contributions to form the right hand side  C     partial update with sum of explicit contributions
426           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)           GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
427       &        + deltaTggl90*(       &        + deltaTggl90*(
428       &        + KappaM*verticalShear       &        + KappaM*verticalShear
429       &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)       &        - KappaH*Nsquare(i,j,k)
430       &        - TKEdissipation       &        - TKEdissipation
431       &        )       &        )
432          ENDDO            ENDDO
433         ENDDO         ENDDO
434        ENDDO      
435  C      #ifdef ALLOW_GGL90_HORIZDIFF
436  C     Implicit time step to update TKE for k=1,Nr; TKE(Nr+1)=0 by default         IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
437  C      C--    Add horiz. diffusion tendency
438            DO j=jMin,jMax
439             DO i=iMin,iMax
440              GGL90TKE(i,j,k,bi,bj) = GGL90TKE(i,j,k,bi,bj)
441         &                          + gTKE(i,j)*deltaTggl90
442             ENDDO
443            ENDDO
444           ENDIF
445    #endif /* ALLOW_GGL90_HORIZDIFF */
446    
447    C--   end of k loop
448          ENDDO
449    
450    C     ============================================
451    C     Implicit time step to update TKE for k=1,Nr;
452    C     TKE(Nr+1)=0 by default
453    C     ============================================
454  C     set up matrix  C     set up matrix
455  C--   Old aLower  C--   Lower diagonal
456        DO j=jMin,jMax        DO j=jMin,jMax
457         DO i=iMin,iMax         DO i=iMin,iMax
458           a(i,j,1) = 0. _d 0           a3d(i,j,1) = 0. _d 0
459         ENDDO         ENDDO
460        ENDDO        ENDDO
461        DO k=2,Nr        DO k=2,Nr
462         km1=MAX(1,k-1)         km1=MAX(2,k-1)
463         DO j=jMin,jMax         DO j=jMin,jMax
464          DO i=iMin,iMax          DO i=iMin,iMax
465           a(i,j,k) = -deltaTggl90  C-    We keep recip_hFacC in the diffusive flux calculation,
466       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)  C-    but no hFacC in TKE volume control
467       &        *.5*(KappaE(i,j, k )+KappaE(i,j,km1))  C-    No need for maskC(k-1) with recip_hFacC(k-1)
468       &        *recip_drC(k)           a3d(i,j,k) = -deltaTggl90
469            IF (recip_hFacI(i,j,km1,bi,bj).EQ.0.) a(i,j,k)=0.       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
470         &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
471         &        *recip_drC(k)*maskC(i,j,k,bi,bj)
472          ENDDO          ENDDO
473         ENDDO         ENDDO
474        ENDDO        ENDDO
475    C--   Upper diagonal
 C--   Old aUpper  
476        DO j=jMin,jMax        DO j=jMin,jMax
477         DO i=iMin,iMax         DO i=iMin,iMax
478           c(i,j,1)  = 0. _d 0           c3d(i,j,1)  = 0. _d 0
          c(i,j,Nr) = 0. _d 0  
479         ENDDO         ENDDO
480        ENDDO        ENDDO
 CML      DO k=1,Nr-1  
481        DO k=2,Nr        DO k=2,Nr
        kp1=min(Nr,k+1)  
482         DO j=jMin,jMax         DO j=jMin,jMax
483          DO i=iMin,iMax          DO i=iMin,iMax
484            c(i,j,k) = -deltaTggl90            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
485       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)  C-    We keep recip_hFacC in the diffusive flux calculation,
486       &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))  C-    but no hFacC in TKE volume control
487       &        *recip_drC(k)  C-    No need for maskC(k) with recip_hFacC(k)
488  C          IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0.) c(i,j,k)=0.            c3d(i,j,k) = -deltaTggl90
489         &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
490         &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
491         &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
492          ENDDO          ENDDO
493         ENDDO         ENDDO
494        ENDDO        ENDDO
495    C--   Center diagonal
 C--   Old aCenter  
496        DO k=1,Nr        DO k=1,Nr
497           km1 = MAX(k-1,1)
498         DO j=jMin,jMax         DO j=jMin,jMax
499          DO i=iMin,iMax          DO i=iMin,iMax
500            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)
501       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + implDissFac*deltaTggl90*GGL90ceps*SQRTTKE(i,j,k)
502       &        /GGL90mixingLength(I,J,K)       &        * rMixingLength(i,j,k)
503         &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
504           ENDDO           ENDDO
505         ENDDO         ENDDO
506        ENDDO        ENDDO
507  C     end set up matrix  C     end set up matrix
508    
 C  
509  C     Apply boundary condition  C     Apply boundary condition
510  C            kp1 = MIN(Nr,kSurf+1)
511        DO J=jMin,jMax        DO j=jMin,jMax
512         DO I=iMin,iMax         DO i=iMin,iMax
513  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
514          uStarSquare = SQRT(          uStarSquare = SQRT(
515       &         ( .5*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(i,  j,  bi,bj)
516       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(i+1,j,  bi,bj) ) )**2
517       &       + ( .5*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(i,  j,  bi,bj)
518       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(i,  j+1,bi,bj) ) )**2
519       &                     )       &                     )
520  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
521          gTKE(I,J,kSurf) = MAX(GGL90TKEmin,GGL90m2*uStarSquare)          GGL90TKE(i,j,kSurf,bi,bj) = maskC(i,j,kSurf,bi,bj)
522       &                     *maskC(I,J,kSurf,bi,bj)       &           *MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
523            GGL90TKE(i,j,kp1,bi,bj) = GGL90TKE(i,j,kp1,bi,bj)
524         &               - a3d(i,j,kp1)*GGL90TKE(i,j,kSurf,bi,bj)
525            a3d(i,j,kp1) = 0. _d 0
526  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
527          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)          kBottom   = MAX(kLowC(i,j,bi,bj),1)
528          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          GGL90TKE(i,j,kBottom,bi,bj) = GGL90TKE(i,j,kBottom,bi,bj)
529       &       - GGL90TKEbottom*c(I,J,kBottom)       &                              - GGL90TKEbottom*c3d(i,j,kBottom)
530         ENDDO          c3d(i,j,kBottom) = 0. _d 0
531        ENDDO             ENDDO
532  C        ENDDO
533  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  
534  C  C     solve tri-diagonal system
535        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
536       I     a, b, c,       I                        a3d, b3d, c3d,
537       U     gTKE,       U                        GGL90TKE(1-OLx,1-OLy,1,bi,bj),
538       I     myThid )       O                        errCode,
539  C       I                        bi, bj, myThid )
540  C     now update TKE  
541  C            DO k=1,Nr
542        DO K=1,Nr         DO j=jMin,jMax
543         DO J=jMin,jMax          DO i=iMin,iMax
         DO I=iMin,iMax  
544  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
545           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
546       &        * maskC(I,J,K,bi,bj)       &                  *MAX( GGL90TKE(i,j,k,bi,bj), GGL90TKEmin )
547  C          ENDDO
548           ENDDO
549          ENDDO
550    
551  C     end of time step  C     end of time step
552  C  C     ===============================
553    
554          DO k=2,Nr
555           DO j=1,sNy
556            DO i=1,sNx
557    #ifdef ALLOW_GGL90_SMOOTH
558             tmpVisc=
559         &  (
560         &   p4 *  GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
561         &  +p8 *( GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
562         &       + GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
563         &       + GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
564         &       + GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
565         &  +p16*( GGL90visctmp(i+1,j+1,k) * mskCor(i+1,j+1,bi,bj)
566         &       + GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj)
567         &       + GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
568         &       + GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj))
569         &  )
570         & /(p4
571         &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
572         &       +       maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
573         &       +       maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
574         &       +       maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
575         &  +p16*(       maskC(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)
576         &       +       maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)
577         &       +       maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
578         &       +       maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))
579         &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
580    #else
581             tmpVisc = GGL90visctmp(i,j,k)
582    #endif
583             tmpVisc = MIN(tmpVisc/TKEPrandtlNumber(i,j,k),GGL90diffMax)
584             GGL90diffKr(i,j,k,bi,bj)= MAX( tmpVisc , diffKrNrT(k) )
585            ENDDO
586           ENDDO
587          ENDDO
588    
589          DO k=2,Nr
590           DO j=1,sNy
591            DO i=1,sNx+1
592    #ifdef ALLOW_GGL90_SMOOTH
593            tmpVisc =
594         & (
595         &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
596         &       +GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj))
597         &  +p8 *(GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
598         &       +GGL90visctmp(i-1,j+1,k) * mskCor(i-1,j+1,bi,bj)
599         &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj)
600         &       +GGL90visctmp(i  ,j+1,k) * mskCor(i  ,j+1,bi,bj))
601         &  )
602         & /(p4 * 2. _d 0
603         &  +p8 *(      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
604         &       +      maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
605         &       +      maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
606         &       +      maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
607         &  )
608         &  *maskC(i  ,j,k,bi,bj)*mskCor(i  ,j,bi,bj)
609         &  *maskC(i-1,j,k,bi,bj)*mskCor(i-1,j,bi,bj)
610    #else
611            tmpVisc = _maskW(i,j,k,bi,bj) *
612         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
613         &                            +GGL90visctmp(i-1,j,k))
614         &                   )
615    #endif
616            tmpVisc = MIN( tmpVisc , GGL90viscMax )
617            GGL90viscArU(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
618          ENDDO          ENDDO
619         ENDDO         ENDDO
620        ENDDO            ENDDO
621  C  
622  C     compute viscosity coefficients        DO k=2,Nr
623  C             DO j=1,sNy+1
624        DO K=2,Nr          DO i=1,sNx
625         DO J=jMin,jMax  #ifdef ALLOW_GGL90_SMOOTH
626          DO I=iMin,iMax          tmpVisc =
627  C     Eq. (11), (18) and (21)       & (
628           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*       &   p4 *(GGL90visctmp(i  ,j  ,k) * mskCor(i  ,j  ,bi,bj)
629       &                  SQRT( GGL90TKE(I,J,K,bi,bj) )       &       +GGL90visctmp(i  ,j-1,k) * mskCor(i  ,j-1,bi,bj))
630           KappaH = KappaM/TKEPrandtlNumber(I,J,K)       &  +p8 *(GGL90visctmp(i-1,j  ,k) * mskCor(i-1,j  ,bi,bj)
631  C     Set a minium (= background) value       &       +GGL90visctmp(i-1,j-1,k) * mskCor(i-1,j-1,bi,bj)
632           KappaM = MAX(KappaM,viscAr)       &       +GGL90visctmp(i+1,j  ,k) * mskCor(i+1,j  ,bi,bj)
633           KappaH = MAX(KappaH,diffKrT)       &       +GGL90visctmp(i+1,j-1,k) * mskCor(i+1,j-1,bi,bj))
634  C     Set a maximum and mask land point       &  )
635           GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)       & /(p4 * 2. _d 0
636       &        * maskC(I,J,K,bi,bj)       &  +p8 *(      maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
637           GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)       &       +      maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj)
638       &        * maskC(I,J,K,bi,bj)       &       +      maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
639          ENDDO       &       +      maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj))
640         ENDDO       &  )
641  C     end third k-loop       &   *maskC(i,j  ,k,bi,bj)*mskCor(i,j  ,bi,bj)
642        ENDDO           &   *maskC(i,j-1,k,bi,bj)*mskCor(i,j-1,bi,bj)
643    #else
644            tmpVisc = _maskS(i,j,k,bi,bj) *
645         &                   (.5 _d 0*(GGL90visctmp(i,j,k)
646         &                            +GGL90visctmp(i,j-1,k))
647         &                   )
648    
649    #endif
650            tmpVisc = MIN( tmpVisc , GGL90viscMax )
651            GGL90viscArV(i,j,k,bi,bj) = MAX( tmpVisc, viscArNr(k) )
652            ENDDO
653           ENDDO
654          ENDDO
655    
656    #ifdef ALLOW_DIAGNOSTICS
657          IF ( useDiagnostics ) THEN
658             CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
659         &                          0,Nr, 1, bi, bj, myThid )
660             CALL DIAGNOSTICS_FILL( GGL90viscArU,'GGL90ArU',
661         &                          0,Nr, 1, bi, bj, myThid )
662             CALL DIAGNOSTICS_FILL( GGL90viscArV,'GGL90ArV',
663         &                          0,Nr, 1, bi, bj, myThid )
664             CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
665         &                          0,Nr, 1, bi, bj, myThid )
666             CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
667         &                          0,Nr, 2, bi, bj, myThid )
668             CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
669         &                          0,Nr, 2, bi, bj, myThid )
670          ENDIF
671    #endif
672    
673  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
674    
675        RETURN        RETURN
676        END        END
   
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.26
  ViewVC Help
Powered by ViewVC 1.1.22