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