Context Navigation

topmons_mod.F90 @ 3026

Last change on this file since 3026 was 2953, checked in by romain.vande, 20 months ago
Mars PCM : Adapt vdifc, co2condens, rocketduststorm and topmons routines to the subslope parametrisation. RV
File size: 53.1 KB

Rev	Line
[2199]	1	MODULE topmons_mod
	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
	7	!=======================================================================
	8	! ENTRAINMENT OF DUST ABOVE THE SUB-GRID SCALE TOPOGRAPHY
	9	!=======================================================================
	10	! calculation of the vertical wind velocity of topdust tracers
	11	! transport scheme of topdust tracers
	12	! detrainement of topdust into background dust
	13	! -----------------------------------------------------------------------
	14	! Authors: M. Vals; C. Wang; T. Bertrand; F. Forget
	15	! Institution: Laboratoire de Meteorologie Dynamique (LMD) Paris, France
	16	! -----------------------------------------------------------------------
	17
	18	SUBROUTINE topmons(ngrid,nlayer,nq,ptime,ptimestep, &
	19	pq,pdq,pt,pdt,pplev,pplay,pzlev, &
	20	pzlay,pdtsw,pdtlw, &
	21	! input for radiative transfer
[2826]	22	icount,zday,zls,tsurf,co2ice,igout, &
	23	aerosol,tauscaling,dust_rad_adjust, &
[2953]	24	IRtoVIScoef,albedo,emis, &
[2199]	25	! input sub-grid scale rocket dust storm
[2246]	26	totstormfract,clearatm, &
[2199]	27	! input sub-grid scale cloud
	28	clearsky,totcloudfrac, &
	29	! input sub-grid scale mountain
[2628]	30	nohmons, &
[2199]	31	! output
[2246]	32	pdqtop,wfin,dsodust,dsords,dsotop, &
[2415]	33	tau_pref_scenario,tau_pref_gcm)
[2199]	34
	35	USE tracer_mod, only: igcm_topdust_mass,igcm_topdust_number &
	36	,igcm_dust_mass,igcm_dust_number &
	37	,rho_dust
	38	USE comcstfi_h, only: r,g,cpp,rcp
[2953]	39	USE dimradmars_mod, only: naerkind
[2199]	40	USE comsaison_h, only: dist_sol,mu0,fract
[2953]	41	USE surfdat_h, only: hmons,summit,alpha_hmons, &
[2628]	42	hsummit,contains_mons
[2226]	43	USE callradite_mod, only: callradite
[2932]	44	use write_output_mod, only: write_output
[2199]	45
	46	IMPLICIT NONE
	47
	48	!--------------------------------------------------------
	49	! Input Variables
	50	!--------------------------------------------------------
	51
	52	INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points
	53	INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points
	54	INTEGER, INTENT(IN) :: nq ! number of tracer species
	55	REAL, INTENT(IN) :: ptime
	56	REAL, INTENT(IN) :: ptimestep
	57
	58	REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq
	59	REAL, INTENT(IN) :: pdq(ngrid,nlayer,nq)! tendancy field pq
	60	REAL, INTENT(IN) :: pt(ngrid,nlayer) ! temperature at mid-layer (K)
	61	REAL, INTENT(IN) :: pdt(ngrid,nlayer) ! tendancy temperature at mid-layer (K)
	62
	63	REAL, INTENT(IN) :: pplay(ngrid,nlayer) ! pressure at middle of the layers
	64	REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure at intermediate levels
	65	REAL, INTENT(IN) :: pzlay(ngrid,nlayer) ! altitude at the middle of the layers
	66	REAL, INTENT(IN) :: pzlev(ngrid,nlayer+1) ! altitude at layer boundaries
	67
	68	REAL, INTENT(IN) :: pdtsw(ngrid,nlayer) ! (K/s) env
	69	REAL, INTENT(IN) :: pdtlw(ngrid,nlayer) ! (K/s) env
	70
	71	! input for second radiative transfer
	72	INTEGER, INTENT(INOUT) :: icount
	73	REAL, INTENT(IN) :: zday
	74	REAL, INTENT(IN) :: zls
	75	REAL, INTENT(IN) :: tsurf(ngrid)
[2953]	76	REAL, INTENT(IN) :: albedo(ngrid,2)
	77	REAL, INTENT(IN) :: emis(ngrid)
[2826]	78	REAL,INTENT(IN) :: co2ice(ngrid) ! co2 ice surface layer (kg.m-2)
[2199]	79	INTEGER, INTENT(IN) :: igout
[2226]	80	REAL, INTENT(INOUT) :: tauscaling(ngrid)
[2634]	81	REAL,INTENT(INOUT) :: dust_rad_adjust(ngrid)
[2643]	82	REAL,INTENT(INOUT) :: IRtoVIScoef(ngrid) ! NB: not modified by this call to callradite,
	83	! the OUT is just here because callradite needs it
	84
[2199]	85	! input sub-grid scale rocket dust storm
	86	LOGICAL, INTENT(IN) :: clearatm
	87	REAL, INTENT(IN) :: totstormfract(ngrid)
	88	! input sub-grid scale water ice clouds
	89	LOGICAL, INTENT(IN) :: clearsky
	90	REAL, INTENT(IN) :: totcloudfrac(ngrid)
	91	! input sub-grid scale mountain
	92	LOGICAL, INTENT(IN) :: nohmons
	93
	94	!--------------------------------------------------------
	95	! Output Variables
	96	!--------------------------------------------------------
	97
	98	REAL, INTENT(OUT) :: pdqtop(ngrid,nlayer,nq) ! tendancy field for dust
	99	REAL, INTENT(OUT) :: wfin(ngrid,nlayer+1) ! final vertical wind velocity: combination of both wup and wrad
	100	! wfin < 0 means upward wind (in pressure levels/s)
	101
	102	! output for second radiative transfer
	103	REAL, INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind)
[2415]	104	REAL, INTENT(OUT) :: dsodust(ngrid,nlayer)
	105	REAL, INTENT(OUT) :: dsords(ngrid,nlayer)
	106	REAL, INTENT(OUT) :: dsotop(ngrid,nlayer)
	107	REAL,INTENT(OUT) :: tau_pref_scenario(ngrid) ! prescribed dust column
	108	! visible opacity at odpref
	109	REAL,INTENT(OUT) :: tau_pref_gcm(ngrid) ! dust column visible opacity at
	110	! odpref in the GCM
[2199]	111
	112	!--------------------------------------------------------
	113	! Local variables
	114	!--------------------------------------------------------
	115	LOGICAL,SAVE :: firstcall=.true.
[2616]	116
	117	!$OMP THREADPRIVATE(firstcall)
	118
[2199]	119	INTEGER l,ig,tsub,iq,ll
	120	REAL zq0(ngrid,nlayer,nq) ! initial tracers
	121	REAL zq(ngrid,nlayer,nq) ! updated tracers
	122	REAL,PARAMETER:: mu0lim=0.001 ! solar zenith angle limit to determine the daytime
	123
	124	! Variables for radiative transfer
	125	REAL fluxsurf_lw1(ngrid)
[2685]	126	REAL fluxsurf_dn_sw1(ngrid,2),fluxsurf_up_sw1(ngrid,2)
[2199]	127	REAL fluxtop_lw1(ngrid)
[2685]	128	REAL fluxtop_dn_sw1(ngrid,2),fluxtop_up_sw1(ngrid,2)
[2199]	129	REAL tau(ngrid,naerkind)
	130	REAL taucloudtes(ngrid)
	131	REAL rdust(ngrid,nlayer)
	132	REAL rstormdust(ngrid,nlayer)
	133	REAL rtopdust(ngrid,nlayer)
	134	REAL rice(ngrid,nlayer)
	135	REAL nuice(ngrid,nlayer)
