Context Navigation

topmons_mod.F90 @ 2740

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

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

Download in other formats:

Original Format