Context Navigation

rocketduststorm_mod.F90 @ 2415

Last change on this file since 2415 was 2415, checked in by emillour, 4 years ago
Mars GCM: Code cleanup: remove "tauref" and replace it with two distinct variables, "tau_pref_gcm" and "tau_pref_scenario" which are repectively the visible dust opacity column (at 610 Pa) in the GCM and prescribed by a dust scenario. EM
File size: 30.8 KB

Line
1	MODULE rocketduststorm_mod
2
3	IMPLICIT NONE
4
5	REAL, SAVE, ALLOCATABLE :: dustliftday(:) ! dust lifting rate (s-1)
6
7	CONTAINS
8
9	!=======================================================================
10	! ROCKET DUST STORM - vertical transport and detrainment
11	!=======================================================================
12	! calculation of the vertical flux
13	! call of van_leer : Van Leer transport scheme of the dust tracers
14	! detrainement of stormdust in the background dust
15	! -----------------------------------------------------------------------
16	! Authors: M. Vals; C. Wang; F. Forget; T. Bertrand
17	! Institution: Laboratoire de Meteorologie Dynamique (LMD) Paris, France
18	! -----------------------------------------------------------------------
19
20	SUBROUTINE rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep, &
21	pq,pdqfi,pt,pdtfi,pplev,pplay,pzlev, &
22	pzlay,pdtsw,pdtlw, &
23	! input for radiative transfer
24	clearatm,icount,zday,zls, &
25	tsurf,igout,totstormfract, &
26	tauscaling, &
27	! input sub-grid scale cloud
28	clearsky,totcloudfrac, &
29	! input sub-grid scale topography
30	nohmons,alpha_hmons, &
31	! output
32	pdqrds,wrad,dsodust,dsords,dsotop, &
33	tau_pref_scenario,tau_pref_gcm)
34
35	USE tracer_mod, only: igcm_stormdust_mass,igcm_stormdust_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,zmea, zstd, zsig, hmons
42	USE callradite_mod, only: callradite
43	IMPLICIT NONE
44
45	include "callkeys.h"
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) :: pdqfi(ngrid,nlayer,nq)! tendancy field pq
59	REAL, INTENT(IN) :: pt(ngrid,nlayer) ! temperature at mid-layer (K)
60	REAL, INTENT(IN) :: pdtfi(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	LOGICAL, INTENT(IN) :: clearatm
72	INTEGER, INTENT(INOUT) :: icount
73	REAL, INTENT(IN) :: zday
74	REAL, INTENT(IN) :: zls
75	REAL, INTENT(IN) :: tsurf(ngrid)
76	INTEGER, INTENT(IN) :: igout
77	REAL, INTENT(IN) :: totstormfract(ngrid)
78	REAL, INTENT(INOUT) :: tauscaling(ngrid)
79
80	! sbgrid scale water ice clouds
81	logical, intent(in) :: clearsky
82	real, intent(in) :: totcloudfrac(ngrid)
83
84	! sbgrid scale topography
85	LOGICAL, INTENT(IN) :: nohmons
86	REAL, INTENT(IN) :: alpha_hmons(ngrid)
87
88	!--------------------------------------------------------
89	! Output Variables
90	!--------------------------------------------------------
91
92	REAL, INTENT(OUT) :: pdqrds(ngrid,nlayer,nq) ! tendancy field for dust when detraining
93	REAL, INTENT(OUT) :: wrad(ngrid,nlayer+1) ! vertical speed within the rocket dust storm
94	REAL, INTENT(OUT) :: dsodust(ngrid,nlayer) ! density scaled opacity of env. dust
95	REAL, INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity of storm dust
96	REAL, INTENT(OUT) :: dsotop(ngrid,nlayer) ! density scaled opacity of topmons dust
97	REAL,INTENT(OUT) :: tau_pref_scenario(ngrid) ! prescribed dust column
98	! visible opacity at odpref
99	REAL,INTENT(OUT) :: tau_pref_gcm(ngrid) ! dust column visible opacity at
100	! odpref in the GCM
101	!--------------------------------------------------------
102	! Local variables
103	!--------------------------------------------------------
104	INTEGER l,ig,iq,ll
105	! local variables from callradite.F
106	REAL zdtlw1(ngrid,nlayer) ! (K/s) storm
107	REAL zdtsw1(ngrid,nlayer) ! (K/s) storm
108	REAL zt(ngrid,nlayer) ! actual temperature at mid-layer (K)
109	REAL zdtvert(ngrid,nlayer) ! dT/dz , lapse rate
110	REAL ztlev(ngrid,nlayer) ! temperature at intermediate levels l+1/2 without last level
111
112	REAL zdtlw1_lev(nlayer),zdtsw1_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2 for stormdust
113	REAL zdtlw_lev(nlayer),zdtsw_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2 for background dust
114
115	REAL zq_stormdust_mass(ngrid,nlayer) ! intermediate tracer stormdust mass
116	REAL zq_stormdust_number(ngrid,nlayer) ! intermediate tracer stormdust number
117	REAL zq_dust_mass(ngrid,nlayer) ! intermediate tracer dust mass
118	REAL zq_dust_number(ngrid,nlayer) ! intermediate tracer dust number
119
120	REAL mr_stormdust_mass(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (stormdust mass)
121	REAL mr_stormdust_number(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (stormdust number)
122	REAL mr_dust_mass(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (dust mass)
123	REAL mr_dust_number(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (sdust number)
124
125	REAL dqvl_stormdust_mass(ngrid,nlayer) ! tendancy of vertical transport (stormdust mass)
126	REAL dqvl_stormdust_number(ngrid,nlayer) ! tendancy of vertical transport (stormdust number)
127	REAL dqvl_dust_mass(ngrid,nlayer) ! tendancy of vertical transport (dust mass)
128	REAL dqvl_dust_number(ngrid,nlayer) ! tendancy of vertical transport (dust number)
129	REAL dqdet_stormdust_mass(ngrid,nlayer) ! tendancy of detrainement (stormdust mass)
130	REAL dqdet_stormdust_number(ngrid,nlayer) ! tendancy of detrainement (stormdust number)
131
132	REAL masse_col(nlayer) ! mass of atmosphere (kg/m2)
133	REAL zq(ngrid,nlayer,nq) ! updated tracers
134
135	REAL w(ngrid,nlayer) ! air mass flux (calculated with the vertical wind velocity profile) used as input in Van Leer (kgair/m2)
136	REAL wqmass(ngrid,nlayer+1) ! tracer (dust_mass) mass flux in Van Leer (kg/m2)
137	REAL wqnumber(ngrid,nlayer+1) ! tracer (dust_number) mass flux in Van Leer (kg/m2)
138
139	LOGICAL storm(ngrid) ! true when there is a dust storm (if the opacity is high): trigger the rocket dust storm scheme
140	REAL detrain(ngrid,nlayer) ! coefficient for detrainment : % of stormdust detrained
141	INTEGER scheme(ngrid) ! triggered scheme
142
143	REAL,PARAMETER:: coefmin =0.025 ! 0<coefmin<1 Minimum fraction of stormdust detrained
144	REAL,PARAMETER:: wmin =0.25 ! stormdust detrainment if wrad < wmin
145	REAL,PARAMETER:: wmax =10. ! maximum vertical velocity of the rocket dust storms (m/s)
146
147	! subtimestep
148	INTEGER tsub
149	INTEGER nsubtimestep !number of subtimestep when calling van_leer
150	REAL subtimestep !ptimestep/nsubtimestep
151	REAL dtmax !considered time needed for dust to cross one layer
152	REAL,PARAMETER:: secu=3.!3. !coefficient on wspeed to avoid dust crossing many layers during subtimestep
153
154	! diagnostics
155	REAL lapserate(ngrid,nlayer)
156	REAL deltahr(ngrid,nlayer+1)
157
158	LOGICAL,SAVE :: firstcall=.true.
159
160	! variables for the radiative transfer
161	REAL fluxsurf_lw1(ngrid)
162	REAL fluxsurf_sw1(ngrid,2)
163	REAL fluxtop_lw1(ngrid)
164	REAL fluxtop_sw1(ngrid,2)
165	REAL tau(ngrid,naerkind)
166	REAL aerosol(ngrid,nlayer,naerkind)
167	REAL taucloudtes(ngrid)
168	REAL rdust(ngrid,nlayer)
169	REAL rstormdust(ngrid,nlayer)
170	REAL rtopdust(ngrid,nlayer)
171	REAL rice(ngrid,nlayer)
172	REAL nuice(ngrid,nlayer)
173
174
175	! **********************************************************************
176	! **********************************************************************
177	! Rocket dust storm parametrization to reproduce the detached dust layers
178	! during the dust storm season:
179	! The radiative warming due to the presence of storm dust is
180	! balanced by the adiabatic cooling. The tracer "storm dust"
181	! is transported by the upward/downward flow.
182	! **********************************************************************
183	! **********************************************************************
184	!! 1. Radiative transfer in storm dust
185	!! 2. Compute vertical velocity for storm dust
186	!! case 1 storm = false: nothing to do
187	!! case 2 rocket dust storm (storm=true)
188	!! 3. Vertical transport (Van Leer)
189	!! 4. Detrainment
190	! **********************************************************************
191	! **********************************************************************
192
193
194	! **********************************************************************
195	! Initializations
196	! **********************************************************************
197	storm(:)=.false.
198	pdqrds(:,:,:) = 0.
199	mr_dust_mass(:,:)=0.
200	mr_dust_number(:,:)=0.
201	mr_stormdust_mass(:,:)=0.
202	mr_stormdust_number(:,:)=0.
203	dqvl_dust_mass(:,:)=0.
204	dqvl_dust_number(:,:)=0.
205	dqvl_stormdust_mass(:,:)=0.
206	dqvl_stormdust_number(:,:)=0.
207	dqdet_stormdust_mass(:,:)=0.
208	dqdet_stormdust_number(:,:)=0.
209	wrad(:,:)=0.
210	w(:,:)=0.
211	wqmass(:,:)=0.
212	wqnumber(:,:)=0.
213	zdtvert(:,:)=0.
214	lapserate(:,:)=0.
215	deltahr(:,:)=0.
216	scheme(:)=0
217	detrain(:,:)=1.
218
219	!! no update for the stormdust tracer and temperature fields
220	!! because previous callradite was for background dust
221	zq(1:ngrid,1:nlayer,1:nq)=pq(1:ngrid,1:nlayer,1:nq)
222	zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer)
223
224
225	zq_dust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_mass)
226	zq_dust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_number)
227	zq_stormdust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_stormdust_mass)
228	zq_stormdust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_stormdust_number)
229
230	! *********************************************************************
231	! 0. Check if there is a storm
232	! *********************************************************************
233	DO ig=1,ngrid
234	storm(ig)=.false.
235	DO l=1,nlayer
236	IF (zq(ig,l,igcm_stormdust_mass) &
237	.gt. zq(ig,l,igcm_dust_mass)*(1.E-4)) THEN
238	storm(ig)=.true.
239	EXIT
240	ENDIF
241	ENDDO
242	ENDDO
243
244	! *********************************************************************
245	! 1. Call the second radiative transfer for stormdust, obtain the extra heating
246	! *********************************************************************
247	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, &
248	emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, &
249	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
250	fluxtop_sw1,tau_pref_scenario,tau_pref_gcm, &
251	tau,aerosol,dsodust,tauscaling, &
252	taucloudtes,rdust,rice,nuice,co2ice,rstormdust,rtopdust, &
253	totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,&
254	clearsky,totcloudfrac)
255
256	! **********************************************************************
257	! 2. Compute vertical velocity for storm dust
258	! **********************************************************************
259	!! **********************************************************************
260	!! 2.1 Nothing to do when no storm
261	!! no storm
262	DO ig=1,ngrid
263	IF (.NOT.(storm(ig))) then
264	scheme(ig)=1
265	cycle
266	ENDIF ! IF (storm(ig))
267	ENDDO ! DO ig=1,ngrid
268
269	!! **********************************************************************
270	!! 2.2 Calculation of the extra heating : computing heating rates
271	!! gradient at boundaries of each layer, start from surface
272	DO ig=1,ngrid
273	IF (storm(ig)) THEN
274
275	scheme(ig)=2
276
277	!! computing heating rates gradient at boundraies of each layer
278	!! start from surface
279	zdtlw1_lev(1)=0.
280	zdtsw1_lev(1)=0.
281	zdtlw_lev(1)=0.
282	zdtsw_lev(1)=0.
283	ztlev(ig,1)=zt(ig,1)
284
285	DO l=1,nlayer-1
286	!! Calculation for the dust storm fraction
287	zdtlw1_lev(l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
288	zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
289	(pzlay(ig,l+1)-pzlay(ig,l))
290
291	zdtsw1_lev(l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
292	zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
293	(pzlay(ig,l+1)-pzlay(ig,l))
294	!! Calculation for the background dust fraction
295	zdtlw_lev(l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
296	pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
297	(pzlay(ig,l+1)-pzlay(ig,l))
298
299	zdtsw_lev(l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
300	pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
301	(pzlay(ig,l+1)-pzlay(ig,l))
302
303	ztlev(ig,l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
304	zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
305	(pzlay(ig,l+1)-pzlay(ig,l))
306	ENDDO ! DO l=1,nlayer-1
307
308	!! This is the env. lapse rate
309	zdtvert(ig,1)=0.
310	DO l=1,nlayer-1
311	zdtvert(ig,l+1)=(ztlev(ig,l+1)-ztlev(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
312	lapserate(ig,l+1)=zdtvert(ig,l+1)
313	ENDDO
314
315	!! **********************************************************************
316	!! 2.3 Calculation of the vertical velocity : extra heating
317	!! balanced by adiabatic cooling
318
319	DO l=1,nlayer
320	deltahr(ig,l)=(zdtlw1_lev(l)+zdtsw1_lev(l)) &
321	-(zdtlw_lev(l)+zdtsw_lev(l))
322	wrad(ig,l)=-deltahr(ig,l)/(g/cpp+ &
323	max(zdtvert(ig,l),-0.99*g/cpp))
324	!! Limit vertical wind in case of lapse rate close to adiabatic
325	wrad(ig,l)=max(wrad(ig,l),-wmax)
326	wrad(ig,l)=min(wrad(ig,l),wmax)
327	ENDDO
328
329	ENDIF ! IF (storm(ig))
330	ENDDO ! DO ig=1,ngrid
331
332	! **********************************************************************
333	! 3. Vertical transport
334	! **********************************************************************
335	!! **********************************************************************
336	!! 3.1 Transport of background dust + storm dust (concentrated)
337	DO ig=1,ngrid
338	IF (storm(ig)) THEN
339	DO l=1,nlayer
340	mr_dust_mass(ig,l) = zq_dust_mass(ig,l)
341	mr_dust_number(ig,l) = zq_dust_number(ig,l)
342	mr_stormdust_mass(ig,l) = zq_dust_mass(ig,l)+ &
343	zq_stormdust_mass(ig,l)/totstormfract(ig)
344	mr_stormdust_number(ig,l) = zq_dust_number(ig,l)+ &
345	zq_stormdust_number(ig,l)/totstormfract(ig)
346	ENDDO ! DO l=1,nlayer
347	ENDIF ! IF (storm(ig))
348	ENDDO ! DO ig=1,ngrid
349
350	DO ig=1,ngrid
351	IF (storm(ig)) THEN
352	!! **********************************************************************
353	!! 3.2 Compute the subtimestep to conserve the mass in the Van Leer transport
354	dtmax=ptimestep
355	DO l=2,nlayer
356	IF (wrad(ig,l).lt.0.) THEN
357	dtmax=min(dtmax,(pzlev(ig,l)-pzlev(ig,l-1))/ &
358	(secu*abs(wrad(ig,l))))
359	ELSE IF (wrad(ig,l).gt.0.) then
360	dtmax=min(dtmax,(pzlev(ig,l+1)-pzlev(ig,l))/ &
361	(secu*abs(wrad(ig,l))))
362	ENDIF
363	ENDDO
364	nsubtimestep= int(ptimestep/dtmax)
365	subtimestep=ptimestep/float(nsubtimestep)
366	!! Mass flux generated by wup in kg/m2
367	DO l=1,nlayer
368	w(ig,l)=wrad(ig,l)pplev(ig,l)/(rztlev(ig,l)) &
369	*subtimestep
370	ENDDO ! l=1,nlayer
371
372	!! **********************************************************************
373	!! 3.3 Vertical transport by a Van Leer scheme
374	!! Mass of atmosphere in the layer
375	DO l=1,nlayer
376	masse_col(l)=(pplev(ig,l)-pplev(ig,l+1))/g
377	ENDDO
378	!! Mass flux in kg/m2 if you are not using the subtimestep
379	!DO l=1,nlayer
380	! w(ig,l)=wrad(ig,l)(pplev(ig,l)/(rztlev(ig,l)))*ptimestep
381	!ENDDO
382	!! Loop over the subtimestep
383	DO tsub=1,nsubtimestep
384	!! Van Leer scheme
385	wqmass(ig,:)=0.
386	wqnumber(ig,:)=0.
387	CALL van_leer(nlayer,mr_stormdust_mass(ig,:),2., &
388	masse_col,w(ig,:),wqmass(ig,:))
389	CALL van_leer(nlayer,mr_stormdust_number(ig,:),2., &
390	masse_col,w(ig,:),wqnumber(ig,:))
391	ENDDO !tsub=...
392
393	ENDIF ! IF storm(ig)
394	ENDDO ! DO ig=1,ngrid
395
396	!! **********************************************************************
397	!! 3.4 Re-calculation of the mixing ratios after vertical transport
398	DO ig=1,ngrid
399	IF (storm(ig)) THEN
400	DO l=1,nlayer
401
402	!! General and "healthy" case
403	IF (mr_stormdust_mass(ig,l).ge.mr_dust_mass(ig,l)) THEN
404	zq_dust_mass(ig,l) = mr_dust_mass(ig,l)
405	zq_dust_number(ig,l) = mr_dust_number(ig,l)
406	zq_stormdust_mass(ig,l) = totstormfract(ig)*(mr_stormdust_mass(ig,l)-mr_dust_mass(ig,l))
407	zq_stormdust_number(ig,l) = totstormfract(ig)*(mr_stormdust_number(ig,l)-mr_dust_number(ig,l))
408	!! Particular case: there is less than initial dust mixing ratio after the vertical transport
409	ELSE
410	zq_dust_mass(ig,l) = (1.-totstormfract(ig))mr_dust_mass(ig,l)+totstormfract(ig)mr_stormdust_mass(ig,l)
411	zq_dust_number(ig,l) = (1.-totstormfract(ig))mr_dust_number(ig,l)+totstormfract(ig)mr_stormdust_number(ig,l)
412	zq_stormdust_mass(ig,l) = 0.
413	zq_stormdust_number(ig,l) = 0.
414	ENDIF
415
416	ENDDO ! DO l=1,nlayer
417	ENDIF ! IF storm(ig)
418	ENDDO ! DO ig=1,ngrid
419
420	!! **********************************************************************
421	!! 3.5 Calculation of the tendencies of the vertical transport
422	DO ig=1,ngrid
423	IF (storm(ig)) THEN
424	DO l=1,nlayer
425	dqvl_stormdust_mass(ig,l) = (zq_stormdust_mass(ig,l)- &
426	zq(ig,l,igcm_stormdust_mass)) /ptimestep
427	dqvl_stormdust_number(ig,l) = (zq_stormdust_number(ig,l)- &
428	zq(ig,l,igcm_stormdust_number)) /ptimestep
429	dqvl_dust_mass(ig,l) = (zq_dust_mass(ig,l)-zq(ig,l,igcm_dust_mass)) /ptimestep
430	dqvl_dust_number(ig,l) = (zq_dust_number(ig,l)-zq(ig,l,igcm_dust_number)) /ptimestep
431	ENDDO
432	ENDIF ! IF storm(ig)
433	ENDDO ! DO ig=1,ngrid
434
435	! **********************************************************************
436	! 4. Detrainment: stormdust is converted to background dust
437	! **********************************************************************
438	!! **********************************************************************
439	!! 4.1 Compute the coefficient of detrainmen
440	DO ig=1,ngrid
441	DO l=1,nlayer
442	IF ((max(abs(wrad(ig,l)),abs(wrad(ig,l+1))) .lt. &
443	wmin) .or. (zq_dust_mass(ig,l) .gt. &
444	10000.*zq_stormdust_mass(ig,l))) THEN
445	detrain(ig,l)=1.
446	ELSE IF (max(abs(wrad(ig,l)),abs(wrad(ig,l+1))) &
447	.le. wmax) THEN
448	detrain(ig,l)=coeff_detrainment* &
449	(((1-coefmin)/(wmin-wmax)*2) &
450	(max(abs(wrad(ig,l)),abs(wrad(ig,l+1)))-wmax)**2 &
451	+coefmin)
452	ELSE IF (max(abs(wrad(ig,l)),abs(wrad(ig,l+1))).gt. wmax ) THEN
453	detrain(ig,l)=coefmin
454	ELSE
455	detrain(ig,l)=coefmin
456	ENDIF
457	ENDDO ! DO l=1,nlayer
458	ENDDO ! DO ig=1,ngrid
459
460	!! **********************************************************************
461	!! 4.2 Calculation of the tendencies of the detrainment
462	DO ig=1,ngrid
463	DO l=1,nlayer
464	dqdet_stormdust_mass(ig,l)=-detrain(ig,l)*zq_stormdust_mass(ig,l) &
465	/ptimestep
466	dqdet_stormdust_number(ig,l)=-detrain(ig,l)*zq_stormdust_number(ig,l) &
467	/ptimestep
468	ENDDO ! DO l=1,nlayer
469	ENDDO ! DO ig=1,ngrid
470
471	! **********************************************************************
472	! 5. Final tendencies: vertical transport + detrainment
473	! **********************************************************************
474	DO ig=1,ngrid
475	DO l=1,nlayer
476	pdqrds(ig,l,igcm_stormdust_mass)=dqdet_stormdust_mass(ig,l) &
477	+dqvl_stormdust_mass(ig,l)
478	pdqrds(ig,l,igcm_stormdust_number)=dqdet_stormdust_number(ig,l) &
479	+dqvl_stormdust_number(ig,l)
480	pdqrds(ig,l,igcm_dust_mass)= -dqdet_stormdust_mass(ig,l) &
481	+dqvl_dust_mass(ig,l)
482	pdqrds(ig,l,igcm_dust_number)= -dqdet_stormdust_number(ig,l) &
483	+dqvl_dust_number(ig,l)
484	ENDDO ! DO l=1,nlayer
485	ENDDO ! DO ig=1,ngrid
486
487	! ! **********************************************************************
488	! ! 6. To prevent from negative values
489	! ! **********************************************************************
490	! DO ig=1,ngrid
491	! DO l=1,nlayer
492	! IF ((pq(ig,l,igcm_stormdust_mass) &
493	! +pdqrds(ig,l,igcm_stormdust_mass)*ptimestep .le. 0.) .or. &
494	! (pq(ig,l,igcm_stormdust_number) &
495	! +pdqrds(ig,l,igcm_stormdust_number)*ptimestep .le. 0.)) THEN
496	! pdqrds(ig,l,igcm_stormdust_mass)=-pq(ig,l,igcm_stormdust_mass)/ptimestep
497	! pdqrds(ig,l,igcm_stormdust_number)=-pq(ig,l,igcm_stormdust_number)/ptimestep
498	! ENDIF
499	! ENDDO ! nlayer
500	! ENDDO ! DO ig=1,ngrid
501	!
502	! DO ig=1,ngrid
503	! DO l=1,nlayer
504	! IF ((pq(ig,l,igcm_dust_mass) &
505	! +pdqrds(ig,l,igcm_dust_mass)*ptimestep .le. 0.) .or. &
506	! (pq(ig,l,igcm_dust_number) &
507	! +pdqrds(ig,l,igcm_dust_number)*ptimestep .le. 0.)) THEN
508	! pdqrds(ig,l,igcm_dust_mass)=-pq(ig,l,igcm_dust_mass)/ptimestep
509	! pdqrds(ig,l,igcm_dust_number)=-pq(ig,l,igcm_dust_number)/ptimestep
510	! ENDIF
511	! ENDDO ! nlayer
512	! ENDDO ! DO ig=1,ngrid
513
514	!=======================================================================
515	! WRITEDIAGFI
516
517	call WRITEDIAGFI(ngrid,'lapserate','lapse rate in the storm', &
518	& 'k/m',3,lapserate)
519	call WRITEDIAGFI(ngrid,'deltahr','extra heating rates', &
520	& 'K/s',3,deltahr)
521	call writediagfi(ngrid,'scheme','which scheme',&
522	' ',2,real(scheme))
523
524	END SUBROUTINE rocketduststorm
525
526	!=======================================================================
527	! **********************************************************************
528	! Subroutine to determine the vertical transport with
529	! a Van Leer advection scheme (copied from the sedimentation scheme --> see vlz_fi.F)
530	!***********************************************************************
531	SUBROUTINE van_leer(nlay,q,pente_max,masse,w,wq)
532
533	IMPLICIT NONE
534
535	!--------------------------------------------------------
536	! Input/Output Variables
537	!--------------------------------------------------------
538	INTEGER,INTENT(IN) :: nlay ! number of atmospheric layers
539	REAL,INTENT(IN) :: masse(nlay) ! mass of the layer Dp/g
540	REAL,INTENT(IN) :: pente_max != 2 conseillee
541	REAL,INTENT(INOUT) :: q(nlay) ! mixing ratio (kg/kg)
542	REAL,INTENT(INOUT) :: w(nlay) ! atmospheric mass "transferred" at each timestep (kg.m-2)
543	REAL,INTENT(INOUT) :: wq(nlay+1)
544
545	!--------------------------------------------------------
546	! Local Variables
547	!--------------------------------------------------------
548
549	INTEGER i,l,j,ii
550	REAL dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax
551	REAL newmasse
552	REAL sigw, Mtot, MQtot
553	INTEGER m
554
555	! **********************************************************************
556	! Mixing ratio vertical gradient at the levels
557	! **********************************************************************
558	do l=2,nlay
559	dzqw(l)=q(l-1)-q(l)
560	adzqw(l)=abs(dzqw(l))
561	enddo
562
563	do l=2,nlay-1
564	if(dzqw(l)*dzqw(l+1).gt.0.) then
565	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
566	else
567	dzq(l)=0.
568	endif
569	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
570	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
571	enddo
572
573	dzq(1)=0.
574	dzq(nlay)=0.
575
576	! **********************************************************************
577	! Vertical advection
578	! **********************************************************************
579
580	!! No flux at the model top:
581	wq(nlay+1)=0.
582
583	!! Surface flux up:
584	if(w(1).lt.0.) wq(1)=0. ! warning : not always valid
585
586	do l = 1,nlay-1
587
588	! 1) Compute wq where w < 0 (up) (UPWARD TRANSPORT)
589	! ===============================
590
591	if(w(l+1).le.0)then
592	! Regular scheme (transfered mass < 1 layer)
593	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
594	if(-w(l+1).le.masse(l))then
595	sigw=w(l+1)/masse(l)
596	wq(l+1)=w(l+1)(q(l)-0.5(1.+sigw)*dzq(l))
597	!!-------------------------------------------------------
598	! The following part should not be needed in the
599	! case of an integration over subtimesteps
600	!!-------------------------------------------------------
601	! Extended scheme (transfered mass > 1 layer)
602	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
603	else
604	m = l-1
605	Mtot = masse(m+1)
606	MQtot = masse(m+1)*q(m+1)
607	if (m.le.0)goto 77
608	do while(-w(l+1).gt.(Mtot+masse(m)))
609	! do while(-w(l+1).gt.Mtot)
610	m=m-1
611	Mtot = Mtot + masse(m+1)
612	MQtot = MQtot + masse(m+1)*q(m+1)
613	if (m.le.0)goto 77
614	end do
615	77 continue
616
617	if (m.gt.0) then
618	sigw=(w(l+1)+Mtot)/masse(m)
619	wq(l+1)= -(MQtot + (-w(l+1)-Mtot)* &
620	(q(m)-0.5(1.+sigw)dzq(m)) )
621	else
622	w(l+1) = -Mtot
623	wq(l+1) = -MQtot
624	end if
625	endif ! (-w(l+1).le.masse(l))
626
627	! 2) Compute wq where w > 0 (down) (DOWNWARD TRANSPORT)
628	! ===============================
629
630	else if(w(l).gt.0.)then
631
632	! Regular scheme (transfered mass < 1 layer)
633	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
634	if(w(l).le.masse(l))then
635	sigw=w(l)/masse(l)
636	wq(l)=w(l)(q(l)+0.5(1.-sigw)*dzq(l))
637	!!-------------------------------------------------------
638	! The following part should not be needed in the
639	! case of an integration over subtimesteps
640	!!-------------------------------------------------------
641	! Extended scheme (transfered mass > 1 layer)
642	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
643	else
644	m=l
645	Mtot = masse(m)
646	MQtot = masse(m)*q(m)
647	if(m.ge.nlay)goto 88
648	do while(w(l).gt.(Mtot+masse(m+1)))
649	m=m+1
650	Mtot = Mtot + masse(m)
651	MQtot = MQtot + masse(m)*q(m)
652	if(m.ge.nlay)goto 88
653	end do
654	88 continue
655	if (m.lt.nlay) then
656	sigw=(w(l)-Mtot)/masse(m+1)
657	wq(l)=(MQtot + (w(l)-Mtot)* &
658	(q(m+1)+0.5(1.-sigw)dzq(m+1)) )
659	else
660	w(l) = Mtot
661	wq(l) = MQtot
662	end if
663	end if
664
665	end if ! w<0 (up)
666
667	enddo ! l = 1,nlay-1
668
669	do l = 1,nlay
670
671	! it cannot entrain more than available mass !
672	if ( (wq(l+1)-wq(l)) .lt. -(masse(l)*q(l)) ) then
673	wq(l+1) = wq(l)-masse(l)*q(l)
674	end if
675
676	q(l)=q(l) + (wq(l+1)-wq(l))/masse(l)
677
678	enddo
679
680	END SUBROUTINE van_leer
681
682	!=======================================================================
683	! Initialization of the module variables
684
685	subroutine ini_rocketduststorm_mod(ngrid)
686
687	implicit none
688
689	integer, intent(in) :: ngrid
690
691	allocate(dustliftday(ngrid))
692
693	end subroutine ini_rocketduststorm_mod
694
695	subroutine end_rocketduststorm_mod
696
697	implicit none
698
699	if (allocated(dustliftday)) deallocate(dustliftday)
700
701	end subroutine end_rocketduststorm_mod
702
703	END MODULE rocketduststorm_mod

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

Download in other formats:

Original Format