Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

lmdz_atke_exchange_coeff.f90 @ 5354

Last change on this file since 5354 was 5268, checked in by abarral, 7 weeks ago
.f90 <-> .F90 depending on cpp key use
File size: 22.7 KB

Rev	Line
[4687]	1	module lmdz_atke_exchange_coeff
[4449]	2
	3	implicit none
	4
	5	contains
	6
[4631]	7	subroutine atke_compute_km_kh(ngrid,nlay,dtime, &
[4653]	8	wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv, &
[5039]	9	tke,eps,tke_shear,tke_buoy,tke_trans,Km_out,Kh_out)
[4449]	10
	11	!========================================================================
	12	! Routine that computes turbulent Km / Kh coefficients with a
	13	! 1.5 order closure scheme (TKE) with or without stationarity assumption
	14	!
	15	! This parameterization has been constructed in the framework of a
	16	! collective and collaborative workshop,
	17	! the so-called 'Atelier TKE (ATKE)' with
[4478]	18	! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, T Dubos,
[4449]	19	! M. Coulon-Decorzens, S. Fromang, G. Riviere, A. Sima, F. Hourdin, E. Vignon
	20	!
	21	! Main assumptions of the model :
[4714]	22	! (1) horizontal homogeneity (Dx=Dy=0.)
[4449]	23	!=======================================================================
	24
[4780]	25	USE lmdz_atke_turbulence_ini, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual, ri0, ri1
[4745]	26	USE lmdz_atke_turbulence_ini, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes, rg, rd, rv, atke_ok_vdiff
[4780]	27	USE lmdz_atke_turbulence_ini, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin, cn
[4449]	28
[4804]	29	!!-------------------------------------------------------------------------------------------------------------
	30	! integer :: iflag_atke ! flag that controls options in atke_compute_km_kh
	31	! integer :: iflag_atke_lmix ! flag that controls the calculation of mixing length in atke
	32	! integer :: iflag_num_atke ! flag that controls the numerical treatment of diffusion coeffiient calculation
	33
	34	! logical :: atke_ok_vdiff ! activate vertical diffusion of TKE or not
	35	! logical :: atke_ok_virtual ! account for vapor for flottability
	36
	37	! real :: kappa = 0.4 ! Von Karman constant
	38	! real :: cn ! Sm value at Ri=0
	39	! real :: cinf ! Sm value at Ri=-Inf
	40	! real :: ri0 ! Richardson number near zero to guarantee continuity in slope of Sm (stability function) at Ri=0
	41	! real :: ri1 ! Richardson number near zero to guarantee continuity in slope of Pr (Prandlt's number) at Ri=0
	42	! real :: lmin ! minimum mixing length corresponding to the Kolmogov dissipation scale in planetary atmospheres (Chen et al 2016, JGR Atmos)
	43	! real :: ctkes ! coefficient for surface TKE
	44	! real :: clmixshear ! coefficient for mixing length depending on local wind shear
	45
	46	! Tunable parameters for the ATKE scheme and their range of values
	47	!!-------------------------------------------------------------------------------------------------------------
	48	! real :: cepsilon ! controls the value of the dissipation length scale, range [1.2 - 10]
	49	! real :: cke ! controls the value of the diffusion coefficient of TKE, range [1 - 5]
	50	! real :: l0 ! asymptotic mixing length far from the ground [m] (Sun et al 2011, JAMC), range [15 - 75]
	51	! real :: clmix ! controls the value of the mixing length in stratified conditions, range [0.1 - 2]
	52	! real :: ric ! critical Richardson number controlling the slope of Sm in stable conditions, range [0.19 - 0.25]
	53	! real :: smmin ! minimum value of Sm in very stable conditions (defined here as minsqrt(Ez/Ek)) at large Ri, range [0.025 - 0.1]
	54	! real :: pr_neut ! neutral value of the Prandtl number (Ri=0), range [0.7 - 1]
	55	! real :: pr_slope ! linear slope of Pr with Ri in the very stable regime, range [3 - 5]
	56	! real :: cinffac ! ratio between cinf and cn controlling the convective limit of Sm, range [1.2 - 5.0]
	57	! real :: pr_asym ! value of Prandlt in the convective limit(Ri=-Inf), range [0.3 - 0.5]
	58	!!-------------------------------------------------------------------------------------------------------------
	59
[4449]	60	implicit none
	61
	62
	63	! Declarations:
	64	!=============
	65
	66	INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
[4631]	67	INTEGER, INTENT(IN) :: nlay ! number of vertical index
[4449]	68
[4631]	69	REAL, INTENT(IN) :: dtime ! physics time step (s)
[4449]	70	REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_u ! zonal velocity (m/s)
	71	REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_v ! meridional velocity (m/s)
[4881]	72	REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K)
	73	REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: qvap ! specific humidity (kg/kg)
	74	REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa)
	75	REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: pinterf ! pressure at interfaces(Pa)
	76	REAL, DIMENSION(ngrid), INTENT(IN) :: cdrag_uv ! surface drag coefficient for momentum
