Context Navigation

topmons_mod.F90 @ 2947

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

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

Download in other formats:

Original Format