/[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.4 by jmc, Sat Dec 4 00:16:49 2004 UTC revision 1.12 by jmc, Thu Oct 8 20:07:18 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)
30  C  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 :: array indices on which to apply calculations
46  c     myTime - Current time in simulation  C     myTime :: Current time in simulation
47    C     myThid :: My Thread Id number
48        INTEGER bi, bj        INTEGER bi, bj
       INTEGER myThid  
49        _RL     myTime        _RL     myTime
50          INTEGER myThid
51    CEOP
52    
53  #ifdef ALLOW_GGL90  #ifdef ALLOW_GGL90
54    
55  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
56  c Local constants  C Local constants
57  C     iMin, iMax, jMin, jMax, I, J - array computation indices  C     iMin, iMax, jMin, jMax, I, J - array computation indices
58  C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices  C     K, Kp1, km1, kSurf, kBottom  - vertical loop indices
59  C     ab15, ab05       - weights for implicit timestepping  C     ab15, ab05       - weights for implicit timestepping
60  C     uStarSquare      - square of friction velocity  C     uStarSquare      - square of friction velocity
61  C     verticalShear    - (squared) vertical shear of horizontal velocity  C     verticalShear    - (squared) vertical shear of horizontal velocity
62  C     Nsquare          - squared buoyancy freqency  C     Nsquare          - squared buoyancy freqency
63  C     RiNumber         - local Richardson number  C     RiNumber         - local Richardson number
64  C     KappaM           - (local) viscosity parameter (eq.10)  C     KappaM           - (local) viscosity parameter (eq.10)
65  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)  C     KappaH           - (local) diffusivity parameter for temperature (eq.11)
66  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)  C     KappaE           - (local) diffusivity parameter for TKE (eq.15)
 C     buoyFreq         - buoyancy freqency  