[4449]	77
[4881]	78	REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers (m2/s2)
[4449]	79
[4881]	80	REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: eps ! output: TKE dissipation rate at interface between layers (m2/s3)
[5039]	81	REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: tke_shear! output: TKE shear production rate (m2/s3)
	82	REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: tke_buoy ! output: TKE buoyancy production rate (m2/s3)
	83	REAL, DIMENSION(ngrid,nlay+1), INTENT(OUT) :: tke_trans! output: TKE transport (diffusion) term (m2/s3)
[4881]	84	REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Km_out ! output: Exchange coefficient for momentum at interface between layers (m2/s)
	85	REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Kh_out ! output: Exchange coefficient for heat flux at interface between layers (m2/s)
[4449]	86
	87	! Local variables
[4478]	88	REAL, DIMENSION(ngrid,nlay+1) :: Km ! Exchange coefficient for momentum at interface between layers
	89	REAL, DIMENSION(ngrid,nlay+1) :: Kh ! Exchange coefficient for heat flux at interface between layers
	90	REAL, DIMENSION(ngrid,nlay) :: theta ! Potential temperature
	91	REAL, DIMENSION(ngrid,nlay+1) :: l_exchange ! Length of exchange (at interface)
	92	REAL, DIMENSION(ngrid,nlay+1) :: z_interf ! Altitude at the interface
	93	REAL, DIMENSION(ngrid,nlay) :: z_lay ! Altitude of layers
	94	REAL, DIMENSION(ngrid,nlay) :: dz_interf ! distance between two consecutive interfaces
	95	REAL, DIMENSION(ngrid,nlay) :: dz_lay ! distance between two layer middles (NB: first and last are half layers)
[4631]	96	REAL, DIMENSION(ngrid,nlay+1) :: N2 ! square of Brunt Vaisala pulsation (at interface)
	97	REAL, DIMENSION(ngrid,nlay+1) :: shear2 ! square of wind shear (at interface)
[4478]	98	REAL, DIMENSION(ngrid,nlay+1) :: Ri ! Richardson's number (at interface)
	99	REAL, DIMENSION(ngrid,nlay+1) :: Prandtl ! Turbulent Prandtl's number (at interface)
	100	REAL, DIMENSION(ngrid,nlay+1) :: Sm ! Stability function for momentum (at interface)
	101	REAL, DIMENSION(ngrid,nlay+1) :: Sh ! Stability function for heat (at interface)
[4449]	102
	103	INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid)
[4478]	104	REAL :: preff ! reference pressure for potential temperature calculations
	105	REAL :: thetam ! mean potential temperature at interface
[4631]	106	REAL :: delta ! discriminant of the second order polynomial
	107	REAL :: qq ! tke=qq**2/2
	108	REAL :: shear ! wind shear
	109	REAL :: lstrat ! mixing length depending on local stratification
	110	REAL :: taustrat ! caracteristic timescale for turbulence in very stable conditions
[4644]	111	REAL :: netloss ! net loss term of tke
	112	REAL :: netsource ! net source term of tke
	113	REAL :: ustar ! friction velocity estimation
[4653]	114	REAL :: invtau
	115	REAL :: rvap
[4449]	116
	117	! Initializations:
	118	!================
	119
	120	DO igrid=1,ngrid
[4478]	121	dz_interf(igrid,1) = 0.0
[4449]	122	z_interf(igrid,1) = 0.0
	123	END DO
	124
[4653]	125	! Calculation of potential temperature: (if vapor -> virtual potential temperature)
[4478]	126	!=====================================
[4449]	127
[4478]	128	preff=100000.
[4653]	129	! results should not depend on the choice of preff
[4478]	130	DO ilay=1,nlay
[4804]	131	DO igrid = 1, ngrid
[4478]	132	theta(igrid,ilay)=temp(igrid,ilay)(preff/play(igrid,ilay))*(rd/rcpd)
[4804]	133	END DO
[4478]	134	END DO
[4449]	135
[4653]	136	! account for water vapor mass for buoyancy calculation
	137	IF (atke_ok_virtual) THEN
