1 | module lmdz_atke_exchange_coeff |
---|
2 | |
---|
3 | implicit none |
---|
4 | |
---|
5 | contains |
---|
6 | |
---|
7 | subroutine atke_compute_km_kh(ngrid,nlay,dtime, & |
---|
8 | wind_u,wind_v,temp,qvap,play,pinterf,cdrag_uv, & |
---|
9 | tke,Km_out,Kh_out) |
---|
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 |
---|
18 | ! K. Arjdal, L. Raillard, C. Dehondt, P. Tiengou, A. Spiga, F. Cheruy, T Dubos, |
---|
19 | ! M. Coulon-Decorzens, S. Fromang, G. Riviere, A. Sima, F. Hourdin, E. Vignon |
---|
20 | ! |
---|
21 | ! Main assumptions of the model : |
---|
22 | ! (1) horizontal homogeneity (Dx=Dy=0.) |
---|
23 | !======================================================================= |
---|
24 | |
---|
25 | USE lmdz_atke_turbulence_ini, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, atke_ok_virtual, ri0, ri1 |
---|
26 | USE lmdz_atke_turbulence_ini, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes, rg, rd, rv, atke_ok_vdiff |
---|
27 | USE lmdz_atke_turbulence_ini, ONLY : viscom, viscoh, clmix, clmixshear, iflag_atke_lmix, lmin, smmin, cn |
---|
28 | |
---|
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 | |
---|
60 | implicit none |
---|
61 | |
---|
62 | |
---|
63 | ! Declarations: |
---|
64 | !============= |
---|
65 | |
---|
66 | INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) |
---|
67 | INTEGER, INTENT(IN) :: nlay ! number of vertical index |
---|
68 | |
---|
69 | REAL, INTENT(IN) :: dtime ! physics time step (s) |
---|
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) |
---|
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 |
---|
77 | |
---|
78 | REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers |
---|
79 | |
---|
80 | REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Km_out ! output: Exchange coefficient for momentum at interface between layers |
---|
81 | REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Kh_out ! output: Exchange coefficient for heat flux at interface between layers |
---|
82 | |
---|
83 | ! Local variables |
---|
84 | REAL, DIMENSION(ngrid,nlay+1) :: Km ! Exchange coefficient for momentum at interface between layers |
---|
85 | REAL, DIMENSION(ngrid,nlay+1) :: Kh ! Exchange coefficient for heat flux at interface between layers |
---|
86 | REAL, DIMENSION(ngrid,nlay) :: theta ! Potential temperature |
---|
87 | REAL, DIMENSION(ngrid,nlay+1) :: l_exchange ! Length of exchange (at interface) |
---|
88 | REAL, DIMENSION(ngrid,nlay+1) :: z_interf ! Altitude at the interface |
---|
89 | REAL, DIMENSION(ngrid,nlay) :: z_lay ! Altitude of layers |
---|
90 | REAL, DIMENSION(ngrid,nlay) :: dz_interf ! distance between two consecutive interfaces |
---|
91 | REAL, DIMENSION(ngrid,nlay) :: dz_lay ! distance between two layer middles (NB: first and last are half layers) |
---|
92 | REAL, DIMENSION(ngrid,nlay+1) :: N2 ! square of Brunt Vaisala pulsation (at interface) |
---|
93 | REAL, DIMENSION(ngrid,nlay+1) :: shear2 ! square of wind shear (at interface) |
---|
94 | REAL, DIMENSION(ngrid,nlay+1) :: Ri ! Richardson's number (at interface) |
---|
95 | REAL, DIMENSION(ngrid,nlay+1) :: Prandtl ! Turbulent Prandtl's number (at interface) |
---|
96 | REAL, DIMENSION(ngrid,nlay+1) :: Sm ! Stability function for momentum (at interface) |
---|
97 | REAL, DIMENSION(ngrid,nlay+1) :: Sh ! Stability function for heat (at interface) |
---|
98 | |
---|
99 | INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid) |
---|
100 | REAL :: preff ! reference pressure for potential temperature calculations |
---|
101 | REAL :: thetam ! mean potential temperature at interface |
---|
102 | REAL :: delta ! discriminant of the second order polynomial |
---|
103 | REAL :: qq ! tke=qq**2/2 |
---|
104 | REAL :: shear ! wind shear |
---|
105 | REAL :: lstrat ! mixing length depending on local stratification |
---|
106 | REAL :: taustrat ! caracteristic timescale for turbulence in very stable conditions |
---|
107 | REAL :: netloss ! net loss term of tke |
---|
108 | REAL :: netsource ! net source term of tke |
---|
109 | REAL :: ustar ! friction velocity estimation |
---|
110 | REAL :: invtau |
---|
111 | REAL :: rvap |
---|
112 | |
---|
113 | ! Initializations: |
---|
114 | !================ |
---|
115 | |
---|
116 | DO igrid=1,ngrid |
---|
117 | dz_interf(igrid,1) = 0.0 |
---|
118 | z_interf(igrid,1) = 0.0 |
---|
119 | END DO |
---|
120 | |
---|
121 | ! Calculation of potential temperature: (if vapor -> virtual potential temperature) |
---|
122 | !===================================== |
---|
123 | |
---|
124 | preff=100000. |
---|
125 | ! results should not depend on the choice of preff |
---|
126 | DO ilay=1,nlay |
---|
127 | DO igrid = 1, ngrid |
---|
128 | theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd) |
---|
129 | END DO |
---|
130 | END DO |
---|
131 | |
---|
132 | ! account for water vapor mass for buoyancy calculation |
---|
133 | IF (atke_ok_virtual) THEN |
---|
134 | DO ilay=1,nlay |
---|
135 | DO igrid = 1, ngrid |
---|
136 | rvap=max(0.,qvap(igrid,ilay)/(1.-qvap(igrid,ilay))) |
---|
137 | theta(igrid,ilay)=theta(igrid,ilay)*(1.+rvap/(RD/RV))/(1.+rvap) |
---|
138 | END DO |
---|
139 | END DO |
---|
140 | ENDIF |
---|
141 | |
---|
142 | |
---|
143 | ! Calculation of altitude of layers' middle and bottom interfaces: |
---|
144 | !================================================================= |
---|
145 | |
---|
146 | DO ilay=2,nlay+1 |
---|
147 | DO igrid=1,ngrid |
---|
148 | dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay)) |
---|
149 | z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1) |
---|
150 | ENDDO |
---|
151 | ENDDO |
---|
152 | |
---|
153 | DO ilay=1,nlay |
---|
154 | DO igrid=1,ngrid |
---|
155 | z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay)) |
---|
156 | ENDDO |
---|
157 | ENDDO |
---|
158 | |
---|
159 | |
---|
160 | ! Computes the gradient Richardson's number and stability functions: |
---|
161 | !=================================================================== |
---|
162 | |
---|
163 | DO ilay=2,nlay |
---|
164 | DO igrid=1,ngrid |
---|
165 | dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1) |
---|
166 | thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1)) |
---|
167 | N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay) |
---|
168 | shear2(igrid,ilay)= (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & |
---|
169 | ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 ) |
---|
170 | Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10) |
---|
171 | |
---|
172 | IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases |
---|
173 | Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn |
---|
174 | Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut |
---|
175 | ELSE ! stable cases |
---|
176 | Sm(igrid,ilay) = max(smmin,cn*(1.-Ri(igrid,ilay)/Ric)) |
---|
177 | ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019 |
---|
178 | Prandtl(igrid,ilay) = pr_neut*exp(-pr_slope/pr_neut*Ri(igrid,ilay)+Ri(igrid,ilay)/pr_neut) & |
---|
179 | + Ri(igrid,ilay) * pr_slope |
---|
180 | IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN |
---|
181 | call abort_physic("atke_compute_km_kh", & |
---|
182 | 'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1) |
---|
183 | ENDIF |
---|
184 | END IF |
---|
185 | |
---|
186 | Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay) |
---|
187 | |
---|
188 | ENDDO |
---|
189 | ENDDO |
---|
190 | |
---|
191 | |
---|
192 | ! Computing the mixing length: |
---|
193 | !============================================================== |
---|
194 | |
---|
195 | |
---|
196 | IF (iflag_atke_lmix .EQ. 1 ) THEN |
---|
197 | ! Blackadar formulation (~kappa l) + buoyancy length scale (Deardoff 1980) for very stable conditions |
---|
198 | DO ilay=2,nlay |
---|
199 | DO igrid=1,ngrid |
---|
200 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
201 | IF (N2(igrid,ilay) .GT. 0.) THEN |
---|
202 | lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)) |
---|
203 | lstrat=max(lstrat,lmin) |
---|
204 | !Inverse interpolation, Van de Wiel et al. 2010 |
---|
205 | l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) |
---|
206 | ENDIF |
---|
207 | ENDDO |
---|
208 | ENDDO |
---|
209 | |
---|
210 | ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN |
---|
211 | ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms |
---|
212 | ! implies 2 tuning coefficients clmix and clmixshear |
---|
213 | DO ilay=2,nlay |
---|
214 | DO igrid=1,ngrid |
---|
215 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
216 | IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN |
---|
217 | lstrat=min(clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)), & |
---|
218 | clmixshear*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))) |
---|
219 | lstrat=max(lstrat,lmin) |
---|
220 | !Inverse interpolation, Van de Wiel et al. 2010 |
---|
221 | l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) |
---|
222 | ENDIF |
---|
223 | ENDDO |
---|
224 | ENDDO |
---|
225 | |
---|
226 | ELSE IF (iflag_atke_lmix .EQ. 3 ) THEN |
---|
227 | ! add effect of wind shear on lstrat following grisogono 2010, qjrms |
---|
228 | ! keeping a single tuning coefficient clmix |
---|
229 | DO ilay=2,nlay |
---|
230 | DO igrid=1,ngrid |
---|
231 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
232 | IF (N2(igrid,ilay) .GT. 0. .AND. shear2(igrid,ilay) .GT. 0.) THEN |
---|
233 | lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(shear2(igrid,ilay))*(1.0+Ri(igrid,ilay)/(2.*Prandtl(igrid,ilay))) |
---|
234 | lstrat=max(lstrat,lmin) |
---|
235 | !Inverse interpolation, Van de Wiel et al. 2010 |
---|
236 | l_exchange(igrid,ilay)=(1./(l_exchange(igrid,ilay))+1./(lstrat))**(-1.0) |
---|
237 | ENDIF |
---|
238 | ENDDO |
---|
239 | ENDDO |
---|
240 | |
---|
241 | ELSE |
---|
242 | ! default Blackadar formulation: neglect effect of local stratification and shear |
---|
243 | DO ilay=2,nlay+1 |
---|
244 | DO igrid=1,ngrid |
---|
245 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
246 | ENDDO |
---|
247 | ENDDO |
---|
248 | ENDIF |
---|
249 | |
---|
250 | |
---|
251 | ! Computing the TKE k>=2: |
---|
252 | !======================== |
---|
253 | IF (iflag_atke == 0) THEN |
---|
254 | |
---|
255 | ! stationary solution (dtke/dt=0) |
---|
256 | |
---|
257 | DO ilay=2,nlay |
---|
258 | DO igrid=1,ngrid |
---|
259 | tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * & |
---|
260 | shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)) |
---|
261 | ENDDO |
---|
262 | ENDDO |
---|
263 | |
---|
264 | ELSE IF (iflag_atke == 1) THEN |
---|
265 | |
---|
266 | ! full implicit scheme resolved with a second order polynomial equation |
---|
267 | ! default solution which shows fair convergence properties |
---|
268 | DO ilay=2,nlay |
---|
269 | DO igrid=1,ngrid |
---|
270 | qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) |
---|
271 | delta=(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime)**2. & |
---|
272 | +4.*(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime*qq + & |
---|
273 | 2.*l_exchange(igrid,ilay)*l_exchange(igrid,ilay)*cepsilon*Sm(igrid,ilay) & |
---|
274 | *shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay))) |
---|
275 | qq=(-2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)/dtime + sqrt(delta))/2. |
---|
276 | qq=max(0.,qq) |
---|
277 | tke(igrid,ilay)=0.5*(qq**2) |
---|
278 | ENDDO |
---|
279 | ENDDO |
---|
280 | |
---|
281 | |
---|
282 | ELSE IF (iflag_atke == 2) THEN |
---|
283 | |
---|
284 | ! semi implicit scheme when l does not depend on tke |
---|
285 | ! positive-guaranteed if pr slope in stable condition >1 |
---|
286 | |
---|
287 | DO ilay=2,nlay |
---|
288 | DO igrid=1,ngrid |
---|
289 | qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) |
---|
290 | qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.) & |
---|
291 | *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) & |
---|
292 | /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) |
---|
293 | tke(igrid,ilay)=0.5*(qq**2) |
---|
294 | ENDDO |
---|
295 | ENDDO |
---|
296 | |
---|
297 | |
---|
298 | ELSE IF (iflag_atke == 3) THEN |
---|
299 | ! numerical resolution adapted from that in MAR (Deleersnijder 1992) |
---|
300 | ! positively defined by construction |
---|
301 | |
---|
302 | DO ilay=2,nlay |
---|
303 | DO igrid=1,ngrid |
---|
304 | qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) |
---|
305 | IF (Ri(igrid,ilay) .LT. 0.) THEN |
---|
306 | netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay)) |
---|
307 | netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) |
---|
308 | ELSE |
---|
309 | netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ & |
---|
310 | l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay) |
---|
311 | netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay) |
---|
312 | ENDIF |
---|
313 | qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss) |
---|
314 | tke(igrid,ilay)=0.5*(qq**2) |
---|
315 | ENDDO |
---|
316 | ENDDO |
---|
317 | |
---|
318 | ELSE IF (iflag_atke == 4) THEN |
---|
319 | ! semi implicit scheme from Arpege (V. Masson methodology with |
---|
320 | ! Taylor expansion of the dissipation term) |
---|
321 | DO ilay=2,nlay |
---|
322 | DO igrid=1,ngrid |
---|
323 | qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10) |
---|
324 | qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) & |
---|
325 | +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) & |
---|
326 | /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.))) |
---|
327 | qq=max(0.,qq) |
---|
328 | tke(igrid,ilay)=0.5*(qq**2) |
---|
329 | ENDDO |
---|
330 | ENDDO |
---|
331 | |
---|
332 | |
---|
333 | ELSE |
---|
334 | call abort_physic("atke_compute_km_kh", & |
---|
335 | 'numerical treatment of TKE not possible yet', 1) |
---|
336 | |
---|
337 | END IF |
---|
338 | |
---|
339 | ! We impose a 0 tke at nlay+1 |
---|
340 | !============================== |
---|
341 | |
---|
342 | DO igrid=1,ngrid |
---|
343 | tke(igrid,nlay+1)=0. |
---|
344 | END DO |
---|
345 | |
---|
346 | |
---|
347 | ! Calculation of surface TKE (k=1) |
---|
348 | !================================= |
---|
349 | ! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note) |
---|
350 | DO igrid=1,ngrid |
---|
351 | ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2)) |
---|
352 | tke(igrid,1)=ctkes*(ustar**2) |
---|
353 | END DO |
---|
354 | |
---|
355 | |
---|
356 | ! vertical diffusion of TKE |
---|
357 | !========================== |
---|
358 | IF (atke_ok_vdiff) THEN |
---|
359 | CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke) |
---|
360 | ENDIF |
---|
361 | |
---|
362 | |
---|
363 | ! Computing eddy diffusivity coefficients: |
---|
364 | !======================================== |
---|
365 | DO ilay=2,nlay ! TODO: also calculate for nlay+1 ? |
---|
366 | DO igrid=1,ngrid |
---|
367 | ! we add the molecular viscosity to Km,h |
---|
368 | Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5 |
---|
369 | Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5 |
---|
370 | END DO |
---|
371 | END DO |
---|
372 | |
---|
373 | ! for output: |
---|
374 | !=========== |
---|
375 | Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay) |
---|
376 | Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay) |
---|
377 | |
---|
378 | end subroutine atke_compute_km_kh |
---|
379 | |
---|
380 | !=============================================================================================== |
---|
381 | subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke) |
---|
382 | |
---|
383 | ! routine that computes the vertical diffusion of TKE by the turbulence |
---|
384 | ! using an implicit resolution (See note by Dufresne and Ghattas (2009)) |
---|
385 | ! E Vignon, July 2023 |
---|
386 | |
---|
387 | USE lmdz_atke_turbulence_ini, ONLY : rd, cke, viscom |
---|
388 | |
---|
389 | |
---|
390 | INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) |
---|
391 | INTEGER, INTENT(IN) :: nlay ! number of vertical index |
---|
392 | |
---|
393 | REAL, INTENT(IN) :: dtime ! physics time step (s) |
---|
394 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: z_lay ! altitude of mid-layers (m) |
---|
395 | REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: z_interf ! altitude of bottom interfaces (m) |
---|
396 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K) |
---|
397 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa) |
---|
398 | REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: l_exchange ! mixing length at interfaces between layers |
---|
399 | REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: Sm ! stability function for eddy diffusivity for momentum at interface between layers |
---|
400 | |
---|
401 | REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers |
---|
402 | |
---|
403 | |
---|
404 | |
---|
405 | INTEGER :: igrid,ilay |
---|
406 | REAL, DIMENSION(ngrid,nlay+1) :: Ke ! eddy diffusivity for TKE |
---|
407 | REAL, DIMENSION(ngrid,nlay+1) :: dtke |
---|
408 | REAL, DIMENSION(ngrid,nlay+1) :: ak, bk, ck, CCK, DDK |
---|
409 | REAL :: gammak,Kem,KKb,KKt |
---|
410 | |
---|
411 | |
---|
412 | ! Few initialisations |
---|
413 | CCK(:,:)=0. |
---|
414 | DDK(:,:)=0. |
---|
415 | dtke(:,:)=0. |
---|
416 | |
---|
417 | |
---|
418 | ! Eddy diffusivity for TKE |
---|
419 | |
---|
420 | DO ilay=2,nlay |
---|
421 | DO igrid=1,ngrid |
---|
422 | Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke |
---|
423 | ENDDO |
---|
424 | ENDDO |
---|
425 | ! at the top of the atmosphere set to 0 |
---|
426 | Ke(:,nlay+1)=0. |
---|
427 | ! at the surface, set it equal to that at the first model level |
---|
428 | Ke(:,1)=Ke(:,2) |
---|
429 | |
---|
430 | |
---|
431 | ! calculate intermediary variables |
---|
432 | |
---|
433 | DO ilay=2,nlay |
---|
434 | DO igrid=1,ngrid |
---|
435 | Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay)) |
---|
436 | KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay)) |
---|
437 | Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1)) |
---|
438 | KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1)) |
---|
439 | gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1)) |
---|
440 | ak(igrid,ilay)=-gammak*dtime*KKb |
---|
441 | ck(igrid,ilay)=-gammak*dtime*KKt |
---|
442 | bk(igrid,ilay)=1.+gammak*dtime*(KKt+KKb) |
---|
443 | ENDDO |
---|
444 | ENDDO |
---|
445 | |
---|
446 | ! calculate CCK and DDK coefficients |
---|
447 | ! downhill phase |
---|
448 | |
---|
449 | DO igrid=1,ngrid |
---|
450 | CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay) |
---|
451 | DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay) |
---|
452 | ENDDO |
---|
453 | |
---|
454 | |
---|
455 | DO ilay=nlay-1,2,-1 |
---|
456 | DO igrid=1,ngrid |
---|
457 | CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) & |
---|
458 | / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1)) |
---|
459 | DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1)) |
---|
460 | ENDDO |
---|
461 | ENDDO |
---|
462 | |
---|
463 | ! calculate TKE |
---|
464 | ! uphill phase |
---|
465 | |
---|
466 | DO ilay=2,nlay+1 |
---|
467 | DO igrid=1,ngrid |
---|
468 | dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay) |
---|
469 | ENDDO |
---|
470 | ENDDO |
---|
471 | |
---|
472 | ! update TKE |
---|
473 | tke(:,:)=tke(:,:)+dtke(:,:) |
---|
474 | |
---|
475 | |
---|
476 | end subroutine atke_vdiff_tke |
---|
477 | |
---|
478 | |
---|
479 | |
---|
480 | end module lmdz_atke_exchange_coeff |
---|