Context Navigation

rocketduststorm_mod.F90 @ 2929

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

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

Download in other formats:

Original Format