1 | module atke_exchange_coeff_mod |
---|
2 | |
---|
3 | implicit none |
---|
4 | |
---|
5 | contains |
---|
6 | |
---|
7 | subroutine atke_compute_km_kh(ngrid,nlay,dtime, & |
---|
8 | wind_u,wind_v,temp,play,pinterf, & |
---|
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) dry atmosphere |
---|
23 | ! (2) horizontal homogeneity (Dx=Dy=0.) |
---|
24 | !======================================================================= |
---|
25 | |
---|
26 | |
---|
27 | |
---|
28 | USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd |
---|
29 | USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, rg, rd |
---|
30 | USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin |
---|
31 | |
---|
32 | implicit none |
---|
33 | |
---|
34 | |
---|
35 | ! Declarations: |
---|
36 | !============= |
---|
37 | |
---|
38 | INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid) |
---|
39 | INTEGER, INTENT(IN) :: nlay ! number of vertical index |
---|
40 | |
---|
41 | REAL, INTENT(IN) :: dtime ! physics time step (s) |
---|
42 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_u ! zonal velocity (m/s) |
---|
43 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: wind_v ! meridional velocity (m/s) |
---|
44 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: temp ! temperature (K) |
---|
45 | REAL, DIMENSION(ngrid,nlay), INTENT(IN) :: play ! pressure (Pa) |
---|
46 | REAL, DIMENSION(ngrid,nlay+1), INTENT(IN) :: pinterf ! pressure at interfaces(Pa) |
---|
47 | |
---|
48 | |
---|
49 | REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT) :: tke ! turbulent kinetic energy at interface between layers |
---|
50 | |
---|
51 | REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Km_out ! output: Exchange coefficient for momentum at interface between layers |
---|
52 | REAL, DIMENSION(ngrid,nlay), INTENT(OUT) :: Kh_out ! output: Exchange coefficient for heat flux at interface between layers |
---|
53 | |
---|
54 | ! Local variables |
---|
55 | REAL, DIMENSION(ngrid,nlay+1) :: Km ! Exchange coefficient for momentum at interface between layers |
---|
56 | REAL, DIMENSION(ngrid,nlay+1) :: Kh ! Exchange coefficient for heat flux at interface between layers |
---|
57 | REAL, DIMENSION(ngrid,nlay) :: theta ! Potential temperature |
---|
58 | REAL, DIMENSION(ngrid,nlay+1) :: l_exchange ! Length of exchange (at interface) |
---|
59 | REAL, DIMENSION(ngrid,nlay+1) :: z_interf ! Altitude at the interface |
---|
60 | REAL, DIMENSION(ngrid,nlay) :: z_lay ! Altitude of layers |
---|
61 | REAL, DIMENSION(ngrid,nlay) :: dz_interf ! distance between two consecutive interfaces |
---|
62 | REAL, DIMENSION(ngrid,nlay) :: dz_lay ! distance between two layer middles (NB: first and last are half layers) |
---|
63 | REAL, DIMENSION(ngrid,nlay+1) :: N2 ! square of Brunt Vaisala pulsation (at interface) |
---|
64 | REAL, DIMENSION(ngrid,nlay+1) :: shear2 ! square of wind shear (at interface) |
---|
65 | REAL, DIMENSION(ngrid,nlay+1) :: Ri ! Richardson's number (at interface) |
---|
66 | REAL, DIMENSION(ngrid,nlay+1) :: Prandtl ! Turbulent Prandtl's number (at interface) |
---|
67 | REAL, DIMENSION(ngrid,nlay+1) :: Sm ! Stability function for momentum (at interface) |
---|
68 | REAL, DIMENSION(ngrid,nlay+1) :: Sh ! Stability function for heat (at interface) |
---|
69 | LOGICAL, DIMENSION(ngrid,nlay+1) :: switch_num ! switch of numerical integration in very stable cases |
---|
70 | |
---|
71 | INTEGER :: igrid,ilay ! horizontal,vertical index (flat grid) |
---|
72 | REAL :: cn,Ri0,Ri1 ! parameter for Sm stability function and Prandlt |
---|
73 | REAL :: preff ! reference pressure for potential temperature calculations |
---|
74 | REAL :: thetam ! mean potential temperature at interface |
---|
75 | REAL :: delta ! discriminant of the second order polynomial |
---|
76 | REAL :: qq ! tke=qq**2/2 |
---|
77 | REAL :: shear ! wind shear |
---|
78 | REAL :: lstrat ! mixing length depending on local stratification |
---|
79 | REAL :: taustrat ! caracteristic timescale for turbulence in very stable conditions |
---|
80 | |
---|
81 | ! Initializations: |
---|
82 | !================ |
---|
83 | |
---|
84 | DO igrid=1,ngrid |
---|
85 | dz_interf(igrid,1) = 0.0 |
---|
86 | z_interf(igrid,1) = 0.0 |
---|
87 | END DO |
---|
88 | |
---|
89 | ! Calculation of potential temperature: (if vapor -> todo virtual potential temperature) |
---|
90 | !===================================== |
---|
91 | |
---|
92 | preff=100000. |
---|
93 | ! The result should not depend on the choice of preff |
---|
94 | DO ilay=1,nlay |
---|
95 | DO igrid = 1, ngrid |
---|
96 | theta(igrid,ilay)=temp(igrid,ilay)*(preff/play(igrid,ilay))**(rd/rcpd) |
---|
97 | END DO |
---|
98 | END DO |
---|
99 | |
---|
100 | |
---|
101 | |
---|
102 | ! Calculation of altitude of layers' middle and bottom interfaces: |
---|
103 | !================================================================= |
---|
104 | |
---|
105 | DO ilay=2,nlay+1 |
---|
106 | DO igrid=1,ngrid |
---|
107 | dz_interf(igrid,ilay-1) = rd*temp(igrid,ilay-1)/rg/play(igrid,ilay-1)*(pinterf(igrid,ilay-1)-pinterf(igrid,ilay)) |
---|
108 | z_interf(igrid,ilay) = z_interf(igrid,ilay-1) + dz_interf(igrid,ilay-1) |
---|
109 | ENDDO |
---|
110 | ENDDO |
---|
111 | |
---|
112 | DO ilay=1,nlay |
---|
113 | DO igrid=1,ngrid |
---|
114 | z_lay(igrid,ilay)=0.5*(z_interf(igrid, ilay+1) + z_interf(igrid, ilay)) |
---|
115 | ENDDO |
---|
116 | ENDDO |
---|
117 | |
---|
118 | |
---|
119 | ! Computes the gradient Richardson's number and stability functions: |
---|
120 | !=================================================================== |
---|
121 | |
---|
122 | ! calculation of cn = Sm value at Ri=0 |
---|
123 | ! direct dependance on cepsilon to guarantee Fm=1 (first-order like stability function) at Ri=0 |
---|
124 | cn=(1./sqrt(cepsilon))**(2/3) |
---|
125 | ! calculation of Ri0 such that continuity in slope of Sm at Ri=0 |
---|
126 | Ri0=2./rpi*(cinf - cn)*ric/cn |
---|
127 | ! calculation of Ri1 to guarantee continuity in slope of Prandlt number at Ri=0 |
---|
128 | Ri1 = -2./rpi * (pr_asym - pr_neut) / pr_slope |
---|
129 | |
---|
130 | |
---|
131 | DO ilay=2,nlay |
---|
132 | DO igrid=1,ngrid |
---|
133 | dz_lay(igrid,ilay)=z_lay(igrid,ilay)-z_lay(igrid,ilay-1) |
---|
134 | thetam=0.5*(theta(igrid,ilay) + theta(igrid,ilay-1)) |
---|
135 | N2(igrid,ilay) = rg * (theta(igrid,ilay) - theta(igrid,ilay-1))/thetam / dz_lay(igrid,ilay) |
---|
136 | shear2(igrid,ilay)= (((wind_u(igrid,ilay) - wind_u(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 + & |
---|
137 | ((wind_v(igrid,ilay) - wind_v(igrid,ilay-1)) / dz_lay(igrid,ilay))**2 ) |
---|
138 | Ri(igrid,ilay) = N2(igrid,ilay) / MAX(shear2(igrid,ilay),1E-10) |
---|
139 | |
---|
140 | IF (Ri(igrid,ilay) < 0.) THEN ! unstable cases |
---|
141 | Sm(igrid,ilay) = 2./rpi * (cinf-cn) * atan(-Ri(igrid,ilay)/Ri0) + cn |
---|
142 | Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut |
---|
143 | ELSE ! stable cases |
---|
144 | Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric)) |
---|
145 | Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope |
---|
146 | IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN |
---|
147 | call abort_physic("atke_compute_km_kh", & |
---|
148 | 'Ri>=Pr in stable conditions -> violates energy conservation principles, change pr_neut or slope', 1) |
---|
149 | ENDIF |
---|
150 | END IF |
---|
151 | |
---|
152 | Sh(igrid,ilay) = Sm(igrid,ilay) / Prandtl(igrid,ilay) |
---|
153 | |
---|
154 | ENDDO |
---|
155 | ENDDO |
---|
156 | |
---|
157 | |
---|
158 | ! Computing the mixing length: |
---|
159 | !============================================================== |
---|
160 | |
---|
161 | switch_num(:,:)=.false. |
---|
162 | |
---|
163 | IF (iflag_atke_lmix .EQ. 1 ) THEN |
---|
164 | |
---|
165 | DO ilay=2,nlay |
---|
166 | DO igrid=1,ngrid |
---|
167 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
168 | IF (N2(igrid,ilay) .GT. 0.) THEN |
---|
169 | lstrat=clmix*sqrt(tke(igrid,ilay))/sqrt(N2(igrid,ilay)) |
---|
170 | IF (lstrat .LT. l_exchange(igrid,ilay)) THEN |
---|
171 | l_exchange(igrid,ilay)=max(lstrat,lmin) |
---|
172 | switch_num(igrid,ilay)=.true. |
---|
173 | ENDIF |
---|
174 | ENDIF |
---|
175 | ENDDO |
---|
176 | ENDDO |
---|
177 | |
---|
178 | ELSE IF (iflag_atke_lmix .EQ. 2 ) THEN |
---|
179 | ! add effect of wind shear on lstrat following grisogono and belusic 2008, qjrms, eq 2 |
---|
180 | DO ilay=2,nlay |
---|
181 | DO igrid=1,ngrid |
---|
182 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
183 | IF (N2(igrid,ilay) .GT. 0.) THEN |
---|
184 | lstrat=clmix*sqrt(tke(igrid,ilay))/(2.*max(sqrt(shear2(igrid,ilay)),1E-10)*(1.+sqrt(Ri(igrid,ilay))/2.)) |
---|
185 | IF (lstrat .LT. l_exchange(igrid,ilay)) THEN |
---|
186 | l_exchange(igrid,ilay)=max(lstrat,lmin) |
---|
187 | ENDIF |
---|
188 | ENDIF |
---|
189 | ENDDO |
---|
190 | ENDDO |
---|
191 | |
---|
192 | |
---|
193 | |
---|
194 | ELSE |
---|
195 | ! default: neglect effect of local stratification and shear |
---|
196 | |
---|
197 | DO ilay=2,nlay+1 |
---|
198 | DO igrid=1,ngrid |
---|
199 | l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0) |
---|
200 | ENDDO |
---|
201 | |
---|
202 | ENDDO |
---|
203 | ENDIF |
---|
204 | |
---|
205 | |
---|
206 | ! Computing the TKE: |
---|
207 | !=================== |
---|
208 | IF (iflag_atke == 0) THEN |
---|
209 | |
---|
210 | DO ilay=2,nlay |
---|
211 | DO igrid=1,ngrid |
---|
212 | tke(igrid,ilay) = cepsilon * l_exchange(igrid,ilay)**2 * Sm(igrid,ilay) * & |
---|
213 | shear2(igrid,ilay) * (1. - Ri(igrid,ilay) / Prandtl(igrid,ilay)) |
---|
214 | ENDDO |
---|
215 | ENDDO |
---|
216 | |
---|
217 | ELSE IF (iflag_atke == 1) THEN |
---|
218 | |
---|
219 | ! full implicit scheme resolved with a second order polynomial equation |
---|
220 | |
---|
221 | DO ilay=2,nlay |
---|
222 | DO igrid=1,ngrid |
---|
223 | qq=sqrt(2.*tke(igrid,ilay)) |
---|
224 | delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * & |
---|
225 | (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*shear2(igrid,ilay) & |
---|
226 | *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) |
---|
227 | qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay) |
---|
228 | tke(igrid,ilay)=0.5*(qq**2) |
---|
229 | ENDDO |
---|
230 | ENDDO |
---|
231 | |
---|
232 | ELSE IF (iflag_atke == 2) THEN |
---|
233 | |
---|
234 | ! full implicit scheme resolved with a second order polynomial equation |
---|
235 | ! + exact resolution for very stable cases (iflag_atke_lmix=1) |
---|
236 | DO ilay=2,nlay |
---|
237 | DO igrid=1,ngrid |
---|
238 | qq=sqrt(2.*tke(igrid,ilay)) |
---|
239 | IF (switch_num(igrid,ilay)) THEN |
---|
240 | taustrat=clmix/2./sqrt(N2(igrid,ilay))*Sm(igrid,ilay)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) & |
---|
241 | - sqrt(N2(igrid,ilay))/2./cepsilon/clmix |
---|
242 | taustrat=-1./taustrat |
---|
243 | qq=qq*exp(-dtime/taustrat) |
---|
244 | ELSE |
---|
245 | delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * & |
---|
246 | (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay) * shear2(igrid,ilay) & |
---|
247 | *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) |
---|
248 | qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay) |
---|
249 | ENDIF |
---|
250 | tke(igrid,ilay)=0.5*(qq**2) |
---|
251 | ENDDO |
---|
252 | ENDDO |
---|
253 | |
---|
254 | |
---|
255 | |
---|
256 | ELSE IF (iflag_atke == 3) THEN |
---|
257 | |
---|
258 | ! semi implicit scheme when l does not depend on tke |
---|
259 | ! positive-guaranteed if pr slope in stable condition >1 |
---|
260 | DO ilay=2,nlay |
---|
261 | DO igrid=1,ngrid |
---|
262 | qq=sqrt(2.*tke(igrid,ilay)) |
---|
263 | qq=(qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) & |
---|
264 | /(1.+dtime*qq/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2))) |
---|
265 | tke(igrid,ilay)=0.5*(qq**2) |
---|
266 | ENDDO |
---|
267 | ENDDO |
---|
268 | |
---|
269 | |
---|
270 | ELSE |
---|
271 | call abort_physic("atke_compute_km_kh", & |
---|
272 | 'numerical treatment of TKE not possible yet', 1) |
---|
273 | |
---|
274 | END IF |
---|
275 | |
---|
276 | |
---|
277 | ! Computing eddy diffusivity coefficients: |
---|
278 | !======================================== |
---|
279 | DO ilay=2,nlay ! TODO: also calculate for nlay+1 ? |
---|
280 | DO igrid=1,ngrid |
---|
281 | ! we add the molecular viscosity to Km,h |
---|
282 | Km(igrid,ilay) = viscom + l_exchange(igrid,ilay) * Sm(igrid,ilay) * tke(igrid,ilay)**0.5 |
---|
283 | Kh(igrid,ilay) = viscoh + l_exchange(igrid,ilay) * Sh(igrid,ilay) * tke(igrid,ilay)**0.5 |
---|
284 | END DO |
---|
285 | END DO |
---|
286 | |
---|
287 | ! for output: |
---|
288 | !=========== |
---|
289 | Km_out(1:ngrid,2:nlay)=Km(1:ngrid,2:nlay) |
---|
290 | Kh_out(1:ngrid,2:nlay)=Kh(1:ngrid,2:nlay) |
---|
291 | |
---|
292 | end subroutine atke_compute_km_kh |
---|
293 | |
---|
294 | |
---|
295 | end module atke_exchange_coeff_mod |
---|