[2459]	136	DOUBLE PRECISION riceco2(ngrid,nlayer)
[2448]	137	REAL nuiceco2(ngrid,nlayer)
[2199]	138
	139	! Temperature profile
	140	REAL zt(ngrid,nlayer) ! actual temperature at mid-layer (K)
	141	REAL ztlev(ngrid,nlayer) ! temperature at intermediate levels l+1/2 without last level
	142	REAL zdtlw1(ngrid,nlayer) ! (K/s) topdust
	143	REAL zdtsw1(ngrid,nlayer) ! (K/s) topdust
	144	REAL zdtlw1_lev(ngrid,nlayer),zdtsw1_lev(ngrid,nlayer) ! rad. heating rate at intermediate levels l+1/2
	145	REAL zdtlw_lev(ngrid,nlayer),zdtsw_lev(ngrid,nlayer) ! rad. heating rate at intermediate levels l+1/2
	146	REAL zdtvert(ngrid,nlayer) ! dT/dz , lapse rate
	147	REAL deltahr(ngrid,nlayer+1) ! extra heating between the concentrated and non-concentrated dust
	148
	149	REAL newzt(ngrid,nlayer) ! recalculated temperature with extra radiative heating
	150	REAL t_top(ngrid,nlayer) ! estimated temperature above the mountain
	151	REAL t_top_lev(ngrid,nlayer)! estimated temperature above the mountain at the levels
	152	REAL dt_top(ngrid) ! temperature difference between the summit and the environment
	153	INTEGER lmons(ngrid) ! layer of the mountain's slope
	154	INTEGER lsummit(ngrid) ! layer of the mountain's summit
	155	REAL t_env(ngrid,nlayer) ! recalculated temperature of the environment air "next" to the mountain
	156
	157	! Vertical velocity profile
	158	REAL wrad(ngrid,nlayer+1) ! vertical wind velocity resulting from radiative heating
	159	REAL wup(ngrid,nlayer+1) ! vertical wind velocity calculated above the sub-grid scale mountain
	160	REAL wup2(ngrid,nlayer+1) ! square of the vertical velocity calculated above the sub-grid scale mountain
	161	REAL w0(ngrid) ! vertical wind velocity convergence at the top of the sub-grid scale mountain
	162	INTEGER lwmax(ngrid) ! level of the maximum vertical velocity
	163
	164	REAL,PARAMETER:: wmax=10. ! maximum vertical velocity resulting from radiative heating (m/s)
	165	REAL,PARAMETER:: secu=3. ! coefficient on wfin to avoid dust crossing many layers during subtimestep
	166	REAL,PARAMETER:: k0=0.25 !max 10m/s: k0=1. !max 3m/s: k0=0.25
	167	REAL,PARAMETER:: k1=5.e-4 !max 10m/s: k1=5.e-4 !max 3m/s: k1=5.e-4
	168	REAL,PARAMETER:: k2=5.e-3 !max 10m/s: k2=30.e-3 !max 3m/s: k2=5.e-3
	169
	170	! Entrainment
	171	REAL entr(ngrid) ! entrainment flux
[2255]	172	REAL rho(ngrid,nlayer),rhobarz(ngrid,nlayer+1) ! density, density at the levels
[2199]	173	REAL masse(ngrid,nlayer) ! air mass
	174	REAL masse_pbl(ngrid) ! total air mass within the PBL
	175	REAL masse_pbl_dust_mass(ngrid),masse_pbl_dust_number(ngrid) ! total air mass + dust mass/number within the PBL (kg/m2)
	176	REAL dqm_mass_pbl(ngrid),dqm_number_pbl(ngrid) ! total mass of the entrained dust mass/number from the PBL (kg/m2)
	177
	178	! Van Leer
	179	REAL zq_dust_mass(ngrid,nlayer),zq_dust_number(ngrid,nlayer) ! mixing ratio background dust (kg/kg)
	180	REAL zq_topdust_mass(ngrid,nlayer),zq_topdust_number(ngrid,nlayer) ! mixing ratio topdust (kg/kg)
	181	REAL mr_dust_mass(ngrid,nlayer),mr_dust_number(ngrid,nlayer) ! mixing ratio background dust (kg/kg)
	182	REAL mr_topdust_mass(ngrid,nlayer),mr_topdust_number(ngrid,nlayer) ! mixing ratio background dust + concentrated topdust (kg/kg)
	183
	184	REAL w(ngrid,nlayer) ! vertical flux above the sub-grid scale mountain
	185	REAL wqmass(ngrid,nlayer+1) ! flux mass
	186	REAL wqnumber(ngrid,nlayer+1) ! flux number
	187	REAL qbar_mass_pbl(ngrid) ! total dust mass mixing ratio within the PBL for Van Leer
	188	REAL qbar_number_pbl(ngrid) ! total dust number mixing ratio within the PBL for Van Leer
	189	REAL masse_col(nlayer) ! masse of the atmosphere in one column
	190
	191	REAL dqvl_dust_mass(ngrid,nlayer) ! tendancy pdq dust mass after vertical transport
	192	REAL dqvl_dust_number(ngrid,nlayer) ! tendancy pdq dust number after vertical transport
	193	REAL dqvl_topdust_mass(ngrid,nlayer) ! tendancy pdq topdust mass after vertical transport
	194	REAL dqvl_topdust_number(ngrid,nlayer) ! tendancy pdq topdust number after vertical transport
	195
	196	INTEGER nsubtimestep(ngrid) ! number of subtimestep when calling Van Leer scheme
	197	REAL subtimestep(ngrid) ! ptimestep/nsubtimestep
	198	REAL dtmax ! considered minimum time interval needed for the dust to cross one vertical layer
	199
	200	! Detrainment
[2628]	201	REAL coefdetrain(ngrid,nlayer) ! coefficient for detrainment : % of topdust detrained
[2199]	202	REAL dqdet_topdust_mass(ngrid,nlayer) ! tendancy pdq topdust mass after detrainment only
	203	REAL dqdet_topdust_number(ngrid,nlayer) ! tendancy pdq topdust number after detrainment only
	204
	205	! Diagnostics
	206	REAL zlaywmax(ngrid)
	207	REAL detrain_rate(ngrid,nlayer)
	208
	209	! **********************************************************************
	210	! **********************************************************************
[2628]	211	! Parametrization of the entrainment of dust by slope winds above the
	212	! converging sub-grid scale topography ("mountain top dust flows")
[2199]	213	! **********************************************************************
	214	! **********************************************************************
	215	!! 1. Radiative transfer in topdust
	216	!! 2. Compute vertical velocity for topdust
	217	!! 3. Vertical transport
	218	!! 4. Detrainment
	219	! **********************************************************************
	220	! **********************************************************************
	221
	222	! **********************************************************************
	223	! 0. Initializations
	224	! **********************************************************************
	225	aerosol(:,:,:)=0.
	226	pdqtop(:,:,:) = 0.
	227	dqvl_topdust_mass(:,:)=0.
	228	dqvl_topdust_number(:,:)=0.
	229	dqvl_dust_mass(:,:)=0.
	230	dqvl_dust_number(:,:)=0.
	231	dqdet_topdust_mass(:,:)=0.
	232	dqdet_topdust_number(:,:)=0.
	233	wup(:,:)=0.
	234	wup2(:,:)=0.
	235	wrad(:,:)=0.
	236	w0(:)=0.
	237	wfin(:,:)=0.
	238	w(:,:)=0.
	239	zq(:,:,:) = 0.
	240	zq0(:,:,:) = 0.
	241	entr(:) = 0.
	242	mr_dust_mass(:,:) = 0.
	243	mr_dust_number(:,:) = 0.
	244	mr_topdust_mass(:,:) = 0.
	245	mr_topdust_number(:,:) = 0.
	246	t_top(:,:) = 0.
	247	dt_top(:) = 0.
	248	newzt(:,:) = 0.
	249	lmons(:) = 1
	250	lsummit(:) =1
	251	lwmax(:) = 1
	252	deltahr(:,:) = 0.
	253	zdtvert(:,:) = 0.
	254	wqmass(:,:) = 0.
	255	wqnumber(:,:) = 0.
	256	coefdetrain(:,:) = 1.
	257	dqm_mass_pbl(:)=0.
	258	dqm_number_pbl(:)=0.
	259	nsubtimestep(:)=1
	260	masse_pbl(:)=0.
	261	detrain_rate(:,:) = 0.
