Context Navigation

rocketduststorm_mod.F90 @ 2934

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

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

Download in other formats:

Original Format