Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

PHY_genTKE_RUN.f90 @ 2100

Last change on this file since 2100 was 2089, checked in by Laurent Fairhead, 10 years ago

Inclusion de la physique de MAR

Integration of MAR physics

File size: 30.6 KB

Line
1	subroutine PHY_genTKE_RUN
2
3	!------------------------------------------------------------------------------+
4	! Sat 22-Jun-2013 MAR \|
5	! MAR PHY_genTKE_RUN \|
6	! subroutine PHY_genTKE_RUN computes Turbulent Kinetic Energy \|
7	! \|
8	! version 3.p.4.1 created by H. Gallee, Fri 15-Mar-2013 \|
9	! Last Modification by H. Gallee, Sat 22-Jun-2013 \|
10	! \|
11	!------------------------------------------------------------------------------+
12	! \|
13	! METHOD: 1. `Standard' Closure of Turbulence: \|
14	! ^^^^^^^ E - epsilon , Duynkerke, JAS 45, 865--880, 1988 \|
15	! 2. E - epsilon , Huang and Raman, BLM 55, 381--407, 1991 \|
16	! 3. K - l , Therry et Lacarrere, BLM 25, 63-- 88, 1983 \|
17	! \|
18	! INPUT : itexpe: Nb of iterations \|
19	! ^^^^^^^^ dt__AT: Local Time Step (s) \|
20	! explIO: Experiment Label (s) \|
21	! \|
22	! INPUT / OUTPUT: The Vertical Turbulent Fluxes are included for: \|
23	! ^^^^^^^^^^^^^^ \|
24	! a) Turbulent Kinetic Energy TKE_AT(kcolq,mzp) (m2/s2) \|
25	! b) Turbulent Kinetic Energy Dissipation eps_AT(kcolq,mzp) (m2/s3) \|
26	! \|
27	! OUTPUT : Kzm_AT(kcolq,mzp): vertic. diffusion coeffic. (momentum) (m2/s) \|
28	! ^^^^^^^^ Kzh_AT(kcolq,mzp): vertic. diffusion coeffic. (scalars ) (m2/s) \|
29	! zi__AT(kcolq) : inversion height (m) \|
30	! \|
31	! OPTIONS: #HY: Latent Heat Exchanges associated with Cloud Microphysics \|
32	! ^^^^^^^^ contribute to Buoyancy \|
33	! #KA: replaces TKE & e below a prescribed height above the surface \|
34	! by their box weighted vertical moving averages \|
35	! #KC: Modification of TKE near the Lower Boundary \|
36	! #LD: Buoyancy includes Loading by the Hydrometeors \|
37	! #RI: Correction of the Prandtl Nb (Kzm/Kzh) \|
38	! by the The Richardson Nb \|
39	! \|
40	!------------------------------------------------------------------------------+
41
42	use Mod_Real
43	use Mod_Phy____dat
44	use Mod_Phy____grd
45	use Mod_PHY_AT_ctr
46	use Mod_PHY_AT_grd
47	use Mod_PHY____kkl
48	use Mod_PHY_AT_kkl
49	use Mod_PHY_DY_kkl
50	use Mod_PHY_CM_kkl
51	use Mod_SISVAT_gpt
52
53
54
55	! Local Variables
56	! ================
57
58	use Mod_genTKE_RUN
59
60
61	IMPLICIT NONE
62
63
64	real(kind=real8) :: z__SBL ! SBL Thickness assumed to be the lowest Model Layer Height [m]
65	real(kind=real8) :: TKE_zi ! TKE at the inversion height [m2/s2]
66	integer :: kTKEzi ! level correponding to TKE_zi
67	real(kind=real8) :: dTKEdk ! TKE difference across model Layers [m2/s2]
68	real(kind=real8) :: kz_inv ! 1 / (k z)
69	real(kind=real8) :: Buoy_F ! Contribution of Buoyancy Force
70	real(kind=real8) :: dTKE_p ! Positive Contribution to TKE
71	real(kind=real8) :: r_eTKE ! e / TKE Ratio (used in the Product./Destruct. Scheme of TKE)
72	real(kind=real8) :: TKE_ds ! Verticaly Integrated TKE
73	real(kind=real8) :: eps_ds ! Verticaly Integrated TKE Dissipation
74	real(kind=real8) :: edt_HR ! Minimum Energy Dissipation Time
75	real(kind=real8) :: TKEnew ! New Value of TKE
76	real(kind=real8) :: TKEsbc ! SBC: TKE [m2/s2]
77	real(kind=real8) :: epssbc ! SBC: TKE Dissipation [m2/s3]
78	real(kind=real8) :: zeta ! SBL: z / LMO [-]
79	real(kind=real8) :: u_star ! SBL: u* (Friction Velocity) [m/s]
80	real(kind=real8) :: kz_max,kz_mix,kz2mix ! Kzh Auxiliary Variables
81	real(kind=real8) :: KzhMAX ! max.vert. Turbulent Diffusion Coefficient (Scalars) [m2/s]
82
83	! #KC real(kind=real8) :: se ! 1 if TKE(mzp-2) > TKE(mzp-1) and 0 otherwise
84	! #KC integer :: ke ! index of the level where TKE is max
85
86	real(kind=real8) :: relHum ! Relative Humidity [-]
87	real(kind=real8) :: QS_mid ! Saturation Sp. Humidity % Liquid Water [kg/kg]
88	real(kind=real8) :: qwsLRT ! Qws * L / (Rd * T)
89	real(kind=real8) :: surSat ! Normalized sur-Saturation (> 0 if relHum > 1)
90	real(kind=real8) :: C_thq ! (Duynkerke & Driedonks 1987 JAS 44(1), Table 1 p.47)
91	real(kind=real8) :: C_q_w ! (Duynkerke & Driedonks 1987)
92	real(kind=real8) :: CX__hi !
93	real(kind=real8) :: Ce1_hi ! ... Ce1 / hi
94	real(kind=real8) :: Ck1_hi ! ... Ck1 / hi
95	real(kind=real8) :: Ck2_hi ! ... Ck2 / hi
96	real(kind=real8) :: Th_Lac ! Therry & Lacarrere (1983) Parameter
97
98	real(kind=real8) :: sgnLMO ! sign(LMO)
99	real(kind=real8) :: absLMO ! \|LMO\|
100	! #RI real(kind=real8) :: fac_Ri,vuzvun,Kz_vun ! Sukorianski Parameterization Parameters
101
102	integer :: ikl
103	integer :: i
104	integer :: j
105	integer :: k
106
107
108
109
110	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
111	! !
112	! ALLOCATION
113	! ==========
114
115	IF (it_RUN.EQ.1 .OR. FlagDALLOC) THEN !
116
117	allocate ( dukkp1(mzpp) ) ! Difference (u(k) - u(k+1)) [m/s]
118	allocate ( dvkkp1(mzpp) ) ! Difference (v(k) - v(k+1)) [m/s]
119	allocate ( kkp1dz(mzpp) ) ! 1 / Difference (Z(k) - Z(k+1)) [1/m]
120	allocate ( zShear(mzp) ) ! Wind Shear Contribution to TKE [m2/s3]
121	allocate ( REq_PT(mzpp) ) ! Reduced (Equivalent) Potential Temperature [K]
122	allocate ( c_Buoy(mzpp) ) ! Buoyancy Coefficient (g/theta) X (dtheta/dz) [1/s2]
123	allocate ( Ri__Nb(mzp) ) ! Richardson Number [-]
124	allocate ( Prandtl(mzp) ) ! Prandtl Number (Kzm/Kzh) [-]
125	allocate ( Ls_inv(mzp) ) ! 1 / Ls (Therry & Lacarrere, 1983) [1/m]
126	allocate ( ML_inv(mzp) ) ! 1 / ML (Mixing Length, Therry & Lacarr, 1983) [1/m]
127	allocate ( DL_inv(mzp) ) ! 1 / DL (Dissipation Length, Therry & Lacarr, 1983) [1/m]
128	allocate ( Dissip(mzp) ) ! Dissipation [m2/s3]
129	allocate ( TKEvav(mzp) ) ! TKE Vertical moving Average [m2/s2]
130	allocate ( epsvav(mzp) ) ! Dissipation Vertical moving Average [m2/s3]
131	allocate ( pkt (mzpp) ) ! Reduced Potential Temperature [X]
132
133	END IF
134	! !
135	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
136
137
138
139
140	! ++++++++++++++
141	DO ikl=1,kcolp
142	! ++++++++++++++
143
144	! i = ii__AP(ikl)
145	! j = jj__AP(ikl)
146
147
148
149
150	! Friction Velocity
151	! =================
152
153	u_star = us__SV_gpt(ikl)
154
155
156
157
158	! Verification: TKE must be Positive Definite
159	! ===========================================
160
161	DO k=1,mzp
162	TKE_AT(ikl,k)=max(eps6,TKE_AT(ikl,k))
163	eps_AT(ikl,k)=max(epsn,eps_AT(ikl,k))
164	END DO
165
166
167
168
169	! Reduced Potential Temperature
170	! =============================
171
172	DO k=1,mzpp
173	pkt (k) = pkt_DY(ikl,k)
174	END DO
175
176
177
178
179	! Inversion Height
180	! ================
181
182	TKE_zi = 0.05*max(max(TKE_AT(ikl,mzp-1), &
183	& TKE_AT(ikl,mzp )),TKEmin)
184	kTKEzi = 1
185
186
187	DO k=1,mzp
188	IF (TKE_AT(ikl,k) .lt. TKE_zi ) THEN
189	kTKEzi = min(mzp-1,k)
190	END IF
191	ENDDO
192
193	k = kTKEzi
194	IF (TKE_AT(ikl,k+1).lt.TKEmin) THEN
195	TKE_zi = ZmidDY(ikl,mzp)-sh__AP(ikl)
196	ELSE
197	dTKEdk = TKE_AT(ikl,k) -TKE_AT(ikl,k+1)
198	TKE_zi = ZmidDY(ikl,k+2) &
199	& +(ZmidDY(ikl,k+1)-ZmidDY(ikl,k+2)) &
200	& *(TKE_zi -TKE_AT(ikl,k+1)) &
201	& /(sign(un_1 , dTKEdk) &
202	& *max(epsn ,abs(dTKEdk))) -sh__AP(ikl)
203	END IF
204
205	zi__AT(ikl) = min(TKE_zi, ZmidDY(ikl, 1)-sh__AP(ikl))
206
207	zi__AT(ikl) = max( ZmidDY(ikl,mzp)-sh__AP(ikl) &
208	& ,zi__AT(ikl))
209	TKE_zi = 0.
210	kTKEzi = 0
211
212
213
214
215	! TKE Production/Destruction by the Vertical Wind Shear
216	! =====================================================
217
218	DO k=1,mzp-1
219	dukkp1(k) = ua__DY(ikl,k) - ua__DY(ikl,k+1)
220	dvkkp1(k) = va__DY(ikl,k) - va__DY(ikl,k+1)
221	END DO
222	k= mzp
223	dukkp1(k) = ua__DY(ikl,k)
224	dvkkp1(k) = va__DY(ikl,k)
225
226	DO k=1,mzp
227	kkp1dz(k) = Z___DY(ikl,k) - Z___DY(ikl,k+1) ! dz(k+1/2)
228	END DO
229
230	DO k=1,mzp-1
231	zShear(k) = &
232	& Kzm_AT(ikl,k)(dukkp1(k) dukkp1(k) + dvkkp1(k) *dvkkp1(k)) &
233	& /(kkp1dz(k) *kkp1dz(k))
234	END DO
235	zShear(mzp) = 0.0
236
237
238
239
240	! Buoyancy
241	! ========
242
243	! Reduced (Equivalent) Potential Temperature
244	! ------------------------------------------
245
246	! Control Martin
247	!PRINT*,'CpdAir=',CpdAir
248	!PRINT*,'minval(Ta__DY(ikl,:))=',minval(Ta__DY(ikl,:))
249	! Control Martin
250
251	DO k=1,mzpp
252	REq_PT(k) = pkt ( k) &
253	! #LD& * (1.0 + 0.715 * ld_H2O(ikl,k) ) &
254	& * exp(LhvH2O * qv__DY(ikl,k) &
255	& / (CpdAir * Ta__DY(ikl,k)) ) &
256	& + 0.0
257	END DO
258
259
260
261	! Buoyancy Coefficients
262	! ---------------------
263
264	DO k=1,mzp
265
266	relHum = 0.50 *(qv__DY(ikl,k ) /qvswCM(ikl,k ) &
267	& +qv__DY(ikl,k+1) /qvswCM(ikl,k+1))
268	QS_mid = 0.50 *(qvswCM(ikl,k ) +qvswCM(ikl,k+1))
269
270	surSat = max(zer0,sign(un_1,relHum +epsp-un_1))
271	qwsLRT = LhvH2OQS_mid /(R_DAir TmidDY(ikl,k))
272
273	! Vectorization of the unsaturated (H<1) and saturated cases (H=1.)
274	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
275	C_thq =1.-surSat & ! H<1.
276	& +surSat*(1.+qwsLRT) & ! H=1.
277	& /(1.+qwsLRT.605LhvH2O & !
278	& /(CpdAir*TmidDY(ikl,k))) !
279	! ... C_thq (Duynkerke & Driedonks 1987 JAS 44(1), Table 1 p.47) !
280
281	C_q_w=(1.-surSat) *(LhvH2O & !
282	& /(CpdAir*TmidDY(ikl,k)) & ! H<1.
283	& - 0.605 ) & !
284	& +surSat ! H=1.
285	! ... C_q_w (Duynkerke and Driedonks 1987) !
286	! with (1-Ra/Rv)/(Ra/Rv) = 0.605 [Ra=287.J/kg/K; !
287	! Rv=461.J/kg/K] !
288
289	! Unsaturated Case
290	! ~~~~~~~~~~~~~~~~
291	! IF(relHum.lt.1.0) THEN !
292	! C_thq = 1.0
293	! C_q_w = LhvH2O/(CpdAir*TmidDY(ikl,k)) & !
294	! & - 0.605 !
295
296	! Saturated Case
297	! ~~~~~~~~~~~~~~~~
298	! ELSE !
299	! qwsLRT= QS_mid LhvH2O/(RDryAiTmidDY(ikl,k)) !
300	! C_thq = (1.0+qwsLRT) & !
301	! & /(1.0+qwsLRT0.605LhvH2O/(CpdAir*TmidDY(ikl,k))) !
302	! C_q_w = 1.0
303	! END IF
304
305
306	! Buoyancy
307	! --------
308
309	IF(k.EQ.mzp) C_q_w = 0.0
310
311	c_Buoy(k) = Grav_F &
312	& kkp1dz(k) ( (REq_PT(k)-REq_PT(k+1)) &
313	& *2./(REq_PT(k)+REq_PT(k+1)) &
314	& * C_thq &
315	& - C_q_w *(qv__DY(ikl,k)-qv__DY(ikl,k+1) &
316	& +qw__CM(ikl,k)-qw__CM(ikl,k+1) &
317	& +qr__CM(ikl,k)-qr__CM(ikl,k+1) &
318	& +qi__CM(ikl,k)-qi__CM(ikl,k+1) &
319	& +qs__CM(ikl,k)-qs__CM(ikl,k+1)) &
320	& )
321	! ... (g/theta) X(dtheta/dz) :
322	! Buoyancy Parameter beta X grad.vert.temp.pot. en k+1/2
323
324	END DO
325
326
327	! Dissipation & Length Scales Parameters (Therry and Lacarrere 1983 Model)
328	! ======================================
329
330	IF (Kl_TherryLac) THEN
331	CX__hi = 1.0 / zi__AT(ikl)
332	Ce1_hi = 15.0 * CX__hi ! ... Ce1/hi
333	Ck1_hi = 5.0 * CX__hi ! ... Ck1/hi
334	Ck2_hi = 11.0 * CX__hi ! ... Ck2/hi
335
336	sgnLMO = sign(un_1,LMO_AT(ikl))
337	absLMO = abs( LMO_AT(ikl))
338	LMO_AT(ikl) = sgnLMO * max(absLMO,epsp)
339	Th_Lac = -min(0.00,sgnLMO) &
340	& /(1.0 -min(0.00,LMO_AT(ikl)) / zi__AT(ikl))
341
342	! Replacement of:
343	! IF (LMO_AT(ikl).lt.0.0) THEN
344	! LMO_AT(ikl) = min(LMO_AT(ikl),-epsp)
345	! Th_Lac = 1.0/(1.d0-LMO_AT(ikl)/zi__AT(ikl))
346	! ELSE
347	! LMO_AT(ikl) = max(LMO_AT(ikl), epsp)
348	! Th_Lac = 0.0
349	! END IF
350	! ... m2
351
352
353	DO k=1,mzp
354	Ls_inv(k) = sqrt( max(zer0,c_Buoy(k)) /TKE_AT(ikl,k) ) ! ... 1/ls
355	END DO
356
357
358	! Dissipation Length
359	! ------------------
360
361	DO k=1,mzp
362	kz_inv=vK_inv/(ZmidDY(ikl,k+1)-sh__AP(ikl)) ! ... 1/kz(i,j,k+1/2)
363	!
364	DL_inv(k)=kz_inv +Ce1_hi &! ... DL_inv=1/Dissipation Length
365	& -(kz_inv+Ck1_hi)Th_Lac/(1.+5.0e-3zi__AT(ikl)*kz_inv) &! (Therry and Lacarrere, 1983 BLM 25 p.75)
366	& +1.5*Ls_inv(k)
367
368
369	! Mixing Length
370	! -------------
371
372	ML_inv(k)=kz_inv +Ce1_hi &! ... ML_inv=1/Mixing Length
373	& -(kz_inv+Ck2_hi)Th_Lac/(1.+2.5e-3zi__AT(ikl)*kz_inv) &! (Therry and Lacarrere, 1983 BLM 25 p.78)
374	& +3.0*Ls_inv(k)
375
376	Dissip(k) = 0.125DL_inv(k)sqrt(TKE_AT(ikl,k))*TKE_AT(ikl,k)
377	ENDDO
378
379
380
381
382	! Dissipation (E-epsilon Models)
383	! ===========
384
385	ELSE
386
387	DO k=1,mzp
388	Dissip(k) = eps_AT(ikl,k)
389	ENDDO
390
391	END IF
392
393
394
395
396	! Contribution of Vertical Wind Shear + Buoyancy + Dissipation to TKE
397	! ===================================================================
398
399	DO k=1,mzp
400	Buoy_F=-Kzh_AT(ikl,k) * c_Buoy(k)
401
402	TKEnew= TKE_AT(ikl,k) * &
403	& (TKE_AT(ikl,k)+dt__AT*(zShear(k) +max(zer0,Buoy_F))) &
404	& /(TKE_AT(ikl,k)+dt__AT*( -min(zer0,Buoy_F) &
405	& + Dissip(k) ))
406	! ... Numerical Scheme : cfr. E. Deleersnijder, 1992 (thesis) pp.59-61
407
408
409
410
411	! Contribution of Vertical Wind Shear + Buoyancy to epsilon
412	! =========================================================
413
414	IF (Ee_Duynkerke) THEN
415	dTKE_p = zShear(k) +max(zer0,Buoy_F)+max(TrT_AT(ikl,k),zer0)
416	ELSE
417	dTKE_p = zShear(k) +max(zer0,Buoy_F)
418	! ... based on standard values of Kitada, 1987, BLM 41, p.220
419	END IF
420
421	! #BH dTKE_p = zShear(k) +max(zer0,Buoy_F) * 1.80 &
422	! #BH& -min(Buoy_F,zer0) * 1.15
423	! ... based on values of Betts et Haroutunian, 1983
424	! can be used by replacing strings `c #KI' (except the previous one)
425	! and `c #BH' by blanks
426	! (cfr. Kitada, 1987, BLM 41, p.220):
427	! buoyancy > 0 (unstability) => (1-ce3) X buoyancy = 1.8 X buoyancy
428	! buoyancy < 0 ( stability) => (1-ce3) X buoyancy =-1.15 X buoyancy
429
430	r_eTKE = eps_AT(ikl,k) /TKE_AT(ikl,k)
431	eps_AT(ikl,k) = &
432	& eps_AT(ikl,k) &
433	& (eps_AT(ikl,k) +dt__AT r_eTKE c1ep dTKE_p) &
434	& /(eps_AT(ikl,k) +dt__AT r_eTKE c2ep *eps_AT(ikl,k))
435	! ... Numerical Scheme : cfr. E. Deleersnijder, 1992 (thesis) pp.59-61
436
437	IF (Kl_TherryLac) &
438	& eps_AT(ikl,k)=Dissip(k)
439
440
441
442
443	! New TKE Value
444	! =============
445
446	TKE_AT(ikl,k)=TKEnew
447
448	END DO
449
450
451
452	! Lower Boundary Conditions
453	! =========================
454
455	TKEsbc = u_star * u_star ! -> TKE SBC
456	z__SBL = Z___DY(ikl,mzp) - sh__AP(ikl) ! z_SBL
457
458	epssbc = TKEsbc * u_star ! -> e SBC
459	TKE_AT(ikl,mzp) = TKEsbc * sqrcmu ! TKE SBC
460
461	zeta = z__SBL / LMO_AT(ikl) ! zeta
462	sgnLMO = max(0.0,sign(un_1,LMO_AT(ikl))) !
463
464	eps_AT(ikl,mzp) = epssbc & !
465	& ( (sgnLMO (1.+A_Stab* zeta) & ! phim Stab.
466	& +(1.0-sgnLMO )/(1.-20.*min(0.,zeta))) & ! phim Inst.
467	& *vK_inv /z__SBL & !
468	& -vK_inv /LMO_AT(ikl))
469	! ... Duynkerke, 1988, JAS (45), (19) p. 869
470
471	! #KI eps_AT(ikl,mzp) = epssbc &
472	! #KI& *vK_inv /z__SBL
473
474
475	! When TKE Closure is Used, TKE is Modified near the Lower Boundary
476	! -----------------------------------------------------------------
477
478	! #KC se = max(0.,sign(un_1,TKE_AT(ikl,mzp-2)-TKE_AT(ikl,mzp-1)))
479	! #KC ke = mzp-1-se
480	! #KC TKE_AT(ikl,mzp-1)= TKE_AT(ikl,ke )
481	! #KC eps_AT(ikl,mzp-1)= eps_AT(ikl,ke )
482	! ... Schayes and Thunis, 1990, Contrib. 60 Inst.Astr.Geoph. p.8, 1.4.4.
483
484	! #KC TKE_AT(ikl,mzp-1)= TKE_AT(ikl,mzp )
485	! #KC eps_AT(ikl,mzp-1)= eps_AT(ikl,mzp )
486
487
488	! Upper Boundary Conditions
489	! =========================
490
491	TKE_AT(ikl,1) = TKE_AT(ikl,2)
492	eps_AT(ikl,1) = eps_AT(ikl,2)
493
494
495
496	! TKE-e Vertical Moving Average
497	! =============================
498
499	DO k= 1,mzp
500	TKEvav(k)= TKE_AT(ikl, k )
501	epsvav(k)= eps_AT(ikl, k )
502	END DO
503	IF (AT_vav) THEN
504	DO k= 2,mzp-1
505	TKEvav(k)=(dsigmi(k1p(k))*TKE_AT(ikl,k1p(k)) &
506	& +dsigmi( k )TKE_AT(ikl, k )2. &
507	& +dsigmi(k1m(k))*TKE_AT(ikl,k1m(k)) ) &
508	& /(dsigmi(k1p(k)) &
509	& +dsigmi( k ) *2. &
510	& +dsigmi(k1m(k)) )
511	epsvav(k)=(dsigmi(k1p(k))*eps_AT(ikl,k1p(k)) &
512	& +dsigmi( k )eps_AT(ikl, k )2. &
513	& +dsigmi(k1m(k))*eps_AT(ikl,k1m(k)) ) &
514	& /(dsigmi(k1p(k)) &
515	& +dsigmi( k ) *2. &
516	& +dsigmi(k1m(k)) )
517	END DO
518	END IF
519
520	! #KA DO k=mz__KA,mzp-1
521	! #KA TKE_AT(ikl,k)= TKEvav(k)
522	! #KA eps_AT(ikl,k)= epsvav(k)
523	! #KA END DO
524
525
526	! Verification: TKE must be Positive Definite
527	! ===========================================
528
529	DO k=1,mzp
530	TKE_AT(ikl,k)=max(eps6,TKE_AT(ikl,k))
531	eps_AT(ikl,k)=max(epsn,eps_AT(ikl,k))
532	TKEvav(k) =max(eps6,TKEvav(k) )
533	epsvav(k) =max(eps6,epsvav(k) )
534	END DO
535
536
537	! Minimum Energy Dissipation Time
538	! ===============================
539
540	IF (Ee_HuangRamn) THEN
541	TKE_ds = 0.0
542	eps_ds = 0.0
543	DO k=1,mzp
544	TKE_ds = TKE_ds + TKE_AT(ikl,k)*dsigma(k)
545	eps_ds = eps_ds + eps_AT(ikl,k)*dsigma(k)
546	END DO
547
548	IF (eps_ds.gt.0.0) THEN
549	edt_HR = betahr * TKE_ds /eps_ds
550	ELSE
551	edt_HR = epsn
552	! ... edt_HR set to an arbitrary small value
553	END IF
554	END IF
555
556
557	! Turbulent Diffusion Coefficients
558	! ================================
559
560	IF (it_EXP.gt.1) THEN
561
562
563	! Richardson Number (contributors)
564	! -----------------
565
566	Prandtl(mzp) = 1.
567
568	DO k=2,mzp-1
569	c_Buoy(k) = 0.0
570	! #RI c_Buoy(k) = (pkt ( k)-pkt ( k+1))*p0_kap &
571	! #RI& / TmidDY(ikl,k)
572	! #RI zShear(k) = max(dukkp1(k)2 +dvkkp1(k) 2 , eps6)
573
574
575	! Richardson Number
576	! -----------------
577
578	Ri__Nb(k) = 0.0
579	! #RI Ri__Nb(k) = (Grav_F/kkp1dz(k)) & ! g * dz (k+1/2)
580	! #RI& *c_Buoy(k) & ! d(theta)/T (k+1/2)
581	! #RI& /zShear(k) ! d\|V\|
582	END DO
583
584
585	! Diffusion Coefficient for Heat
586	! ------------------------------
587
588	DO k=2,mzp
589
590	IF (Kl_TherryLac) THEN
591	Kzh_AT(ikl,k)= 0.50 * sqrt(TKEvav(k ))/ML_inv(k)
592	! ... nu_t =c_mu X ECT X ECT / eps
593
594	ELSE IF (Ee_HuangRamn) THEN
595	Kzh_AT(ikl,k)= &
596	& cmu * TKE_AT(ikl,k)* TKE_AT(ikl,k)/(eps_AT(ikl,k) &
597	& +TKE_AT(ikl,k)/ edt_HR )
598	! ... nu_t =c_mu X ECT X ECT / eps
599
600	ELSE ! (Ee_Duynkerke / Ee_Kitada)
601	Kzh_AT(ikl,k)= &
602	& cmu * TKEvav(k )* TKEvav(k )/(epsvav(k ))
603	! ... nu_t =c_mu X ECT X ECT / eps
604
605	END IF
606
607	kz_max= vonKrm * Z___DY(ikl,k+1)-Z___DY(ikl,mzpp)
608	kz_mix= kz_max / (1. +kz_max *0.1)
609	kz2mix= kz_mix * kz_mix
610	KzhMAX= max( 5000. , 100. &
611	& * kz2mix *abs((WindDY (ikl,k)-WindDY (ikl,k1p(k))) &
612	& *kkp1dz(k) ))
613	Kzh_AT(ikl,k)= min( KzhMAX , Kzh_AT(ikl,k))
614	Kzh_AT(ikl,k)= max( A_MolV , Kzh_AT(ikl,k))
615	Kzh0AT(ikl,k)= Kzh_AT(ikl,k)
616	END DO
617
618
619
620	! Flux Limitor (Limitation of pkt Flux)
621	! ------------
622
623	IF (AT_LIM .AND. mzp.GE.3) THEN
624	DO k=min(3,mzp),mzp
625	IF (pkt_DY(ikl,k+1) .GT. pkt_DY(ikl,k )) THEN
626	Kz_MAX=( (Z___DY(ikl,k ) - Z___DY(ikl,k+1)) &
627	& /(pkt_DY(ikl,k+1) - pkt_DY(ikl,k ))) &
628	& (Dpkt_X 0.5 *(Z___DY(ikl,k-1) - Z___DY(ikl,k+1)) &
629	& -Kzh_AT(ikl,k-1)*(pkt_DY(ikl,k-1) - pkt_DY(ikl,k )) &
630	& /(Z___DY(ikl,k-1) - Z___DY(ikl,k )))
631	ELSE IF (pkt_DY(ikl,k+1) .LT. pkt_DY(ikl,k )) THEN
632	Kz_MAX=( (Z___DY(ikl,k ) - Z___DY(ikl,k+1)) &
633	& /(pkt_DY(ikl,k ) - pkt_DY(ikl,k+1))) &
634	& (Dpkt_X 0.5 *(Z___DY(ikl,k-1) - Z___DY(ikl,k+1)) &
635	& +Kzh_AT(ikl,k-1)*(pkt_DY(ikl,k-1) - pkt_DY(ikl,k )) &
636	& /(Z___DY(ikl,k-1) - Z___DY(ikl,k )))
637	END IF
638	Kzh_AT(ikl,k) = min(Kzh_AT(ikl,k),Kz_MAX)
639	Kzh_AT(ikl,k) = max(Kzh_AT(ikl,k),A_MolV)
640	ENDDO
641	ENDIF
642
643
644
645	! Prandtl Number (Sukoriansky et al., 2005,
646	! -------------- BLM 117: 231-257, Eq.15, 19, 20 & Fig.2)
647
648	DO k=2,mzp-1
649	! #RI fac_Ri= 5.0 * max(Ri__Nb(i,j,k), eps6)
650	! #RI vuzvun= 0.4 *(1.-(fac_Ri-1./fac_Ri)/(fac_Ri+1./fac_Ri)) + 0.2
651	! #RI fac_Ri= 4.2 * max(Ri__Nb(i,j,k), eps6)
652	! #RI Kz_vun= 0.7 *(1.-(fac_Ri-1./fac_Ri)/(fac_Ri+1./fac_Ri))
653	Prandtl(k) = 1.
654	! #RI Prandtl(k) = max(0. ,sign(1. , Kzh_AT(ikl,k)-0.20)) &
655	! #RI& - min(0. ,sign(1. , Kzh_AT(ikl,k)-0.20)) &
656	! #RI& * min(vuzvun / max(eps6,Kz_vun), 20.00)
657	END DO
658
659
660	! Diffusion Coefficient for Momentum
661	! ----------------------------------
662
663	DO k=2,mzp
664	IF (Kl_TherryLac) THEN
665	Kzm_AT(ikl,k)= 0.7 * Kzh_AT(ikl,k)
666
667	ELSE
668	Kzm_AT(ikl,k)= Kzh_AT(ikl,k)
669	! ... cfr Machiels, 1992, TFE (FSA/UCL) (3.21) p.21
670
671	! #RI Kzm_AT(ikl,k)= Prandtl(k) * Kzh_AT(ikl,k)
672	END IF
673
674	END DO
675
676	END IF
677
678
679
680
681	! Work Arrays Reset
682	! =================
683
684	DO k=1,mzp
685	TrT_AT(ikl,k) = 0.0
686	END DO
687
688
689
690
691	! +++++++++++++++++++
692	ENDDO ! ikl=1,kcolp
693	! +++++++++++++++++++
694
695
696
697
698	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
699	! !
700	! DE-ALLOCATION
701	! =============
702
703	IF (FlagDALLOC) THEN !
704	deallocate ( dukkp1 ) ! Difference (u(k) - u(k+1)) [m/s]
705	deallocate ( dvkkp1 ) ! Difference (v(k) - v(k+1)) [m/s]
706	deallocate ( kkp1dz ) ! 1 / Difference (Z(k) - Z(k+1))
707	deallocate ( zShear ) ! Wind Shear Contribution to TKE [m2/s3]
708	deallocate ( REq_PT ) ! Reduced (Equivalent) Potential Temperature [K]
709	deallocate ( c_Buoy ) ! Buoyancy Coefficient (g/theta) X (dtheta/dz) [1/s2]
710	deallocate ( Ri__Nb ) ! Richardson Number [-]
711	deallocate ( Prandtl) ! Prandtl Number (Kzm/Kzh) [-]
712	deallocate ( Ls_inv ) ! 1 / Ls (Therry & Lacarrere, 1983) [1/m]
713	deallocate ( ML_inv ) ! 1 / ML (Mixing Length, Therry & Lacarr, 1983) [1/m]
714	deallocate ( DL_inv ) ! 1 / DL (Dissipation Length, Therry & Lacarr, 1983) [1/m]
715	deallocate ( Dissip ) ! Dissipation [m2/s3]
716	deallocate ( TKEvav ) ! TKE Vertical moving Average [m2/s2]
717	deallocate ( epsvav ) ! Dissipation Vertical moving Average [m2/s3]
718	deallocate ( pkt ) ! Reduced Potential Temperature [X]
719	END IF !
720	! !
721	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
722
723
724
725	return
726
727	end subroutine PHY_genTKE_RUN

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format