/[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.11 by dfer, Fri Jan 30 02:23:56 2009 UTC revision 1.15 by dfer, Sat Nov 14 14:27:56 2009 UTC
# Line 11  C !INTERFACE: ========================== Line 11  C !INTERFACE: ==========================
11       I     bi, bj, myTime, myThid )       I     bi, bj, myTime, 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 :: array indices on which to apply calculations
45  c     myTime - Current time in simulation  C     myTime :: Current time in simulation
46    C     myThid :: My Thread Id number
47        INTEGER bi, bj        INTEGER bi, bj
       INTEGER myThid  
48        _RL     myTime        _RL     myTime
49          INTEGER myThid
50    CEOP
51    
52  #ifdef ALLOW_GGL90  #ifdef ALLOW_GGL90
53    
54  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
55  c Local constants  C Local constants
56  C     iMin, iMax, jMin, jMax, I, J - array computation indices  C     iMin, iMax, jMin, jMax, I, J - array computation indices
57  C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices  C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices
58  C     ab15, ab05       - weights for implicit timestepping  C     ab15, ab05       - weights for implicit timestepping
# Line 86  c     _RL     SQRTTKE Line 84  c     _RL     SQRTTKE
84        _RL     RiNumber        _RL     RiNumber
85        _RL     TKEdissipation        _RL     TKEdissipation
86        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
87        _RL     MaxLength        _RL     MaxLength, tmpmlx
88        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
90        _RL         rMixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     rMixingLength    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
91          _RL     mxLength_Dn      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
92        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     KappaE           (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
93        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoK             (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     rhoKm1           (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     totalDepth       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     gTKE             (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
97  C     tri-diagonal matrix  C-    tri-diagonal matrix
98        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     a(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
99        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     b(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
100        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     c(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101          INTEGER errCode
102  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
103  C     xA, yA   - area of lateral faces  C-    xA, yA   - area of lateral faces
104  C     dfx, dfy - diffusive flux across lateral faces  C-    dfx, dfy - diffusive flux across lateral faces
105        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 113  C     dfx, dfy - diffusive flux across l Line 113  C     dfx, dfy - diffusive flux across l
113        p8=0.125 _d 0        p8=0.125 _d 0
114        p16=0.0625 _d 0        p16=0.0625 _d 0
115  #endif  #endif
 CEOP  
116        iMin = 2-OLx        iMin = 2-OLx
117        iMax = sNx+OLx-1        iMax = sNx+OLx-1
118        jMin = 2-OLy        jMin = 2-OLy
# Line 121  CEOP Line 120  CEOP
120    
121  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
122        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
123  C  
124        kSurf = 1        kSurf = 1
125  C     implicit timestepping weights for dissipation  C     implicit timestepping weights for dissipation
126        ab15 =  1.5 _d 0        ab15 =  1.5 _d 0
# Line 135  C     Initialize local fields Line 134  C     Initialize local fields
134          DO I=1-Olx,sNx+Olx          DO I=1-Olx,sNx+Olx
135           gTKE(I,J,K)              = 0. _d 0           gTKE(I,J,K)              = 0. _d 0
136           KappaE(I,J,K)            = 0. _d 0           KappaE(I,J,K)            = 0. _d 0
137           TKEPrandtlNumber(I,J,K)  = 0. _d 0           TKEPrandtlNumber(I,J,K)  = 1. _d 0
138           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin
              rMixingLength(I,J,K) = 0. _d 0  
139          ENDDO          ENDDO
140         ENDDO         ENDDO
141        ENDDO        ENDDO
142        DO J=1-Oly,sNy+Oly        DO J=1-Oly,sNy+Oly
143         DO I=1-Olx,sNx+Olx         DO I=1-Olx,sNx+Olx
144          rhoK    (I,J) = 0. _d 0          rhoK(I,J)          = 0. _d 0
145          rhoKm1  (I,J) = 0. _d 0          rhoKm1(I,J)        = 0. _d 0
146          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)
147            rMixingLength(i,j,1) = 0. _d 0
148            mxLength_Dn(I,J,1) = GGL90mixingLengthMin
149            SQRTTKE(i,j,1) = SQRT( GGL90TKE(i,j,1,bi,bj) )
150         ENDDO         ENDDO
151        ENDDO        ENDDO
152    
# Line 167  c      Kp1 = MIN(Nr,K+1) Line 168  c      Kp1 = MIN(Nr,K+1)
168         DO J=jMin,jMax         DO J=jMin,jMax
169          DO I=iMin,iMax          DO I=iMin,iMax
170           SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )
171  C  
172  C     buoyancy frequency  C     buoyancy frequency
 C  
173           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)
174       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)
175  cC     vertical shear term (dU/dz)^2+(dV/dz)^2  cC     vertical shear term (dU/dz)^2+(dV/dz)^2
# Line 192  C     mixing length Line 192  C     mixing length
192         ENDDO         ENDDO
193        ENDDO        ENDDO
194    
195  C- Impose upper bound for mixing length (total depth)  C- Impose upper and lower bound for mixing length
196        IF ( mxlMaxFlag .EQ. 0 ) THEN        IF ( mxlMaxFlag .EQ. 0 ) THEN
197    C-
198         DO k=2,Nr         DO k=2,Nr
199          DO J=jMin,jMax          DO J=jMin,jMax
200           DO I=iMin,iMax           DO I=iMin,iMax
201            MaxLength=totalDepth(I,J)            MaxLength=totalDepth(I,J)
202            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
203       &        MaxLength)       &                                   MaxLength)
204           ENDDO           ENDDO
205          ENDDO          ENDDO
206         ENDDO         ENDDO
207    
208           DO k=2,Nr
209            DO J=jMin,jMax
210             DO I=iMin,iMax
211              GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
212         &                                   GGL90mixingLengthMin)
213              rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)
214             ENDDO
215            ENDDO
216           ENDDO
217    
218        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
219    C-
220         DO k=2,Nr         DO k=2,Nr
221          DO J=jMin,jMax          DO J=jMin,jMax
222           DO I=iMin,iMax           DO I=iMin,iMax
223            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))
224  c         MaxLength=MAX(MaxLength,20. _d 0)  c         MaxLength=MAX(MaxLength,20. _d 0)
225            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
226       &        MaxLength)       &                                   MaxLength)
227             ENDDO
228            ENDDO
229           ENDDO
230    
231           DO k=2,Nr
232            DO J=jMin,jMax
233             DO I=iMin,iMax
234              GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
235         &                                   GGL90mixingLengthMin)
236              rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)
237           ENDDO           ENDDO
238          ENDDO          ENDDO
239         ENDDO         ENDDO
240    
241        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN        ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
242    C-
243         DO k=2,Nr         DO k=2,Nr
244          DO J=jMin,jMax          DO J=jMin,jMax
245           DO I=iMin,iMax           DO I=iMin,iMax
# Line 234  c         MaxLength=MAX(MaxLength,20. _d Line 259  c         MaxLength=MAX(MaxLength,20. _d
259           DO I=iMin,iMax           DO I=iMin,iMax
260            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),            GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
261       &        GGL90mixingLength(I,J,K+1)+drF(k))       &        GGL90mixingLength(I,J,K+1)+drF(k))
262           ENDDO           ENDDO
263          ENDDO            ENDDO
264         ENDDO         ENDDO
       ELSE  
        STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing lenght limit)'  
       ENDIF  
