Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

vdif_cd_mod.F90

Last change on this file was 3333, checked in by llange, 6 months ago
Mars PCM Fixing a bug in vdif_cd: a "residual", used as criterion to enter an iterative loop, was wrongly initialized. Hence, for some points, the algorithm does not go into the loop, and a wrong value of Cd, Ch was computed. Also some cleaning/small fixing with save variables with OMP. LL
File size: 14.3 KB

Rev	Line
[3219]	1	MODULE vdif_cd_mod
[3150]	2
	3	!======================================================================================================================!
	4	! Subject:
	5	!---------
	6	! Module used to compute the exchange coefficient in the surface layers Cd; Ch
	7	!----------------------------------------------------------------------------------------------------------------------!
	8	! Reference:
	9	!-----------
	10	! Colaïtis, A., Spiga, A., Hourdin, F., Rio, C., Forget, F., and Millour, E. (2013), A thermal plume model for the Martian convective boundary layer, J. Geophys. Res. Planets, 118, 1468–1487, doi:10.1002/jgre.20104.
	11	!
	12	!======================================================================================================================!
	13
	14	IMPLICIT NONE
	15
	16	CONTAINS
	17
[3325]	18	SUBROUTINE vdif_cd(ngrid,nlay,nslope,pz0,pg,pz,pp,pu,pv,wstar,pts,ph,mumean,pqvap,pqsurf,write_outputs,pcdv,pcdh)
[3150]	19
	20	use comcstfi_h
	21	use turb_mod, only: turb_resolved
[3153]	22	use watersat_mod, only: watersat
[3219]	23	use lmdz_atke_turbulence_ini, only : smmin, ric, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi
[3325]	24	use paleoclimate_mod, only: include_waterbuoyancy
	25	use write_output_mod, only: write_output
	26	use comslope_mod, ONLY: iflat
[3150]	27	IMPLICIT NONE
[3167]	28	include "callkeys.h"
[3325]	29
	30
[3150]	31	!=======================================================================
	32	!
	33	! Subject: computation of the surface drag coefficient using the
	34	! ------- approch developed by Loui for ECMWF.
	35	!
	36	! Author: Frederic Hourdin 15 /10 /93
	37	! Modified by : Arnaud Colaitis 03/08/11; rewritten by LL to switch to F90
	38	! -------
	39	!
	40	! Arguments:
	41	! ----------
	42	!
	43	! inputs:
	44	! ------
[3153]	45	! ngrid size of the horizontal grid
	46	! pg gravity (m s -2)
	47	! pz(ngrid,nlay) height of layers
	48	! pp(ngrid,nlay) pressure of layers
	49	! pu(ngrid,nlay) u component of the wind
	50	! pv(ngrid,nlay) v component of the wind
	51	! pts(ngrid) surface temperature
	52	! ph(ngrid) potential temperature T*(p/ps)^kappa
	53	! mumean Molecular mass of the atmosphere (kg/mol)
	54	! pqvap(ngrid,nlay) Water vapor mixing ratio (kg/kg) to account for vapor flottability
	55	! pqsurf(ngrid) Water ice frost on the surface(kg/m^2) to account for vapor flottability
[3150]	56	!
	57	! outputs:
	58	! --------
	59	! pcdv(ngrid) Cd for the wind
	60	! pcdh(ngrid) Cd for potential temperature
	61	!
	62	!=======================================================================
	63
	64
	65	!-----------------------------------------------------------------------
	66	! Declarations:
	67	! -------------
	68
	69
	70	! Arguments:
	71
	72	! Inputs:
	73	! ----------
[3165]	74	INTEGER, INTENT(IN) :: ngrid,nlay,nslope ! Number of points in the physical grid and atmospheric grid
[3150]	75	REAL, INTENT(IN) :: pz0(ngrid) ! Surface Roughness (m)
[3153]	76	REAL, INTENT(IN) :: pg ! Mars gravity (m/s^2)
	77	REAL, INTENT(IN) :: pz(ngrid,nlay),pp(ngrid,nlay) ! Layers pseudo altitudes (m) and pressure (Pa)
[3150]	78	REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) ! Zonal and Meriditionnal winds (m/s)
[3165]	79	REAL, INTENT(IN) :: pts(ngrid,nslope),ph(ngrid,nlay) ! Surface Temperature and atmospheric temperature (K)
[3150]	80	REAL, INTENT(IN) :: wstar(ngrid) ! Vertical velocity due to turbulence (m/s)
[3153]	81	REAL, INTENT(IN) :: mumean(ngrid) ! Molecular mass of the atmosphere (kg/mol)
	82	REAL, INTENT(IN) :: pqvap(ngrid,nlay) ! Water vapor mixing ratio (kg/kg) to account for vapor flottability