[4804]	138	DO ilay=1,nlay
	139	DO igrid = 1, ngrid
	140	rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay)))
	141	theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap)
	142	END DO
	143	END DO
[4653]	144	ENDIF
[4449]	145
[4478]	146
	147	! Calculation of altitude of layers' middle and bottom interfaces:
	148	!=================================================================
	149
[4449]	150	DO ilay=2,nlay+1
	151	DO igrid=1,ngrid
[4478]	152	dz_interf(igrid,ilay-1) = rdtemp(igrid,ilay-1)/rg/play(igrid,ilay-1)(pinterf(igrid,ilay-1)-pinterf(igrid,ilay))
	153	z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1)
[4449]	154	ENDDO
	155	ENDDO
	156
[4478]	157	DO ilay=1,nlay
[4804]	158	DO igrid=1,ngrid
	159	z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay))
	160	ENDDO
[4478]	161	ENDDO
[4449]	162
[4478]	163
[4449]	164	! Computes the gradient Richardson's number and stability functions:
	165	!===================================================================
	166
[4478]	167	DO ilay=2,nlay
[4449]	168	DO igrid=1,ngrid
[4478]	169	dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1)
	170	thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1))
[4631]	171	N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay)
	172	shear2(igrid,ilay)= (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + &
	173	((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 )
	174	Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10)
[4481]	175
	176	IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases
[4478]	177	Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn
	178	Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut
[4481]	179	ELSE ! stable cases
[4644]	180	Sm(igrid,ilay) = max(smmin,cn*(1.-Ri(igrid,ilay)/Ric))
[4688]	181	! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019
	182	Prandtl(igrid,ilay) = pr_neutexp(-pr_slope/pr_neutRi(igrid,ilay)+Ri(igrid,ilay)/pr_neut) &
	183	+ Ri(igrid,ilay) * pr_slope
[4481]	184	IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN
[4804]	185	call abort_physic("atke_compute_km_kh", &
	186	'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1)
[4481]	187	ENDIF
[4449]	188	END IF
	189
	190	Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay)
	191
	192	ENDDO
	193	ENDDO
	194
[4631]	195
	196	! Computing the mixing length:
	197	!==============================================================
	198
	199
	200	IF (iflag_atke_lmix .EQ. 1 ) THEN
[4804]	201	! Blackadar formulation (~kappa l) + buoyancy length scale (Deardoff 1980) for very stable conditions
	202	DO ilay=2,nlay
	203	DO igrid=1,ngrid
	204	l_exchange(igrid,ilay) = kappal0z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
	205	IF (N2(igrid,ilay) .GT. 0.) THEN
	206	lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay))
	207	lstrat=max(lstrat,lmin)
	208	!Inverse interpolation, Van de Wiel et al. 2010
	209	l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
	210	ENDIF
	211	ENDDO
	212	ENDDO
[4631]	213
[4804]	214	ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN
	215	! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms
	216	! implies 2 tuning coefficients clmix and clmixshear
	217	DO ilay=2,nlay
	218	DO igrid=1,ngrid
	219	l_exchange(igrid,ilay) = kappal0z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
	220	IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN
	221	lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), &
[4663]	222	clmixshear*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay)))
[4804]	223	lstrat=max(lstrat,lmin)
	224	!Inverse interpolation, Van de Wiel et al. 2010
	225	l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
	226	ENDIF
	227	ENDDO
	228	ENDDO
[4632]	229
[4804]	230	ELSE IF (iflag_atke_lmix .EQ. 3 ) THEN
	231	! add effect of wind shear on lstrat following grisogono 2010, qjrms
	232	! keeping a single tuning coefficient clmix
	233	DO ilay=2,nlay
	234	DO igrid=1,ngrid
	235	l_exchange(igrid,ilay) = kappal0z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
	236	IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN
	237	lstrat=clmixsqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))(1.0+Ri(igrid,ilay)/(2.*Prandtl(igrid,ilay)))
	238	lstrat=max(lstrat,lmin)
	239	!Inverse interpolation, Van de Wiel et al. 2010
	240	l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0)
	241	ENDIF
	242	ENDDO
	243	ENDDO
