Context Navigation

rocketduststorm_mod.F90 @ 2673

Last change on this file since 2673 was 2643, checked in by abierjon, 3 years ago

Mars GCM:
Implementation of the IRabs-to-VISext dust scenario conversion, depending on the GCM effective radius :

new flag 'reff_driven_IRtoVIS_scenario', false by default, must be set to true to use this dynamic conversion (otherwise, the coefficient takes the constant value of 2.6, eqv to reff=1.5um, as before) A specific line is added in deftank/callphys.def.GCM6
this flag requires the 'dustiropacity' to be set to 'tes' to match the IR scenario's wavelength. 'mcs'-like dso diagnostics can only be produced by the use of the post-processing tool util/aeroptical.F90
the variable IRtoVIScoef is computed in aeropacity via the GCM CDOD ratio : tau_pref_gcm_VIS/tau_pref_gcm_IR (only at the first call to callradite during each physical timestep)
change read_dust_scenario into module read_dust_scenario_mod

File size: 31.3 KB

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

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

Download in other formats:

Original Format