[3165]	83	REAL, INTENT(IN) :: pqsurf(ngrid,nslope) ! Water ice frost on the surface (kg/m^2) to account for vapor flottability
[3325]	84	LOGICAL, INTENT(IN) :: write_outputs ! write_output in xios or not.
[3150]	85
	86	! Outputs:
	87	! ----------
[3165]	88	REAL, INTENT(OUT) :: pcdv(ngrid,nslope) ! momentum drag coefficient (1)
	89	REAL, INTENT(OUT) :: pcdh(ngrid,nslope) ! heat drag coefficient (1)
[3150]	90
	91	! Local:
	92	! ------
	93
[3165]	94	INTEGER ig,islope ! Loop variable
[3150]	95
[3165]	96	REAL karman,nu ! Von Karman constant and fluid kinematic viscosity
[3150]	97
	98	LOGICAL firstcal
	99	DATA karman,nu/.41,0.001/
	100	DATA firstcal/.true./
	101	SAVE karman,nu
	102
	103	!$OMP THREADPRIVATE(karman,nu)
	104
[3325]	105	REAL rib(ngrid,nslope) ! Bulk Richardson number [virtual or dry] (1)
	106	REAL rib_dry(ngrid,nslope) ! Dry bulk Richardson number [virtual or dry] (1)
	107	REAL fm(ngrid,nslope) ! stability function for momentum (1)
	108	REAL fh(ngrid,nslope) ! stability function for heat (1)
[3150]	109
	110	! Formalism for stability functions fm= 1/phim^2; fh = 1/(phim*phih) where
	111	! phim = 1+betamzeta or (1-bmzeta)**am
	112	! phih = alphah + betahzeta or alphah(1.-bhzeta)**ah
	113	! ah and am are assumed to be -0.25 and -0.5 respectively
	114	! lambda is an intermediate variable to simplify the expression
	115	REAL betam, betah, alphah, bm, bh, lambda
	116
