Context Navigation

topmons_mod.F90 @ 2682

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

Mars GCM:
Implementation of the IRabs-to-VISext dust scenario conversion, depending on the GCM effective radius :

new flag 'reff_driven_IRtoVIS_scenario', false by default, must be set to true to use this dynamic conversion (otherwise, the coefficient takes the constant value of 2.6, eqv to reff=1.5um, as before) A specific line is added in deftank/callphys.def.GCM6
this flag requires the 'dustiropacity' to be set to 'tes' to match the IR scenario's wavelength. 'mcs'-like dso diagnostics can only be produced by the use of the post-processing tool util/aeroptical.F90
the variable IRtoVIScoef is computed in aeropacity via the GCM CDOD ratio : tau_pref_gcm_VIS/tau_pref_gcm_IR (only at the first call to callradite during each physical timestep)
change read_dust_scenario into module read_dust_scenario_mod

File size: 53.2 KB

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

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

Download in other formats:

Original Format