Context Navigation

cdrag_mod.F90 @ 5449

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

Rev	Line
[4725]	1	!$Id: cdrag_mod.F90 5144 2024-07-29 21:01:04Z fhourdin $
[5099]	2
	3
[3817]	4	MODULE cdrag_mod
[5099]	5
[5144]	6	! This module contains some procedures for calculation of the cdrag
	7	! coefficients for turbulent diffusion at surface
[5099]	8
[3817]	9	IMPLICIT NONE
	10
	11	CONTAINS
[5099]	12
[5144]	13	!****************************************************************************************
[5099]	14
[5144]	15	!r original routine svn3623
[5099]	16
[5144]	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)
[3817]	22
[5144]	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
[3817]	33
[5144]	34	IMPLICIT NONE
	35	! ================================================================= c
[5099]	36
[5144]	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
[5099]	42
[5144]	43	! Modified history:
	44	! writting on the 20/05/2016
	45	! modified on the 13/12/2016 to be adapted to LMDZ6
[5099]	46
[5144]	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
[5099]	66
[5144]	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
[5099]	71
[5144]	72	! iflag_corr_insta:
	73	! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h)
	74	! 2: Louis 1982
	75	! 3: Laurent Li
[5099]	76
[5144]	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)
[5099]	83
	84
[5144]	85	!*****************************************************************
	86	! Parametres d'entree
	87	!*****************************************************************
[3817]	88
[5144]	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
[3817]	100
[5144]	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
[3817]	106
[4722]	107
[3817]	108
[5144]	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
[3817]	115
[5144]	116	! Variables Locales
	117	!******************************************************************
[3817]	118
[5144]	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
[3817]	136
[5144]	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
[3817]	148
[5144]	149	!------------------RAJOUT POUR ECUME -------------------
[4722]	150
[5144]	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
[4722]	162
[5144]	163	LOGICAL :: OPRECIP
	164	LOGICAL :: OPWEBB
	165	LOGICAL :: OPERTFLUX
	166	LOGICAL :: LPRECIP
	167	LOGICAL :: LPWG
[4722]	168
[5144]	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)
[4722]	177
[5144]	178	!===================================================================c
	179	! Valeurs numeriques des constantes
	180	!===================================================================c
[4722]	181
	182
[5144]	183	! Minimum du carre du vent
[3817]	184
[5144]	185	CEPDU2 = (0.1)**2
[3817]	186
[5144]	187	! Louis 1982
[3817]	188
[5144]	189	CB = 5.0
	190	CC = 5.0
	191	CD = 5.0
[3817]	192
	193
[5144]	194	! King 2001
[3817]	195
[5144]	196	C2 = 0.25
	197	C3 = 0.0625
[3817]	198
	199
[5144]	200	! Louis 1979
[3817]	201
[5144]	202	BPRIME = 4.7
	203	B = 9.4
[3817]	204
	205
[5144]	206	!MO
[3817]	207
[5144]	208	ALPHA = 5.0
[3817]	209
[5144]	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
[3817]	214
[4478]	215
[5144]	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
[4481]	241
[3817]	242
	243
[5144]	244	!=============================================================================c
	245	! Calcul du cdrag
	246	!=============================================================================c
[3817]	247
[5144]	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.
[3817]	254
[5144]	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)
[3817]	258
[5144]	259	firstCALL = .FALSE.
	260	ENDIF
[3817]	261
[5144]	262	!------------------ Rajout (OT2018) --------------------
	263	!--------- Rajout pour itération sur rugosité ----------------
	264	rugos_itm(:, 2) = z0m
	265	rugos_itm(:, 1) = 3. * tol_it_rugos * z0m
[4722]	266
[5144]	267	rugos_ith(:, 2) = z0h !cp nouveau rugos_it
	268	rugos_ith(:, 1) = 3. * tol_it_rugos * z0h
	269	!--------------------------------------------------------------------
[4722]	270
[5144]	271	!xxxxxxxxxxxxxxxxxxxxxxx
	272	DO i = 1, knon ! Boucle sur l'horizontal
	273	!xxxxxxxxxxxxxxxxxxxxxxx
[3817]	274
	275
[5144]	276	! calculs preliminaires:
	277	!***********************
[3817]	278
[5144]	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
[5103]	299	mixte = .TRUE.
[5144]	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)
[3817]	310
[5144]	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)
[3817]	326
[5144]	327	cdmm(i) = coeffs(1)
	328	cdhh(i) = coeffs(2)
	329	cdm(i) = cdmm(i)
	330	cdh(i) = cdhh(i)
[3817]	331
[5144]	332	WRITE(, ) "ecume cd ch", cdmm(i), cdhh(i)
[5060]	333
[5144]	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)
[3817]	349
[5144]	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)
[3817]	356
[5144]	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
[4481]	381
	382
[5144]	383	!----------------------- Rajout des itérations --------------
	384	DO k = 1, nit_bulk
[3817]	385
[5144]	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))))
[3817]	391
[5144]	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
[3817]	445
[5144]	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
[3817]	489
[5144]	490	! Calcul des drags
[3817]	491
[5144]	492	cdmm(i) = cdmn(i) * FM(i)
	493	cdhh(i) = f_cdrag_ter * cdhn(i) * FH(i)
[3817]	494
[5144]	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
[3817]	499
[5144]	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)
[4725]	510
[5144]	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
[3817]	535
[5144]	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
[3817]	637	END MODULE cdrag_mod

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

Download in other formats:

Original Format