[3325]	117	REAL pz0t ! initial thermal roughness length z0t for the iterative scheme (m)
	118	REAL ric_colaitis ! critical richardson number (1)
	119	REAL reynolds(ngrid,nslope) ! Reynolds number (1)
	120	REAL prandtl(ngrid) ! Prandtl number (1)
	121	REAL pz0tcomp(ngrid) ! computed z0t during the iterative scheme (m)
	122	REAL z0t(ngrid,nslope) ! computed z0t at the last step the iterative scheme (m)
	123	REAL ite ! Number of iteration (1)
	124	REAL itemax ! Maximum number of iteration (1)
	125	REAL residual ! Residual between pz0t and pz0tcomp (m)
	126	REAL tol_iter ! Tolerance for the residual to ensure convergence (1=
[3150]	127
[3325]	128	REAL zu2(ngrid) ! Near surface winds (m/s)
[3150]	129
[3325]	130	REAL cdn(ngrid) ! neutral momentum drag coefficient (1)
	131	REAL chn(ngrid) ! neutral heat drag coefficient (1)
	132	REAL z1z0,z1z0t ! ratios z1/z0 and z1/z0T
	133	REAL z1,zcd0 ! Neutral roughness (m) and Cd/Ch coefficient when call richls is not called
	134	REAL tsurf_v(ngrid,nslope) ! Virtual surface temperature, accounting for vapor flottability
	135	REAL temp_v(ngrid) ! Potential virtual air temperature in the frist layer, accounting for vapor flottability
	136	REAL mu_h2o ! Molecular mass of water vapor (kg/mol)
	137	REAL tol_frost ! Tolerance to consider the effect of frost on the surface
	138	REAL qsat(ngrid) ! saturated mixing ratio of water vapor
[3153]	139
[3325]	140	REAL sm ! Stability function computed with the ATKE scheme
	141	REAL f_ri_cd_min ! minimum of the stability function with the ATKE scheme
[3167]	142
[3150]	143	! Code:
	144	! --------
	145
	146	! Initialisation
	147	itemax= 10
	148	tol_iter = 0.01
[3153]	149	mu_h2o = 18e-3
[3154]	150	tol_frost = 1e-4
[3325]	151	reynolds(:,:) = 0.
[3150]	152	pz0t = 0.
	153	pz0tcomp(:) = 0.1*pz0(:)
[3325]	154	rib(:,:) = 0.
	155	rib_dry(:,:) = 0.
[3165]	156	pcdv(:,:) = 0.
	157	pcdh(:,:) = 0.
[3325]	158	z0t(:,:) = 0.
[3333]	159	fm(:,:) = 0.
	160	fh(:,:) = 0.
[3167]	161	f_ri_cd_min = 0.01
[3150]	162	! this formulation assumes alphah=1., implying betah=betam
	163	! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
	164	bm = 16.
	165	bh = 16.
	166	alphah = 1.
	167	betam = 5.
	168	betah = 5.
	169	lambda = (sqrt(bh/bm))/alphah
	170
[3219]	171	ric_colaitis = betah/(betam**2)
[3153]	172
[3325]	173	IF(include_waterbuoyancy) then
[3165]	174	temp_v(:) = ph(:,1)*(1.+pqvap(:,1)/(mu_h2o/mumean(:)))/(1.+pqvap(:,1))
	175	DO islope = 1,nslope
	176	DO ig = 1,ngrid
	177	IF(pqsurf(ig,islope).gt.tol_frost) then
	178	call watersat(1,pts(ig,islope),pp(ig,1),qsat(ig))
	179	tsurf_v(ig,islope) = pts(ig,islope)*(1.+qsat(ig)/(mu_h2o/mumean(ig)))/(1.+qsat(ig))
	180	ELSE
	181	tsurf_v(ig,islope) = pts(ig,islope)*(1.+pqvap(ig,1)/(mu_h2o/mumean(ig)))/(1.+pqvap(ig,1))
	182	ENDIF
	183	ENDDO
[3153]	184	ENDDO
	185	ELSE
[3165]	186	tsurf_v(:,:) = pts(:,:)
[3153]	187	temp_v(:) = ph(:,1)
	188	ENDIF
[3150]	189
	190	IF(.not.callrichsl) then
	191	! Original formulation as in LMDZ Earth: Cd = Ch = (kappa/(ln(1+z1/z0)))^2
	192	! NB: In forget et al., 1999, it's Cd = Ch = (kappa/(ln(z1/z0)))^2
[3165]	193	DO ig = 1,ngrid
	194	z1 = 1.E+0 + pz(ig,1)/pz0(ig)
	195	zcd0 = karman/log(z1)
	196	zcd0 = zcd0*zcd0
	197	DO islope = 1,nslope
	198	pcdv(ig,islope) = zcd0
	199	pcdh(ig,islope) = zcd0
	200	ENDDO
[3150]	201	ENDDO
	202	ELSE
	203	! We follow the parameterization from Colaitis et al., 2013 (supplementary material)
[3165]	204	DO islope = 1,nslope
	205	DO ig=1,ngrid
	206	ite = 0.
[3333]	207	residual = 100tol_iterpz0(ig)
[3165]	208	z1z0=pz(ig,1)/pz0(ig)
	209	cdn(ig)=karman/log(z1z0)
	210	cdn(ig)=cdn(ig)*cdn(ig)
	211	DO WHILE((residual .gt. tol_iter*pz0(ig)) .and. (ite .lt. itemax))
[3150]	212	! Computations of z0T; iterated until z0T converges
[3165]	213	pz0t = pz0tcomp(ig)
	214	z1z0t=pz(ig,1)/pz0t
	215	chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t)
	216	IF ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then
	217	zu2(ig) = pu(ig,1)pu(ig,1)+pv(ig,1)pv(ig,1) + (log(1.+0.7wstar(ig) + 2.3wstar(ig)2))2 ! Near surface winds account for buoyancy
	218	IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.)
[3150]	219	! Richardson number formulation proposed by D.E. England et al. (1995)
[3325]	220	rib(ig,islope) = (pg/tsurf_v(ig,islope))sqrt(pz(ig,1)pz0(ig))(((log(pz(ig,1)/pz0(ig)))2)/(log(pz(ig,1)/pz0t)))(temp_v(ig)-tsurf_v(ig,islope))/zu2(ig)
	221	rib_dry(ig,islope) = (pg/pts(ig,islope))sqrt(pz(ig,1)pz0(ig))(((log(pz(ig,1)/pz0(ig)))2)/(log(pz(ig,1)/pz0t)))(ph(ig,1)-pts(ig,islope))/zu2(ig)
[3165]	222	ELSE
	223	print*,'warning, infinite Richardson at surface'
	224	print*,pu(ig,1),pv(ig,1)
[3325]	225	rib(ig,islope) = ric_colaitis
	226	rib_dry(ig,islope) = ric_colaitis
[3165]	227	ENDIF
[3167]	228	! Compute the stability functions fm; fh depending on the stability of the surface layer
[3150]	229
[3167]	230	IF(callatke) THEN
	231	! Computation following parameterizaiton from ATKE
[3325]	232	IF(rib(ig,islope) .gt. 0) THEN
	233	! STABLE BOUNDARY LAYER :include_waterbuoyancy
	234	sm = MAX(smmin,cn*(1.-rib(ig,islope)/ric))
