Context Navigation

topmons_mod.F90 @ 2930

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

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

Download in other formats:

Original Format