Context Navigation

rocketduststorm_mod.F90 @ 3741

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

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

Download in other formats:

Original Format