Context Navigation

rocketduststorm_mod.F90 @ 2282

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

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

Download in other formats:

Original Format