Context Navigation

rocketduststorm_mod.F90 @ 2616

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

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

Download in other formats:

Original Format