Context Navigation

topmons_mod.F90 @ 2825

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

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

Download in other formats:

Original Format