Context Navigation

topmons_mod.F90 @ 2628

Last change on this file since 2628 was 2628, checked in by abierjon, 3 years ago

Mars GCM:
Big changes on mountain top dust flows for GCM6:

the scheme now activates only in grid meshes that contain a mountain among a hard-written list, instead of every meshes. This is done to prevent strong artificial reinjections of dust in places that don't present huge converging slopes enabling the concentration of dust (ex: Valles Marineris, Hellas). Topdust is now also detrained as soon as it leaves the column it originated from.
the list of the 19 allowed mountains is used by the subroutine topmons_setup in module topmons_mod, to compute a logical array contains_mons(ngrid). alpha_hmons and hsummit are also set up once and for all by this subroutine, which is called in physiq_mod's firstcall.
contains_mons, alpha_hmons and hsummit are now saved variables of the module surfdat_h, and are called as such and not as arguments in the subroutines using them.
the logical flag "slpwind" and the comments in the code have also been updated to the new terminology "mountain top dust flows", accordingly to ticket #71. The new flag read in callphys.def is "topflows".

File size: 52.8 KB

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

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

Download in other formats:

Original Format