Directory: | ./ |
---|---|
File: | phys/cdrag_mod.f90 |
Date: | 2022-01-11 19:19:34 |
Exec | Total | Coverage | |
---|---|---|---|
Lines: | 106 | 179 | 59.2% |
Branches: | 44 | 82 | 53.7% |
Line | Branch | Exec | Source |
---|---|---|---|
1 | ! | ||
2 | MODULE cdrag_mod | ||
3 | ! | ||
4 | ! This module contains some procedures for calculation of the cdrag | ||
5 | ! coefficients for turbulent diffusion at surface | ||
6 | ! | ||
7 | IMPLICIT NONE | ||
8 | |||
9 | CONTAINS | ||
10 | ! | ||
11 | !**************************************************************************************** | ||
12 | ! | ||
13 | !r original routine svn3623 | ||
14 | ! | ||
15 | 28939168 | SUBROUTINE cdrag(knon, nsrf, & | |
16 | 4320 | speed, t1, q1, zgeop1, & | |
17 | psol, tsurf, qsurf, z0m, z0h, & | ||
18 | cdm, cdh, zri, pref) | ||
19 | |||
20 | USE dimphy | ||
21 | USE indice_sol_mod | ||
22 | USE print_control_mod, ONLY: lunout, prt_level | ||
23 | USE ioipsl_getin_p_mod, ONLY : getin_p | ||
24 | |||
25 | IMPLICIT NONE | ||
26 | ! ================================================================= c | ||
27 | ! | ||
28 | ! Objet : calcul des cdrags pour le moment (pcfm) et | ||
29 | ! les flux de chaleur sensible et latente (pcfh) d'apr??s | ||
30 | ! Louis 1982, Louis 1979, King et al 2001 | ||
31 | ! ou Zilitinkevich et al 2002 pour les cas stables, Louis 1979 | ||
32 | ! et 1982 pour les cas instables | ||
33 | ! | ||
34 | ! Modified history: | ||
35 | ! writting on the 20/05/2016 | ||
36 | ! modified on the 13/12/2016 to be adapted to LMDZ6 | ||
37 | ! | ||
38 | ! References: | ||
39 | ! Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the | ||
40 | ! atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202. | ||
41 | ! Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the | ||
42 | ! operational PBL parametrization at ECMWF'. Workshop on boundary layer | ||
43 | ! parametrization, November 1981, ECMWF, Reading, England. | ||
44 | ! Page: 19. Equations in Table 1. | ||
45 | ! Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS | ||
46 | ! IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM | ||
47 | ! Boundary-Layer Meteorology 72: 331-344 | ||
48 | ! Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. | ||
49 | ! European Centre for Medium-Range Weather Forecasts. | ||
50 | ! Equations: 110-113. Page 40. | ||
51 | ! Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF | ||
52 | ! model to the parameterization of evaporation from the tropical oceans. J. | ||
53 | ! Climate, 5:418-434. | ||
54 | ! King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate | ||
55 | ! to surface and boundary-layer flux parametrizations | ||
56 | ! QJRMS, 127, pp 779-794 | ||
57 | ! | ||
58 | ! ================================================================= c | ||
59 | ! ================================================================= c | ||
60 | ! On choisit le couple de fonctions de correction avec deux flags: | ||
61 | ! Un pour les cas instables, un autre pour les cas stables | ||
62 | ! | ||
63 | ! iflag_corr_insta: | ||
64 | ! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h) | ||
65 | ! 2: Louis 1982 | ||
66 | ! 3: Laurent Li | ||
67 | ! | ||
68 | ! iflag_corr_sta: | ||
69 | ! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h) | ||
70 | ! 2: Louis 1982 | ||
71 | ! 3: Laurent Li | ||
72 | ! 4: King 2001 (SHARP) | ||
73 | ! 5: MO 1st order theory (allow collapse of turbulence) | ||
74 | ! | ||
75 | ! | ||
76 | !***************************************************************** | ||
77 | ! Parametres d'entree | ||
78 | !***************************************************************** | ||
79 | |||
80 | INTEGER, INTENT(IN) :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface | ||
81 | REAL, DIMENSION(klon), INTENT(IN) :: speed ! module du vent au 1er niveau du modele | ||
82 | REAL, DIMENSION(klon), INTENT(IN) :: zgeop1! geopotentiel au 1er niveau du modele | ||
83 | REAL, DIMENSION(klon), INTENT(IN) :: tsurf ! Surface temperature (K) | ||
84 | REAL, DIMENSION(klon), INTENT(IN) :: qsurf ! Surface humidity (Kg/Kg) | ||
85 | REAL, DIMENSION(klon), INTENT(IN) :: z0m, z0h ! Rugosity at surface (m) | ||
86 | REAL, DIMENSION(klon), INTENT(IN) :: t1 ! Temperature au premier niveau (K) | ||
87 | REAL, DIMENSION(klon), INTENT(IN) :: q1 ! humidite specifique au premier niveau (kg/kg) | ||
88 | REAL, DIMENSION(klon), INTENT(IN) :: psol ! pression au sol | ||
89 | |||
90 | |||
91 | |||
92 | ! Parametres de sortie | ||
93 | !****************************************************************** | ||
94 | REAL, DIMENSION(klon), INTENT(OUT) :: cdm ! Drag coefficient for momentum | ||
95 | REAL, DIMENSION(klon), INTENT(OUT) :: cdh ! Drag coefficient for heat flux | ||
96 | REAL, DIMENSION(klon), INTENT(OUT) :: zri ! Richardson number | ||
97 | REAL, DIMENSION(klon), INTENT(OUT) :: pref ! Pression au niveau zgeop/RG | ||
98 | |||
99 | ! Variables Locales | ||
100 | !****************************************************************** | ||
101 | |||
102 | |||
103 | INCLUDE "YOMCST.h" | ||
104 | INCLUDE "YOETHF.h" | ||
105 | INCLUDE "clesphys.h" | ||
106 | |||
107 | |||
108 | REAL, PARAMETER :: CKAP=0.40, CKAPT=0.42 | ||
109 | REAL CEPDU2 | ||
110 | REAL ALPHA | ||
111 | REAL CB,CC,CD,C2,C3 | ||
112 | REAL MU, CM, CH, B, CMstar, CHstar | ||
113 | REAL PM, PH, BPRIME | ||
114 | REAL C | ||
115 | INTEGER ng_q1 ! Number of grids that q1 < 0.0 | ||
116 | INTEGER ng_qsurf ! Number of grids that qsurf < 0.0 | ||
117 | INTEGER i | ||
118 | REAL zdu2, ztsolv | ||
119 | REAL ztvd, zscf | ||
120 | REAL zucf, zcr | ||
121 | REAL friv, frih | ||
122 | 8640 | REAL, DIMENSION(klon) :: FM, FH ! stability functions | |
123 | 8640 | REAL, DIMENSION(klon) :: cdmn, cdhn ! Drag coefficient in neutral conditions | |
124 | REAL zzzcd | ||
125 | |||
126 | LOGICAL, SAVE :: firstcall = .TRUE. | ||
127 | !$OMP THREADPRIVATE(firstcall) | ||
128 | INTEGER, SAVE :: iflag_corr_sta | ||
129 | !$OMP THREADPRIVATE(iflag_corr_sta) | ||
130 | INTEGER, SAVE :: iflag_corr_insta | ||
131 | !$OMP THREADPRIVATE(iflag_corr_insta) | ||
132 | |||
133 | !===================================================================c | ||
134 | ! Valeurs numeriques des constantes | ||
135 | !===================================================================c | ||
136 | |||
137 | |||
138 | ! Minimum du carre du vent | ||
139 | |||
140 | CEPDU2 = (0.1)**2 | ||
141 | |||
142 | ! Louis 1982 | ||
143 | |||
144 | CB=5.0 | ||
145 | CC=5.0 | ||
146 | CD=5.0 | ||
147 | |||
148 | |||
149 | ! King 2001 | ||
150 | |||
151 | C2=0.25 | ||
152 | C3=0.0625 | ||
153 | |||
154 | |||
155 | ! Louis 1979 | ||
156 | |||
157 | BPRIME=4.7 | ||
158 | B=9.4 | ||
159 | |||
160 | |||
161 | !MO | ||
162 | |||
163 | ALPHA=5.0 | ||
164 | |||
165 | |||
166 | ! ================================================================= c | ||
167 | ! Tests avant de commencer | ||
168 | ! Fuxing WANG, 04/03/2015 | ||
169 | ! To check if there are negative q1, qsurf values. | ||
170 | !====================================================================c | ||
171 | 4320 | ng_q1 = 0 ! Initialization | |
172 | 4320 | ng_qsurf = 0 ! Initialization | |
173 |
2/2✓ Branch 0 taken 1824584 times.
✓ Branch 1 taken 4320 times.
|
1828904 | DO i = 1, knon |
174 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1824584 times.
|
1824584 | IF (q1(i).LT.0.0) ng_q1 = ng_q1 + 1 |
175 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1824584 times.
|
1828904 | IF (qsurf(i).LT.0.0) ng_qsurf = ng_qsurf + 1 |
176 | ENDDO | ||
177 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 4320 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
4320 | IF (ng_q1.GT.0 .and. prt_level > 5) THEN |
178 | ✗ | WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !" | |
179 | ✗ | WRITE(lunout,*)" The total number of the grids is: ", ng_q1 | |
180 | ✗ | WRITE(lunout,*)" The negative q1 is set to zero " | |
181 | ! abort_message="voir ci-dessus" | ||
182 | ! CALL abort_physic(modname,abort_message,1) | ||
183 | ENDIF | ||
184 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 4320 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
4320 | IF (ng_qsurf.GT.0 .and. prt_level > 5) THEN |
185 | ✗ | WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !" | |
186 | ✗ | WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf | |
187 | ✗ | WRITE(lunout,*)" The negative qsurf is set to zero " | |
188 | ! abort_message="voir ci-dessus" | ||
189 | ! CALL abort_physic(modname,abort_message,1) | ||
190 | ENDIF | ||
191 | |||
192 | |||
193 | |||
194 | !=============================================================================c | ||
195 | ! Calcul du cdrag | ||
196 | !=============================================================================c | ||
197 | |||
198 | ! On choisit les fonctions de stabilite utilisees au premier appel | ||
199 | !************************************************************************** | ||
200 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 4319 times.
|
4320 | IF (firstcall) THEN |
201 | 1 | iflag_corr_sta=2 | |
202 | 1 | iflag_corr_insta=2 | |
203 | |||
204 | 1 | CALL getin_p('iflag_corr_sta',iflag_corr_sta) | |
205 | 1 | CALL getin_p('iflag_corr_insta',iflag_corr_insta) | |
206 | |||
207 | 1 | firstcall = .FALSE. | |
208 | ENDIF | ||
209 | |||
210 | !xxxxxxxxxxxxxxxxxxxxxxx | ||
211 |
2/2✓ Branch 0 taken 1824584 times.
✓ Branch 1 taken 4320 times.
|
1828904 | DO i = 1, knon ! Boucle sur l'horizontal |
212 | !xxxxxxxxxxxxxxxxxxxxxxx | ||
213 | |||
214 | |||
215 | ! calculs preliminaires: | ||
216 | !*********************** | ||
217 | |||
218 | |||
219 | 1824584 | zdu2 = MAX(CEPDU2, speed(i)**2) | |
220 | pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* & | ||
221 | 1824584 | (1.+ RETV * max(q1(i),0.0)))) ! negative q1 set to zero | |
222 | 1824584 | ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0)) ! negative qsurf set to zero | |
223 | ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*max(q1(i),0.0))) & | ||
224 | 1824584 | *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero | |
225 | 1824584 | zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd) | |
226 | |||
227 | |||
228 | ! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)): | ||
229 | !******************************************************************** | ||
230 | |||
231 | 1824584 | zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i))) | |
232 | 1824584 | cdmn(i) = zzzcd*zzzcd | |
233 | 1824584 | cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i)))) | |
234 | |||
235 | |||
236 | ! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi : | ||
237 | !******************************************************* | ||
238 | |||
239 | !'''''''''''''' | ||
240 | ! Cas instables | ||
241 | !'''''''''''''' | ||
242 | |||
243 | 4320 | IF (zri(i) .LT. 0.) THEN | |
244 | |||
245 | |||
246 | ✗ | SELECT CASE (iflag_corr_insta) | |
247 | |||
248 | CASE (1) ! Louis 1979 + Mascart 1995 | ||
249 | |||
250 | ✗ | MU=LOG(MAX(z0m(i)/z0h(i),0.01)) | |
251 | ✗ | CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3) | |
252 | ✗ | PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3) | |
253 | ✗ | CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3) | |
254 | ✗ | PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3) | |
255 | CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & | ||
256 | & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i))) & | ||
257 | ✗ | & * ((zgeop1(i)/(RG*z0h(i)))**PH) | |
258 | CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & | ||
259 | & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & | ||
260 | ✗ | & * ((zgeop1(i)/(RG*z0m(i)))**PM) | |
261 | |||
262 | |||
263 | |||
264 | |||
265 | ✗ | FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i)))) | |
266 | ✗ | FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i)))) | |
267 | |||
268 | CASE (2) ! Louis 1982 | ||
269 | |||
270 | zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & | ||
271 | 1198238 | *(1.0+zgeop1(i)/(RG*z0m(i))))) | |
272 | 1198238 | FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) | |
273 | 1198238 | FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) | |
274 | |||
275 | |||
276 | CASE (3) ! Laurent Li | ||
277 | |||
278 | |||
279 | ✗ | FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) | |
280 | ✗ | FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) | |
281 | |||
282 | |||
283 | |||
284 | CASE default ! Louis 1982 | ||
285 | |||
286 | zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & | ||
287 | ✗ | *(1.0+zgeop1(i)/(RG*z0m(i))))) | |
288 | ✗ | FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) | |
289 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 1198238 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
1198238 | FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) |
290 | |||
291 | |||
292 | END SELECT | ||
293 | |||
294 | |||
295 | |||
296 | ! Calcul des drags | ||
297 | |||
298 | |||
299 | 1198238 | cdm(i)=cdmn(i)*FM(i) | |
300 | 1198238 | cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) | |
301 | |||
302 | |||
303 | ! Traitement particulier des cas oceaniques | ||
304 | ! on applique Miller et al 1992 en l'absence de gustiness | ||
305 | |||
306 | 1198238 | IF (nsrf == is_oce) THEN | |
307 | ! cdh(i)=f_cdrag_oce*cdhn(i)*FH(i) | ||
308 | |||
309 | 629952 | IF(iflag_gusts==0) THEN | |
310 | ✗ | zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.) | |
311 | ✗ | cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25) | |
312 | ENDIF | ||
313 | |||
314 | |||
315 | 629952 | cdm(i)=MIN(cdm(i),cdmmax) | |
316 | 629952 | cdh(i)=MIN(cdh(i),cdhmax) | |
317 | |||
318 | END IF | ||
319 | |||
320 | |||
321 | |||
322 | ELSE | ||
323 | |||
324 | !''''''''''''''' | ||
325 | ! Cas stables : | ||
326 | !''''''''''''''' | ||
327 | 626346 | zri(i) = MIN(20.,zri(i)) | |
328 | |||
329 | ✗ | SELECT CASE (iflag_corr_sta) | |
330 | |||
331 | CASE (1) ! Louis 1979 + Mascart 1995 | ||
332 | |||
333 | ✗ | FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min) | |
334 | ✗ | FH(i)=FM(i) | |
335 | |||
336 | |||
337 | CASE (2) ! Louis 1982 | ||
338 | |||
339 | ✗ | zscf = SQRT(1.+CD*ABS(zri(i))) | |
340 | ✗ | FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) | |
341 | ✗ | FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) | |
342 | |||
343 | |||
344 | CASE (3) ! Laurent Li | ||
345 | |||
346 | ✗ | FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min) | |
347 | 626346 | FH(i)=FM(i) | |
348 | |||
349 | |||
350 | CASE (4) ! King 2001 | ||
351 | |||
352 |
2/2✓ Branch 0 taken 318325 times.
✓ Branch 1 taken 308021 times.
|
626346 | if (zri(i) .LT. C2/2.) then |
353 | 318325 | FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min) | |
354 | 318325 | FH(i)= FM(i) | |
355 | |||
356 | |||
357 | else | ||
358 | 308021 | FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min) | |
359 | 308021 | FH(i)= FM(i) | |
360 | endif | ||
361 | |||
362 | |||
363 | CASE (5) ! MO | ||
364 | |||
365 | ✗ | if (zri(i) .LT. 1./alpha) then | |
366 | |||
367 | ✗ | FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min) | |
368 | ✗ | FH(i)=FM(i) | |
369 | |||
370 | else | ||
371 | |||
372 | |||
373 | ✗ | FM(i)=MAX(1E-7,f_ri_cd_min) | |
374 | ✗ | FH(i)=FM(i) | |
375 | |||
376 | endif | ||
377 | |||
378 | |||
379 | |||
380 | |||
381 | |||
382 | CASE default ! Louis 1982 | ||
383 | |||
384 | ✗ | zscf = SQRT(1.+CD*ABS(zri(i))) | |
385 | ✗ | FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) | |
386 |
1/6✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 626346 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
626346 | FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) |
387 | |||
388 | |||
389 | |||
390 | END SELECT | ||
391 | |||
392 | ! Calcul des drags | ||
393 | |||
394 | |||
395 | 626346 | cdm(i)=cdmn(i)*FM(i) | |
396 | 626346 | cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) | |
397 | |||
398 | 626346 | IF(nsrf.EQ.is_oce) THEN | |
399 | |||
400 | 96002 | cdh(i)=f_cdrag_oce*cdhn(i)*FH(i) | |
401 | 96002 | cdm(i)=MIN(cdm(i),cdmmax) | |
402 | 96002 | cdh(i)=MIN(cdh(i),cdhmax) | |
403 | |||
404 | ENDIF | ||
405 | |||
406 | |||
407 | ENDIF | ||
408 | |||
409 | !xxxxxxxxxxx | ||
410 | END DO ! Fin de la boucle sur l'horizontal | ||
411 | !xxxxxxxxxxx | ||
412 | ! ================================================================= c | ||
413 | |||
414 | 4320 | END SUBROUTINE cdrag | |
415 | |||
416 | ! | ||
417 | 14400 | SUBROUTINE cdragn_ri1(knon, nsrf, & | |
418 | 14400 | speed, t1, q1, zgeop1, & | |
419 | psol, tsurf, qsurf, z0m, z0h, & | ||
420 | ri1, iri1, & | ||
421 | cdm, cdh, zri, pref) | ||
422 | |||
423 |
7/8✓ Branch 0 taken 1198238 times.
✓ Branch 1 taken 626346 times.
✓ Branch 2 taken 629952 times.
✓ Branch 3 taken 568286 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 629952 times.
✓ Branch 6 taken 96002 times.
✓ Branch 7 taken 530344 times.
|
3649168 | USE dimphy |
424 | USE indice_sol_mod | ||
425 | USE print_control_mod, ONLY: lunout, prt_level | ||
426 | USE ioipsl_getin_p_mod, ONLY : getin_p | ||
427 | |||
428 | IMPLICIT NONE | ||
429 | ! ================================================================= c | ||
430 | ! | ||
431 | ! Objet : calcul des cdrags pour le moment (pcfm) et | ||
432 | ! les flux de chaleur sensible et latente (pcfh) d'apr??s | ||
433 | ! Louis 1982, Louis 1979, King et al 2001 | ||
434 | ! ou Zilitinkevich et al 2002 pour les cas stables, Louis 1979 | ||
435 | ! et 1982 pour les cas instables | ||
436 | ! | ||
437 | ! Modified history: | ||
438 | ! writting on the 20/05/2016 | ||
439 | ! modified on the 13/12/2016 to be adapted to LMDZ6 | ||
440 | ! | ||
441 | ! References: | ||
442 | ! Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the | ||
443 | ! atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202. | ||
444 | ! Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the | ||
445 | ! operational PBL parametrization at ECMWF'. Workshop on boundary layer | ||
446 | ! parametrization, November 1981, ECMWF, Reading, England. | ||
447 | ! Page: 19. Equations in Table 1. | ||
448 | ! Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS | ||
449 | ! IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM | ||
450 | ! Boundary-Layer Meteorology 72: 331-344 | ||
451 | ! Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. | ||
452 | ! European Centre for Medium-Range Weather Forecasts. | ||
453 | ! Equations: 110-113. Page 40. | ||
454 | ! Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF | ||
455 | ! model to the parameterization of evaporation from the tropical oceans. J. | ||
456 | ! Climate, 5:418-434. | ||
457 | ! King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate | ||
458 | ! to surface and boundary-layer flux parametrizations | ||
459 | ! QJRMS, 127, pp 779-794 | ||
460 | ! | ||
461 | ! ================================================================= c | ||
462 | ! ================================================================= c | ||
463 | ! On choisit le couple de fonctions de correction avec deux flags: | ||
464 | ! Un pour les cas instables, un autre pour les cas stables | ||
465 | ! | ||
466 | ! iflag_corr_insta: | ||
467 | ! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h) | ||
468 | ! 2: Louis 1982 | ||
469 | ! 3: Laurent Li | ||
470 | ! | ||
471 | ! iflag_corr_sta: | ||
472 | ! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h) | ||
473 | ! 2: Louis 1982 | ||
474 | ! 3: Laurent Li | ||
475 | ! 4: King 2001 (SHARP) | ||
476 | ! 5: MO 1st order theory (allow collapse of turbulence) | ||
477 | ! | ||
478 | ! | ||
479 | !***************************************************************** | ||
480 | ! Parametres d'entree | ||
481 | !***************************************************************** | ||
482 | |||
483 | INTEGER, INTENT(IN) :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface | ||
484 | REAL, DIMENSION(klon), INTENT(IN) :: speed ! module du vent au 1er niveau du modele | ||
485 | REAL, DIMENSION(klon), INTENT(IN) :: zgeop1! geopotentiel au 1er niveau du modele | ||
486 | REAL, DIMENSION(klon), INTENT(IN) :: tsurf ! Surface temperature (K) | ||
487 | REAL, DIMENSION(klon), INTENT(IN) :: qsurf ! Surface humidity (Kg/Kg) | ||
488 | REAL, DIMENSION(klon), INTENT(IN) :: z0m, z0h ! Rugosity at surface (m) | ||
489 | REAL, DIMENSION(klon), INTENT(IN) :: ri1 ! Richardson 1ere couche | ||
490 | INTEGER, INTENT(IN) :: iri1 ! | ||
491 | REAL, DIMENSION(klon), INTENT(IN) :: t1 ! Temperature au premier niveau (K) | ||
492 | REAL, DIMENSION(klon), INTENT(IN) :: q1 ! humidite specifique au premier niveau (kg/kg) | ||
493 | REAL, DIMENSION(klon), INTENT(IN) :: psol ! pression au sol | ||
494 | |||
495 | |||
496 | |||
497 | ! Parametres de sortie | ||
498 | !****************************************************************** | ||
499 | REAL, DIMENSION(klon), INTENT(OUT) :: cdm ! Drag coefficient for heat flux | ||
500 | REAL, DIMENSION(klon), INTENT(OUT) :: cdh ! Drag coefficient for momentum | ||
501 | REAL, DIMENSION(klon), INTENT(OUT) :: zri ! Richardson number | ||
502 | REAL, DIMENSION(klon), INTENT(OUT) :: pref ! Pression au niveau zgeop/RG | ||
503 | ! Variables Locales | ||
504 | !****************************************************************** | ||
505 | |||
506 | |||
507 | INCLUDE "YOMCST.h" | ||
508 | INCLUDE "YOETHF.h" | ||
509 | INCLUDE "clesphys.h" | ||
510 | |||
511 | |||
512 | REAL, PARAMETER :: CKAP=0.40, CKAPT=0.42 | ||
513 | REAL CEPDU2 | ||
514 | REAL ALPHA | ||
515 | REAL CB,CC,CD,C2,C3 | ||
516 | REAL MU, CM, CH, B, CMstar, CHstar | ||
517 | REAL PM, PH, BPRIME | ||
518 | REAL C | ||
519 | INTEGER ng_q1 ! Number of grids that q1 < 0.0 | ||
520 | INTEGER ng_qsurf ! Number of grids that qsurf < 0.0 | ||
521 | INTEGER i | ||
522 | REAL zdu2, ztsolv | ||
523 | REAL ztvd, zscf | ||
524 | REAL zucf, zcr | ||
525 | REAL friv, frih | ||
526 | 28800 | REAL, DIMENSION(klon) :: FM, FH ! stability functions | |
527 | 28800 | REAL, DIMENSION(klon) :: cdmn, cdhn ! Drag coefficient in neutral conditions | |
528 | REAL zzzcd | ||
529 | |||
530 | LOGICAL, SAVE :: firstcall = .TRUE. | ||
531 | !$OMP THREADPRIVATE(firstcall) | ||
532 | INTEGER, SAVE :: iflag_corr_sta | ||
533 | !$OMP THREADPRIVATE(iflag_corr_sta) | ||
534 | INTEGER, SAVE :: iflag_corr_insta | ||
535 | !$OMP THREADPRIVATE(iflag_corr_insta) | ||
536 | |||
537 | !===================================================================c | ||
538 | ! Valeurs numeriques des constantes | ||
539 | !===================================================================c | ||
540 | |||
541 | |||
542 | ! Minimum du carre du vent | ||
543 | |||
544 | CEPDU2 = (0.1)**2 | ||
545 | ! Louis 1982 | ||
546 | |||
547 | CB=5.0 | ||
548 | CC=5.0 | ||
549 | CD=5.0 | ||
550 | |||
551 | |||
552 | ! King 2001 | ||
553 | |||
554 | C2=0.25 | ||
555 | C3=0.0625 | ||
556 | |||
557 | |||
558 | ! Louis 1979 | ||
559 | |||
560 | BPRIME=4.7 | ||
561 | B=9.4 | ||
562 | |||
563 | |||
564 | !MO | ||
565 | |||
566 | ALPHA=5.0 | ||
567 | |||
568 | |||
569 | ! ================================================================= c | ||
570 | ! Tests avant de commencer | ||
571 | ! Fuxing WANG, 04/03/2015 | ||
572 | ! To check if there are negative q1, qsurf values. | ||
573 | !====================================================================c | ||
574 | 14400 | ng_q1 = 0 ! Initialization | |
575 | 14400 | ng_qsurf = 0 ! Initialization | |
576 |
2/2✓ Branch 0 taken 6216792 times.
✓ Branch 1 taken 14400 times.
|
6231192 | DO i = 1, knon |
577 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6216792 times.
|
6216792 | IF (q1(i).LT.0.0) ng_q1 = ng_q1 + 1 |
578 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6216792 times.
|
6231192 | IF (qsurf(i).LT.0.0) ng_qsurf = ng_qsurf + 1 |
579 | ENDDO | ||
580 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 14400 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
14400 | IF (ng_q1.GT.0 .and. prt_level > 5) THEN |
581 | ✗ | WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !" | |
582 | ✗ | WRITE(lunout,*)" The total number of the grids is: ", ng_q1 | |
583 | ✗ | WRITE(lunout,*)" The negative q1 is set to zero " | |
584 | ! abort_message="voir ci-dessus" | ||
585 | ! CALL abort_physic(modname,abort_message,1) | ||
586 | ENDIF | ||
587 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 14400 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
14400 | IF (ng_qsurf.GT.0 .and. prt_level > 5) THEN |
588 | ✗ | WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !" | |
589 | ✗ | WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf | |
590 | ✗ | WRITE(lunout,*)" The negative qsurf is set to zero " | |
591 | ! abort_message="voir ci-dessus" | ||
592 | ! CALL abort_physic(modname,abort_message,1) | ||
593 | ENDIF | ||
594 | |||
595 | |||
596 | |||
597 | !=============================================================================c | ||
598 | ! Calcul du cdrag | ||
599 | !=============================================================================c | ||
600 | |||
601 | ! On choisit les fonctions de stabilite utilisees au premier appel | ||
602 | !************************************************************************** | ||
603 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 14399 times.
|
14400 | IF (firstcall) THEN |
604 | 1 | iflag_corr_sta=2 | |
605 | 1 | iflag_corr_insta=2 | |
606 | |||
607 | 1 | CALL getin_p('iflag_corr_sta',iflag_corr_sta) | |
608 | 1 | CALL getin_p('iflag_corr_insta',iflag_corr_insta) | |
609 | |||
610 | 1 | firstcall = .FALSE. | |
611 | ENDIF | ||
612 | |||
613 | !xxxxxxxxxxxxxxxxxxxxxxx | ||
614 |
2/2✓ Branch 0 taken 6216792 times.
✓ Branch 1 taken 14400 times.
|
6231192 | DO i = 1, knon ! Boucle sur l'horizontal |
615 | !xxxxxxxxxxxxxxxxxxxxxxx | ||
616 | |||
617 | |||
618 | ! calculs preliminaires: | ||
619 | !*********************** | ||
620 | |||
621 | |||
622 | 6216792 | zdu2 = MAX(CEPDU2, speed(i)**2) | |
623 | pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* & | ||
624 | 6216792 | (1.+ RETV * max(q1(i),0.0)))) ! negative q1 set to zero | |
625 | 6216792 | ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0)) ! negative qsurf set to zero | |
626 | ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*max(q1(i),0.0))) & | ||
627 | 6216792 | *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero | |
628 | 6216792 | zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd) | |
629 | |||
630 |
2/2✓ Branch 0 taken 2072264 times.
✓ Branch 1 taken 4144528 times.
|
6216792 | IF (iri1.EQ.1) THEN |
631 | 2072264 | zri(i) = ri1(i) | |
632 | ENDIF | ||
633 | |||
634 | ! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)): | ||
635 | !******************************************************************** | ||
636 | |||
637 | 6216792 | zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i))) | |
638 | 6216792 | cdmn(i) = zzzcd*zzzcd | |
639 | 6216792 | cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i)))) | |
640 | |||
641 | |||
642 | ! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi : | ||
643 | !******************************************************* | ||
644 | |||
645 | !'''''''''''''' | ||
646 | ! Cas instables | ||
647 | !'''''''''''''' | ||
648 | |||
649 |
2/2✓ Branch 0 taken 4033605 times.
✓ Branch 1 taken 2183187 times.
|
6231192 | IF (zri(i) .LT. 0.) THEN |
650 | |||
651 | |||
652 | ✗ | SELECT CASE (iflag_corr_insta) | |
653 | |||
654 | CASE (1) ! Louis 1979 + Mascart 1995 | ||
655 | |||
656 | ✗ | MU=LOG(MAX(z0m(i)/z0h(i),0.01)) | |
657 | ✗ | CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3) | |
658 | ✗ | PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3) | |
659 | ✗ | CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3) | |
660 | ✗ | PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3) | |
661 | CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & | ||
662 | & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i))) & | ||
663 | ✗ | & * ((zgeop1(i)/(RG*z0h(i)))**PH) | |
664 | CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & | ||
665 | & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & | ||
666 | ✗ | & * ((zgeop1(i)/(RG*z0m(i)))**PM) | |
667 | |||
668 | |||
669 | |||
670 | |||
671 | ✗ | FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i)))) | |
672 | ✗ | FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i)))) | |
673 | |||
674 | CASE (2) ! Louis 1982 | ||
675 | |||
676 | zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & | ||
677 | 4033605 | *(1.0+zgeop1(i)/(RG*z0m(i))))) | |
678 | 4033605 | FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) | |
679 | 4033605 | FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) | |
680 | |||
681 | |||
682 | CASE (3) ! Laurent Li | ||
683 | |||
684 | |||
685 | ✗ | FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) | |
686 | ✗ | FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) | |
687 | |||
688 | |||
689 | |||
690 | CASE default ! Louis 1982 | ||
691 | |||
692 | zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & | ||
693 | ✗ | *(1.0+zgeop1(i)/(RG*z0m(i))))) | |
694 | ✗ | FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) | |
695 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 4033605 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
4033605 | FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) |
696 | |||
697 | |||
698 | END SELECT | ||
699 | |||
700 | |||
701 | |||
702 | ! Calcul des drags | ||
703 | |||
704 | |||
705 | 4033605 | cdm(i)=cdmn(i)*FM(i) | |
706 | 4033605 | cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) | |
707 | |||
708 | |||
709 | ! Traitement particulier des cas oceaniques | ||
710 | ! on applique Miller et al 1992 en l'absence de gustiness | ||
711 | |||
712 |
2/2✓ Branch 0 taken 1889818 times.
✓ Branch 1 taken 2143787 times.
|
4033605 | IF (nsrf == is_oce) THEN |
713 | ! cdh(i)=f_cdrag_oce*cdhn(i)*FH(i) | ||
714 | |||
715 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1889818 times.
|
1889818 | IF(iflag_gusts==0) THEN |
716 | ✗ | zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.) | |
717 | ✗ | cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25) | |
718 | ENDIF | ||
719 | |||
720 | |||
721 | 1889818 | cdm(i)=MIN(cdm(i),cdmmax) | |
722 | 1889818 | cdh(i)=MIN(cdh(i),cdhmax) | |
723 | |||
724 | END IF | ||
725 | |||
726 | |||
727 | |||
728 | ELSE | ||
729 | |||
730 | !''''''''''''''' | ||
731 | ! Cas stables : | ||
732 | !''''''''''''''' | ||
733 | 2183187 | zri(i) = MIN(20.,zri(i)) | |
734 | |||
735 | ✗ | SELECT CASE (iflag_corr_sta) | |
736 | |||
737 | CASE (1) ! Louis 1979 + Mascart 1995 | ||
738 | |||
739 | ✗ | FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min) | |
740 | ✗ | FH(i)=FM(i) | |
741 | |||
742 | |||
743 | CASE (2) ! Louis 1982 | ||
744 | |||
745 | ✗ | zscf = SQRT(1.+CD*ABS(zri(i))) | |
746 | ✗ | FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) | |
747 | ✗ | FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) | |
748 | |||
749 | |||
750 | CASE (3) ! Laurent Li | ||
751 | |||
752 | ✗ | FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min) | |
753 | ✗ | FH(i)=FM(i) | |
754 | |||
755 | |||
756 | CASE (4) ! King 2001 | ||
757 | |||
758 |
2/2✓ Branch 0 taken 1357227 times.
✓ Branch 1 taken 825960 times.
|
2183187 | if (zri(i) .LT. C2/2.) then |
759 | 1357227 | FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min) | |
760 | 1357227 | FH(i)= FM(i) | |
761 | |||
762 | |||
763 | else | ||
764 | 825960 | FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min) | |
765 | 825960 | FH(i)= FM(i) | |
766 | endif | ||
767 | |||
768 | |||
769 | CASE (5) ! MO | ||
770 | |||
771 | ✗ | if (zri(i) .LT. 1./alpha) then | |
772 | |||
773 | ✗ | FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min) | |
774 | ✗ | FH(i)=FM(i) | |
775 | |||
776 | else | ||
777 | |||
778 | |||
779 | ✗ | FM(i)=MAX(1E-7,f_ri_cd_min) | |
780 | ✗ | FH(i)=FM(i) | |
781 | |||
782 | endif | ||
783 | |||
784 | |||
785 | |||
786 | |||
787 | |||
788 | CASE default ! Louis 1982 | ||
789 | |||
790 | ✗ | zscf = SQRT(1.+CD*ABS(zri(i))) | |
791 | ✗ | FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) | |
792 |
1/6✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2183187 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
2183187 | FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) |
793 | |||
794 | |||
795 | |||
796 | END SELECT | ||
797 | |||
798 | ! Calcul des drags | ||
799 | |||
800 | |||
801 | 2183187 | cdm(i)=cdmn(i)*FM(i) | |
802 | 2183187 | cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) | |
803 | |||
804 |
2/2✓ Branch 0 taken 288044 times.
✓ Branch 1 taken 1895143 times.
|
2183187 | IF(nsrf.EQ.is_oce) THEN |
805 | |||
806 | 288044 | cdh(i)=f_cdrag_oce*cdhn(i)*FH(i) | |
807 | 288044 | cdm(i)=MIN(cdm(i),cdmmax) | |
808 | 288044 | cdh(i)=MIN(cdh(i),cdhmax) | |
809 | |||
810 | ENDIF | ||
811 | |||
812 | |||
813 | ENDIF | ||
814 | |||
815 | !xxxxxxxxxxx | ||
816 | END DO ! Fin de la boucle sur l'horizontal | ||
817 | !xxxxxxxxxxx | ||
818 | ! ================================================================= c | ||
819 | |||
820 | 14400 | END SUBROUTINE cdragn_ri1 | |
821 | |||
822 | END MODULE cdrag_mod | ||
823 |