Context Navigation

topmons_mod.F90 @ 3118

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

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

Download in other formats:

Original Format