67  C     TKEdissipation   - dissipation of TKE  C     TKEdissipation   - dissipation of TKE
68  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse  C     GGL90mixingLength- mixing length of scheme following Banke+Delecuse
69    C         rMixingLength- inverse of mixing length
70  C     totalDepth       - thickness of water column (inverse of recip_Rcol)  C     totalDepth       - thickness of water column (inverse of recip_Rcol)
71  C     TKEPrandtlNumber - here, an empirical function of the Richardson number  C     TKEPrandtlNumber - here, an empirical function of the Richardson number
72  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)  C     rhoK, rhoKm1     - density at layer K and Km1 (relative to K)
# Line 78  C     gTKE             - right hand side Line 77  C     gTKE             - right hand side
77        _RL     uStarSquare        _RL     uStarSquare
78        _RL     verticalShear        _RL     verticalShear
79        _RL     KappaM, KappaH        _RL     KappaM, KappaH
80        _RL     Nsquare  c     _RL     Nsquare
81          _RL     Nsquare(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
82        _RL     deltaTggl90        _RL     deltaTggl90
83        _RL     SQRTTKE  c     _RL     SQRTTKE
84          _RL     SQRTTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
85        _RL     RiNumber        _RL     RiNumber
86        _RL     TKEdissipation        _RL     TKEdissipation
87        _RL     tempU, tempV, prTemp        _RL     tempU, tempV, prTemp
88          _RL     MaxLength
89        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     TKEPrandtlNumber (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
90        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL     GGL90mixingLength(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
91          _RL         rMixingLength(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)
# Line 103  C     dfx, dfy - diffusive flux across l Line 106  C     dfx, dfy - diffusive flux across l
106        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL     dfy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
108  #endif /* ALLOW_GGL90_HORIZDIFF */  #endif /* ALLOW_GGL90_HORIZDIFF */
109  CEOP  #ifdef ALLOW_GGL90_SMOOTH
110          _RL p4, p8, p16, tmpdiffKrS
111          p4=0.25 _d 0
112          p8=0.125 _d 0
113          p16=0.0625 _d 0
114    #endif
115        iMin = 2-OLx        iMin = 2-OLx
116        iMax = sNx+OLx-1        iMax = sNx+OLx-1
117        jMin = 2-OLy        jMin = 2-OLy
# Line 111  CEOP Line 119  CEOP
119    
120  C     set separate time step (should be deltaTtracer)  C     set separate time step (should be deltaTtracer)
121        deltaTggl90 = dTtracerLev(1)        deltaTggl90 = dTtracerLev(1)
122  C      
123        kSurf = 1        kSurf = 1
124  C     implicit timestepping weights for dissipation  C     implicit timestepping weights for dissipation
125        ab15 =  1.5 _d 0        ab15 =  1.5 _d 0
# Line 126  C     Initialize local fields Line 134  C     Initialize local fields
134           gTKE(I,J,K)              = 0. _d 0           gTKE(I,J,K)              = 0. _d 0
135           KappaE(I,J,K)            = 0. _d 0           KappaE(I,J,K)            = 0. _d 0
136           TKEPrandtlNumber(I,J,K)  = 0. _d 0           TKEPrandtlNumber(I,J,K)  = 0. _d 0
137           GGL90mixingLength(I,J,K) = 0. _d 0           GGL90mixingLength(I,J,K) = GGL90mixingLengthMin
138                 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)   = 0. _d 0          totalDepth(I,J) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
         IF ( recip_Rcol(I,J,bi,bj) .NE. 0. )  
      &       totalDepth(I,J) = 1./recip_Rcol(I,J,bi,bj)  
147         ENDDO         ENDDO
148        ENDDO        ENDDO
149    
150  C     start k-loop  C     start k-loop
151        DO K = 2, Nr        DO K = 2, Nr
152         Km1 = K-1         Km1 = K-1
153         Kp1 = MIN(Nr,K+1)  c      Kp1 = MIN(Nr,K+1)
154         CALL FIND_RHO(         CALL FIND_RHO_2D(
155       I      bi, bj, iMin, iMax, jMin, jMax, Km1, K,       I      iMin, iMax, jMin, jMax, K,
156       I      theta, salt,       I      theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),
157       O      rhoKm1,       O      rhoKm1,
158       I      myThid )       I      Km1, bi, bj, myThid )
159         CALL FIND_RHO(  
160       I      bi, bj, iMin, iMax, jMin, jMax, K, K,         CALL FIND_RHO_2D(
161       I      theta, salt,       I      iMin, iMax, jMin, jMax, K,
162         I      theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),
163       O      rhoK,       O      rhoK,
164       I      myThid )       I      K, bi, bj, myThid )
165         DO J=jMin,jMax         DO J=jMin,jMax
166          DO I=iMin,iMax          DO I=iMin,iMax
167           SQRTTKE=SQRT( GGL90TKE(I,J,K,bi,bj) )           SQRTTKE(i,j,k)=SQRT( GGL90TKE(I,J,K,bi,bj) )
168  C  C
169  C     buoyancy frequency  C     buoyancy frequency
170  C  C
171           Nsquare = - gravity*recip_rhoConst*recip_drC(K)           Nsquare(i,j,k) = - gravity*recip_rhoConst*recip_drC(K)
172       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)       &        * ( rhoKm1(I,J) - rhoK(I,J) )*maskC(I,J,K,bi,bj)
173    cC     vertical shear term (dU/dz)^2+(dV/dz)^2
174    c         tempU= .5 _d 0*( uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)
175    c     &                 -( uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )
176    c     &        *recip_drC(K)
177    c         tempV= .5 _d 0*( vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)
178    c     &                 -( vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )
179    c     &        *recip_drC(K)
180    c         verticalShear = tempU*tempU + tempV*tempV
181    c         RiNumber   = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
182    cC     compute Prandtl number (always greater than 0)
183    c         prTemp = 1. _d 0
184    c         IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
185    c         TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)
186    C     mixing length
187             GGL90mixingLength(I,J,K) = SQRTTWO *
188         &        SQRTTKE(i,j,k)/SQRT( MAX(Nsquare(i,j,k),GGL90eps) )
189            ENDDO
190           ENDDO
191          ENDDO
192    
193    C- Impose upper bound for mixing length (total depth)
194          IF ( mxlMaxFlag .EQ. 0 ) THEN
195           DO k=2,Nr
196            DO J=jMin,jMax
197             DO I=iMin,iMax
198              MaxLength=totalDepth(I,J)
199              GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
200         &        MaxLength)
201             ENDDO
202            ENDDO
203           ENDDO
204          ELSEIF ( mxlMaxFlag .EQ. 1 ) THEN
205           DO k=2,Nr
206            DO J=jMin,jMax
207             DO I=iMin,iMax
208              MaxLength=MIN(Ro_surf(I,J,bi,bj)-rF(k),rF(k)-R_low(I,J,bi,bj))
209    c         MaxLength=MAX(MaxLength,20. _d 0)
210              GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
211         &        MaxLength)
212             ENDDO
213            ENDDO
214           ENDDO
215          ELSEIF ( mxlMaxFlag .EQ. 2 ) THEN
216           DO k=2,Nr
217            DO J=jMin,jMax
218             DO I=iMin,iMax
219              GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
220         &        GGL90mixingLength(I,J,K-1)+drF(k-1))
221             ENDDO
222            ENDDO
223           ENDDO
224           DO J=jMin,jMax
225            DO I=iMin,iMax
226              GGL90mixingLength(I,J,Nr) = MIN(GGL90mixingLength(I,J,Nr),
227         &       GGL90mixingLengthMin+drF(Nr))
228            ENDDO
229           ENDDO
230           DO k=Nr-1,2,-1
231            DO J=jMin,jMax
232             DO I=iMin,iMax
233              GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),
234         &        GGL90mixingLength(I,J,K+1)+drF(k))
235             ENDDO
236            ENDDO
237           ENDDO
238          ELSE
239           STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing lenght limit)'
240          ENDIF
241    
242    C- Impose minimum mixing length (to avoid division by zero)
243          DO k=2,Nr
244           DO J=jMin,jMax
245            DO I=iMin,iMax
246             GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),
247         &        GGL90mixingLengthMin)
248             rMixingLength(I,J,K) = 1. _d 0 /GGL90mixingLength(I,J,K)
249            ENDDO
250           ENDDO
251          ENDDO
252    
253          DO k=2,Nr
254           Km1 = K-1
255           DO J=jMin,jMax
256            DO I=iMin,iMax
257  C     vertical shear term (dU/dz)^2+(dV/dz)^2  C     vertical shear term (dU/dz)^2+(dV/dz)^2
258           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)
259       &             - (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)) )
260       &        *recip_drC(K)       &        *recip_drC(K)
261           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)
262       &             - (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)) )
263       &        *recip_drC(K)       &        *recip_drC(K)
264           verticalShear = tempU*tempU + tempV*tempV           verticalShear = tempU*tempU + tempV*tempV
265           RiNumber   = MAX(Nsquare,0. _d 0)/(verticalShear+GGL90eps)           RiNumber = MAX(Nsquare(i,j,k),0. _d 0)/(verticalShear+GGL90eps)
266  C     compute Prandtl number (always greater than 0)  C     compute Prandtl number (always greater than 0)
267           prTemp = 1. _d 0           prTemp = 1. _d 0
268           IF ( RiNumber .GE. 0.2 ) prTemp = 5.0 * RiNumber           IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber
269           TKEPrandtlNumber(I,J,K) = MIN(10.,prTemp)           TKEPrandtlNumber(I,J,K) = MIN(10. _d 0,prTemp)
270  C     mixing length  
271           GGL90mixingLength(I,J,K) =  C     viscosity and diffusivity
272       &        SQRTTKE/SQRT( MAX(Nsquare,GGL90eps) )           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE(i,j,k)
273  C     impose upper bound for mixing length (total depth)           KappaH = KappaM/TKEPrandtlNumber(I,J,K)
274           GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K),  
275       &        totalDepth(I,J))  C     Set a minium (= background) and maximum value
276  C     impose minimum mixing length (to avoid division by zero)           KappaM = MAX(KappaM,viscArNr(k))
277           GGL90mixingLength(I,J,K) = MAX(GGL90mixingLength(I,J,K),           KappaH = MAX(KappaH,diffKrNrT(k))
278       &        GGL90mixingLengthMin)           KappaM = MIN(KappaM,GGL90viscMax)
279  C     viscosity of last timestep           KappaH = MIN(KappaH,GGL90diffMax)
280           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*SQRTTKE  
281           KappaE(I,J,K) = KappaM*GGL90alpha  C     Mask land points and storage
282             GGL90viscAr(I,J,K,bi,bj) = KappaM * maskC(I,J,K,bi,bj)
283             GGL90diffKr(I,J,K,bi,bj) = KappaH * maskC(I,J,K,bi,bj)
284             KappaE(I,J,K) = GGL90alpha * GGL90viscAr(I,J,K,bi,bj)
285    
286  C     dissipation term  C     dissipation term
287           TKEdissipation = ab05*GGL90ceps           TKEdissipation = ab05*GGL90ceps
288       &        *SQRTTKE/GGL90mixingLength(I,J,K)       &        *SQRTTKE(i,j,k)*rMixingLength(I,J,K)
289       &        *GGL90TKE(I,J,K,bi,bj)             &        *GGL90TKE(I,J,K,bi,bj)
290  C     sum up contributions to form the right hand side  C     sum up contributions to form the right hand side
291           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)           gTKE(I,J,K) = GGL90TKE(I,J,K,bi,bj)
292       &        + deltaTggl90*(       &        + deltaTggl90*(
293       &        + KappaM*verticalShear       &        + KappaM*verticalShear
294       &        - KappaM*Nsquare/TKEPrandtlNumber(I,J,K)       &        - KappaH*Nsquare(i,j,k)
295       &        - TKEdissipation       &        - TKEdissipation
296       &        )       &        )
297          ENDDO            ENDDO
298         ENDDO         ENDDO
299        ENDDO        ENDDO
300  C     horizontal diffusion of TKE (requires an exchange in  
301  C       do_fields_blocking_exchanges)  C     horizontal diffusion of TKE (requires an exchange in
302    C      do_fields_blocking_exchanges)
303  #ifdef ALLOW_GGL90_HORIZDIFF  #ifdef ALLOW_GGL90_HORIZDIFF
304        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN        IF ( GGL90diffTKEh .GT. 0. _d 0 ) THEN
305         DO K=2,Nr         DO K=2,Nr
# Line 214  C     common factors Line 311  C     common factors
311            yA(i,j) = _dxG(i,j,bi,bj)            yA(i,j) = _dxG(i,j,bi,bj)
312       &         *drF(k)*_hFacS(i,j,k,bi,bj)       &         *drF(k)*_hFacS(i,j,k,bi,bj)
313           ENDDO           ENDDO
314          ENDDO            ENDDO
315  C     Compute diffusive fluxes  C     Compute diffusive fluxes
316  C     ... across x-faces  C     ... across x-faces
317          DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
318           dfx(1-Olx,j)=0.           dfx(1-Olx,j)=0. _d 0
319           DO i=1-Olx+1,sNx+Olx           DO i=1-Olx+1,sNx+Olx
320            dfx(i,j) = -GGL90diffTKEh*xA(i,j)            dfx(i,j) = -GGL90diffTKEh*xA(i,j)
321       &      *_recip_dxC(i,j,bi,bj)       &      *_recip_dxC(i,j,bi,bj)
# Line 228  C     ... across x-faces Line 325  C     ... across x-faces
325          ENDDO          ENDDO
326  C     ... across y-faces  C     ... across y-faces
327          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
328           dfy(i,1-Oly)=0.           dfy(i,1-Oly)=0. _d 0
329          ENDDO          ENDDO
330          DO j=1-Oly+1,sNy+Oly          DO j=1-Oly+1,sNy+Oly
331           DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
# Line 239  C     ... across y-faces Line 336  C     ... across y-faces
336       &      *CosFacV(j,bi,bj)       &      *CosFacV(j,bi,bj)
337  #endif /* ISOTROPIC_COS_SCALING */  #endif /* ISOTROPIC_COS_SCALING */
338           ENDDO           ENDDO
339          ENDDO            ENDDO
340  C     Compute divergence of fluxes  C     Compute divergence of fluxes
341          DO j=1-Oly,sNy+Oly-1          DO j=1-Oly,sNy+Oly-1
342           DO i=1-Olx,sNx+Olx-1           DO i=1-Olx,sNx+Olx-1
343            gTKE(i,j,k)=gTKE(i,j,k)            gTKE(i,j,k)=gTKE(i,j,k)
344       &   -_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)
345       &         *( (dfx(i+1,j)-dfx(i,j))       &         *( (dfx(i+1,j)-dfx(i,j))
346       &           +(dfy(i,j+1)-dfy(i,j))       &           +(dfy(i,j+1)-dfy(i,j))
347       &           )       &           )
348           ENDDO             ENDDO
349          ENDDO          ENDDO
350  C       end of k-loop  C       end of k-loop
351         ENDDO         ENDDO
352  C     end if GGL90diffTKEh .eq. 0.  C     end if GGL90diffTKEh .eq. 0.
353        ENDIF        ENDIF
# Line 264  C     set up matrix Line 361  C     set up matrix
361  C--   Lower diagonal  C--   Lower diagonal
362        DO j=jMin,jMax        DO j=jMin,jMax
363         DO i=iMin,iMax         DO i=iMin,iMax
364           a(i,j,1) = 0. _d 0           a(i,j,1) = 0. _d 0
365         ENDDO         ENDDO
366        ENDDO        ENDDO
367        DO k=2,Nr        DO k=2,Nr
368         km1=MAX(1,k-1)         km1=max(2,k-1)
369         DO j=jMin,jMax         DO j=jMin,jMax
370          DO i=iMin,iMax          DO i=iMin,iMax
371           a(i,j,k) = -deltaTggl90           a(i,j,k) = -deltaTggl90
372       &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)  c     &        *recip_drF(km1)*recip_hFacI(i,j,k,bi,bj)
373       &        *.5*(KappaE(i,j, k )+KappaE(i,j,km1))       &        *recip_drF(k-1)*recip_hFacC(i,j,k-1,bi,bj)
374       &        *recip_drC(k)       &        *.5 _d 0*(KappaE(i,j, k )+KappaE(i,j,km1))
375            IF (recip_hFacI(i,j,km1,bi,bj).EQ.0.) a(i,j,k)=0.       &        *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
376          ENDDO          ENDDO
377         ENDDO         ENDDO
378        ENDDO        ENDDO
# Line 283  C--   Upper diagonal Line 380  C--   Upper diagonal
380        DO j=jMin,jMax        DO j=jMin,jMax
381         DO i=iMin,iMax         DO i=iMin,iMax
382           c(i,j,1)  = 0. _d 0           c(i,j,1)  = 0. _d 0
          c(i,j,Nr) = 0. _d 0  
