Context Navigation

rocketduststorm_mod.F90 @ 2079

Last change on this file since 2079 was 2079, checked in by mvals, 6 years ago

Mars GCM:

Update of rocketduststorm_mod.F90 :

We want to separate both parametrizations related to the formation of the detached dust layers. Therefore, rocketduststorm_mod.F90 now only comprises the rocket dust storm scheme, whereas it contained
also before the calculation of the vertical velocity induced by the presence of the sub-grid scale topography. This latter part is under development and will be integrated as a fully independant
parametrization: the aim is to simulate the entrainment of dust by slope winds, from the boundary layer up to the top of the sub-grid scale topography.

Addition of initial surface parameters "summit" and "base" to prepare the previously described slope wind parametrization

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

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

Download in other formats:

Original Format