Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

PHY_genTKE_RUN.f90 @ 5412

Last change on this file since 5412 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

Rev	Line
[2089]	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