[4632]	244
[4631]	245	ELSE
[4804]	246	! default Blackadar formulation: neglect effect of local stratification and shear
	247	DO ilay=2,nlay+1
	248	DO igrid=1,ngrid
	249	l_exchange(igrid,ilay) = kappal0z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
	250	ENDDO
	251	ENDDO
[4631]	252	ENDIF
	253
	254
[4644]	255	! Computing the TKE k>=2:
	256	!========================
[4449]	257	IF (iflag_atke == 0) THEN
	258
[4804]	259	! stationary solution (dtke/dt=0)
[4644]	260
[4804]	261	DO ilay=2,nlay
	262	DO igrid=1,ngrid
	263	tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)*2 Sm(igrid,ilay) * &
	264	shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))
[4881]	265	eps(igrid,ilay) = (tke(igrid,ilay)*(3./2))/(cepsilonl_exchange(igrid,ilay))
[5039]	266	tke_shear(igrid,ilay)=l_exchange(igrid,ilay)Sm(igrid,ilay)sqrt(tke(igrid,ilay))*shear2(igrid,ilay)
	267	tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)Sm(igrid,ilay)sqrt(tke(igrid,ilay))*shear2(igrid,ilay) &
	268	*(Ri(igrid,ilay) / Prandtl(igrid,ilay))
[4804]	269	ENDDO
[4449]	270	ENDDO
	271
[4631]	272	ELSE IF (iflag_atke == 1) THEN
	273
[4804]	274	! full implicit scheme resolved with a second order polynomial equation
	275	! default solution which shows fair convergence properties
	276	DO ilay=2,nlay
	277	DO igrid=1,ngrid
	278	qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
	279	delta=(2.sqrt(2.)cepsilonl_exchange(igrid,ilay)/dtime)*2. &
	280	+4.(2.sqrt(2.)cepsilonl_exchange(igrid,ilay)/dtime*qq + &
	281	2.l_exchange(igrid,ilay)l_exchange(igrid,ilay)cepsilonSm(igrid,ilay) &
	282	shear2(igrid,ilay) (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)))
	283	qq=(-2.sqrt(2.)cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2.
	284	qq=max(0.,qq)
[4881]	285	tke(igrid,ilay)=0.5(qq*2)
[5039]	286	eps(igrid,ilay) = (tke(igrid,ilay)*(3./2))/(cepsilonl_exchange(igrid,ilay))
	287	tke_shear(igrid,ilay)=l_exchange(igrid,ilay)Sm(igrid,ilay)sqrt(tke(igrid,ilay))*shear2(igrid,ilay)
	288	tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)Sm(igrid,ilay)sqrt(tke(igrid,ilay))*shear2(igrid,ilay) &
	289	*(Ri(igrid,ilay) / Prandtl(igrid,ilay))
[4804]	290	ENDDO
[4631]	291	ENDDO
	292
[4644]	293
[4631]	294	ELSE IF (iflag_atke == 2) THEN
	295
[4804]	296	! semi implicit scheme when l does not depend on tke
	297	! positive-guaranteed if pr slope in stable condition >1
[4644]	298
[4804]	299	DO ilay=2,nlay
	300	DO igrid=1,ngrid
	301	qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
	302	qq=(qq+l_exchange(igrid,ilay)Sm(igrid,ilay)dtime/sqrt(2.) &
	303	shear2(igrid,ilay)(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
[5039]	304	/(1.+qqdtime/(cepsilonl_exchange(igrid,ilay)2.sqrt(2.)))
	305	tke_shear(igrid,ilay)=l_exchange(igrid,ilay)Sm(igrid,ilay)qq/sqrt(2.)*shear2(igrid,ilay)
	306	tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)Sm(igrid,ilay)qq/sqrt(2.)*shear2(igrid,ilay) &
	307	*(Ri(igrid,ilay) / Prandtl(igrid,ilay))
[4804]	308	tke(igrid,ilay)=0.5(qq*2)
[4881]	309	eps(igrid,ilay) = (tke(igrid,ilay)*(3./2))/(cepsilonl_exchange(igrid,ilay))
[4804]	310	ENDDO
[4644]	311	ENDDO
	312
	313
	314	ELSE IF (iflag_atke == 3) THEN
[4804]	315	! numerical resolution adapted from that in MAR (Deleersnijder 1992)
	316	! positively defined by construction
[4644]	317
[4804]	318	DO ilay=2,nlay
	319	DO igrid=1,ngrid
[4881]	320	eps(igrid,ilay) = (tke(igrid,ilay)*(3./2))/(cepsilonl_exchange(igrid,ilay))
[4804]	321	qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
[5039]	322	tke_shear(igrid,ilay)=l_exchange(igrid,ilay)Sm(igrid,ilay)qq/sqrt(2.)*shear2(igrid,ilay)
	323	tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)Sm(igrid,ilay)qq/sqrt(2.)*shear2(igrid,ilay) &
	324	*(Ri(igrid,ilay) / Prandtl(igrid,ilay))