265    
266  C- Impose minimum mixing length (to avoid division by zero)         DO k=2,Nr
267        DO k=2,Nr          DO J=jMin,jMax
268             DO I=iMin,iMax
269              GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
270         &                                   GGL90mixingLengthMin)
271              rMixingLength(I,J,K) = 1. _d 0 / GGL90mixingLength(I,J,K)
272             ENDDO
273            ENDDO
274           ENDDO
275    
276          ELSEIF ( mxlMaxFlag .EQ. 3 ) THEN
277    C-
278           DO k=2,Nr
279            DO J=jMin,jMax
280             DO I=iMin,iMax
281              mxLength_Dn(I,J,K) = MIN(GGL90mixingLength(I,J,K),
282         &        mxLength_Dn(I,J,K-1)+drF(k-1))
283             ENDDO
284            ENDDO
285           ENDDO
286         DO J=jMin,jMax         DO J=jMin,jMax
287          DO I=iMin,iMax          DO I=iMin,iMax
288           GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),            GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),
289       &        GGL90mixingLengthMin)       &       GGL90mixingLengthMin+drF(Nr))
290           rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)          ENDDO
291           ENDDO
292           DO k=Nr-1,2,-1
293            DO J=jMin,jMax
294             DO I=iMin,iMax
295              GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
296         &        GGL90mixingLength(I,J,K+1)+drF(k))
297             ENDDO
298            ENDDO
299           ENDDO
300    
301           DO k=2,Nr
302            DO J=jMin,jMax
303             DO I=iMin,iMax
304              GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
305         &                                  mxLength_Dn(I,J,K))
306              tmpmlx = SQRT( GGL90mixingLength(I,J,K)*mxLength_Dn(I,J,K) )
307              tmpmlx = MAX( tmpmlx, GGL90mixingLengthMin)
308              rMixingLength(I,J,K) = 1. _d 0 / tmpmlx
309             ENDDO
310          ENDDO          ENDDO
311         ENDDO         ENDDO
312        ENDDO  
313          ELSE
314           STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing length limit)'
315          ENDIF
316    
317    C- Impose minimum mixing length (to avoid division by zero)
318    c      DO k=2,Nr
319    c      DO J=jMin,jMax
320    c       DO I=iMin,iMax
321    c        GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
322    c    &        GGL90mixingLengthMin)
323    c        rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)
324    c       ENDDO
325    c      ENDDO
326    c     ENDDO
327    
328        DO k=2,Nr        DO k=2,Nr
329         Km1 = K-1         Km1 = K-1
# Line 264  C     vertical shear term (dU/dz)^2+(dV/ Line 337  C     vertical shear term (dU/dz)^2+(dV/
337       &                 -( 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)) )
338       &        *recip_drC(K)       &        *recip_drC(K)
339           verticalShear = tempU*tempU + tempV*tempV           verticalShear = tempU*tempU + tempV*tempV
340           RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)           RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
341  C     compute Prandtl number (always greater than 0)  C     compute Prandtl number (always greater than 0)
342           prTemp = 1. _d 0           prTemp = 1. _d 0
343           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
344           TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)           TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)
345    c         TKEPrandtlNumber(I,J,K) = 1. _d 0
346    
347  C     viscosity and diffusivity  C     viscosity and diffusivity
348           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)
349           KappaH = KappaM/TKEPrandtlNumber(I,J,K)           KappaH = KappaM/TKEPrandtlNumber(I,J,K)
350    
351  C     Set a minium (= background) and maximum value  C     Set a minium (= background) and maximum value
352           KappaM = MAX(KappaM,viscAr)           KappaM = MAX(KappaM,viscArNr(k))
353           KappaH = MAX(KappaH,diffKrNrT(k))           KappaH = MAX(KappaH,diffKrNrT(k))
354           KappaM = MIN(KappaM,GGL90viscMax)           KappaM = MIN(KappaM,GGL90viscMax)
355           KappaH = MIN(KappaH,GGL90diffMax)           KappaH = MIN(KappaH,GGL90diffMax)
# Line 346  C     Compute divergence of fluxes Line 420  C     Compute divergence of fluxes
420       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
421       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
422       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
423       &           )       &           )*deltaTggl90
424           ENDDO           ENDDO
425          ENDDO          ENDDO
426  C       end of k-loop  C       end of k-loop
# Line 367  C--   Lower diagonal Line 441  C--   Lower diagonal
441         ENDDO         ENDDO
442        ENDDO        ENDDO
443        DO k=2,Nr        DO k=2,Nr
444         km1=max(2,k-1)         km1=MAX(2,k-1)
445         DO j=jMin,jMax         DO j=jMin,jMax
446          DO i=iMin,iMax          DO i=iMin,iMax
447    C-    We keep recip_hFacC in the diffusive flux calculation,
448    C-    but no hFacC in TKE volume control
449    C-    No need for maskC(k-1) with recip_hFacC(k-1)
450           a(i,j,k) = -deltaTggl90           a(i,j,k) = -deltaTggl90
 c     &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)  