[2825]	262	zlaywmax(:)=-1.
[2199]	263
	264	!! Update the initial temperature
	265	zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer)+pdt(1:ngrid,1:nlayer)*ptimestep
	266	newzt(1:ngrid,1:nlayer)=zt(1:ngrid,1:nlayer)
	267	t_env(1:ngrid,1:nlayer)=zt(1:ngrid,1:nlayer)
	268	t_top(1:ngrid,1:nlayer)=zt(1:ngrid,1:nlayer)
	269
	270	!! Update the initial mixing ratios
[2248]	271	zq0(1:ngrid,1:nlayer,1:nq)=pq(1:ngrid,1:nlayer,1:nq)+pdq(1:ngrid,1:nlayer,1:nq)*ptimestep ! update of the background dust after rocket dust storm scheme
[2199]	272	zq(1:ngrid,1:nlayer,1:nq)=zq0(1:ngrid,1:nlayer,1:nq)
	273	zq_dust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_mass)
	274	zq_dust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_number)
	275	zq_topdust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_topdust_mass)
	276	zq_topdust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_topdust_number)
	277
	278	! ----------------------------------------------------------------------
	279	! 1. Radiative heating
	280	! ----------------------------------------------------------------------
	281
	282	! *********************************************************************
	283	! 1.1. Call the second radiative transfer for topdust, obtain the extra heating
	284	! *********************************************************************
[2628]	285
	286	! NB: since the only grid meshes that matter for the radiative transfer
	287	! are the ones that contain a mount (contains_mons(ig)=.true.),
	288	! it could be relevant to optimize the code by calling the RT
	289	! only for those grid meshes.
	290	! See Ticket #92 on trac.lmd.jussieu.fr/Planeto
[2199]	291	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,zq,albedo, &
	292	emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, &
[2685]	293	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_dn_sw1,fluxsurf_up_sw1, &
	294	fluxtop_lw1,fluxtop_dn_sw1,fluxtop_up_sw1, &
	295	tau_pref_scenario,tau_pref_gcm, &
[2643]	296	tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef, &
[2448]	297	taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,rstormdust,rtopdust, &
[2628]	298	totstormfract,clearatm,dsords,dsotop,nohmons,&
[2199]	299	clearsky,totcloudfrac)
	300	! **********************************************************************
	301	! 1.2. Compute radiative vertical velocity for entrained topdust
	302	! **********************************************************************
	303	DO ig=1,ngrid
[2628]	304	IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) ) THEN
	305	!! mu0(ig)>mu0lim ensures that it is daytime
	306	!! contains_mons=True ensures that there is a mount in the mesh and alpha_hmons>0
	307
[2199]	308	!! **********************************************************************
	309	!! Temperature profile above the mountain and in the close environment
	310	call t_topmons(nlayer,summit(ig),hsummit(ig),hmons(ig),zt(ig,:),pzlay(ig,:), &
	311	t_top(ig,:),dt_top(ig),t_env(ig,:),lsummit(ig),lmons(ig))
	312	!! **********************************************************************
	313	!! Calculation of the extra heating : computing heating rates
	314	!! gradient at boundaries of each layer
	315	zdtlw1_lev(ig,1)=0.
	316	zdtsw1_lev(ig,1)=0.
	317	zdtlw_lev(ig,1)=0.
	318	zdtsw_lev(ig,1)=0.
	319	ztlev(ig,1)=zt(ig,1)
[2226]	320	t_top_lev(ig,1)=tsurf(ig)
	321