[4804]	325	IF (Ri(igrid,ilay) .LT. 0.) THEN
	326	netloss=qq/(2.sqrt(2.)cepsilon*l_exchange(igrid,ilay))
	327	netsource=l_exchange(igrid,ilay)Sm(igrid,ilay)/sqrt(2.)shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))
	328	ELSE
	329	netloss=qq/(2.sqrt(2.)cepsilon*l_exchange(igrid,ilay))+ &
	330	l_exchange(igrid,ilay)Sm(igrid,ilay)/sqrt(2.)N2(igrid,ilay)/Prandtl(igrid,ilay)
	331	netsource=l_exchange(igrid,ilay)Sm(igrid,ilay)/sqrt(2.)shear2(igrid,ilay)
	332	ENDIF
	333	qq=((qq*2)/dtime+qqnetsource)/(qq/dtime+netloss)
	334	tke(igrid,ilay)=0.5(qq*2)
	335	ENDDO
[4631]	336	ENDDO
	337
[4644]	338	ELSE IF (iflag_atke == 4) THEN
[4804]	339	! semi implicit scheme from Arpege (V. Masson methodology with
	340	! Taylor expansion of the dissipation term)
	341	DO ilay=2,nlay
	342	DO igrid=1,ngrid
	343	qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
[5039]	344	tke_shear(igrid,ilay)=l_exchange(igrid,ilay)Sm(igrid,ilay)qq/sqrt(2.)*shear2(igrid,ilay)
	345	tke_buoy(igrid,ilay)=-l_exchange(igrid,ilay)Sm(igrid,ilay)qq/sqrt(2.)*shear2(igrid,ilay) &
	346	*(Ri(igrid,ilay) / Prandtl(igrid,ilay))
[4804]	347	qq=(l_exchange(igrid,ilay)Sm(igrid,ilay)/sqrt(2.)shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
	348	+qq(1.+dtimeqq/(cepsilonl_exchange(igrid,ilay)2.*sqrt(2.)))) &
	349	/(1.+2.qqdtime/(cepsilonl_exchange(igrid,ilay)2.*sqrt(2.)))
	350	qq=max(0.,qq)
	351	tke(igrid,ilay)=0.5(qq*2)
[4881]	352	eps(igrid,ilay) = (tke(igrid,ilay)*(3./2))/(cepsilonl_exchange(igrid,ilay))
[4804]	353	ENDDO
[4631]	354	ENDDO
	355
	356
	357	ELSE
[4804]	358	call abort_physic("atke_compute_km_kh", &
[4631]	359	'numerical treatment of TKE not possible yet', 1)
[4449]	360
	361	END IF
	362
[4644]	363	! We impose a 0 tke at nlay+1
	364	!==============================
[4449]	365
[4644]	366	DO igrid=1,ngrid
[4804]	367	tke(igrid,nlay+1)=0.
[4884]	368	eps(igrid,nlay+1)=0.
[5039]	369	tke_shear(igrid,nlay+1)=0.
	370	tke_buoy(igrid,nlay+1)=0.
[4644]	371	END DO
	372
	373
	374	! Calculation of surface TKE (k=1)
	375	!=================================
	376	! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note)
	377	DO igrid=1,ngrid
[4804]	378	ustar=sqrt(cdrag_uv(igrid)(wind_u(igrid,1)2+wind_v(igrid,1)*2))
	379	tke(igrid,1)=ctkes(ustar*2)
[4884]	380	eps(igrid,1)=0. ! arbitrary as TKE is not properly defined at the surface
[5039]	381	tke_shear(igrid,1)=0.
	382	tke_buoy(igrid,1)=0.
[4644]	383	END DO
	384
	385
	386	! vertical diffusion of TKE
	387	!==========================
