Context Navigation

rocketduststorm_mod.F90 @ 2566

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

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

Download in other formats:

Original Format