[2199]	322	DO l=1,nlayer-1
	323	!! Calculation for the background dust AND the concentrated topdust
	324	zdtlw1_lev(ig,l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
	325	zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
	326	(pzlay(ig,l+1)-pzlay(ig,l))
	327
	328	zdtsw1_lev(ig,l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
	329	zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
	330	(pzlay(ig,l+1)-pzlay(ig,l))
	331	!! Calculation for the background dust only
	332	zdtlw_lev(ig,l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
	333	pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
	334	(pzlay(ig,l+1)-pzlay(ig,l))
	335
	336	zdtsw_lev(ig,l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
	337	pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
	338	(pzlay(ig,l+1)-pzlay(ig,l))
	339	!! Temperature at the levels
	340	ztlev(ig,l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
	341	zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
	342	(pzlay(ig,l+1)-pzlay(ig,l))
	343	!! Temperature above the mountain at the levels
	344	t_top_lev(ig,l+1)=(t_top(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
	345	t_top(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
	346	(pzlay(ig,l+1)-pzlay(ig,l))
	347	ENDDO ! DO l=1,nlayer-1
	348	!! **********************************************************************
	349	!! Vertical velocity profile for extra heating balanced by adiabatic cooling
	350	!! Gradient at boundaries of each layer, start from surface
	351	zdtvert(ig,1)=0.
	352	DO l=1,nlayer-1
	353	zdtvert(ig,l+1)=(t_top_lev(ig,l+1)-t_top_lev(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
	354	ENDDO
	355	!! Extra heating balanced by adiabatic cooling
	356	DO l=1,nlayer
	357	deltahr(ig,l)=(zdtlw1_lev(ig,l)+zdtsw1_lev(ig,l)) &
	358	-(zdtlw_lev(ig,l)+zdtsw_lev(ig,l))
	359	wrad(ig,l)=-deltahr(ig,l)/(g/cpp+ &
	360	max(zdtvert(ig,l),-0.99*g/cpp))
	361	!! Limit of the vertical wind velocity in case of lapse rate close to adiabatic
	362	wrad(ig,l)=max(wrad(ig,l),-wmax)
	363	wrad(ig,l)=min(wrad(ig,l),wmax)
	364	!! ATTENTION !! to simplify here, no downward velocity calculated
	365	if (wrad(ig,l).gt.0.) then
	366	wrad(ig,l)=0.
	367	endif
	368	ENDDO
[2628]	369	ENDIF ! IF ((mu0(ig) .gt. mu0lim) .and. contains_mons(ig))
[2199]	370	ENDDO ! DO ig=1,ngrid
	371
	372	! ----------------------------------------------------------------------
	373	! 2. Vertical velocity profile
	374	! ----------------------------------------------------------------------
	375
	376	! **********************************************************************
	377	! 2.1 Compute total vertical velocity for entrained topdust: buoyancy + radiative
	378	! **********************************************************************
	379	DO ig=1,ngrid
[2628]	380	IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) ) THEN
[2199]	381	!! Positive buoyancy: negative vertical velocity entrains UP
	382	IF (dt_top(ig) .gt. 0.) THEN
	383	!! Convergence of vertical winds at the top of the mountain
	384	w0(ig)=-k0gsqrt((dt_top(ig)/t_env(ig,lsummit(ig))))
	385	!! Case the vertical velocity at the top of the mountain > wrad
	386	IF (w0(ig).lt.wrad(ig,lsummit(ig)+1)) then
	387	wup(ig,lsummit(ig)+1)=w0(ig)
	388	wfin(ig,lsummit(ig)+1)=w0(ig)
	389	newzt(ig,lsummit(ig))= t_top(ig,lsummit(ig))
	390	!! Temperature profile from the top of the summit to the top of the atmosphere
	391	DO l=lsummit(ig)+1,nlayer-1
	392	!ATTENTION: test without wrad, T is only decreasing because of the adiabatic gradient
	393	!newzt(ig,l)=newzt(ig,l-1)-(g/cpp)* &
	394	! (pzlay(ig,l)-pzlay(ig,l-1))
	395	!! Case of superadiabatic layer
	396	IF (zdtvert(ig,l) .lt. -g/cpp) then
	397	newzt(ig,l)=t_top(ig,l)
	398	!! New temperature due to radiative heating balanced by adiabatic cooling
	399	ELSE
	400	newzt(ig,l)=newzt(ig,l-1)+ &
[2226]	401	( (deltahr(ig,l)/(-wfin(ig,l)))-g/cpp )* &
[2199]	402	( pzlay(ig,l)-pzlay(ig,l-1) )
	403	ENDIF ! (zdtvert(ig,l) .lt. -g/cpp)
	404	!! Vertical velocity profile due to the presence of the mountain with the new temperature
	405	!! Square of the vertical velocity profile with the new temperature
	406	wup2(ig,l+1)=(1.-2k1(pzlev(ig,l+1)-pzlev(ig,l)))*&
	407	(wup(ig,l)*2)+2.k2g(pzlev(ig,l+1) &
	408	-pzlev(ig,l))*((newzt(ig,l) &
	409	-t_env(ig,l))/t_env(ig,l))
	410	!! Vertical velocity profile if W2 > 0
	411	IF (wup2(ig,l+1) .gt. 0.) THEN
	412	wup(ig,l+1)=-sqrt(wup2(ig,l+1))
	413	wfin(ig,l+1)=wup(ig,l+1)
	414	IF (wup(ig,l+1) .gt. wrad(ig,l+1)) then
	415	DO ll=l,nlayer-1
	416	newzt(ig,ll)=zt(ig,ll)
	417	wfin(ig,ll+1)=wrad(ig,ll+1)
	418	ENDDO
	419	EXIT
	420	ENDIF
	421	!! Vertical velocity profile if W2 < 0
	422	ELSE
	423	DO ll=l,nlayer-1
	424	newzt(ig,ll)=zt(ig,ll)
	425	wfin(ig,ll+1)=wrad(ig,ll+1)
	426	ENDDO
	427	EXIT
	428	ENDIF ! IF (wup2(ig,l+1) .gt. 0.)
	429	ENDDO ! DO l=lsummit(ig)+1,nlayer-1
	430	!! Case the vertical velocity at the top of the mountain < wrad
	431	ELSE
	432	DO ll=lsummit(ig),nlayer-1
	433	newzt(ig,ll)=zt(ig,ll)
	434	wfin(ig,ll+1)=wrad(ig,ll+1)
	435	ENDDO
	436	ENDIF !(w0(ig).lt.wrad(ig,lsummit(ig)+1))
	437	!! Positive buoyancy: positive vertical velocity entrains DOWN
	438	ELSE
	439	DO l=lsummit(ig)+1,nlayer
	440	wfin(ig,l)=wrad(ig,l)
	441	ENDDO
	442	ENDIF ! IF (dt_top(ig) .gt. 0.)
	443	! **********************************************************************
	444	! 2.2 Find the level lwmax of the maximum vertical velocity caused by the
	445	! mountain presence (wup), where to entrain the dust from the PBL
	446	! **********************************************************************
	447	IF (minloc(wup(ig,:),dim=1).lt.nlayer) THEN
	448	lwmax(ig)=minloc(wup(ig,:),dim=1)
	449	ENDIF
	450	zlaywmax(ig)=pzlay(ig,lwmax(ig)) ! wup(lwmax) acts on level lwmax, going to layer lwmax
	451	! **********************************************************************
	452	! 2.3 Compute the subtimestep to conserve the mass in the Van Leer transport
	453	! **********************************************************************
	454	!! Calculation of the subtimesteps
	455	dtmax=ptimestep
	456	DO l=lwmax(ig),nlayer
	457	IF (wfin(ig,l).lt.0.) THEN
	458	dtmax=min(dtmax,(pzlev(ig,l)-pzlev(ig,l-1))/ &
	459	(secu*abs(wfin(ig,l))))
	460	ELSE IF (wfin(ig,l).gt.0.) then
	461	dtmax=min(dtmax,(pzlev(ig,l+1)-pzlev(ig,l))/ &
	462	(secu*abs(wfin(ig,l))))
	463	ENDIF
	464	ENDDO
	465	nsubtimestep(ig)= int(ptimestep/dtmax)
	466	subtimestep(ig)=ptimestep/float(nsubtimestep(ig))
	467	!! Mass flux generated by wup in kg/m2
	468	DO l=1,nlayer!+1
	469	w(ig,l)=wfin(ig,l)pplev(ig,l)/(rztlev(ig,l)) &
	470	*subtimestep(ig)
	471	ENDDO ! l=1,nlayer
	472
	473	ENDIF ! ((mu0(ig) .gt. mu0lim)...
	474	ENDDO ! DO ig=1,ngrid
	475
	476	! ----------------------------------------------------------------------
	477	! 3. Dust entrainment by slope winds above the mountain
	478	! ----------------------------------------------------------------------
	479	!! rho on the layer
	480	rho(:,:)=pplay(:,:)/(r*pt(:,:))
	481	!! rhobarz(l) = rho(l+1/2): rho on the levels
	482	rhobarz(:,1)=rho(:,1)
	483	DO l=2,nlayer
	484	rhobarz(:,l)=pplev(:,l)/(r0.5(pt(:,l)+pt(:,l-1)))
	485	ENDDO
[2255]	486	rhobarz(:,nlayer+1)=0. !top of the model
[2199]	487	!! Mass computation
	488	DO l=1,nlayer
	489	masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
	490	ENDDO
	491
	492	DO ig=1,ngrid
[2628]	493	IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) ) THEN
[2199]	494	!! Total air mass within the PBL before entrainment (=> by PBL we mean between the surface and the layer where the vertical wind is maximum)
	495	masse_pbl(ig)=0.
	496	DO l=1,lwmax(ig)-1
	497	masse_pbl(ig)=masse_pbl(ig)+(pplev(ig,l)-pplev(ig,l+1))/g
	498	ENDDO
	499	ENDIF ! ((mu0(ig) .gt. mu0lim)...
	500	ENDDO ! ig=1,ngrid
	501	! **********************************************************************
	502	! 3.1. Entrainment
	503	! **********************************************************************
	504	DO ig=1,ngrid
[2628]	505	IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) .and. (masse_pbl(ig) .gt. 0.) ) THEN
[2199]	506	!! Transport of background dust + concentrated topdust above lwmax
	507	DO l=lwmax(ig),nlayer
	508	mr_dust_mass(ig,l) = zq_dust_mass(ig,l)
	509	mr_dust_number(ig,l) = zq_dust_number(ig,l)
	510	mr_topdust_mass(ig,l) = zq_dust_mass(ig,l) &
	511	+zq_topdust_mass(ig,l)/alpha_hmons(ig)
	512	mr_topdust_number(ig,l) = zq_dust_number(ig,l) &
	513	+zq_topdust_number(ig,l)/alpha_hmons(ig)
	514	ENDDO ! l=lwmax(ig),nlayer
	515	!! Entrainment flux from the PBL
	516	entr(ig) = alpha_hmons(ig)wup(ig,lwmax(ig))rhobarz(ig,lwmax(ig))/masse_pbl(ig)
	517	!! If entrains more than available masse_pbl (abs(entr) x mass_pbl x ptimestep > masse_pbl)
	518	if (abs(entr(ig)) .gt. 1./ptimestep) then
	519	entr(ig) = -(1./ptimestep)
	520	end if
	521	! **********************************************************************
	522	! 3.2. Vertical transport: Van Leer
	523	! **********************************************************************
	524	!! Loop over the subtimesteps
	525	DO tsub=1,nsubtimestep(ig)
	526	!! qbar_pbl: total mixing ratio within the PBL
	527	masse_pbl_dust_mass(ig)=0.
	528	masse_pbl_dust_number(ig)=0.
	529	DO l=1,lwmax(ig)-1
	530	masse_pbl_dust_mass(ig)=masse_pbl_dust_mass(ig)+zq_dust_mass(ig,l)*(pplev(ig,l)-pplev(ig,l+1))/g
	531	masse_pbl_dust_number(ig)=masse_pbl_dust_number(ig)+zq_dust_number(ig,l)*(pplev(ig,l)-pplev(ig,l+1))/g
	532	ENDDO
	533	qbar_mass_pbl(ig)=masse_pbl_dust_mass(ig)/masse_pbl(ig)
	534	qbar_number_pbl(ig)=masse_pbl_dust_mass(ig)/masse_pbl(ig)
	535	!! integration of the equation dq/dt = -(entr)q, with entr constant --> w0 < 0 when up so the equation is dq/dt = (entr)q
	536	dqm_mass_pbl(:)=0.
	537	dqm_number_pbl(:)=0.
	538	DO l=1,lwmax(ig)-1 ! this time we take the dust from the surface up to the level, where the velocity is maximum
	539	!! Total dust entrained from the PBL:
	540	dqm_mass_pbl(ig)=dqm_mass_pbl(ig)+(1.-exp(entr(ig)subtimestep(ig)))zq_dust_mass(ig,l)*(pplev(ig,l)-pplev(ig,l+1))/g
	541	dqm_number_pbl(ig)=dqm_number_pbl(ig)+(1.-exp(entr(ig)subtimestep(ig)))zq_dust_number(ig,l)*(pplev(ig,l)-pplev(ig,l+1))/g
	542	!! Background dust after entrainment (explicit: qdust(ig,l)=qdust(ig,l)-entr(ig)qdust(ig,l)ptimestep)
	543	zq_dust_mass(ig,l)=zq_dust_mass(ig,l)-(1.-exp(entr(ig)subtimestep(ig)))zq_dust_mass(ig,l)
	544	zq_dust_number(ig,l)=zq_dust_number(ig,l)-(1.-exp(entr(ig)subtimestep(ig)))zq_dust_number(ig,l)
	545	ENDDO
	546	!! Van Leer scheme (from lwmax to nlayer)
	547	wqmass(ig,:)=0.
	548	wqnumber(ig,:)=0.
	549	CALL van_leer(nlayer,mr_topdust_mass(ig,:),2.,lwmax(ig), &
	550	masse(ig,:),w(ig,:),wqmass(ig,:),masse_pbl(ig),dqm_mass_pbl(ig),alpha_hmons(ig),qbar_mass_pbl(ig))
	551	CALL van_leer(nlayer,mr_topdust_number(ig,:),2.,lwmax(ig), &
	552	masse(ig,:),w(ig,:),wqnumber(ig,:),masse_pbl(ig),dqm_number_pbl(ig),alpha_hmons(ig),qbar_number_pbl(ig))
	553	ENDDO !tsub=...
	554	! **********************************************************************
	555	! 3.3. Re-calculation of the mixing ratios after vertical transport
	556	! **********************************************************************
	557	!! Redistribution from lwmax to nlayer
	558	DO l=lwmax(ig),nlayer
	559	!! General and "healthy" case
	560	IF (mr_topdust_mass(ig,l).ge.mr_dust_mass(ig,l)) THEN
	561	zq_dust_mass(ig,l) = mr_dust_mass(ig,l)
	562	zq_dust_number(ig,l) = mr_dust_number(ig,l)
	563	zq_topdust_mass(ig,l) = alpha_hmons(ig)*(mr_topdust_mass(ig,l)-mr_dust_mass(ig,l))
	564	zq_topdust_number(ig,l) = alpha_hmons(ig)*(mr_topdust_number(ig,l)-mr_dust_number(ig,l))
	565	!! Particular case: there is less than initial dust mixing ratio after the vertical transport
	566	ELSE
	567	zq_dust_mass(ig,l) = (1.-alpha_hmons(ig))mr_dust_mass(ig,l)+alpha_hmons(ig)mr_topdust_mass(ig,l)
	568	zq_dust_number(ig,l) = (1.-alpha_hmons(ig))mr_dust_number(ig,l)+alpha_hmons(ig)mr_topdust_number(ig,l)
	569	zq_topdust_mass(ig,l) = 0.
	570	zq_topdust_number(ig,l) = 0.
	571	ENDIF
	572	ENDDO ! DO l=lwmax(ig),nlayer
	573	! **********************************************************************
	574	! 3.4. Calculation of the tendencies of the vertical transport from the surface up to the top of the atmosphere
	575	! **********************************************************************
	576	DO l=1,nlayer
	577	dqvl_topdust_mass(ig,l) = (zq_topdust_mass(ig,l)- &
	578	zq0(ig,l,igcm_topdust_mass)) /ptimestep
	579	dqvl_topdust_number(ig,l) = (zq_topdust_number(ig,l)- &
	580	zq0(ig,l,igcm_topdust_number)) /ptimestep
	581	dqvl_dust_mass(ig,l) = (zq_dust_mass(ig,l)-zq0(ig,l,igcm_dust_mass)) /ptimestep
	582	dqvl_dust_number(ig,l) = (zq_dust_number(ig,l)-zq0(ig,l,igcm_dust_number)) /ptimestep
	583	ENDDO
	584
	585	ENDIF ! ((mu0(ig) .gt. mu0lim)...
	586	ENDDO ! ig=1,ngrid
	587
	588	! ----------------------------------------------------------------------
	589	! 4. Detrainment: topdust is converted to background dust
	590	! ----------------------------------------------------------------------
	591
	592	!! **********************************************************************
	593	!! 4.1 Compute the coefficient of detrainment
	594	!! **********************************************************************
	595	DO ig=1,ngrid
	596	!DO l=lwmax(ig),nlayer-1
	597	DO l=1,nlayer!-1
	598	!! Detrainment during the day
[2628]	599	IF ( (mu0(ig) .gt. mu0lim) .and. (zq_topdust_mass(ig,l) .gt. zq_dust_mass(ig,l)*0.01) .and. (contains_mons(ig)) ) THEN
[2199]	600	coefdetrain(ig,l)=1.( rhobarz(ig,l+1)abs(wfin(ig,l+1)) - rhobarz(ig,l)*abs(wfin(ig,l)) ) / masse(ig,l)
	601	!! Detrainment when abs(w(l)) > abs(w(l+1)), i.e. coefdetrain < 0
	602	if ( coefdetrain(ig,l).lt.0.) then
	603	dqdet_topdust_mass(ig,l)=-(1.-exp(coefdetrain(ig,l)ptimestep))zq_topdust_mass(ig,l)/ptimestep
	604	dqdet_topdust_number(ig,l)=-(1.-exp(coefdetrain(ig,l)ptimestep))zq_topdust_number(ig,l)/ptimestep
	605	! for diagnostics
	606	detrain_rate(ig,l)=(1.-exp(coefdetrain(ig,l)*ptimestep))
	607	!! One cannot detrain more than available topdust
	608	if ( (abs(dqdet_topdust_mass(ig,l)).gt.zq_topdust_mass(ig,l)/ptimestep) .or. (abs(dqdet_topdust_number(ig,l)).gt.zq_topdust_number(ig,l)/ptimestep) ) then
	609	dqdet_topdust_mass(ig,l)=-zq_topdust_mass(ig,l)/ptimestep
	610	dqdet_topdust_number(ig,l)=-zq_topdust_number(ig,l)/ptimestep
	611	detrain_rate(ig,l)=1.
	612	endif
	613	!else!entrainment
	614	! dqdet_topdust_mass(ig,l)=-(1.-exp(coefdetrain(ig,l)ptimestep))zq_dust_mass(ig,l)/ptimestep
	615	! dqdet_topdust_number(ig,l)=-(1.-exp(coefdetrain(ig,l)ptimestep))zq_dust_number(ig,l)/ptimestep
	616	endif
[2628]	617	!! Full detrainment imposed during the night or when topdust leaves its origin column (contains_mons=False)
[2199]	618	ELSE
	619	dqdet_topdust_mass(ig,l)=-zq_topdust_mass(ig,l)/ptimestep
	620	dqdet_topdust_number(ig,l)=-zq_topdust_number(ig,l)/ptimestep
	621	detrain_rate(ig,l)=1.
	622	ENDIF
	623	ENDDO ! DO l=1,nlayer
	624	ENDDO ! DO ig=1,ngrid
	625
	626	! ----------------------------------------------------------------------
	627	! 5. Final tendencies: vertical transport + detrainment
	628	! ----------------------------------------------------------------------
	629
	630	DO ig=1,ngrid
	631	DO l=1,nlayer
	632	pdqtop(ig,l,igcm_topdust_mass)=dqdet_topdust_mass(ig,l) &
	633	+dqvl_topdust_mass(ig,l)
	634	pdqtop(ig,l,igcm_topdust_number)=dqdet_topdust_number(ig,l) &
	635	+dqvl_topdust_number(ig,l)
	636	pdqtop(ig,l,igcm_dust_mass)= -dqdet_topdust_mass(ig,l) &
	637	+dqvl_dust_mass(ig,l)
	638	pdqtop(ig,l,igcm_dust_number)= -dqdet_topdust_number(ig,l) &
	639	+dqvl_dust_number(ig,l)
	640	ENDDO ! DO l=1,nlayer
	641	ENDDO ! DO ig=1,ngrid
	642
	643
	644	! **********************************************************************
	645	! WRITEDIAGFI
	646	! **********************************************************************
[2934]	647	! CALL write_output('wup_top', &
	648	! 'wup_top','',wup(:,:))
	649	! CALL write_output('wlwmax_top', &
	650	! 'wlwmax_top','',wup(:,lwmax(1)))
	651	! CALL write_output('w0_top', &
	652	! 'w0_top','',wup(:,lsummit(1)+1))
	653	! CALL write_output('wup2_top', &
	654	! 'wup2_top','',wup2(:,:))
	655	! CALL write_output('wrad_top', &
	656	! 'wrad_top','',wrad(:,:))
	657	! CALL write_output('wfin_top', &
	658	! 'wfin_top','',wfin(:,:))
	659	! CALL write_output('hmons', &
	660	! 'hmons','',hmons(:))
	661	! CALL write_output('hsummit', &
	662	! 'hsummit','',hsummit(:))
	663	! CALL write_output('alpha_hmons', &
	664	! 'alpha_hmons','',alpha_hmons(:))
	665	! CALL write_output('zt_top', &
	666	! 'zt_top','',t_top(:,:))
	667	! CALL write_output('dt_top', &
	668	! 'dt_top','',dt_top(:))
	669	! CALL write_output('ztemp', &
	670	! 'envtemp','',zt(:,:))
	671	! CALL write_output('zt_env', &
	672	! 't_env','',t_env(:,:))
	673	! CALL write_output('zdtvert_top', &
	674	! 'zdtvert_top','',zdtvert(:,:))
	675	! CALL write_output('newzt_top', &
	676	! 'newzt_top','',newzt(:,:))
	677	! CALL write_output('deltahr_top', &
	678	! 'deltahr_top','',deltahr(:,:))
	679	! CALL write_output('rhobarz', &
	680	! 'rhobarz','',rhobarz(:,:))
	681	! CALL write_output('zlwmax', &
	682	! 'zlwmax','',zlaywmax(:))
	683	! CALL write_output('coefdetrain', &
	684	! 'coefdetrain','',coefdetrain(:,:))
	685	! CALL write_output('detrain_rate', &
	686	! 'detrain_rate', &
	687	! '',detrain_rate(:,:))
	688	! CALL write_output('rholev', &
	689	! 'rholev','',rho(:,:))
	690	! CALL write_output('zlsummit', &
	691	! 'zlsummit','',pzlay(:,lsummit(1)))
	692	! CALL write_output('pzlev_top', &
	693	! 'pzlev_top','',pzlev(:,:))
	694	! CALL write_output('entr', &
	695	! 'entr','',entr(:))
[2199]	696
	697	END SUBROUTINE topmons
	698
	699	!=======================================================================
	700
	701	! **********************************************************************
	702	! Subroutine to determine both temperature profiles above and near
	703	! a sub-grid scale mountain
	704	!***********************************************************************
	705	SUBROUTINE t_topmons(nlayer,summit,hsummit,hmons,zt,zlay,t_top,dt_top,t_env,lsummit,lmons)
	706
	707	USE comcstfi_h, only: g,cpp
	708
	709	IMPLICIT NONE
	710
	711	!--------------------------------------------------------
	712	! Input Variables
	713	!--------------------------------------------------------
	714	integer,intent(in) :: nlayer
	715	real,intent(in) :: summit ! Height of the mountain summit
	716	real,intent(in) :: hmons ! Height of the mountain slope
	717	real,intent(in) :: hsummit ! Height of the mountain summit from the GCM surface
	718	real,intent(in) :: zt(nlayer) ! large scale temperature profile
	719	real,intent(in) :: zlay(nlayer) ! height at the middle of each layer
	720	!--------------------------------------------------------
	721	! Output Variables
	722	!--------------------------------------------------------
	723	real,intent(out) :: t_top(nlayer) ! temperature above the mountain
	724	real,intent(out) :: dt_top ! temperature difference between the slope and
	725	! the environment
	726	real,intent(out) :: t_env(nlayer) ! temperature of the environment next to the mountain
	727	integer,intent(out) :: lsummit ! layer reached by the summit height in the GCM
	728	integer,intent(out) :: lmons ! layer of the real temperature to be considered near to the mountain top
	729
	730	!--------------------------------------------------------
	731	! Local Variables
	732	!--------------------------------------------------------
	733	integer l,ll
	734	real zlmons
	735
	736	! **********************************************************************
	737	t_env(:)=zt(:)
	738	t_top(:)=zt(:)
	739	dt_top=0.
	740	lsummit=1
	741	lmons=1
	742
	743	!! Layer where the dust is injected
	744	do while (hsummit .ge. zlay(lsummit))
	745	!! highest layer reached by the mountain (we don't interpolate and define zlay(lsummit) as the top of the mountain)
	746	lsummit=lsummit+1
	747	enddo
	748	!! temperature above the mountain is set to the surface temperature
	749	t_top(lsummit)=zt(1)
	750	do l=lsummit,nlayer-1
	751	!! temperature above the mountain follows the adiabatic
	752	t_top(l+1)=t_top(l)-(zlay(l+1)-zlay(l))*g/cpp
	753	enddo
	754
	755	!! Layer where to take the temperature above the mountain
	756	do while (hmons .ge. zlay(lmons))
	757	lmons=lmons+1
	758	enddo
	759	!! Real temperature of the environment near to the mountain is t_env(lsummit)=zt(lmons)
	760	!! (More simple and makes no real difference is to impose t_env(:)=zt(:))
	761	if (lmons.gt.lsummit) then
	762	t_env(lsummit)=zt(lmons)
	763	zlmons=zlay(lmons)
	764	do l=0,nlayer-lsummit-1!2
	765	zlmons=zlmons+(zlay(lsummit+l+1)-zlay(lsummit+l))
	766	ll=lmons
	767	do while (zlmons .ge. zlay(ll) .and. ll.lt.nlayer )
	768	ll=ll+1
	769	enddo
	770	ll=ll-1
	771	!! Interpolation above lsummit
	772	t_env(lsummit+l+1)= zt(ll) + ((zlmons-zlay(ll))*(zt(ll+1)-zt(ll))/(zlay(ll+1)-zlay(ll)))
	773	enddo
	774	t_env(nlayer)=zt(nlayer)
	775	endif ! (lmons.gt.lsummit)
	776	do l=1,nlayer
	777	!! temperature above the mountain follows the adiabatic up to reach the environment temperature
	778	if (t_top(l).le.t_env(l)) then
	779	t_top(l)=t_env(l)
	780	endif
	781	enddo
	782	!! Temperature delta created at the top of the mountain
	783	dt_top = t_top(lsummit)-t_env(lsummit)
	784
	785	END SUBROUTINE t_topmons
	786
	787	!=======================================================================
	788	! **********************************************************************
	789	! Subroutine to determine the vertical transport with
	790	! a Van Leer advection scheme (copied from the sedimentation scheme --> see vlz_fi.F)
	791	!***********************************************************************
	792	SUBROUTINE van_leer(nlay,q,pente_max,lwmax,masse,w,wq,masse_pbl,dqm_pbl,alpha_hmons,qbar_pbl)
	793
	794	IMPLICIT NONE
	795
	796	!--------------------------------------------------------
	797	! Input/Output Variables
	798	!--------------------------------------------------------
	799	INTEGER,INTENT(IN) :: nlay ! number of atmospheric layers
	800	REAL,INTENT(IN) :: masse(nlay) ! mass of the layer Dp/g
	801	REAL,INTENT(IN) :: pente_max != 2 conseillee
	802	INTEGER,INTENT(IN) :: lwmax ! layer of dust injection above the mountain, where the vertical velocity is maximum
	803	REAL,INTENT(INOUT) :: w(nlay) ! atmospheric mass "transferred" at each timestep (kg.m-2)
	804	REAL,INTENT(INOUT) :: wq(nlay+1)
	805	REAL,INTENT(INOUT) :: q(nlay) ! mixing ratio (kg/kg)
	806	REAL,INTENT(IN) :: masse_pbl
	807	REAL,INTENT(IN) :: dqm_pbl
	808	REAL,INTENT(IN) :: alpha_hmons
	809	REAL,INTENT(IN) :: qbar_pbl
	810
	811	!--------------------------------------------------------
	812	! Local Variables
	813	!--------------------------------------------------------
	814
	815	INTEGER i,l,j,ii
	816	REAL dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax
	817	REAL sigw, Mtot, MQtot
	818	INTEGER m
	819
	820	! **********************************************************************
	821	! Mixing ratio vertical gradient at the levels
	822	! **********************************************************************
	823	do l=lwmax+1,nlay !l=2,nlay
	824	dzqw(l)=q(l-1)-q(l)
	825	adzqw(l)=abs(dzqw(l))
	826	enddo
	827
	828	do l=lwmax+1,nlay-1 !l=2,nlay-1
	829	if(dzqw(l)*dzqw(l+1).gt.0.) then
	830	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
	831	else
	832	dzq(l)=0.
	833	endif
	834	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
	835	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
	836	enddo
	837
	838	dzq(lwmax)=0. !! dzq(1)=0. => Attention: in the case of van leer only above the mountain it is very important to initialize the transport from lwmax !!
	839	dzq(nlay)=0.
	840	! **********************************************************************
	841	! Vertical advection
	842	! **********************************************************************
	843
	844	!! No flux at the model top:
	845	wq(nlay+1)=0.
	846
	847	!! Mass flux injected at lwmax
	848	wq(lwmax) = -dqm_pbl/alpha_hmons ! dqm_pbl = alpha_hmons x wup(wmax) x rho(lwmax) x ptimestep x qbar_pbl
	849	! /alpha_hmons because the variable considered is mr_topdust
	850	do l = lwmax,nlay-1
	851
	852	! 1) Compute wq where w < 0 (up) (UPWARD TRANSPORT)
	853	! ===============================
	854
	855	if (w(l+1).le.0) then
	856	! Regular scheme (transfered mass < 1 layer)
	857	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	858	if(-w(l+1).le.masse(l))then
	859	sigw=w(l+1)/masse(l)
	860	wq(l+1)=w(l+1)(q(l)-0.5(1.+sigw)*dzq(l))
	861	!!-------------------------------------------------------
	862	! The following part should not be needed in the
	863	! case of an integration over subtimesteps
	864	!!-------------------------------------------------------
	865	!! Extended scheme (transfered mass > 1 layer)
	866	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	867	!! else
	868	!! m = l-1
	869	!! Mtot = masse(m+1)
	870	!! MQtot = masse(m+1)*q(m+1)
	871	!! if (m.lt.lwmax)goto 77!if (m.le.0)goto 77
	872	!! do while(-w(l+1).gt.(Mtot+masse(m)))
	873	!! m=m-1
	874	!! Mtot = Mtot + masse(m+1)
	875	!! MQtot = MQtot + masse(m+1)*q(m+1)
	876	!! if (m.lt.lwmax)goto 77!if (m.le.0)goto 77
	877	!! end do
	878	!! 77 continue
	879	!!
	880	!! if (m.gt.lwmax) then!if (m.gt.0) then
	881	!! sigw=(w(l+1)+Mtot)/masse(m)
	882	!! wq(l+1)= -(MQtot + (-w(l+1)-Mtot)* &
	883	!! (q(m)-0.5(1.+sigw)dzq(m)))
	884	!! else if ((-w(l+1).gt.(Mtot))) then
	885	!! Mtot=Mtot+masse_pbl!*alpha_hmons
	886	!! MQtot=MQtot+masse_pblqbar_pbl!alpha_hmons
	887	!! sigw=(w(l+1)+Mtot)/masse(m)
	888	!! wq(l+1)= -(MQtot + (-w(l+1)-Mtot)*qbar_pbl)
	889	!! else
	890	!! w(l+1) = -Mtot
	891	!! wq(l+1) = -MQtot
	892	!! end if
	893	!!
	894	endif
	895
	896	! 2) Compute wq where w > 0 (down) (DOWNWARD TRANSPORT)
	897	! ===============================
	898
	899	else if (w(l).gt.0.) then
	900	! Regular scheme (transfered mass < 1 layer)
	901	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	902	if(w(l).le.masse(l))then
	903	sigw=w(l)/masse(l)
	904	wq(l)=w(l)(q(l)+0.5(1.-sigw)*dzq(l))
	905	!!-------------------------------------------------------
	906	! The following part should not be needed in the
	907	! case of an integration over subtimesteps
	908	!!-------------------------------------------------------
	909	!! Extended scheme (transfered mass > 1 layer)
	910	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	911	!! else
	912	!! m=l
	913	!! Mtot = masse(m)
	914	!! MQtot = masse(m)*q(m)
	915	!! if(m.ge.nlay)goto 88
	916	!! do while(w(l).gt.(Mtot+masse(m+1)))
	917	!! m=m+1
	918	!! Mtot = Mtot + masse(m)
	919	!! MQtot = MQtot + masse(m)*q(m)
	920	!! if(m.ge.nlay)goto 88
	921	!! end do
	922	!! 88 continue
	923	!! if (m.lt.nlay) then
	924	!! sigw=(w(l)-Mtot)/masse(m+1)
	925	!! wq(l)=(MQtot + (w(l)-Mtot)* &
	926	!! (q(m+1)+0.5(1.-sigw)dzq(m+1)) )
	927	!! else
	928	!! w(l) = Mtot
	929	!! wq(l) = MQtot
	930	!! end if
	931	end if ! (w(l).le.masse(l))
	932
	933	end if ! (w(l+1).le.0)
	934
	935	enddo ! l = lwmax+1,nlay-1
	936
	937	! **********************************************************************
	938	! Mixing ratio update after the vertical transport
	939	! **********************************************************************
	940	do l = lwmax,nlay-1 !l = 1,nlay-1 ! loop different than when w>0
	941
	942	if ( (wq(l+1)-wq(l)) .lt. -(masse(l)*q(l)) .and. (abs(wq(l+1)).gt. 0.) .and. (abs(wq(l)).gt. 0.) ) then
	943	wq(l+1) = wq(l)-masse(l)*q(l)
	944	end if
	945
	946	q(l)=q(l) + (wq(l+1)-wq(l))/masse(l)
	947
	948	enddo
	949
	950
[2628]	951	end subroutine van_leer
	952
[2199]	953	!********************************************************************************
[2628]	954	subroutine topmons_setup(ngrid)
	955	! Purpose:
	956	! 1) Fill the logical array contains_mons(:),
	957	! with contains_mons(ig)=True if there is
	958	! a mountain in the mesh ig
	959	! 2) Compute alpha_hmons(:) and hsummit(:)
	960	use surfdat_h,only: phisfi,hmons,summit,base,&
	961	alpha_hmons,hsummit,contains_mons
	962	use comcstfi_h,only: pi,g
	963	use planetwide_mod,only: planetwide_maxval, planetwide_minval
[2930]	964	use mod_grid_phy_lmdz, only: nvertex !number of encompassing vertex
[2628]	965	use geometry_mod,only: longitude_deg,latitude_deg,&
[2930]	966	boundslon,boundslat ! if nvertex==4 then we have
	967	! boundslon/lat(ngrid,4) :
[2628]	968	! \|------------------------------\|
	969	! \|north_west=2 north_east=1\|
	970	! \| \|
	971	! \| (ig) \|
	972	! \| \|
	973	! \|south_west=3 south_east=4\|
	974	! \|------------------------------\|
[2930]	975	! if nvertex==6 (isocahedral grid)
	976	! there is no dedicated ordering
[2628]	977
	978	implicit none
	979	integer,intent(in) :: ngrid
	980
	981	! Local variables
	982	integer,parameter :: ntop_max = 19 ! total number of mounts written in the hmons list
	983	real lon_top(ntop_max),lat_top(ntop_max) ! coordinates of the mounts (in deg)
	984	! Mountains list :
	985	! Olympus Mons,Ascraeus Mons,Elysium Mons,Arsia Mons,Pavonis Mons,
	986	! Hecates Tholus,Tharsis Tholus,Ceraunius Tholus,Alba Mons,Apollinaris Mons,
	987	! Albor Tholus,Biblis Tholus,Anseris Mons,Ulysses Tholus,Aeolis Mons,
	988	! Euripus Mons,Hadriacus Mons,Tyrrhenus Mons,Uranius Mons
	989	!
	990	! NB: in 64x48 horiz. resolution, Biblis Tholus & Ulysses Tholus fall in the
	991	! same mesh, hence only Biblis Tholus is kept by the alpha_hmons computation
	992	data lon_top/-134.,-104.5,146.9,-121.1,-113.4,&
	993	150.2,-90.8,-97.4,-109.6,174.4,&
	994	150.4,-124.6,86.6,-121.6,137.8,&
	995	105.,91.8,106.5,-92.2/
	996	data lat_top/18.4,11.8,24.8,-8.4,-0.8,&
	997	31.8,13.4,24.,40.4,-9.3,&
	998	18.8,2.6,-29.8,2.9,-5.4,&
	999	-44.8,-32.1,-21.1,26.8/
	1000	integer,parameter :: ntop = 19 ! the topmons scheme is limited to the first ntop mounts
[2930]	1001	real :: boundslon_deg(nvertex),boundslat_deg(nvertex)
[2628]	1002	real :: hmax,hmin
[2930]	1003	real :: minlon,maxlon ! min/max encompassing vertex longitude (deg)
	1004	real :: minlat,maxlat ! min/max encompassing vertex latitude (deg)
[2628]	1005	integer :: ig,itop
	1006
	1007
	1008
	1009	IF (ngrid.gt.1) THEN
	1010	! Sanity check
	1011	if (ntop.gt.ntop_max) then
	1012	call abort_physic("topmons_setup","Number of mountains ntop greater than ntop_max",1)
	1013	endif
[2199]	1014
[2628]	1015	! Determine contains_mons
	1016	contains_mons(:)=.false.
[2199]	1017
[2628]	1018	do ig=1,ngrid
	1019	boundslon_deg(:)=boundslon(ig,:)/pi*180.
	1020	boundslat_deg(:)=boundslat(ig,:)/pi*180.
[2930]	1021
	1022	minlon=minval(boundslon_deg)
	1023	maxlon=maxval(boundslon_deg)
	1024	minlat=minval(boundslat_deg)
	1025	maxlat=maxval(boundslat_deg)
[2628]	1026
	1027	do itop=1,ntop
[2930]	1028	if ((lon_top(itop).ge.minlon).and.(lon_top(itop).le.maxlon).and. &
	1029	(lat_top(itop).ge.minlat).and.(lat_top(itop).le.maxlat)) then
[2628]	1030	contains_mons(ig)=.true.
[2930]	1031	write(,) "topmons_setup: Found a mount for ig:",ig
	1032	write(,) " at : lon=",longitude_deg(ig)," lat=",latitude_deg(ig)
	1033	if (nvertex==4) then
	1034	write(,) "(",boundslon_deg(2),",",boundslat_deg(2),") (",boundslon_deg(1),",",boundslat_deg(1),")"
	1035	write(,) " ((",lon_top(itop),",",lat_top(itop),"))"
	1036	write(,) "(",boundslon_deg(3),",",boundslat_deg(3),") (",boundslon_deg(4),",",boundslat_deg(4),")"
	1037	endif
[2628]	1038	endif
[2930]	1039	enddo ! of do itop=1,ntop
	1040	enddo ! of do ig=1,ngrid
[2628]	1041
	1042	! Compute alpha_hmons
	1043	call planetwide_maxval(hmons,hmax)
	1044	call planetwide_minval(hmons,hmin)
[2930]	1045
	1046	! Add a sanity check in case there is no topography:
	1047	if (hmin==hmax) then
	1048	write(,) "topmons_setup: hmin=",hmin," hmax=",hmax
	1049	write(,) "topmons_setup: thus no topflows !!!"
	1050	contains_mons(:)=.false.
	1051	endif
	1052
[2628]	1053	do ig=1,ngrid
	1054	if (contains_mons(ig)) then
	1055	! the mesh ig contains a mountain
	1056	alpha_hmons(ig)= 0.5*(hmons(ig)-hmin)/(hmax-hmin)
	1057	! Sanity check
	1058	if (alpha_hmons(ig).le.0) then
	1059	call abort_physic("topmons_setup","ERROR: alpha_hmons cannot be <0 "//&
	1060	"if the mesh contains a mountain. Please check your "//&
	1061	"formula or your start files.",1)
	1062	endif
	1063	else
	1064	! the mesh ig doesn't contain a mountain
	1065	alpha_hmons(ig)= 0
	1066	endif
[2930]	1067	enddo ! of do ig=1,ngrid
[2628]	1068
	1069	! Compute hsummit
	1070	hsummit(:)=summit(:)-phisfi(:)/g
[2199]	1071
[2628]	1072	ELSE ! 1D case
	1073	hmin=0.
	1074	hmax=23162.1 !set here the height of the sub-grid scale topography
	1075	do ig=1,ngrid
	1076	alpha_hmons(ig)= (hmons(ig)-hmin)/(hmax-hmin) !0.1*(hmons(ig)-hmin)/(hmax-hmin)
	1077	print*,"1D, hmons=",hmons(ig),"alpha=",alpha_hmons(ig)
	1078	enddo
	1079
	1080	hsummit(:)=14000.
	1081	ENDIF ! (ngrid.gt.1)
	1082
	1083	end subroutine topmons_setup
	1084	!********************************************************************************
	1085
[2199]	1086	END MODULE topmons_mod

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

Download in other formats:

Original Format