Context Navigation

rocketduststorm_mod.F90 @ 2167

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

Mars GCM:
Set adpatable parameters for the rocket dust storm scheme (parameters included in callkeys.h, and adaptable according to the callphys.def with the function "call getin" in conf_phys.F):

ti_injection, tf_injection (by default: ti_injection=10. and tf_injection=12., impacted files: compute_dtau_mod.F90, vdifc_mod.F)
coeff_detrainment (by default: coeff_detrainment=0., impacted files: rocketduststorm_mod.F90)

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

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

Download in other formats:

Original Format