Context Navigation

rocketduststorm_mod.F90 @ 1974

Last change on this file since 1974 was 1974, checked in by mvals, 6 years ago
Mars GCM: Integration of the detached dust layer parametrizations (rocket dust storm, slope wind lifting, CW, and dust injection scheme, DB). Still experimental, default behaviour (rdstorm=.false., dustinjection=0) identical to previous revision. NB: Updated newstart requires an updated "surface.nc" containing the "hmons" field. EM+MV
File size: 44.4 KB

Line
1	MODULE rocketduststorm_mod
2
3	IMPLICIT NONE
4
5	REAL, SAVE, ALLOCATABLE :: dustliftday(:) ! dust lifting rate (s-1)
6	REAL, SAVE, ALLOCATABLE :: alpha_hmons(:) ! slope winds lifting mesh fraction
7
8	CONTAINS
9
10	!=======================================================================
11	! ROCKET DUST STORM - vertical transport and detrainment
12	!=======================================================================
13	! calculating the vertical flux
14	! calling vl_storm : transport scheme of stormdust tracers
15	! detrainement of stormdust into nomal background dust
16	! -----------------------------------------------------------------------
17	! Authors: C. Wang; F. Forget; T. Bertrand
18	! Institution: Laboratoire de Meteorologie Dynamique (LMD) Paris, France
19	! -----------------------------------------------------------------------
20
21	SUBROUTINE rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep, &
22	pq,pdqfi,pt,pdtfi,pplev,pplay,pzlev, &
23	pzlay,pdtsw,pdtlw, &
24	! input for radiative transfer
25	clearatm,icount,zday,zls, &
26	tsurf,igout,totstormfract, &
27	! input sub-grid scale cloud
28	clearsky,totcloudfrac, &
29	! output
30	pdqrds,wspeed,dsodust,dsords)
31
32	use tracer_mod, only: igcm_stormdust_mass,igcm_stormdust_number &
33	,igcm_dust_mass,igcm_dust_number &
34	,rho_dust
35	USE comcstfi_h, only: r,g,cpp,rcp
36	use dimradmars_mod, only: albedo,naerkind
37	use comsaison_h, only: dist_sol,mu0,fract
38	use surfdat_h, only: emis,co2ice,zmea, zstd, zsig, hmons
39	use planetwide_mod, only: planetwide_maxval,planetwide_minval
40	! use rocketstorm_h, only: rdsinject
41	use callradite_mod
42	implicit none
43
44	!--------------------------------------------------------
45	! Input Variables
46	!--------------------------------------------------------
47
48	INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points
49	INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points
50	INTEGER, INTENT(IN) :: nq ! number of tracer species
51	REAL, INTENT(IN) :: ptime
52	REAL, INTENT(IN) :: ptimestep
53
54	REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq
55	REAL, INTENT(IN) :: pdqfi(ngrid,nlayer,nq)! tendancy field pq
56	REAL, INTENT(IN) :: pt(ngrid,nlayer) ! temperature at mid-layer (K)
57	REAL, INTENT(IN) :: pdtfi(ngrid,nlayer) ! tendancy temperature at mid-layer (K)
58
59	REAL, INTENT(IN) :: pplay(ngrid,nlayer) ! pressure at middle of the layers
60	REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure at intermediate levels
61	REAL, INTENT(IN) :: pzlay(ngrid,nlayer) ! altitude at the middle of the layers
62	REAL, INTENT(IN) :: pzlev(ngrid,nlayer+1) ! altitude at layer boundaries
63
64	REAL, INTENT(IN) :: pdtsw(ngrid,nlayer) ! (K/s) env
65	REAL, INTENT(IN) :: pdtlw(ngrid,nlayer) ! (K/s) env
66
67	! input for second radiative transfer
68	logical, INTENT(IN) :: clearatm
69	INTEGER, INTENT(INOUT) :: icount
70	real, intent(in) :: zday
71	real, intent(in) :: zls
72	real, intent(in) :: tsurf(ngrid)
73	integer, intent(in) :: igout
74	real, intent(in) :: totstormfract(ngrid)
75	! sbgrid scale water ice clouds
76	logical, intent(in) :: clearsky
77	real, intent(in) :: totcloudfrac(ngrid)
78
79	!--------------------------------------------------------
80	! Output Variables
81	!--------------------------------------------------------
82
83	REAL, INTENT(OUT) :: pdqrds(ngrid,nlayer,nq) ! tendancy field for dust when detraining
84	REAL, INTENT(OUT) :: wspeed(ngrid,nlayer+1) ! vertical speed within the rocket dust storm
85	REAL, INTENT(OUT) :: dsodust(ngrid,nlayer) ! density scaled opacity of env. dust
86	REAL, INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity of storm dust
87
88	!--------------------------------------------------------
89	! Local variables
90	!--------------------------------------------------------
91	INTEGER l,ig,tsub,iq,ll
92	! chao local variables from callradite.F
93	REAL zdtlw1(ngrid,nlayer) ! (K/s) storm
94	REAL zdtsw1(ngrid,nlayer) ! (K/s) storm
95	REAL zt(ngrid,nlayer) ! actual temperature at mid-layer (K)
96	REAL zdtvert(nlayer) ! dT/dz , lapse rate
97	REAL ztlev(nlayer) ! temperature at intermediate levels l+1/2 without last level
98
99	REAL zdtlw1_lev(nlayer),zdtsw1_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2
100	REAL zdtlw_lev(nlayer),zdtsw_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2
101
102	REAL zqstorm_mass(ngrid,nlayer) ! tracer pq mass intermediate
103	REAL zqstorm_mass_col(nlayer)
104	REAL zqstorm_number(ngrid,nlayer) ! tracer field pq number intermediate
105	REAL zqstorm_number_col(nlayer)
106
107	REAL zqi_mass(ngrid,nlayer) ! tracer pq mass intermediate
108	REAL zqi_number(ngrid,nlayer) ! tracer pq mass intermediate
109	REAL zdqvlstorm_mass(ngrid,nlayer) ! tendancy pdq mass after vertical transport
110	REAL zdqvlstorm_number(ngrid,nlayer) ! tendancy pdq number after vertical transport
111	REAL zdqdetstorm_mass(ngrid,nlayer) ! tendancy field pq mass after detrainment only
112	REAL zdqdetstorm_number(ngrid,nlayer) ! tendancy field pq number after detrainment only
113
114	REAL zdqenv_mass(ngrid,nlayer) ! tendancy pdq mass from dust->
115	! stormdust in slp
116	REAL zdqenv_number(ngrid,nlayer) ! tendancy pdq number from dust->
117	! stormdust in slp
118
119	REAL masse(nlayer)
120	REAL zq(ngrid,nlayer,nq) ! updated tracers
121
122	REAL wrds(nlayer) ! vertical flux within the rocket dust storm
123	REAL wqrdsmass(nlayer+1) ! flux mass from vl_storm
124	REAL wqrdsnumber(nlayer+1) ! flux number from vl_storm
125
126	INTEGER nsubtimestep !number of subtimestep when calling vl_storm
127	REAL subtimestep !ptimestep/nsubtimestep
128	REAL dtmax !considered time needed for dust to cross one layer
129	!minimal value over a column
130	logical storm(ngrid) !logical : true if you have some storm dust in the column
131	! real slope(ngrid) !logical : true if you don't have storm and have
132	!a slope
133	! real wslplev(ngrid,nlayer)
134	! real wslp(ngrid,nlayer)
135	REAL coefdetrain !coefficient for detrainment : % of stormdust detrained
136
137	REAL,PARAMETER:: coefmin =0.025 !C 0<c<1 Minimum fraction of stormdust detrained
138	REAL,PARAMETER:: detrainlim =0.1!0.25 !L stormdust detrained if wspeed < detrainlim
139	REAL,PARAMETER:: wlim =10. ! maximum vertical speed of rocket storms (m/s)
140	REAL,PARAMETER:: secu=3. !coefficient on wspeed to avoid dust crossing many layers during subtimestep
141
142	! terms for buoyancy and W^2 in equation:
143	! wdw/dz = k1g(T'-T)/T - k2w^2
144	real,parameter:: k1=0.00001
145	real,parameter:: k2=0.25
146	real,parameter:: mu0lim=0.1
147
148	! diagnose
149	REAL detrainment(ngrid,nlayer)
150	real lapserate(ngrid,nlayer)
151	real deltahr(ngrid,nlayer+1)
152	real stormdust_m0(ngrid,nlayer)
153	real stormdust_m1(ngrid,nlayer)
154	real stormdust_m2(ngrid,nlayer)
155
156	real hmax,hmin
157	! real zh(ngrid,nlayer)
158
159	logical,save :: firstcall=.true.
160	real alpha(ngrid) ! scale of the vertical motion (applicable for
161	! rds and slp
162	! variables for radiative transfer
163	REAL fluxsurf_lw1(ngrid)
164	REAL fluxsurf_sw1(ngrid,2)
165	REAL fluxtop_lw1(ngrid)
166	REAL fluxtop_sw1(ngrid,2)
167	REAL tauref(ngrid)
168	REAL tau(ngrid,naerkind)
169	REAL aerosol(ngrid,nlayer,naerkind)
170	REAL tauscaling(ngrid)
171	REAL taucloudtes(ngrid)
172	REAL rdust(ngrid,nlayer)
173	REAL rstormdust(ngrid,nlayer)
174	REAL rice(ngrid,nlayer)
175	REAL nuice(ngrid,nlayer)
176
177	!variables related to slope,reference layer...
178	! integer lref(ngrid),llref ! the reference layer of slopewind
179	! real buoyt(nlayer) ! buoyancy term when there is a subgrid mountain
180	real slpbg(ngrid) ! temperature difference at half height of a mountain
181
182	real zqdustslp(ngrid,nlayer),zndustslp(ngrid,nlayer)
183	real zqstormdustslp(ngrid,nlayer),znstormdustslp(ngrid,nlayer)
184
185	real rdsdustqvl0(ngrid,nlayer)
186	real rdsdustqvl1(ngrid,nlayer)
187	! real q2rds(ngrid,nlayer)
188
189	!for second formule of wslp
190	real wtemp(ngrid,nlayer) ! a intermediate variable for wspeed
191
192	!merge rds and slp
193	real newzt(ngrid,nlayer) !temperature with perturbation (integrated from
194	! vetical motion)
195	real w0(ngrid) !prescribed slope winds at the first layer of atmosphere
196	! real w1(ngrid) !prescribed slope winds at the first layer of atmosphere
197	! real ztb1(ngrid)
198	real wadiabatic(ngrid,nlayer) !for diagnosis
199	! real czt(nlayer),czlay(nlayer),czlev(nlayer+1) !temporary variables
200	real tnew(nlayer) !interpolated temperature profile above the top of Mons
201	real envtemp(nlayer) !interpolated background temperature profile
202	! as the same height as tnew
203	real envt(ngrid,nlayer) ! output,for diagnosing
204	integer scheme(ngrid) ! diagnose
205
206	! Chao: for checking conservation of dust
207	! real totdust0(ngrid)
208	! real totdust1(ngrid)
209
210
211	! **********************************************************************
212	! **********************************************************************
213	! Detached dust layers parametrization: two processes are included
214	! 1) rocket dust storm
215	! The radiative warming due to the presence of storm dust is
216	! balanced by the adiabatic cooling. The tracer "storm dust"
217	! is transported by the upward/downward flow.
218	! 2) daytime slope winds
219	! The daytime thermally driven upslope wind blows dust from the
220	! bottom to the top of the mountain, upward flow keeps rising untill
221	! the velocity becomes zero. Both the storm dust and environment dust
222	! will be transported.
223	! **********************************************************************
224	! **********************************************************************
225	!! 1. Radiative transfer in storm dust
226	!! 2. Compute vertical velocity for storm dust
227	!! case 1 storm = false and nighttime: nothing to do
228	!! case 2 daytime slope wind scheme: (mu0(ig) > mu0lim and with storm=false)
229	!! case 3 rocket dust storm (storm=true)
230	!! 3. Vertical transport
231	!! 4. Detrainment
232	! **********************************************************************
233	! **********************************************************************
234
235
236	!! **********************************************************************
237	!! Firstcall: Evaluate slope wind mesh fraction
238	IF (firstcall) then
239	call planetwide_maxval(hmons,hmax )
240	call planetwide_minval(hmons,hmin )
241	do ig=1,ngrid
242	! It's hard to know the exact the scale of upward flow,
243	! we assume that the maximun is 10% of the mesh area.
244	alpha_hmons(ig)= 0.1*(hmons(ig)-hmin)/(hmax-hmin)
245	enddo
246	firstcall = .false.
247	ENDIF !firstcall
248
249	! **********************************************************************
250	! 0. Initializations
251	! **********************************************************************
252	storm(:)=.false.
253	pdqrds(:,:,:) = 0.
254	zdqdetstorm_mass(:,:)=0.
255	zdqdetstorm_number(:,:)=0.
256	wspeed(:,:)=0.
257	detrainment(:,:)=0.
258	zqstorm_mass_col(:)=0.
259	zqstorm_number_col(:)=0.
260	lapserate(:,:)=0.
261	deltahr(:,:)=0.
262	rdsdustqvl0(:,:)=0.
263	rdsdustqvl1(:,:)=0.
264	zqstormdustslp(:,:)=0.
265	znstormdustslp(:,:)=0.
266	zqdustslp(:,:)=0.
267	zndustslp(:,:)=0.
268	zq(:,:,:) = 0.
269
270	w0(:)=0.
271	! w1(:)=0.
272	! ztb1(:)=0.
273	newzt(:,:)=0
274	wtemp(:,:)=0.
275	wadiabatic(:,:)=0.
276	slpbg(:)=0.
277	! buoyt(:)=0.
278	tnew(:)=0.
279	envtemp(:)=0.
280	envt(:,:)=0.
281	scheme(:)=0
282	alpha(:)=0.
283	stormdust_m0(:,:)=0.
284	stormdust_m1(:,:)=0.
285	stormdust_m2(:,:)=0.
286	! totdust0(:)=0.
287	! totdust1(:)=0.
288
289	!! no update for the stormdust tracer and temperature fields
290	!! because previous callradite was for background dust
291	zq(1:ngrid,1:nlayer,1:nq)=pq(1:ngrid,1:nlayer,1:nq)
292	zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer)
293
294	!! get actual q for stormdust and dust tracers
295	do l=1,nlayer
296	do ig=1, ngrid
297	zqi_mass(ig,l)=zq(ig,l,igcm_dust_mass)
298	zqi_number(ig,l)=zq(ig,l,igcm_dust_number)
299	zqstorm_mass(ig,l)=zq(ig,l,igcm_stormdust_mass)
300	zqstorm_number(ig,l)=zq(ig,l,igcm_stormdust_number)
301	!for diagnostics:
302	stormdust_m0(ig,l)=zqstorm_mass(ig,l)
303	enddo
304	enddo ! of do l=1,nlayer
305
306	!! Check if there is a rocket dust storm
307	do ig=1,ngrid
308	storm(ig)=.false.
309	do l=1,nlayer
310	if (zqstorm_mass(ig,l)/zqi_mass(ig,l) .gt. 1.E-4) then
311	storm(ig)=.true.
312	exit
313	endif
314	enddo
315	enddo
316
317	! *********************************************************************
318	! 1. Call the second radiative transfer for stormdust, obtain the extra heating
319	! *********************************************************************
320	CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, &
321	emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, &
322	zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
323	fluxtop_sw1,tauref,tau,aerosol,dsodust,tauscaling, &
324	taucloudtes,rdust,rice,nuice,co2ice,rstormdust, &
325	totstormfract,clearatm,dsords, &
326	clearsky,totcloudfrac)
327
328	! **********************************************************************
329	! 2. Compute vertical velocity for storm dust
330	! **********************************************************************
331	DO ig=1,ngrid
332	!! **********************************************************************
333	!! 2.1 case 1: Nothing to do when no storm and no slope
334	! IF ((mu0(ig) .LE. mu0lim) .AND. .NOT.(storm(ig)) ) then
335	! scheme(ig)=1
336	! cycle
337	! endif
338	IF ((alpha_hmons(ig) .EQ. 0.) .AND. .NOT.(storm(ig))) then
339	scheme(ig)=1
340	cycle !!no slope
341	endif
342
343
344	! whatever the situation is, we need the vertical velocity computed by
345	! the rds scheme
346	zdtvert(1)=0. !This is the env. lapse rate
347	DO l=1,nlayer-1
348	zdtvert(l+1)=(zt(ig,l+1)-zt(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
349	lapserate(ig,l+1)=zdtvert(l+1) !for diagnosing
350	ENDDO
351
352	! computing heating rates gradient at boundraies of each layer
353	! start from surface
354	zdtlw1_lev(1)=0.
355	zdtsw1_lev(1)=0.
356	zdtlw_lev(1)=0.
357	zdtsw_lev(1)=0.
358	ztlev(1)=zt(ig,1)
359
360	DO l=1,nlayer-1
361	! Calculation for the dust storm fraction
362	zdtlw1_lev(l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
363	zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
364	(pzlay(ig,l+1)-pzlay(ig,l))
365
366	zdtsw1_lev(l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
367	zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
368	(pzlay(ig,l+1)-pzlay(ig,l))
369	!MV18: calculation for the background dust fraction
370	zdtlw_lev(l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
371	pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
372	(pzlay(ig,l+1)-pzlay(ig,l))
373
374	zdtsw_lev(l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
375	pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
376	(pzlay(ig,l+1)-pzlay(ig,l))
377
378	ztlev(l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
379	zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / &
380	(pzlay(ig,l+1)-pzlay(ig,l))
381	ENDDO
382
383	DO l=1,nlayer
384	deltahr(ig,l)=(zdtlw1_lev(l)+zdtsw1_lev(l)) &
385	-(zdtlw_lev(l)+zdtsw_lev(l))
386	wadiabatic(ig,l)=-deltahr(ig,l)/(g/cpp+ &
387	max(zdtvert(l),-0.99*g/cpp))
388
389	!limit vertical wind in case of lapse rate close to adiabat
390	wadiabatic(ig,l)=max(wadiabatic(ig,l),-wlim)
391	wadiabatic(ig,l)=min(wadiabatic(ig,l),wlim)
392	ENDDO
393
394	!! **********************************************************************
395	!! 2.2 case 2: daytime slope wind scheme
396	IF ((mu0(ig) .gt. mu0lim) .and. .not. storm(ig)) then
397	scheme(ig)=2
398	alpha(ig) = alpha_hmons(ig)
399	! interpolate the env. temperature above the mountain top
400	call intep_vtemp(nlayer,hmons(ig),zt(ig,:),pzlay(ig,:), &
401	envtemp,slpbg(ig))
402	envt(ig,:)=envtemp(:) !for diagnosis
403
404	! second: estimate the vertical velocity at boundraies of each layer
405	!wspeed(ig,1)=0. ! at surface, already initialized
406
407	!!!!!!!the first layer of atmosphere!!!!!!!!!!
408	IF (slpbg(ig) .gt. 0.) THEN !only postive buoyancy
409	! if slpbg(ig) lt 0, means the slope flow is colder than env. (night or
410	! early morning ?) !!!!!!!!!!!
411	! ideal method
412	!w1(ig)=-sqrt(gslpbg(ig)/zt(ig,1)hmons(ig))*sin(zsig(ig))
413	! new scheme, simply proportional to temperature and
414	! mountain height
415	w0(ig)=-1.5e-4gslpbg(ig)/zt(ig,1)*hmons(ig)
416	! otherwise, w0(ig) =0.
417	wspeed(ig,2)=w0(ig)
418	ELSE
419	wspeed(ig,2)=wadiabatic(ig,2) !! for early morning ?
420	ENDIF
421
422	! prepare the integration, NOTE: if w is too small, may have artificials
423	IF (abs(wspeed(ig,2)) .lt. 0.01 ) &
424	wspeed(ig,2)=sign(0.01,wspeed(ig,2))
425
426	newzt(ig,1)= zt(ig,1) !temperature of the first layer atmosphere
427	! above the mountain top (radiative
428	! equilibrium on Mars)
429
430	! estimate the vertical velocities
431	! if w0(ig) >= 0, means downward motion, no upslope winds, we switch to
432	! assume that the extra heating integrally convert to
433	! vertical motion.
434	if ( w0(ig) .ge. 0 ) then !! normal, it is impossible,
435	!! because mu(ig)>0.1 here
436	do l=3,nlayer
437	wspeed(ig,l)=wadiabatic(ig,l)
438	enddo
439	else
440	! estimate the velocities by taking into account the heating due
441	! to storm dust, the cooling due to vertical motion ...
442	!!!!!!!!!!!the simple scheme!!!!!!!!!
443	do l=2,nlayer-1
444	!if superadiabatic layer
445	if ( zdtvert(l) .lt. -g/cpp) then
446	!case 1
447	! test, also decrease adiabatically ?
448	!newzt(ig,l)= &
449	! zt(ig,l-1)-g/cpp*(pzlay(ig,l)-pzlay(ig,l-1))
450	newzt(ig,l)=zt(ig,l)
451	!wspeed(ig,l+1)=wspeed(ig,l)
452	else
453	!not superadiabatic
454	newzt(ig,l)=newzt(ig,l-1)+(deltahr(ig,l)/ &
455	(-wspeed(ig,l))-g/cpp)* &
456	(pzlay(ig,l)-pzlay(ig,l-1))
457	endif ! end of if superadiabatic or not
458
459	!wtemp(ig,l+1)=wspeed(ig,l)*2+2.g*(pzlev(ig,l+1) &
460	! -pzlev(ig,l))(k1(newzt(ig,l) &
461	! -envtemp(l))/envtemp(l))
462	wtemp(ig,l+1)=(1.-2.k1(pzlev(ig,l+1)-pzlev(ig,l)))*&
463	wspeed(ig,l)*2+2.k2g(pzlev(ig,l+1) &
464	-pzlev(ig,l))*((newzt(ig,l) &
465	-envtemp(l))/envtemp(l))
466
467	if (wtemp(ig,l+1) .gt. 0.) then
468	!case 2
469	wspeed(ig,l+1)=-sqrt(wtemp(ig,l+1))
470
471	! if \|wspeed\| < \|wadiabatic\| then go to wadiabatic
472	if (wspeed(ig,l+1) .gt. wadiabatic(ig,l+1)) then
473	do ll=l,nlayer-1
474	newzt(ig,ll)=envtemp(ll)
475	wspeed(ig,ll+1)=wadiabatic(ig,ll+1)
476	enddo
477	exit
478	endif
479
480	! avoid artificials
481	if (abs(wspeed(ig,l+1)) .lt. 0.01 ) &
482	wspeed(ig,l+1)=sign(0.01,wspeed(ig,l+1))
483
484	else if (l .lt. nlayer) then
485	!case 3
486	do ll=l,nlayer-1
487	newzt(ig,ll)=envtemp(ll)
488	wspeed(ig,ll+1)=wadiabatic(ig,ll+1)
489	enddo !overshot
490	exit
491
492	else
493	wspeed(ig,l+1)=0.
494	endif
495
496	enddo
497
498	endif !w0
499
500	ELSE
501	!! **********************************************************************
502	!! 2.3 case 3: storm=true
503	if (storm(ig)) then
504	scheme(ig)=3
505	alpha(ig) = totstormfract(ig)
506	do l=1,nlayer
507	wspeed(ig,l)=wadiabatic(ig,l)
508	enddo
509	endif ! storm=1
510
511	ENDIF ! rds or slp
512
513
514	!!!!!!!! estimate the amount of dust for diagnostics
515	DO l=1,nlayer
516	! transfer background dust + storm dust (concentrated)
517	zqstormdustslp(ig,l) =zqi_mass(ig,l)+ &
518	zqstorm_mass(ig,l)/alpha(ig)
519	znstormdustslp(ig,l) =zqi_number(ig,l)+ &
520	zqstorm_number(ig,l)/alpha(ig)
521	zqdustslp(ig,l) = zqi_mass(ig,l)
522	zndustslp(ig,l) = zqi_number(ig,l)
523
524	rdsdustqvl0(ig,l)=zqstormdustslp(ig,l) !for diagnosis
525	ENDDO
526
527	! **********************************************************************
528	! 3. Vertical transport
529	! **********************************************************************
530	do l=1,nlayer
531	masse(l)=(pplev(ig,l)-pplev(ig,l+1))/g
532	enddo
533	!Estimation of "dtmax" (s) to be used for vertical transport
534	dtmax=ptimestep
535	!secu is a margin on subtimstep to avoid dust crossing many layers
536	do l=2,nlayer
537	if (wspeed(ig,l).lt.0.) then ! case up
538	dtmax=min(dtmax,(pzlev(ig,l)-pzlev(ig,l-1))/ &
539	(secu*abs(wspeed(ig,l))))
540	else if (wspeed(ig,l).gt.0.) then
541	dtmax=min(dtmax,(pzlev(ig,l+1)-pzlev(ig,l))/ &
542	(secu*abs(wspeed(ig,l))))
543	endif
544	enddo
545
546	nsubtimestep= int(ptimestep/dtmax)
547	subtimestep=ptimestep/float(nsubtimestep)
548
549	do l=1,nlayer
550	wrds(l)=wspeed(ig,l)pplev(ig,l)/(rztlev(l)) &
551	*subtimestep
552	enddo
553
554	do l=1,nlayer
555	zqstorm_mass_col(l)= zqstormdustslp(ig,l) !zqstorm_mass(ig,l)
556	zqstorm_number_col(l)=znstormdustslp(ig,l) ! zqstorm_number(ig,l)
557	enddo
558
559	do tsub=1,nsubtimestep
560	wqrdsmass(:)=0.
561	wqrdsnumber(:)=0.
562	CALL vl_storm(nlayer,zqstorm_mass_col,2., &
563	masse,wrds ,wqrdsmass)
564	CALL vl_storm(nlayer,zqstorm_number_col,2., &
565	masse,wrds ,wqrdsnumber)
566	enddo
567
568	!!!!!generate the "extra" dust
569	do l=1,nlayer
570	rdsdustqvl1(ig,l)=zqstorm_mass_col(l) ! for diagnosis
571
572	! extra dust = storm dust
573	!zqdustslp(ig,l)=zqi_mass(ig,l) !(1.-alpha(ig))*zqi_mass(ig,l)
574	!zndustslp(ig,l)=zqi_number(ig,l) !(1.-alpha(ig))*zqi_number(ig,l)
575	!zqstorm_mass_col(l)=alpha(ig)*zqstorm_mass_col(l)
576	!zqstorm_number_col(l)=alpha(ig)*zqstorm_number_col(l)
577
578	!with compensation
579	if (zqstorm_mass_col(l) .lt. zqi_mass(ig,l) ) then
580	! the following two equations are easier to understand
581	zqdustslp(ig,l)=(1.-alpha(ig))zqi_mass(ig,l)+alpha(ig) &
582	zqstorm_mass_col(l)
583	zndustslp(ig,l)=(1.-alpha(ig))*zqi_number(ig,l)+alpha(ig)&
584	*zqstorm_number_col(l)
585	!with a bug, should be zqi+alpha****
586	!zqdustslp(ig,l)=zqi_mass(ig,l)-alpha(ig)* &
587	! (zqstorm_mass_col(l)-zqi_mass(ig,l))
588	!zndustslp(ig,l)=zqi_number(ig,l)-alpha(ig)* &
589	! (zqstorm_number_col(l)-zqi_number(ig,l))
590	zqstorm_mass_col(l)=0.
591	zqstorm_number_col(l)=0.
592	else
593	zqstorm_mass_col(l)=alpha(ig)* &
594	(zqstorm_mass_col(l)-zqi_mass(ig,l))
595	zqstorm_number_col(l)=alpha(ig)* &
596	(zqstorm_number_col(l)-zqi_number(ig,l))
597	! the mass mixing ratio of environmental dust doesn't change.
598	endif
599	!diagnose
600	stormdust_m1(ig,l)=zqstorm_mass_col(l)
601	enddo
602
603	!=======================================================================
604	! calculate the tendencies due to vertical transport
605	do l=1,nlayer
606	! tendencies due to vertical transport
607	zdqvlstorm_mass(ig,l)= (zqstorm_mass_col(l)- &
608	zqstorm_mass(ig,l)) /ptimestep
609	zdqvlstorm_number(ig,l)= (zqstorm_number_col(l)- &
610	zqstorm_number(ig,l)) /ptimestep
611
612	zdqenv_mass(ig,l)=(zqdustslp(ig,l)-zqi_mass(ig,l))/ptimestep
613	zdqenv_number(ig,l)=(zndustslp(ig,l)-zqi_number(ig,l)) &
614	/ptimestep
615
616	! chao for output only
617	!qstormdustvl1(ig,l)=zqstorm_mass_col(l)
618	!nstormdustvl1(ig,l)=zqstorm_number_col(l)
619	!stormdust_m_col1(ig)=stormdust_m_col1(ig)+zqstorm_mass_col(l) &
620	! *(pplev(ig,l)-pplev(ig,l+1))/g
621	!rdsdustqvl1(ig,l)=zqstorm_mass_col(l)
622	enddo
623
624	! **********************************************************************
625	! 4. Detrainment: convert dust storm to background dust
626	! **********************************************************************
627	do l=1,nlayer
628	! compute the coefficient of detrainment
629	if ((max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))) .lt. &
630	detrainlim) .or. (zqdustslp(ig,l) .gt. &
631	10000.*zqstorm_mass_col(l))) then
632	coefdetrain=1.
633	else if (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))) &
634	.le. wlim) then
635	! case where detrainment depends on vertical wind
636	! coefdetrain=0.5(((1-coefmin)/(detrainlim-3.)2) &
637	! (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1)))-3.)**2 &
638	! +coefmin)
639	coefdetrain=1.(((1-coefmin)/(detrainlim-wlim)2) &
640	(max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1)))-wlim)**2 &
641	+coefmin)
642	!coefdetrain=0.5
643	else if (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))).gt. 10. )&
644	then
645	coefdetrain=0.025
646	else
647	coefdetrain=coefmin
648	endif
649
650	detrainment(ig,l)=coefdetrain !for diagnose
651
652	! Calculate tendancies corresponding to pdq after detrainement
653	! pdqdet = tendancy corresponding to detrainment only
654	zdqdetstorm_mass(ig,l)=-coefdetrain*zqstorm_mass_col(l) &
655	/ptimestep
656	zdqdetstorm_number(ig,l)=-coefdetrain*zqstorm_number_col(l) &
657	/ptimestep
658
659	! pdqrds ( tendancy corresponding to vertical transport and
660	! detrainment) = zdqvlstorm + pdqdet
661	pdqrds(ig,l,igcm_stormdust_mass)=zdqdetstorm_mass(ig,l) &
662	+zdqvlstorm_mass(ig,l)
663	pdqrds(ig,l,igcm_stormdust_number)=zdqdetstorm_number(ig,l) &
664	+zdqvlstorm_number(ig,l)
665	pdqrds(ig,l,igcm_dust_mass)= zdqenv_mass(ig,l) &
666	-zdqdetstorm_mass(ig,l)
667	pdqrds(ig,l,igcm_dust_number)= zdqenv_number(ig,l) &
668	-zdqdetstorm_number(ig,l)
669
670	!diagnose
671	stormdust_m2(ig,l)=zqstorm_mass_col(l)-coefdetrain*zqstorm_mass_col(l)
672	enddo ! nlayer
673	! endif
674	!=======================================================================
675	enddo ! end do ig=1,ngrid
676
677	! !chao check conservation here
678	! do l=1,nlayer
679	! do ig=1,ngrid
680	! totdust0(ig)=totdust0(ig)+ &
681	! zq(ig,l,igcm_stormdust_mass)* &
682	! ((pplev(ig,l) - pplev(ig,l+1)) / g) &
683	! + zq(ig,l,igcm_dust_mass)* &
684	! ((pplev(ig,l) - pplev(ig,l+1)) / g)
685
686	! totdust1(ig)=totdust1(ig)+ &
687	! (zq(ig,l,igcm_stormdust_mass) + &
688	! pdqrds(ig,l,igcm_stormdust_mass)ptimestep) &
689	! ((pplev(ig,l) - pplev(ig,l+1)) / g) &
690	! + ( zq(ig,l,igcm_dust_mass)+ &
691	! pdqrds(ig,l,igcm_dust_mass)ptimestep) &
692	! ((pplev(ig,l) - pplev(ig,l+1)) / g)
693	! enddo
694	! enddo
695
696	! call writediagfi(ngrid,'totdust0','total dust before rds', &
697	! ' ',2,totdust0)
698	! call writediagfi(ngrid,'totdust1','total dust after rds', &
699	! ' ',2,totdust1)
700	!output for diagnosis
701	call WRITEDIAGFI(ngrid,'detrainment', &
702	'coefficient of detrainment',' ',3,detrainment)
703	!call WRITEDIAGFI(ngrid,'qstormvl1','mmr of stormdust after rds_vl', &
704	! & 'kg/kg',3,qstormdustvl1)
705	call WRITEDIAGFI(ngrid,'lapserate','lapse rate in the storm', &
706	& 'k/m',3,lapserate)
707	call WRITEDIAGFI(ngrid,'deltahr','extra heating rates', &
708	& 'K/s',3,deltahr)
709	call WRITEDIAGFI(ngrid,'wold', &
710	'wind generated from adiabatic colling ', &
711	& 'm/s',3,wadiabatic)
712	call WRITEDIAGFI(ngrid,'newzt','perturbated temperature', &
713	& 'K/s',3,newzt)
714	call WRITEDIAGFI(ngrid,'zt','unperturbated temperature', &
715	& 'K/s',3,zt)
716	call WRITEDIAGFI(ngrid,'wtemp','under square root', &
717	& 'K/s',3,wtemp)
718	!call WRITEDIAGFI(ngrid,'stormdust_m_col1','mass of stormdust after rds_vl', &
719	! & 'kg/kg',2,stormdust_m_col1)
720	!call WRITEDIAGFI(ngrid,'temprds','temp for calculating zdtvert', &
721	! & 'k',3,temprds)
722	call WRITEDIAGFI(ngrid,'stormdust_m0','mass col of stormdust before rds_vl', &
723	& 'kg/kg',3,stormdust_m0)
724	call WRITEDIAGFI(ngrid,'stormdust_m1','mass col of stormdust after rds_vl', &
725	& 'kg/kg',3,stormdust_m1)
726	call WRITEDIAGFI(ngrid,'stormdust_m2','mass col of stormdust after rds_vl', &
727	& 'kg/kg',3,stormdust_m2)
728
729	! call writediagfi(ngrid,'wslp','estimated slope winds','m/s',3,wslp)
730	! call writediagfi(ngrid,'wslp2','estimated slope winds2','m/s',3,wslp2)
731	! call writediagfi(ngrid,'zhb','estimated slope winds2','m/s',3,zhb)
732	! call writediagfi(ngrid,'bruntf','bouyancy frequency',' ',3,bruntf)
733	! call writediagfi(ngrid,'slpdepth','slope depth','m',2,slpdepth)
734	! call writediagfi(ngrid,'slpu','perbulation u','m/s',3,slpu)
735	! call writediagfi(ngrid,'slpzh','perbulation zh',' ',3,slpzh)
736	! call writediagfi(ngrid,'zqslp','zq in rocketduststorm','ikg/kg',3, &
737	! zq(:,:,igcm_dust_mass))
738	! call writediagfi(ngrid,'zrdsqslp','zq in rocketduststorm','ikg/kg',3, &
739	! zq(:,:,igcm_stormdust_mass))
740	! call writediagfi(ngrid,'wslplev','estimated slope winds','m/s',3,wslplev)
741	! call writediagfi(ngrid,'slope','identified slope wind effect',' ',2,slope)
742	call writediagfi(ngrid,'w0','max of slope wind',' ',2,w0)
743	! call writediagfi(ngrid,'w1','max of slope wind',' ',2,w1)
744	call writediagfi(ngrid,'mu0','cosine of solar incidence angle',&
745	' ',2,mu0)
746	! call writediagfi(ngrid,'storm','identified rocket dust storm',&
747	! ' ',2,real(storm))
748	call writediagfi(ngrid,'scheme','which scheme',&
749	' ',2,real(scheme))
750	call writediagfi(ngrid,'alpha','coefficient alpha',' ',2,alpha)
751	! call writediagfi(ngrid,'q2rds','alpha zq',' ', &
752	! 3,q2rds)
753	call writediagfi(ngrid,'rdsdustqvl0','not vl storm slp', &
754	'kg/kg',3,zqstormdustslp)
755	call writediagfi(ngrid,'rdsdustqvl1','vled storm slp', &
756	'kg/kg',3,rdsdustqvl1)
757	call writediagfi(ngrid,'dustqvl0','not vl slp', &
758	'kg/kg',3,zqi_mass)
759	call writediagfi(ngrid,'dustqvl1','vled slp', &
760	'kg/kg',3,zqdustslp)
761	! call WRITEDIAGFI(ngrid,'lmax_th2', &
762	! 'hauteur du thermique','point', &
763	! 2,real(lmax_th(:)))
764	! call WRITEDIAGFI(ngrid,'zmax_th2', &
765	! 'hauteur du thermique','m', &
766	! 2,zmax_th)
767	! call writediagfi(ngrid,'lslope','lenght of slope',' ',2,lslope)
768	! call writediagfi(ngrid,'hmon','identified slope wind effect',' ',2,hmon)
769	call writediagfi(ngrid,'envt','interpolated env. temp.', &
770	'K',3,envt)
771	call writediagfi(ngrid,'hmons','identified slope wind effect', &
772	' ',2,hmons)
773	! call writediagfi(ngrid,'slpbg','temp. diff along slope', &
774	! ' ',2,slpbg)
775	! call writediagfi(ngrid,'zmea','identified slope wind effect', &
776	! ' ',2,zmea)
777	! call writediagfi(ngrid,'zsig','identified slope wind effect', &
778	! ' ',2,zsig)
779	! call writediagfi(ngrid,'zhslpenv','difference of zh above mons',' ',2,zhslpenv)
780	! call writediagfi(ngrid,'lref','identified slope wind effect',' ',2,real(lref))
781	END SUBROUTINE rocketduststorm
782
783	!********************************************************************************
784	!********************************************************************************
785	SUBROUTINE vl_storm(nlay,q,pente_max,masse,w,wq)
786	!
787	! Auteurs: P.Le Van, F.Hourdin, F.Forget
788	!
789	! ********************************************************************
790	! Shema d'advection " pseudo amont " dans la verticale
791	! pour appel dans la physique (sedimentation)
792	! ********************************************************************
793	! q rapport de melange (kg/kg)...
794	! masse : masse de la couche Dp/g
795	! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
796	! pente_max = 2 conseillee
797	!
798	!
799	! --------------------------------------------------------------------
800	IMPLICIT NONE
801	!
802
803	! Arguments:
804	! ----------
805	integer,intent(in) :: nlay ! number of atmospheric layers
806	real masse(nlay),pente_max
807	REAL q(nlay)
808	REAL w(nlay)
809	REAL wq(nlay+1)
810	!
811	! Local
812	! ---------
813	!
814	INTEGER i,l,j,ii
815	!
816
817	real dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax
818	real newmasse
819	real sigw, Mtot, MQtot
820	integer m
821
822	REAL SSUM,CVMGP,CVMGT
823	integer ismax,ismin
824
825
826	! On oriente tout dans le sens de la pression c'est a dire dans le
827	! sens de W
828
829	do l=2,nlay
830	dzqw(l)=q(l-1)-q(l)
831	adzqw(l)=abs(dzqw(l))
832	enddo
833
834	do l=2,nlay-1
835	#ifdef CRAY
836	dzq(l)=0.5*
837	, cvmgp(dzqw(l)+dzqw(l+1),0.,dzqw(l)*dzqw(l+1))
838	#else
839	if(dzqw(l)*dzqw(l+1).gt.0.) then
840	dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
841	else
842	dzq(l)=0.
843	endif
844	#endif
845	dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
846	dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
847	enddo
848
849	dzq(1)=0.
850	dzq(nlay)=0.
851
852	! write(,),'TB14 wq before up',wq(1,:)
853	! write(,),'TB14 q before up',q(1,:)
854	! ---------------------------------------------------------------
855	! .... calcul des termes d'advection verticale .......
856	! ---------------------------------------------------------------
857
858	! calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq
859	!
860	! No flux at the model top:
861	wq(nlay+1)=0.
862
863	! 1) Compute wq where w < 0 (up) (NOT USEFUL FOR SEDIMENTATION)
864	! ===============================
865
866	! Surface flux up:
867	if(w(1).lt.0.) wq(1)=0. ! warning : not always valid
868
869	do l = 1,nlay-1 ! loop different than when w>0
870	if(w(l+1).le.0)then
871
872	! Regular scheme (transfered mass < 1 layer)
873	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
874	if(-w(l+1).le.masse(l))then
875	sigw=w(l+1)/masse(l)
876	wq(l+1)=w(l+1)(q(l)-0.5(1.+sigw)*dzq(l))
877	! Extended scheme (transfered mass > 1 layer)
878	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
879	else
880	m = l-1
881	Mtot = masse(m+1)
882	MQtot = masse(m+1)*q(m+1)
883	if (m.le.0)goto 77
884	do while(-w(l+1).gt.(Mtot+masse(m)))
885	! do while(-w(l+1).gt.Mtot)
886	m=m-1
887	Mtot = Mtot + masse(m+1)
888	MQtot = MQtot + masse(m+1)*q(m+1)
889	if (m.le.0)goto 77
890	end do
891	77 continue
892
893	if (m.gt.0) then
894	sigw=(w(l+1)+Mtot)/masse(m)
895	wq(l+1)= (MQtot + (-w(l+1)-Mtot)* &
896	(q(m)-0.5(1.+sigw)dzq(m)) )
897	else
898	! new
899	w(l+1) = -Mtot
900	wq(l+1) = -MQtot
901	! write(,) 'TB14 MQtot = ',MQtot,l
902
903	! old
904	! wq(l+1)= (MQtot + (-w(l+1)-Mtot)*qm(1))
905	! write(,) 'a rather weird situation in vlz_fi !'
906	! stop
907	end if
908	endif
909	endif ! w<0 (up)
910	enddo
911
912	do l = 1,nlay-1 ! loop different than when w>0
913
914	q(l)=q(l) + (wq(l+1)-wq(l))/masse(l)
915
916	enddo
917
918	! write(,),'TB14 masse',masse(1,:)
919	! write(,),'TB14 wq before down after up',wq(1,:)
920	! write(,),'TB14 q before down',q(1,:)
921
922	! 2) Compute wq where w > 0 (down) (ALWAYS FOR SEDIMENTATION)
923	! ===============================
924
925	! Initialisation wq = 0 to consider now only downward flux
926	wq(:)=0. !
927
928	do l = 1,nlay ! loop different than when w<0
929
930	if(w(l).gt.0.)then
931
932	! Regular scheme (transfered mass < 1 layer)
933	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
934	if(w(l).le.masse(l))then
935	sigw=w(l)/masse(l)
936	wq(l)=w(l)(q(l)+0.5(1.-sigw)*dzq(l))
937	! write(,),'TB14 wq after up',wq(1,:)
938
939
940	! Extended scheme (transfered mass > 1 layer)
941	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
942	else
943	m=l
944	Mtot = masse(m)
945	MQtot = masse(m)*q(m)
946	if(m.ge.nlay)goto 88
947	do while(w(l).gt.(Mtot+masse(m+1)))
948	m=m+1
949	Mtot = Mtot + masse(m)
950	MQtot = MQtot + masse(m)*q(m)
951	if(m.ge.nlay)goto 88
952	end do
953	88 continue
954	if (m.lt.nlay) then
955	sigw=(w(l)-Mtot)/masse(m+1)
956	wq(l)=(MQtot + (w(l)-Mtot)* &
957	(q(m+1)+0.5(1.-sigw)dzq(m+1)) )
958	else
959	w(l) = Mtot
960	wq(l) = MQtot
961	end if
962	end if
963	end if ! w>0 (down)
964	enddo
965
966	do l = 1,nlay ! loop different than when w<0
967
968	q(l)=q(l) + (wq(l+1)-wq(l))/masse(l)
969
970	enddo
971	end subroutine vl_storm
972	!********************************************************************************
973	SUBROUTINE intep_vtemp(nlayer,hm,temp,zlay,envtemp,slpb)
974
975	USE comcstfi_h, only: g,cpp
976
977	implicit none
978	! this subroutine is using for obtaining the environmental
979	! temperature profile when a subgrid mountain exists.
980
981
982	! input:
983	integer,intent(in) :: nlayer
984	real,intent(in) :: hm ! the height of mountain
985	real,intent(in) :: temp(nlayer) !large scale temp. profile of one mesh
986	real,intent(in) :: zlay(nlayer) ! height at the middle of each layer
987
988	! output:
989	real,intent(out) :: envtemp(nlayer) ! environment temperature
990	real,intent(out) :: slpb !the temperature difference between slope and
991	! env. at the half height of mountain
992
993	! local variables
994	integer l,il,tmpl
995	integer lnew !the layer of atmosphere above the mountain
996	! corresponding to the env. (for buoyancy
997	! calc. )
998	real newh(nlayer) !height at the middle of each layer
999	! account for the exist of mountain
1000	! real g,cpp
1001	real halfh ! half the height of a mountain
1002
1003	!initilize the array
1004	lnew=0
1005	newh(:)=0.
1006	envtemp(:)=0.
1007	tmpl=1
1008
1009	do l=1,nlayer
1010	newh(l)=hm+zlay(l)
1011
1012	do il=tmpl,nlayer-1 !MV18: added the -1
1013	if (newh(l) .ge. zlay(il) .and. newh(l) .lt. zlay(il+1))then
1014	! find the corresponding layer
1015	lnew=il
1016
1017	! interpolate
1018	envtemp(l) = temp(il)+(newh(l)-zlay(lnew))*&
1019	(temp(il+1)-temp(il))/(zlay(il+1)-zlay(il))
1020
1021	exit !go to the next layer
1022	else if (newh(l) .ge. zlay(nlayer)) then
1023	! higher than the highest layer
1024	lnew=nlayer
1025	envtemp(l)=temp(nlayer)
1026
1027	endif
1028	enddo
1029	! this can accelerate the loop
1030	tmpl=lnew
1031
1032	enddo
1033
1034	halfh=0.5*hm
1035	if (halfh .le. zlay(1) ) then
1036	slpb=0.
1037	else if (halfh .gt. zlay(nlayer)) then
1038	!normally, impossible for halfh gt zlay(l), anyway...
1039	tmpl=nlayer
1040	!difference between surface and atmosphere (env.)
1041	slpb=temp(1)-(temp(nlayer-1)+((halfh-zlay(nlayer-1))* &
1042	(temp(tmpl)-temp(tmpl-1))/(zlay(tmpl)-zlay(tmpl-1))))
1043	else
1044	do l=1,nlayer-1
1045	if ((halfh .gt. zlay(l)) .and. (halfh .le. zlay(l+1)))then
1046	tmpl= l
1047	slpb=temp(1)-(temp(tmpl)+(halfh-zlay(tmpl))* &
1048	(temp(tmpl+1)-temp(tmpl))/(zlay(tmpl+1)-zlay(tmpl)))
1049	endif
1050	enddo
1051	endif
1052	end subroutine intep_vtemp
1053
1054	! initialization module variables
1055	subroutine ini_rocketduststorm_mod(ngrid)
1056
1057	implicit none
1058
1059	integer, intent(in) :: ngrid
1060
1061	allocate(dustliftday(ngrid))
1062	allocate(alpha_hmons(ngrid))
1063
1064	end subroutine ini_rocketduststorm_mod
1065
1066	subroutine end_rocketduststorm_mod
1067
1068	implicit none
1069
1070	if (allocated(dustliftday)) deallocate(dustliftday)
1071	if (allocated(alpha_hmons)) deallocate(alpha_hmons)
1072
1073	end subroutine end_rocketduststorm_mod
1074
1075	END MODULE rocketduststorm_mod

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

Download in other formats:

Original Format