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" |
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 |
112 |
p8=0.125 _d 0 |
p8=0.125 _d 0 |
113 |
p16=0.0625 _d 0 |
p16=0.0625 _d 0 |
114 |
#endif |
#endif |
|
CEOP |
|
115 |
iMin = 2-OLx |
iMin = 2-OLx |
116 |
iMax = sNx+OLx-1 |
iMax = sNx+OLx-1 |
117 |
jMin = 2-OLy |
jMin = 2-OLy |
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 |
232 |
DO I=iMin,iMax |
DO I=iMin,iMax |
233 |
GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K), |
GGL90mixingLength(I,J,K) = MIN(GGL90mixingLength(I,J,K), |
234 |
& GGL90mixingLength(I,J,K+1)+drF(k)) |
& GGL90mixingLength(I,J,K+1)+drF(k)) |
235 |
ENDDO |
ENDDO |
236 |
ENDDO |
ENDDO |
237 |
ENDDO |
ENDDO |
238 |
ELSE |
ELSE |
239 |
STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing lenght limit)' |
STOP 'GGL90_CALC: Wrong mxlMaxFlag (mixing lenght limit)' |
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(i,j,k),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 _d 0 ) prTemp = 5. _d 0 * RiNumber |
IF ( RiNumber .GE. 0.2 _d 0 ) prTemp = 5. _d 0 * RiNumber |
273 |
KappaH = KappaM/TKEPrandtlNumber(I,J,K) |
KappaH = KappaM/TKEPrandtlNumber(I,J,K) |
274 |
|
|
275 |
C Set a minium (= background) and maximum value |
C Set a minium (= background) and maximum value |
276 |
KappaM = MAX(KappaM,viscAr) |
KappaM = MAX(KappaM,viscArNr(k)) |
277 |
KappaH = MAX(KappaH,diffKrNrT(k)) |
KappaH = MAX(KappaH,diffKrNrT(k)) |
278 |
KappaM = MIN(KappaM,GGL90viscMax) |
KappaM = MIN(KappaM,GGL90viscMax) |
279 |
KappaH = MIN(KappaH,GGL90diffMax) |
KappaH = MIN(KappaH,GGL90diffMax) |
503 |
|
|
504 |
RETURN |
RETURN |
505 |
END |
END |
|
|
|