Context Navigation

rocketduststorm_mod.F90 @ 2628

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

Mars GCM:
Big changes on mountain top dust flows for GCM6:

the scheme now activates only in grid meshes that contain a mountain among a hard-written list, instead of every meshes. This is done to prevent strong artificial reinjections of dust in places that don't present huge converging slopes enabling the concentration of dust (ex: Valles Marineris, Hellas). Topdust is now also detrained as soon as it leaves the column it originated from.
the list of the 19 allowed mountains is used by the subroutine topmons_setup in module topmons_mod, to compute a logical array contains_mons(ngrid). alpha_hmons and hsummit are also set up once and for all by this subroutine, which is called in physiq_mod's firstcall.
contains_mons, alpha_hmons and hsummit are now saved variables of the module surfdat_h, and are called as such and not as arguments in the subroutines using them.
the logical flag "slpwind" and the comments in the code have also been updated to the new terminology "mountain top dust flows", accordingly to ticket #71. The new flag read in callphys.def is "topflows".

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

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

Download in other formats:

Original Format