[3167]	235	! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019
[3325]	236	prandtl(ig) = pr_neutexp(-pr_slope/pr_neutrib(ig,islope)+rib(ig,islope)/pr_neut) + rib(ig,islope) * pr_slope
	237	fm(ig,islope) = max((sm*(3./2.) sqrt(cepsilon) * (1 - rib(ig,islope) / prandtl(ig))**(1./2.)),f_ri_cd_min)
	238	fh(ig,islope) = max((fm(ig,islope) / prandtl(ig)), f_ri_cd_min)
[3167]	239	ELSE
	240	! UNSTABLE BOUNDARY LAYER :
[3325]	241	sm = 2./rpi * (cinf - cn) * atan(-rib(ig,islope)/ri0) + cn
	242	prandtl(ig) = -2./rpi * (pr_asym - pr_neut) * atan(rib(ig,islope)/ri1) + pr_neut
	243	fm(ig,islope) = MAX((sm*(3./2.) sqrt(cepsilon) * (1 - rib(ig,islope) / prandtl(ig))**(1./2.)),f_ri_cd_min)
	244	fh(ig,islope) = MAX((fm(ig,islope) / prandtl(ig)), f_ri_cd_min)
[3167]	245	ENDIF ! Rib < 0
	246	ELSE
	247	! Computation following parameterizaiton from from D.E. England et al. (95)
[3325]	248	IF (rib(ig,islope) .gt. 0.) THEN
[3167]	249	! STABLE BOUNDARY LAYER :
	250	prandtl(ig) = 1.
[3325]	251	IF(rib(ig,islope) .lt. ric_colaitis) THEN
[3150]	252	! Assuming alphah=1. and bh=bm for stable conditions :
[3325]	253	fm(ig,islope) = ((ric_colaitis-rib(ig,islope))/ric_colaitis)**2
	254	fh(ig,islope) = ((ric_colaitis-rib(ig,islope))/ric_colaitis)**2
[3167]	255	ELSE
	256	! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
[3325]	257	fm(ig,islope) = 1.
	258	fh(ig,islope) = 1.
[3167]	259	ENDIF
[3165]	260	ELSE
[3167]	261	! UNSTABLE BOUNDARY LAYER :
[3325]	262	fm(ig,islope) = sqrt(1.-lambdabmrib(ig,islope))
	263	fh(ig,islope) = (1./alphah)((1.-lambdabhrib(ig,islope))0.5)(1.-lambdabmrib(ig,islope))**0.25
	264	prandtl(ig) = alphah((1.-lambdabmrib(ig,islope))0.25)/((1.-lambdabhrib(ig,islope))*0.5)
[3167]	265	ENDIF ! Rib < 0
	266	ENDIF ! atke
[3150]	267	! Recompute the reynolds and z0t
[3325]	268	reynolds(ig,islope) = karmansqrt(fm(ig,islope))sqrt(zu2(ig))pz0(ig)/(log(z1z0)nu)
	269	pz0tcomp(ig) = pz0(ig)exp(-karman7.3(reynolds(ig,islope)0.25)(prandtl(ig)*0.5)+5karman)
[3165]	270	residual = abs(pz0t-pz0tcomp(ig))
	271	ite = ite+1
	272	ENDDO ! of while
[3150]	273
	274	! Compute the coefficients Cdv; Cdh : neutral coefficient x stability functions
[3325]	275	pcdv(ig,islope)=cdn(ig)*fm(ig,islope)
	276	pcdh(ig,islope)=chn(ig)*fh(ig,islope)
	277	z0t(ig,islope) = pz0tcomp(ig)
[3165]	278	pz0t = 0. ! for next grid point
	279	ENDDO ! of ngrid
	280	enddo
[3150]	281	ENDIF !of if call richardson surface layer
[3325]	282
	283	IF(write_outputs) then
[3150]	284
[3325]	285	call write_output("rib_dry_vdif_cd","Dry Richardson number in vdif_cd","1",rib_dry(:,iflat))
	286
	287	call write_output("rib_vdif_cd","Richardson number in vdif_cd","1",rib(:,iflat))
	288
	289	call write_output("fm_vdif_cd","Momentum stability function in vdif_cd","1",fm(:,iflat))
	290
	291	call write_output("fh_vdif_cd","Heat stability function in vdif_cd","1",fh(:,iflat))
	292
	293	call write_output("z0t_vdif_cd","Thermal roughness length in vdif_cd","m",z0t(:,iflat))
	294
	295	call write_output("z0_vdif_cd","Momentum roughness length in vdif_cd","m",pz0(:))
	296
	297	call write_output("Reynolds_vdif_cd","Reynolds number in vdif_cd","1",reynolds(:,iflat))
	298
	299	ENDIF
	300
[3150]	301	RETURN
	302
	303	END SUBROUTINE vdif_cd
	304
	305	END MODULE vdif_cd_mod
	306
	307
	308
	309
	310
	311
	312
	313
[3219]	314
	315
[3150]	316
	317

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

Download in other formats:

Original Format