Context Navigation

cdrag_mod.F90 @ 4722

Last change on this file since 4722 was 4722, checked in by Laurent Fairhead, 7 months ago
Modification by O. Torres to the cdrag routines to include different bulk formulae to calculate cdrag coefficients over ocean as well as an iteration of that calculation. The iteration is controlled by flag ok_cdrag_iter which if set to FALSE by default to converge with previous results. The choice of bulk formulae is set with the choix_bulk parameter The number of iterations to run is set with nit_bulk OT, PB, CD, LF
Property svn:executable set to ``*
File size: 25.1 KB

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

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

Download in other formats:

Original Format