451       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
452       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
453       &        *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k,bi,bj)
454          ENDDO          ENDDO
455         ENDDO         ENDDO
456        ENDDO        ENDDO
# Line 387  C--   Upper diagonal Line 463  C--   Upper diagonal
463        DO k=2,Nr        DO k=2,Nr
464         DO j=jMin,jMax         DO j=jMin,jMax
465          DO i=iMin,iMax          DO i=iMin,iMax
466            kp1=min(klowC(i,j,bi,bj),k+1)            kp1=MAX(1,MIN(klowC(i,j,bi,bj),k+1))
467    C-    We keep recip_hFacC in the diffusive flux calculation,
468    C-    but no hFacC in TKE volume control
469    C-    No need for maskC(k) with recip_hFacC(k)
470            c(i,j,k) = -deltaTggl90            c(i,j,k) = -deltaTggl90
 c     &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)  
471       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
472       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
473       &        *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)       &        *recip_drC(k)*maskC(i,j,k-1,bi,bj)
474          ENDDO          ENDDO
475         ENDDO         ENDDO
476        ENDDO        ENDDO
477  C--   Center diagonal  C--   Center diagonal
478        DO k=1,Nr        DO k=1,Nr
479           km1 = MAX(k-1,1)
480         DO j=jMin,jMax         DO j=jMin,jMax
481          DO i=iMin,iMax          DO i=iMin,iMax
482            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)            b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
483       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + ab15*deltaTggl90*GGL90ceps*SQRTTKE(I,J,K)
484       &        *rMixingLength(I,J,K)*maskC(i,j,k,bi,bj)       &        * rMixingLength(I,J,K)
485         &        * maskC(i,j,k,bi,bj)*maskC(i,j,km1,bi,bj)
486           ENDDO           ENDDO
487         ENDDO         ENDDO
488        ENDDO        ENDDO
489  C     end set up matrix  C     end set up matrix
490    
 C  
