Context Navigation

rocketduststorm_mod.F90 @ 2613

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

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

Download in other formats:

Original Format