Context Navigation

cdrag_mod.F90 @ 5409

Last change on this file since 5409 was 5144, checked in by abarral, 5 months ago
Put YOMCST.h into modules
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 26.9 KB

Line
1	!$Id: cdrag_mod.F90 5144 2024-07-29 21:01:04Z evignon $
2
3
4	MODULE cdrag_mod
5
6	! This module contains some procedures for calculation of the cdrag
7	! coefficients for turbulent diffusion at surface
8
9	IMPLICIT NONE
10
11	CONTAINS
12
13	!****************************************************************************************
14
15	!r original routine svn3623
16
17	SUBROUTINE cdrag(knon, nsrf, &
18	speed, t1, q1, zgeop1, &
19	psol, pblh, tsurf, qsurf, z0m, z0h, &
20	ri_in, iri_in, &
21	cdm, cdh, zri, pref, prain, tsol, pat1)
22
23	USE dimphy
24	USE coare_cp_mod, ONLY: coare_cp
25	USE coare30_flux_cnrm_mod, ONLY: coare30_flux_cnrm
26	USE indice_sol_mod
27	USE lmdz_print_control, ONLY: lunout, prt_level
28	USE lmdz_ioipsl_getin_p, ONLY: getin_p
29	USE lmdz_atke_turbulence_ini, ONLY: smmin, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut
30	USE lmdz_clesphys
31	USE lmdz_yoethf
32	USE lmdz_yomcst
33
34	IMPLICIT NONE
35	! ================================================================= c
36
37	! Objet : calcul des cdrags pour le moment (pcfm) et
38	! les flux de chaleur sensible et latente (pcfh) d'apr??s
39	! Louis 1982, Louis 1979, King et al 2001
40	! ou Zilitinkevich et al 2002 pour les cas stables, Louis 1979
41	! et 1982 pour les cas instables
42
43	! Modified history:
44	! writting on the 20/05/2016
45	! modified on the 13/12/2016 to be adapted to LMDZ6
46
47	! References:
48	! Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the
49	! atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202.
50	! Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the
51	! operational PBL parametrization at ECMWF'. Workshop on boundary layer
52	! parametrization, November 1981, ECMWF, Reading, England.
53	! Page: 19. Equations in Table 1.
54	! Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS
55	! IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM
56	! Boundary-Layer Meteorology 72: 331-344
57	! Anton Beljaars. May 1992. The parametrization of the planetary boundary layer.
58	! European Centre for Medium-Range Weather Forecasts.
59	! Equations: 110-113. Page 40.
60	! Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF
61	! model to the parameterization of evaporation from the tropical oceans. J.
62	! Climate, 5:418-434.
63	! King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate
64	! to surface and boundary-layer flux parametrizations
65	! QJRMS, 127, pp 779-794
66
67	! ================================================================= c
68	! ================================================================= c
69	! On choisit le couple de fonctions de correction avec deux flags:
70	! Un pour les cas instables, un autre pour les cas stables
71
72	! iflag_corr_insta:
73	! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
74	! 2: Louis 1982
75	! 3: Laurent Li
76
77	! iflag_corr_sta:
78	! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
79	! 2: Louis 1982
80	! 3: Laurent Li
81	! 4: King 2001 (SHARP)
82	! 5: MO 1st order theory (allow collapse of turbulence)
83
84
85	!*****************************************************************
86	! Parametres d'entree
87	!*****************************************************************
88
89	INTEGER, INTENT(IN) :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface
90	REAL, DIMENSION(klon), INTENT(IN) :: speed ! module du vent au 1er niveau du modele
91	REAL, DIMENSION(klon), INTENT(IN) :: zgeop1 ! geopotentiel au 1er niveau du modele
92	REAL, DIMENSION(klon), INTENT(IN) :: tsurf ! Surface temperature (K)
93	REAL, DIMENSION(klon), INTENT(IN) :: qsurf ! Surface humidity (Kg/Kg)
94	REAL, DIMENSION(klon), INTENT(INOUT) :: z0m, z0h ! Rugosity at surface (m)
95	REAL, DIMENSION(klon), INTENT(IN) :: ri_in ! Input Richardson 1st layer for first guess calculations of screen var.
96	INTEGER, INTENT(IN) :: iri_in! iflag to activate cdrag calculation using ri1
97	REAL, DIMENSION(klon), INTENT(IN) :: t1 ! Temperature au premier niveau (K)
98	REAL, DIMENSION(klon), INTENT(IN) :: q1 ! humidite specifique au premier niveau (kg/kg)
99	REAL, DIMENSION(klon), INTENT(IN) :: psol ! pression au sol
100
101	!------------------ Rajout pour COARE (OT2018) --------------------
102	REAL, DIMENSION(klon), INTENT(INOUT) :: pblh !hauteur de CL
103	REAL, DIMENSION(klon), INTENT(IN) :: prain !rapport au precipitation
104	REAL, DIMENSION(klon), INTENT(IN) :: tsol !SST imposé sur la surface oceanique
105	REAL, DIMENSION(klon), INTENT(IN) :: pat1 !pression premier lev
106
107
108
109	! Parametres de sortie
110	!******************************************************************
111	REAL, DIMENSION(klon), INTENT(OUT) :: cdm ! Drag coefficient for momentum
112	REAL, DIMENSION(klon), INTENT(OUT) :: cdh ! Drag coefficient for heat flux
113	REAL, DIMENSION(klon), INTENT(OUT) :: zri ! Richardson number
114	REAL, DIMENSION(klon), INTENT(INOUT) :: pref ! Pression au niveau zgeop/RG
115
116	! Variables Locales
117	!******************************************************************
118
119	REAL, PARAMETER :: CKAP = 0.40, CKAPT = 0.42
120	REAL CEPDU2
121	REAL ALPHA
122	REAL CB, CC, CD, C2, C3
123	REAL MU, CM, CH, B, CMstar, CHstar
124	REAL PM, PH, BPRIME
125	INTEGER ng_q1 ! Number of grids that q1 < 0.0
126	INTEGER ng_qsurf ! Number of grids that qsurf < 0.0
127	INTEGER i, k
128	REAL zdu2, ztsolv
129	REAL ztvd, zscf
130	REAL zucf, zcr
131	REAL, DIMENSION(klon) :: FM, FH ! stability functions
132	REAL, DIMENSION(klon) :: cdmn, cdhn ! Drag coefficient in neutral conditions
133	REAL zzzcd
134	REAL, DIMENSION(klon) :: sm, prandtl ! Stability function from atke turbulence scheme
135	REAL ri0, ri1, cn ! to have dimensionless term in sm and prandtl
136
137	!------------------ Rajout (OT2018) --------------------
138	!------------------ Rajout pour les appelles BULK (OT) --------------------
139	REAL, DIMENSION(klon, 2) :: rugos_itm
140	REAL, DIMENSION(klon, 2) :: rugos_ith
141	REAL, PARAMETER :: tol_it_rugos = 1.e-4
142	REAL, DIMENSION(3) :: coeffs
143	LOGICAL :: mixte
144	REAL :: z_0m
145	REAL :: z_0h
146	REAL, DIMENSION(klon) :: cdmm
147	REAL, DIMENSION(klon) :: cdhh
148
149	!------------------RAJOUT POUR ECUME -------------------
150
151	REAL, DIMENSION(klon) :: PQSAT
152	REAL, DIMENSION(klon) :: PSFTH
153	REAL, DIMENSION(klon) :: PFSTQ
154	REAL, DIMENSION(klon) :: PUSTAR
155	REAL, DIMENSION(klon) :: PCD ! Drag coefficient for momentum
156	REAL, DIMENSION(klon) :: PCDN ! Drag coefficient for momentum
157	REAL, DIMENSION(klon) :: PCH ! Drag coefficient for momentum
158	REAL, DIMENSION(klon) :: PCE ! Drag coefficient for momentum
159	REAL, DIMENSION(klon) :: PRI
160	REAL, DIMENSION(klon) :: PRESA
161	REAL, DIMENSION(klon) :: PSSS
162
163	LOGICAL :: OPRECIP
164	LOGICAL :: OPWEBB
165	LOGICAL :: OPERTFLUX
166	LOGICAL :: LPRECIP
167	LOGICAL :: LPWG
168
169	LOGICAL, SAVE :: firstCALL = .TRUE.
170	!$OMP THREADPRIVATE(firstcall)
171	INTEGER, SAVE :: iflag_corr_sta
172	!$OMP THREADPRIVATE(iflag_corr_sta)
173	INTEGER, SAVE :: iflag_corr_insta
174	!$OMP THREADPRIVATE(iflag_corr_insta)
175	LOGICAL, SAVE :: ok_cdrag_iter
176	!$OMP THREADPRIVATE(ok_cdrag_iter)
177
178	!===================================================================c
179	! Valeurs numeriques des constantes
180	!===================================================================c
181
182
183	! Minimum du carre du vent
184
185	CEPDU2 = (0.1)**2
186
187	! Louis 1982
188
189	CB = 5.0
190	CC = 5.0
191	CD = 5.0
192
193
194	! King 2001
195
196	C2 = 0.25
197	C3 = 0.0625
198
199
200	! Louis 1979
201
202	BPRIME = 4.7
203	B = 9.4
204
205
206	!MO
207
208	ALPHA = 5.0
209
210	! Consistent with atke scheme
211	cn = (1. / sqrt(cepsilon))**(2. / 3.)
212	ri0 = 2. / rpi * (cinf - cn) * ric / cn
213	ri1 = -2. / rpi * (pr_asym - pr_neut) / pr_slope
214
215
216	! ================================================================= c
217	! Tests avant de commencer
218	! Fuxing WANG, 04/03/2015
219	! To check if there are negative q1, qsurf values.
220	!====================================================================c
221	ng_q1 = 0 ! Initialization
222	ng_qsurf = 0 ! Initialization
223	DO i = 1, knon
224	IF (q1(i)<0.0) ng_q1 = ng_q1 + 1
225	IF (qsurf(i)<0.0) ng_qsurf = ng_qsurf + 1
226	ENDDO
227	IF (ng_q1>0 .AND. prt_level > 5) THEN
228	WRITE(lunout, )" ** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !"
229	WRITE(lunout, *)" The total number of the grids is: ", ng_q1
230	WRITE(lunout, *)" The negative q1 is set to zero "
231	! abort_message="voir ci-dessus"
232	! CALL abort_physic(modname,abort_message,1)
233	ENDIF
234	IF (ng_qsurf>0 .AND. prt_level > 5) THEN
235	WRITE(lunout, )" ** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !"
236	WRITE(lunout, *)" The total number of the grids is: ", ng_qsurf
237	WRITE(lunout, *)" The negative qsurf is set to zero "
238	! abort_message="voir ci-dessus"
239	! CALL abort_physic(modname,abort_message,1)
240	ENDIF
241
242
243
244	!=============================================================================c
245	! Calcul du cdrag
246	!=============================================================================c
247
248	! On choisit les fonctions de stabilite utilisees au premier appel
249	!**************************************************************************
250	IF (firstcall) THEN
251	iflag_corr_sta = 2
252	iflag_corr_insta = 2
253	ok_cdrag_iter = .FALSE.
254
255	CALL getin_p('iflag_corr_sta', iflag_corr_sta)
256	CALL getin_p('iflag_corr_insta', iflag_corr_insta)
257	CALL getin_p('ok_cdrag_iter', ok_cdrag_iter)
258
259	firstCALL = .FALSE.
260	ENDIF
261
262	!------------------ Rajout (OT2018) --------------------
263	!--------- Rajout pour itération sur rugosité ----------------
264	rugos_itm(:, 2) = z0m
265	rugos_itm(:, 1) = 3. * tol_it_rugos * z0m
266
267	rugos_ith(:, 2) = z0h !cp nouveau rugos_it
268	rugos_ith(:, 1) = 3. * tol_it_rugos * z0h
269	!--------------------------------------------------------------------
270
271	!xxxxxxxxxxxxxxxxxxxxxxx
272	DO i = 1, knon ! Boucle sur l'horizontal
273	!xxxxxxxxxxxxxxxxxxxxxxx
274
275
276	! calculs preliminaires:
277	!***********************
278
279	zdu2 = MAX(CEPDU2, speed(i)**2)
280	pref(i) = EXP(LOG(psol(i)) - zgeop1(i) / (RD * t1(i) * &
281	(1. + RETV * max(q1(i), 0.0)))) ! negative q1 set to zero
282	ztsolv = tsurf(i) * (1.0 + RETV * max(qsurf(i), 0.0)) ! negative qsurf set to zero
283	ztvd = (t1(i) + zgeop1(i) / RCPD / (1. + RVTMP2 * q1(i))) &
284	* (1. + RETV * max(q1(i), 0.0)) ! negative q1 set to zero
285
286	!------------------ Rajout (OT2018) --------------------
287	!-------------- ON DUPLIQUE POUR METTRE ITERATION SUR OCEAN ------------------------
288	IF (iri_in==1) THEN
289	zri(i) = ri_in(i)
290	ENDIF
291
292	IF (nsrf == is_oce) THEN
293
294	!------------------ Pour Core 2 choix Core Pur ou Core Mixte --------------------
295	IF ((choix_bulk > 1 .AND. choix_bulk < 4) .AND. (nsrf == is_oce)) THEN
296	IF(choix_bulk == 2) THEN
297	mixte = .FALSE.
298	ELSE
299	mixte = .TRUE.
300	ENDIF
301	CALL clc_core_cp (sqrt(zdu2), t1(i) - tsurf(i), q1(i) - qsurf(i), t1(i), q1(i), &
302	zgeop1(i) / RG, zgeop1(i) / RG, zgeop1(i) / RG, &
303	psol(i), nit_bulk, mixte, &
304	coeffs, z_0m, z_0h)
305	cdmm(i) = coeffs(1)
306	cdhh(i) = coeffs(2)
307	cdm(i) = cdmm(i)
308	cdh(i) = cdhh(i)
309	WRITE(, ) "clc_core cd ch", cdmm(i), cdhh(i)
310
311	!------------------ Pour ECUME --------------------
312	ELSE IF ((choix_bulk == 4) .AND. (nsrf == is_oce)) THEN
313	OPRECIP = .FALSE.
314	OPWEBB = .FALSE.
315	OPERTFLUX = .FALSE.
316	IF (nsrf == is_oce) THEN
317	PSSS = 0.0
318	ENDIF
319	CALL ini_csts
320	CALL ecumev6_flux(z_0m, t1(i), tsurf(i), &
321	q1(i), qsurf(i), sqrt(zdu2), zgeop1(i) / RG, PSSS, zgeop1(i) / RG, &
322	psol(i), pat1(i), OPRECIP, OPWEBB, &
323	PSFTH, PFSTQ, PUSTAR, PCD, PCDN, &
324	PCH, PCE, PRI, PRESA, prain, &
325	z_0h, OPERTFLUX, coeffs)
326
327	cdmm(i) = coeffs(1)
328	cdhh(i) = coeffs(2)
329	cdm(i) = cdmm(i)
330	cdh(i) = cdhh(i)
331
332	WRITE(, ) "ecume cd ch", cdmm(i), cdhh(i)
333
334	!------------------ Pour COARE CNRM --------------------
335	ELSE IF ((choix_bulk == 5) .AND. (nsrf == is_oce)) THEN
336	LPRECIP = .FALSE.
337	LPWG = .FALSE.
338	CALL ini_csts
339	block
340	REAL, DIMENSION(1) :: z0m_1d, z_0h_1d, sqrt_zdu2_1d, zgeop1_rg_1d ! convert scalar to 1D for call
341	z0m_1d = z0m
342	z_0h_1d = z0h
343	sqrt_zdu2_1d = sqrt(zdu2)
344	zgeop1_rg_1d = zgeop1(i) / RG
345	CALL coare30_flux_cnrm(z0m_1d, t1(i), tsurf(i), q1(i), &
346	sqrt_zdu2_1d, zgeop1_rg_1d, zgeop1_rg_1d, psol(i), qsurf(i), PQSAT, &
347	PSFTH, PFSTQ, PUSTAR, PCD, PCDN, PCH, PCE, PRI, &
348	PRESA, prain, pat1(i), z_0h_1d, LPRECIP, LPWG, coeffs)
349
350	end block
351	cdmm(i) = coeffs(1)
352	cdhh(i) = coeffs(2)
353	cdm(i) = cdmm(i)
354	cdh(i) = cdhh(i)
355	WRITE(, ) "coare CNRM cd ch", cdmm(i), cdhh(i)
356
357	!------------------ Pour COARE Maison --------------------
358	ELSE IF ((choix_bulk == 1) .AND. (nsrf == is_oce)) THEN
359	IF (pblh(i) == 0.) THEN
360	pblh(i) = 1500.
361	ENDIF
362	WRITE(, ) "debug size", size(coeffs)
363	CALL coare_cp(sqrt(zdu2), t1(i) - tsurf(i), q1(i) - qsurf(i), &
364	t1(i), q1(i), &
365	zgeop1(i) / RG, zgeop1(i) / RG, zgeop1(i) / RG, &
366	psol(i), pblh(i), &
367	nit_bulk, &
368	coeffs, z_0m, z_0h)
369	cdmm(i) = coeffs(1)
370	cdhh(i) = coeffs(3)
371	cdm(i) = cdmm(i)
372	cdh(i) = cdhh(i)
373	WRITE(, ) "coare cd ch", cdmm(i), cdhh(i)
374	ELSE
375	!------------------ Pour La param LMDZ (ocean) --------------------
376	zri(i) = zgeop1(i) * (ztvd - ztsolv) / (zdu2 * ztvd)
377	IF (iri_in==1) THEN
378	zri(i) = ri_in(i)
379	ENDIF
380	ENDIF
381
382
383	!----------------------- Rajout des itérations --------------
384	DO k = 1, nit_bulk
385
386	! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
387	!********************************************************************
388	zzzcd = CKAP / LOG(1. + zgeop1(i) / (RG * rugos_itm(i, 2)))
389	cdmn(i) = zzzcd * zzzcd
390	cdhn(i) = zzzcd * (CKAP / LOG(1. + zgeop1(i) / (RG * rugos_ith(i, 2))))
391
392	! Calcul des fonctions de stabilite FMs, FHs, FMi, FHi :
393	!*******************************************************
394	IF (zri(i) < 0.) THEN
395	SELECT CASE (iflag_corr_insta)
396	CASE (1) ! Louis 1979 + Mascart 1995
397	MU = LOG(MAX(z0m(i) / z0h(i), 0.01))
398	CMstar = 6.8741 + 2.6933 * MU - 0.3601 * (MU*2) + 0.0154 (MU**3)
399	PM = 0.5233 - 0.0815 * MU + 0.0135 * (MU*2) - 0.001 (MU**3)
400	CHstar = 3.2165 + 4.3431 * MU + 0.536 * (MU*2) - 0.0781 (MU**3)
401	PH = 0.5802 - 0.1571 * MU + 0.0327 * (MU*2) - 0.0026 (MU**3)
402	CH = CHstar * B * CKAP / LOG(z0m(i) + zgeop1(i) / (RG * z0m(i))) &
403	* CKAPT / LOG(z0h(i) + zgeop1(i) / (RG * z0h(i))) &
404	* ((zgeop1(i) / (RG * z0h(i)))**PH)
405	CM = CMstar * B * CKAP / LOG(z0m(i) + zgeop1(i) / (RG * z0m(i))) &
406	* CKAP / LOG(z0m(i) + zgeop1(i) / (RG * z0m(i))) &
407	* ((zgeop1(i) / (RG * z0m(i)))**PM)
408	FM(i) = 1. - B * zri(i) / (1. + CM * SQRT(ABS(zri(i))))
409	FH(i) = 1. - B * zri(i) / (1. + CH * SQRT(ABS(zri(i))))
410	CASE (2) ! Louis 1982
411	zucf = 1. / (1. + 3.0 * CB * CC * cdmn(i) * SQRT(ABS(zri(i)) &
412	* (1.0 + zgeop1(i) / (RG * z0m(i)))))
413	FM(i) = AMAX1((1. - 2.0 * CB * zri(i) * zucf), f_ri_cd_min)
414	FH(i) = AMAX1((1. - 3.0 * CB * zri(i) * zucf), f_ri_cd_min)
415	CASE (3) ! Laurent Li
416	FM(i) = MAX(SQRT(1.0 - 18.0 * zri(i)), f_ri_cd_min)
417	FH(i) = MAX(SQRT(1.0 - 18.0 * zri(i)), f_ri_cd_min)
418	CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
419	sm(i) = 2. / rpi * (cinf - cn) * atan(-zri(i) / ri0) + cn
420	prandtl(i) = -2. / rpi * (pr_asym - pr_neut) * atan(zri(i) / ri1) + pr_neut
421	FM(i) = MAX((sm(i)*(3. / 2.) sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1. / 2.)), f_ri_cd_min)
422	FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
423	CASE default ! Louis 1982
424	zucf = 1. / (1. + 3.0 * CB * CC * cdmn(i) * SQRT(ABS(zri(i)) &
425	* (1.0 + zgeop1(i) / (RG * z0m(i)))))
426	FM(i) = AMAX1((1. - 2.0 * CB * zri(i) * zucf), f_ri_cd_min)
427	FH(i) = AMAX1((1. - 3.0 * CB * zri(i) * zucf), f_ri_cd_min)
428	END SELECT
429	! Calcul des drags
430	cdmm(i) = cdmn(i) * FM(i)
431	cdhh(i) = f_cdrag_ter * cdhn(i) * FH(i)
432	! Traitement particulier des cas oceaniques
433	! on applique Miller et al 1992 en l'absence de gustiness
434	IF (nsrf == is_oce) THEN
435	! cdh(i)=f_cdrag_ocecdhn(i)FH(i)
436	IF (iflag_gusts==0) THEN
437	zcr = (0.0016 / (cdmn(i) * SQRT(zdu2))) * ABS(ztvd - ztsolv)**(1. / 3.)
438	cdhh(i) = f_cdrag_oce * cdhn(i) * (1.0 + zcr1.25)(1. / 1.25)
439	ENDIF
440	cdmm(i) = MIN(cdmm(i), cdmmax)
441	cdhh(i) = MIN(cdhh(i), cdhmax)
442	! WRITE(,) "LMDZ cd ch",cdmm(i),cdhh(i)
443	END IF
444	ELSE
445
446	!'''''''''''''''
447	! Cas stables :
448	!'''''''''''''''
449	zri(i) = MIN(20., zri(i))
450	SELECT CASE (iflag_corr_sta)
451	CASE (1) ! Louis 1979 + Mascart 1995
452	FM(i) = MAX(1. / ((1 + BPRIME * zri(i))**2), f_ri_cd_min)
453	FH(i) = FM(i)
454	CASE (2) ! Louis 1982
455	zscf = SQRT(1. + CD * ABS(zri(i)))
456	FM(i) = AMAX1(1. / (1. + 2. * CB * zri(i) / zscf), f_ri_cd_min)
457	FH(i) = AMAX1(1. / (1. + 3. * CB * zri(i) * zscf), f_ri_cd_min)
458	CASE (3) ! Laurent Li
459	FM(i) = MAX(1.0 / (1.0 + 10.0 * zri(i) * (1 + 8.0 * zri(i))), f_ri_cd_min)
460	FH(i) = FM(i)
461	CASE (4) ! King 2001
462	IF (zri(i) < C2 / 2.) THEN
463	FM(i) = MAX((1. - zri(i) / C2)**2, f_ri_cd_min)
464	FH(i) = FM(i)
465	ELSE
466	FM(i) = MAX(C3 * ((C2 / zri(i))**2), f_ri_cd_min)
467	FH(i) = FM(i)
468	ENDIF
469	CASE (5) ! MO
470	IF (zri(i) < 1. / alpha) THEN
471	FM(i) = MAX((1. - alpha * zri(i))**2, f_ri_cd_min)
472	FH(i) = FM(i)
473	else
474	FM(i) = MAX(1E-7, f_ri_cd_min)
475	FH(i) = FM(i)
476	endif
477	CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
478	sm(i) = MAX(smmin, cn * (1. - zri(i) / ric))
479	! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019
480	prandtl(i) = pr_neut * exp(-pr_slope / pr_neut * zri(i) + zri(i) / pr_neut) &
481	+ zri(i) * pr_slope
482	FM(i) = MAX((sm(i)*(3. / 2.) sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1. / 2.)), f_ri_cd_min)
483	FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
484	CASE default ! Louis 1982
485	zscf = SQRT(1. + CD * ABS(zri(i)))
486	FM(i) = AMAX1(1. / (1. + 2. * CB * zri(i) / zscf), f_ri_cd_min)
487	FH(i) = AMAX1(1. / (1. + 3. * CB * zri(i) * zscf), f_ri_cd_min)
488	END SELECT
489
490	! Calcul des drags
491
492	cdmm(i) = cdmn(i) * FM(i)
493	cdhh(i) = f_cdrag_ter * cdhn(i) * FH(i)
494
495	IF (choix_bulk == 0) THEN
496	cdm(i) = cdmn(i) * FM(i)
497	cdh(i) = f_cdrag_ter * cdhn(i) * FH(i)
498	ENDIF
499
500	IF (nsrf==is_oce) THEN
501	cdhh(i) = f_cdrag_oce * cdhn(i) * FH(i)
502	cdmm(i) = MIN(cdmm(i), cdmmax)
503	cdhh(i) = MIN(cdhh(i), cdhmax)
504	ENDIF
505	IF (ok_cdrag_iter) THEN
506	rugos_itm(i, 1) = rugos_itm(i, 2)
507	rugos_ith(i, 1) = rugos_ith(i, 2)
508	rugos_itm(i, 2) = 0.018 * cdmm(i) * (speed(i)) / RG &
509	+ 0.11 * 14e-6 / SQRT(cdmm(i) * zdu2)
510
511	!---------- Version SEPARATION DES Z0 ----------------------
512	IF (iflag_z0_oce==0) THEN
513	rugos_ith(i, 2) = rugos_itm(i, 2)
514	ELSE IF (iflag_z0_oce==1) THEN
515	rugos_ith(i, 2) = 0.40 * 14e-6 / SQRT(cdmm(i) * zdu2)
516	ENDIF
517	ENDIF
518	ENDIF
519	IF (ok_cdrag_iter) THEN
520	rugos_itm(i, 2) = MAX(1.5e-05, rugos_itm(i, 2))
521	rugos_ith(i, 2) = MAX(1.5e-05, rugos_ith(i, 2))
522	ENDIF
523	ENDDO
524	IF (nsrf==is_oce) THEN
525	cdm(i) = MIN(cdmm(i), cdmmax)
526	cdh(i) = MIN(cdhh(i), cdhmax)
527	ENDIF
528	z0m = rugos_itm(:, 2)
529	z0h = rugos_ith(:, 2)
530	ELSE ! (nsrf == is_oce)
531	zri(i) = zgeop1(i) * (ztvd - ztsolv) / (zdu2 * ztvd)
532	IF (iri_in==1) THEN
533	zri(i) = ri_in(i)
534	ENDIF
535
536	! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)):
537	!********************************************************************
538	zzzcd = CKAP / LOG(1. + zgeop1(i) / (RG * z0m(i)))
539	cdmn(i) = zzzcd * zzzcd
540	cdhn(i) = zzzcd * (CKAP / LOG(1. + zgeop1(i) / (RG * z0h(i))))
541
542
543	! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi :
544	!*******************************************************
545	!''''''''''''''
546	! Cas instables
547	!''''''''''''''
548	IF (zri(i) < 0.) THEN
549	SELECT CASE (iflag_corr_insta)
550	CASE (1) ! Louis 1979 + Mascart 1995
551	MU = LOG(MAX(z0m(i) / z0h(i), 0.01))
552	CMstar = 6.8741 + 2.6933 * MU - 0.3601 * (MU*2) + 0.0154 (MU**3)
553	PM = 0.5233 - 0.0815 * MU + 0.0135 * (MU*2) - 0.001 (MU**3)
554	CHstar = 3.2165 + 4.3431 * MU + 0.536 * (MU*2) - 0.0781 (MU**3)
555	PH = 0.5802 - 0.1571 * MU + 0.0327 * (MU*2) - 0.0026 (MU**3)
556	CH = CHstar * B * CKAP / LOG(z0m(i) + zgeop1(i) / (RG * z0m(i))) &
557	* CKAPT / LOG(z0h(i) + zgeop1(i) / (RG * z0h(i))) &
558	* ((zgeop1(i) / (RG * z0h(i)))**PH)
559	CM = CMstar * B * CKAP / LOG(z0m(i) + zgeop1(i) / (RG * z0m(i))) &
560	* CKAP / LOG(z0m(i) + zgeop1(i) / (RG * z0m(i))) &
561	* ((zgeop1(i) / (RG * z0m(i)))**PM)
562	FM(i) = 1. - B * zri(i) / (1. + CM * SQRT(ABS(zri(i))))
563	FH(i) = 1. - B * zri(i) / (1. + CH * SQRT(ABS(zri(i))))
564	CASE (2) ! Louis 1982
565	zucf = 1. / (1. + 3.0 * CB * CC * cdmn(i) * SQRT(ABS(zri(i)) &
566	* (1.0 + zgeop1(i) / (RG * z0m(i)))))
567	FM(i) = AMAX1((1. - 2.0 * CB * zri(i) * zucf), f_ri_cd_min)
568	FH(i) = AMAX1((1. - 3.0 * CB * zri(i) * zucf), f_ri_cd_min)
569	CASE (3) ! Laurent Li
570	FM(i) = MAX(SQRT(1.0 - 18.0 * zri(i)), f_ri_cd_min)
571	FH(i) = MAX(SQRT(1.0 - 18.0 * zri(i)), f_ri_cd_min)
572	CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
573	sm(i) = 2. / rpi * (cinf - cn) * atan(-zri(i) / ri0) + cn
574	prandtl(i) = -2. / rpi * (pr_asym - pr_neut) * atan(zri(i) / ri1) + pr_neut
575	FM(i) = MAX((sm(i)*(3. / 2.) sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1. / 2.)), f_ri_cd_min)
576	FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
577	CASE default ! Louis 1982
578	zucf = 1. / (1. + 3.0 * CB * CC * cdmn(i) * SQRT(ABS(zri(i)) &
579	* (1.0 + zgeop1(i) / (RG * z0m(i)))))
580	FM(i) = AMAX1((1. - 2.0 * CB * zri(i) * zucf), f_ri_cd_min)
581	FH(i) = AMAX1((1. - 3.0 * CB * zri(i) * zucf), f_ri_cd_min)
582	END SELECT
583	! Calcul des drags
584	cdm(i) = cdmn(i) * FM(i)
585	cdh(i) = f_cdrag_ter * cdhn(i) * FH(i)
586	ELSE
587	!'''''''''''''''
588	! Cas stables :
589	!'''''''''''''''
590	zri(i) = MIN(20., zri(i))
591	SELECT CASE (iflag_corr_sta)
592	CASE (1) ! Louis 1979 + Mascart 1995
593	FM(i) = MAX(1. / ((1 + BPRIME * zri(i))**2), f_ri_cd_min)
594	FH(i) = FM(i)
595	CASE (2) ! Louis 1982
596	zscf = SQRT(1. + CD * ABS(zri(i)))
597	FM(i) = AMAX1(1. / (1. + 2. * CB * zri(i) / zscf), f_ri_cd_min)
598	FH(i) = AMAX1(1. / (1. + 3. * CB * zri(i) * zscf), f_ri_cd_min)
599	CASE (3) ! Laurent Li
600	FM(i) = MAX(1.0 / (1.0 + 10.0 * zri(i) * (1 + 8.0 * zri(i))), f_ri_cd_min)
601	FH(i) = FM(i)
602	CASE (4) ! King 2001
603	IF (zri(i) < C2 / 2.) THEN
604	FM(i) = MAX((1. - zri(i) / C2)**2, f_ri_cd_min)
605	FH(i) = FM(i)
606	else
607	FM(i) = MAX(C3 * ((C2 / zri(i))**2), f_ri_cd_min)
608	FH(i) = FM(i)
609	endif
610	CASE (5) ! MO
611	IF (zri(i) < 1. / alpha) THEN
612	FM(i) = MAX((1. - alpha * zri(i))**2, f_ri_cd_min)
613	FH(i) = FM(i)
614	else
615	FM(i) = MAX(1E-7, f_ri_cd_min)
616	FH(i) = FM(i)
617	endif
618	CASE (6) ! Consistent with turbulence scheme (in stationary case) derived in atke (2023)
619	sm(i) = MAX(0., cn * (1. - zri(i) / ric))
620	prandtl(i) = pr_neut + zri(i) * pr_slope
621	FM(i) = MAX((sm(i)*(3. / 2.) sqrt(cepsilon) * (1 - zri(i) / prandtl(i))**(1. / 2.)), f_ri_cd_min)
622	FH(i) = MAX((FM(i) / prandtl(i)), f_ri_cd_min)
623	CASE default ! Louis 1982
624	zscf = SQRT(1. + CD * ABS(zri(i)))
625	FM(i) = AMAX1(1. / (1. + 2. * CB * zri(i) / zscf), f_ri_cd_min)
626	FH(i) = AMAX1(1. / (1. + 3. * CB * zri(i) * zscf), f_ri_cd_min)
627	END SELECT
628	! Calcul des drags
629	cdm(i) = cdmn(i) * FM(i)
630	cdh(i) = f_cdrag_ter * cdhn(i) * FH(i)
631	ENDIF
632	ENDIF ! fin du if (nsrf == is_oce)
633	END DO ! Fin de la boucle sur l'horizontal
634
635	END SUBROUTINE cdrag
636
637	END MODULE cdrag_mod

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

Download in other formats:

Original Format