491  C     Apply boundary condition  C     Apply boundary condition
492  C        kp1 = MIN(Nr,kSurf+1)
493        DO J=jMin,jMax        DO J=jMin,jMax
494         DO I=iMin,iMax         DO I=iMin,iMax
495  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
# Line 423  C     estimate friction velocity uStar f Line 502  C     estimate friction velocity uStar f
502  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
503          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
504       &                     *maskC(I,J,kSurf,bi,bj)       &                     *maskC(I,J,kSurf,bi,bj)
505            gTKE(i,j,kp1) = gTKE(i,j,kp1)
506         &                - a(i,j,kp1)*gTKE(i,j,kSurf)
507            a(i,j,kp1) = 0. _d 0
508  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
509          kBottom   = MAX(kLowC(I,J,bi,bj),1)          kBottom   = MAX(kLowC(I,J,bi,bj),1)
510          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)
# Line 430  C     Dirichlet bottom boundary conditio Line 512  C     Dirichlet bottom boundary conditio
512          c(I,J,kBottom) = 0. _d 0          c(I,J,kBottom) = 0. _d 0
513         ENDDO         ENDDO
514        ENDDO        ENDDO
515  C  
516  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)
517  C        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
518        CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,       I                        a, b, c,
519       I     a, b, c,       U                        gTKE,
520       U     gTKE,       O                        errCode,
521       I     myThid )       I                        bi, bj, myThid )
522  C  c     CALL GGL90_SOLVE( bi, bj, iMin, iMax, jMin, jMax,
523    c    I     a, b, c,
524    c    U     gTKE,
525    c    I     myThid )
526    
527  C     now update TKE  C     now update TKE
 C  
528        DO K=1,Nr        DO K=1,Nr
529         DO J=jMin,jMax         DO J=jMin,jMax
530          DO I=iMin,iMax          DO I=iMin,iMax
# Line 505  C     =============================== Line 590  C     ===============================
590    
591        RETURN        RETURN
592        END        END
   
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22