383         ENDDO         ENDDO
384        ENDDO        ENDDO
385  CML      DO k=1,Nr-1        DO k=2,Nr
       DO k=2,Nr-1  
        kp1=min(Nr,k+1)  
386         DO j=jMin,jMax         DO j=jMin,jMax
387          DO i=iMin,iMax          DO i=iMin,iMax
388              kp1=min(klowC(i,j,bi,bj),k+1)
389            c(i,j,k) = -deltaTggl90            c(i,j,k) = -deltaTggl90
390       &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)  c     &        *recip_drF( k )*recip_hFacI(i,j,k,bi,bj)
391       &               *.5*(KappaE(i,j,k)+KappaE(i,j,kp1))       &        *recip_drF( k ) * recip_hFacC(i,j,k,bi,bj)
392       &        *recip_drC(k)       &        *.5 _d 0*(KappaE(i,j,k)+KappaE(i,j,kp1))
393            IF (recip_hFacI(i,j,kp1,bi,bj).EQ.0.) c(i,j,k)=0.       &        *recip_drC(k)*maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
394          ENDDO          ENDDO
395         ENDDO         ENDDO
396        ENDDO        ENDDO
# Line 305  C--   Center diagonal Line 400  C--   Center diagonal
400          DO i=iMin,iMax          DO i=iMin,iMax
401            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)
402       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))       &        + ab15*deltaTggl90*GGL90ceps*SQRT(GGL90TKE(I,J,K,bi,bj))
403       &        /GGL90mixingLength(I,J,K)       &        *rMixingLength(I,J,K)*maskC(i,j,k,bi,bj)
404           ENDDO           ENDDO
405         ENDDO         ENDDO
406        ENDDO        ENDDO
# Line 313  C     end set up matrix Line 408  C     end set up matrix
408    
409  C  C
410  C     Apply boundary condition  C     Apply boundary condition
411  C      C
412        DO J=jMin,jMax        DO J=jMin,jMax
413         DO I=iMin,iMax         DO I=iMin,iMax
414  C     estimate friction velocity uStar from surface forcing  C     estimate friction velocity uStar from surface forcing
415          uStarSquare = SQRT(          uStarSquare = SQRT(
416       &         ( .5*( surfaceForcingU(I,  J,  bi,bj)       &    ( .5 _d 0*( surfaceForcingU(I,  J,  bi,bj)
417       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2       &              + surfaceForcingU(I+1,J,  bi,bj) ) )**2
418       &       + ( .5*( surfaceForcingV(I,  J,  bi,bj)       &  + ( .5 _d 0*( surfaceForcingV(I,  J,  bi,bj)
419       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2       &              + surfaceForcingV(I,  J+1,bi,bj) ) )**2
420       &                     )       &                     )
421  C     Dirichlet surface boundary condition for TKE  C     Dirichlet surface boundary condition for TKE
422          gTKE(I,J,kSurf) = MAX(GGL90TKEmin,GGL90m2*uStarSquare)          gTKE(I,J,kSurf) = MAX(GGL90TKEsurfMin,GGL90m2*uStarSquare)
423       &                     *maskC(I,J,kSurf,bi,bj)       &                     *maskC(I,J,kSurf,bi,bj)
424  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom  C     Dirichlet bottom boundary condition for TKE = GGL90TKEbottom
425          kBottom   = MIN(MAX(kLowC(I,J,bi,bj),1),Nr)          kBottom   = MAX(kLowC(I,J,bi,bj),1)
426          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)          gTKE(I,J,kBottom) = gTKE(I,J,kBottom)
427       &       - GGL90TKEbottom*c(I,J,kBottom)       &                    - GGL90TKEbottom*c(I,J,kBottom)
428            c(I,J,kBottom) = 0. _d 0
429         ENDDO         ENDDO
430        ENDDO            ENDDO
431  C  C
432  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)  C     solve tri-diagonal system, and store solution on gTKE (previously rhs)
433  C  C
# Line 341  C Line 437  C
437       I     myThid )       I     myThid )
438  C  C
439  C     now update TKE  C     now update TKE
440  C      C
441        DO K=1,Nr        DO K=1,Nr
442         DO J=jMin,jMax         DO J=jMin,jMax
443          DO I=iMin,iMax          DO I=iMin,iMax
444  C     impose minimum TKE to avoid numerical undershoots below zero  C     impose minimum TKE to avoid numerical undershoots below zero
445           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )           GGL90TKE(I,J,K,bi,bj) = MAX( gTKE(I,J,K), GGL90TKEmin )
446       &        * maskC(I,J,K,bi,bj)       &        * maskC(I,J,K,bi,bj)
447          ENDDO          ENDDO
448         ENDDO         ENDDO
449        ENDDO            ENDDO
450  C  
451  C     end of time step  C     end of time step
452  C     ===============================  C     ===============================
453  C     compute viscosity coefficients  
454  C      #ifdef ALLOW_GGL90_SMOOTH
455        DO K=2,Nr        DO K=1,Nr
456         DO J=jMin,jMax         DO J=jMin,jMax
457          DO I=iMin,iMax          DO I=iMin,iMax
458  C     Eq. (11), (18) and (21)           tmpdiffKrS=
459           KappaM = GGL90ck*GGL90mixingLength(I,J,K)*       &  (
460       &                  SQRT( GGL90TKE(I,J,K,bi,bj) )       &   p4 *  GGL90viscAr(i  ,j  ,k,bi,bj) * mskCor(i  ,j  ,bi,bj)
461           KappaH = KappaM/TKEPrandtlNumber(I,J,K)       &  +p8 *( GGL90viscAr(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
462  C     Set a minium (= background) value       &       + GGL90viscAr(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
463           KappaM = MAX(KappaM,viscAr)       &       + GGL90viscAr(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
464           KappaH = MAX(KappaH,diffKrNrT(k))       &       + GGL90viscAr(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
465  C     Set a maximum and mask land point       &  +p16*( GGL90viscAr(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)
466           GGL90viscAr(I,J,K,bi,bj) = MIN(KappaM,GGL90viscMax)       &       + GGL90viscAr(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)
467       &        * maskC(I,J,K,bi,bj)       &       + GGL90viscAr(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
468           GGL90diffKr(I,J,K,bi,bj) = MIN(KappaH,GGL90diffMax)       &       + GGL90viscAr(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))
469       &        * maskC(I,J,K,bi,bj)       &  )
470         & /(p4
471         &  +p8 *(       maskC(i-1,j  ,k,bi,bj) * mskCor(i-1,j  ,bi,bj)
472         &       +       maskC(i  ,j-1,k,bi,bj) * mskCor(i  ,j-1,bi,bj)
473         &       +       maskC(i+1,j  ,k,bi,bj) * mskCor(i+1,j  ,bi,bj)
474         &       +       maskC(i  ,j+1,k,bi,bj) * mskCor(i  ,j+1,bi,bj))
475         &  +p16*(       maskC(i+1,j+1,k,bi,bj) * mskCor(i+1,j+1,bi,bj)
476         &       +       maskC(i+1,j-1,k,bi,bj) * mskCor(i+1,j-1,bi,bj)
477         &       +       maskC(i-1,j+1,k,bi,bj) * mskCor(i-1,j+1,bi,bj)
478         &       +       maskC(i-1,j-1,k,bi,bj) * mskCor(i-1,j-1,bi,bj))
479         &  )*maskC(i,j,k,bi,bj)*mskCor(i,j,bi,bj)
480         &   /TKEPrandtlNumber(i,j,k)
481             GGL90diffKrS(I,J,K,bi,bj)= MAX( tmpdiffKrS , diffKrNrT(k) )
482          ENDDO          ENDDO
483         ENDDO         ENDDO
484  C     end third k-loop        ENDDO
485        ENDDO      #endif
486    
487    #ifdef ALLOW_DIAGNOSTICS
488          IF ( useDiagnostics ) THEN
489             CALL DIAGNOSTICS_FILL( GGL90TKE   ,'GGL90TKE',
490         &                          0,Nr, 1, bi, bj, myThid )
491             CALL DIAGNOSTICS_FILL( GGL90viscAr,'GGL90Ar ',
492         &                          0,Nr, 1, bi, bj, myThid )
493             CALL DIAGNOSTICS_FILL( GGL90diffKr,'GGL90Kr ',
494         &                          0,Nr, 1, bi, bj, myThid )
495             CALL DIAGNOSTICS_FILL( TKEPrandtlNumber ,'GGL90Prl',
496         &                          0,Nr, 2, bi, bj, myThid )
497             CALL DIAGNOSTICS_FILL( GGL90mixingLength,'GGL90Lmx',
498         &                          0,Nr, 2, bi, bj, myThid )
499          ENDIF
500    #endif
501    
502  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
503    
504        RETURN        RETURN
505        END        END
   
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22