[5039]	388	tke_trans(:,:)=0.
[4653]	389	IF (atke_ok_vdiff) THEN
[5039]	390	CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke,tke_trans)
[4644]	391	ENDIF
	392
	393
[4449]	394	! Computing eddy diffusivity coefficients:
	395	!========================================
[4478]	396	DO ilay=2,nlay ! TODO: also calculate for nlay+1 ?
[4449]	397	DO igrid=1,ngrid
[4478]	398	! we add the molecular viscosity to Km,h
	399	Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5
	400	Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5
[4449]	401	END DO
	402	END DO
	403
[4478]	404	! for output:
	405	!===========
	406	Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay)
	407	Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay)
[4449]	408
	409	end subroutine atke_compute_km_kh
	410
[4644]	411	!===============================================================================================
[5039]	412	subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke,tke_trans)
[4449]	413
[4644]	414	! routine that computes the vertical diffusion of TKE by the turbulence
[4653]	415	! using an implicit resolution (See note by Dufresne and Ghattas (2009))
[4644]	416	! E Vignon, July 2023
	417
[4687]	418	USE lmdz_atke_turbulence_ini, ONLY : rd, cke, viscom
[4644]	419
	420
	421	INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
	422	INTEGER, INTENT(IN) :: nlay ! number of vertical index
	423
	424	REAL, INTENT(IN) :: dtime ! physics time step (s)
	425	REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: z_lay ! altitude of mid-layers (m)
	426	REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: z_interf ! altitude of bottom interfaces (m)
	427	REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K)
	428	REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa)
	429	REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: l_exchange ! mixing length at interfaces between layers
	430	REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: Sm ! stability function for eddy diffusivity for momentum at interface between layers
	431
	432	REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers
[5039]	433	REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke_trans ! turbulent kinetic energy transport term (m2/s3)
[4644]	434
	435
	436	INTEGER :: igrid,ilay
	437	REAL, DIMENSION(ngrid,nlay+1) :: Ke ! eddy diffusivity for TKE
	438	REAL, DIMENSION(ngrid,nlay+1) :: dtke
	439	REAL, DIMENSION(ngrid,nlay+1) :: ak, bk, ck, CCK, DDK
	440	REAL :: gammak,Kem,KKb,KKt
	441
	442
	443	! Few initialisations
	444	CCK(:,:)=0.
	445	DDK(:,:)=0.
	446	dtke(:,:)=0.
	447
	448
	449	! Eddy diffusivity for TKE
	450
	451	DO ilay=2,nlay
	452	DO igrid=1,ngrid
[4804]	453	Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)Sm(igrid,ilay)sqrt(tke(igrid,ilay)))*cke
[4644]	454	ENDDO
	455	ENDDO
	456	! at the top of the atmosphere set to 0
	457	Ke(:,nlay+1)=0.
	458	! at the surface, set it equal to that at the first model level
	459	Ke(:,1)=Ke(:,2)
	460
	461
	462	! calculate intermediary variables
	463
	464	DO ilay=2,nlay
	465	DO igrid=1,ngrid
	466	Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay))
	467	KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay))
	468	Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1))
	469	KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1))
	470	gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1))
	471	ak(igrid,ilay)=-gammakdtimeKKb
	472	ck(igrid,ilay)=-gammakdtimeKKt
	473	bk(igrid,ilay)=1.+gammakdtime(KKt+KKb)
	474	ENDDO
	475	ENDDO
	476
	477	! calculate CCK and DDK coefficients
	478	! downhill phase
	479
	480	DO igrid=1,ngrid
[4804]	481	CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay)
	482	DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay)
[4644]	483	ENDDO
	484
	485
	486	DO ilay=nlay-1,2,-1
	487	DO igrid=1,ngrid
	488	CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) &
[4804]	489	/ (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
[4644]	490	DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
	491	ENDDO
	492	ENDDO
	493
	494	! calculate TKE
	495	! uphill phase
	496
	497	DO ilay=2,nlay+1
	498	DO igrid=1,ngrid
	499	dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay)
	500	ENDDO
	501	ENDDO
	502
	503	! update TKE
	504	tke(:,:)=tke(:,:)+dtke(:,:)
[5039]	505	tke_trans(:,:)=dtke(:,:)/dtime
[4644]	506
	507
	508	end subroutine atke_vdiff_tke
	509
	510
	511
[4687]	512	end module lmdz_atke_exchange_coeff

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

Download in other formats:

Original Format