Context Navigation

rocketduststorm_mod.F90 @ 2214

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

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

Download in other formats:

Original Format