Context Navigation

rocketduststorm_mod.F90 @ 2437

Last change on this file since 2437 was 2417, checked in by emillour, 5 years ago

Mars GCM:
Add a new scheme to handle the dust and its radiative impact. Triggered using
a new flag dustscaling_mode=2 (dustscaling_mod=0: no rescaling at all, and
dustscaling_mode=1: full rescaling using tauscaling, GCMv5.3 style). Rescaling
is then only done on the radiative impact (see dust_scaling_mod.F90) of dust.
Moreover the scaling factor "dust_rad_adjust" is evaluated using the target dust
scenario opacity for the next sol and left to evolve linearly until then to not
impose the diurnal evolution of dust.
In practice, main changes or additions in the code are:

renamed flag "tauscaling_mode" as "dustscaling_mode"
moved parameter "t_scenario_sol" to "dust_param_mod"
adapted "compute_dustscaling" routine in "dust_scaling_mod"
added module "dust_rad_adjust_mod"
2D fields "dust_rad_adjust_prev" and "dust_rad_adjust_next" required to compute coefficient "dust_rad_adjust" need to be stored in (re)startfi files

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

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

Download in other formats:

Original Format