Context Navigation

rocketduststorm_mod.F90 @ 2953

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

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

Download in other formats:

Original Format