Context Navigation

callsedim_mod.F @ 2323

Last change on this file since 2323 was 2323, checked in by mvals, 5 years ago
Mars GCM: Follow-up of the last commit for the transport of the isotopic ratio: forgot the initialisation of zq0() in callsedim_mod.F ! MV
File size: 25.0 KB

Line
1	MODULE callsedim_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	SUBROUTINE callsedim(ngrid,nlay,ptimestep,
8	& pplev,zlev,zlay,pt,pdt,
9	& rdust,rstormdust,rtopdust,
10	& rice,rsedcloud,rhocloud,
11	& pq,pdqfi,pdqsed,pdqs_sed,nq,
12	& tau,tauscaling)
13
14	USE ioipsl_getin_p_mod, only: getin_p
15	USE updaterad, only: updaterdust,updaterice_micro,updaterice_typ
16	USE tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
17	& rho_dust, rho_q, radius, varian,
18	& igcm_ccn_mass, igcm_ccn_number,
19	& igcm_h2o_ice, igcm_hdo_ice,
20	& nuice_sed, nuice_ref,
21	& igcm_ccnco2_mass,igcm_ccnco2_number,
22	& igcm_co2_ice, igcm_stormdust_mass,
23	& igcm_stormdust_number,igcm_topdust_mass,
24	& igcm_topdust_number,
25	& iqfils,nqfils,qperemin,masseqmin ! MVals: variables isotopes
26	USE newsedim_mod, ONLY: newsedim
27	USE comcstfi_h, ONLY: g
28	USE dimradmars_mod, only: naerkind
29	IMPLICIT NONE
30
31	c=======================================================================
32	c Sedimentation of the Martian aerosols
33	c depending on their density and radius
34	c
35	c F.Forget 1999
36	c
37	c Modified by J.-B. Madeleine 2010: Now includes the doubleq
38	c technique in order to have only one call to callsedim in
39	c physiq.F.
40	c
41	c Modified by J. Audouard 09/16: Now includes the co2clouds case
42	c If the co2 microphysics is on, then co2 theice & ccn tracers
43	c are being sedimented in the microtimestep (co2cloud.F), not
44	c in this routine.
45	c
46	c=======================================================================
47
48	c-----------------------------------------------------------------------
49	c declarations:
50	c -------------
51
52	include "callkeys.h"
53
54	c
55	c arguments:
56	c ----------
57
58	integer,intent(in) :: ngrid ! number of horizontal grid points
59	integer,intent(in) :: nlay ! number of atmospheric layers
60	real,intent(in) :: ptimestep ! physics time step (s)
61	real,intent(in) :: pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
62	real,intent(in) :: zlev(ngrid,nlay+1) ! altitude at layer boundaries
63	real,intent(in) :: zlay(ngrid,nlay) ! altitude at the middle of the layers
64	real,intent(in) :: pt(ngrid,nlay) ! temperature at mid-layer (K)
65	real,intent(in) :: pdt(ngrid,nlay) ! tendency on temperature, from
66	! previous processes (K/s)
67	c Aerosol radius provided by the water ice microphysical scheme:
68	real,intent(out) :: rdust(ngrid,nlay) ! Dust geometric mean radius (m)
69	real,intent(out) :: rstormdust(ngrid,nlay) ! Stormdust geometric mean radius (m)
70	real,intent(out) :: rtopdust(ngrid,nlay) ! topdust geometric mean radius (m)
71	real,intent(out) :: rice(ngrid,nlay) ! H2O Ice geometric mean radius (m)
72	c Sedimentation radius of water ice
73	real,intent(in) :: rsedcloud(ngrid,nlay)
74	c Cloud density (kg.m-3)
75	real,intent(inout) :: rhocloud(ngrid,nlay)
76	c Traceurs :
77	real,intent(in) :: pq(ngrid,nlay,nq) ! tracers (kg/kg)
78	real,intent(in) :: pdqfi(ngrid,nlay,nq) ! tendency before sedimentation (kg/kg.s-1)
79	real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
80	real,intent(out) :: pdqs_sed(ngrid,nq) ! flux at surface (kg.m-2.s-1)
81	integer,intent(in) :: nq ! number of tracers
82	real,intent(in) :: tau(ngrid,naerkind) ! dust opacity
83	real,intent(in) :: tauscaling(ngrid)
84
85	c local:
86	c ------
87
88	INTEGER l,ig, iq
89
90	real zqi(ngrid,nlay,nq) ! to locally store tracers
91	real zt(ngrid,nlay) ! to locally store temperature
92	real masse (ngrid,nlay) ! Layer mass (kg.m-2)
93	real epaisseur (ngrid,nlay) ! Layer thickness (m)
94	real wq(ngrid,nlay+1),w(ngrid,nlay) ! displaced tracer mass wq (kg.m-2), MVals: displaced "pere" tracer mass w (kg.m-2)
95	real r0(ngrid,nlay) ! geometric mean radius used for
96	! sedimentation (m)
97	real r0dust(ngrid,nlay) ! geometric mean radius used for
98	! dust (m)
99	real r0stormdust(ngrid,nlay) ! Geometric mean radius used for stormdust (m)
100	! ! CCNs (m)
101	real r0topdust(ngrid,nlay) ! Geometric mean radius used for topdust (m)
102	! ! CCNs (m)
103	real,save :: beta ! correction for the shape of the ice particles (cf. newsedim)
104	c for ice radius computation
105	REAL Mo,No
106	REAl ccntyp
107	character(len=20),parameter :: modname="callsedim"
108	c MVals: transport of the isotopic ratio
109	REAL Ratio0(ngrid,nlay),Ratio(ngrid,nlay)
110	REAL masseq(ngrid,nlay)
111	REAL newmasse
112	REAL zq0(ngrid,nlay,nq)
113	INTEGER ifils,iq2
114
115
116	c Discrete size distributions (doubleq)
117	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
118	c 1) Parameters used to represent the changes in fall
119	c velocity as a function of particle size;
120	integer ir
121	integer,parameter :: nr=12 !(nr=7) ! number of bins
122	real,save :: rd(nr)
123	real qr(ngrid,nlay,nr)
124	real,save :: rdi(nr+1) ! extreme and intermediate radii
125	real Sq(ngrid,nlay)
126	real,parameter :: rdmin=1.e-8
127	real,parameter :: rdmax=30.e-6
128	real,parameter :: rdimin=1.e-8 ! 1.e-7
129	real,parameter :: rdimax=1.e-4
130
131	c 2) Second size distribution for the log-normal integration
132	c (the mass mixing ratio is computed for each radius)
133
134	integer iint
135	integer,parameter :: ninter=4 ! number of points between each rdi radii
136	real,save :: rr(ninter,nr)
137	integer radpower
138	real sigma0
139
140	c 3) Other local variables used in doubleq
141
142	INTEGER,SAVE :: idust_mass ! index of tracer containing dust mass
143	! mix. ratio
144	INTEGER,SAVE :: idust_number ! index of tracer containing dust number
145	! mix. ratio
146	INTEGER,SAVE :: iccn_mass ! index of tracer containing CCN mass
147	! mix. ratio
148	INTEGER,SAVE :: iccn_number ! index of tracer containing CCN number
149	! mix. ratio
150	INTEGER,SAVE :: istormdust_mass ! index of tracer containing
151	!stormdust mass mix. ratio
152	INTEGER,SAVE :: istormdust_number ! index of tracer containing
153	!stormdust number mix. ratio
154	INTEGER,SAVE :: itopdust_mass ! index of tracer containing
155	!topdust mass mix. ratio
156	INTEGER,SAVE :: itopdust_number ! index of tracer containing
157	!topdust number mix. ratio
158	INTEGER,SAVE :: iccnco2_number ! index of tracer containing CCN number
159	INTEGER,SAVE :: iccnco2_mass ! index of tracer containing CCN number
160	INTEGER,SAVE :: ico2_ice ! index of tracer containing CCN number
161
162
163	LOGICAL,SAVE :: firstcall=.true.
164
165
166
167	c ** un petit test de coherence
168	c --------------------------
169	! AS: firstcall OK absolute
170	IF (firstcall) THEN
171
172	c Doubleq: initialization
173	IF (doubleq) THEN
174	do ir=1,nr
175	rd(ir)= rdmin(rdmax/rdmin)*(float(ir-1)/float(nr-1))
176	end do
177	rdi(1)=rdimin
178	do ir=2,nr
179	rdi(ir)= sqrt(rd(ir-1)*rd(ir))
180	end do
181	rdi(nr+1)=rdimax
182
183	do ir=1,nr
184	do iint=1,ninter
185	rr(iint,ir)=
186	& rdi(ir)*
187	& (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
188	c write(,) rr(iint,ir)
189	end do
190	end do
191
192	! identify tracers corresponding to mass mixing ratio and
193	! number mixing ratio
194	idust_mass=0 ! dummy initialization
195	idust_number=0 ! dummy initialization
196
197	do iq=1,nq
198	if (noms(iq).eq."dust_mass") then
199	idust_mass=iq
200	write(,)"callsedim: idust_mass=",idust_mass
201	endif
202	if (noms(iq).eq."dust_number") then
203	idust_number=iq
204	write(,)"callsedim: idust_number=",idust_number
205	endif
206	enddo
207
208	! check that we did find the tracers
209	if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
210	write(,) 'callsedim: error! could not identify'
211	write(,) ' tracers for dust mass and number mixing'
212	write(,) ' ratio and doubleq is activated!'
213	call abort_physic(modname,"missing dust tracers",1)
214	endif
215	ENDIF !of if (doubleq)
216
217	IF (microphys) THEN
218	iccn_mass=0
219	iccn_number=0
220	do iq=1,nq
221	if (noms(iq).eq."ccn_mass") then
222	iccn_mass=iq
223	write(,)"callsedim: iccn_mass=",iccn_mass
224	endif
225	if (noms(iq).eq."ccn_number") then
226	iccn_number=iq
227	write(,)"callsedim: iccn_number=",iccn_number
228	endif
229	enddo
230	! check that we did find the tracers
231	if ((iccn_mass.eq.0).or.(iccn_number.eq.0)) then
232	write(,) 'callsedim: error! could not identify'
233	write(,) ' tracers for ccn mass and number mixing'
234	write(,) ' ratio and microphys is activated!'
235	call abort_physic(modname,"missing ccn tracers",1)
236	endif
237	ENDIF !of if (microphys)
238
239	IF (co2clouds) THEN
240	iccnco2_mass=0
241	iccnco2_number=0
242	ico2_ice=0
243	do iq=1,nq
244	if (noms(iq).eq."ccnco2_mass") then
245	iccnco2_mass=iq
246	write(,)"callsedim: iccnco2_mass=",iccnco2_mass
247	endif
248	if (noms(iq).eq."co2_ice") then
249	ico2_ice=iq
250	write(,)"callsedim: ico2_ice=",ico2_ice
251	endif
252	if (noms(iq).eq."ccnco2_number") then
253	iccnco2_number=iq
254	write(,)"callsedim: iccnco2_number=",iccnco2_number
255	endif
256	enddo
257	! check that we did find the tracers
258	if ((iccnco2_mass.eq.0).or.(iccnco2_number.eq.0)) then
259	write(,) 'callsedim: error! could not identify'
260	write(,) ' tracers for ccn co2 mass and number mixing'
261	write(,) ' ratio and co2clouds are activated!'
262	call abort_physic(modname,"missing co2 ccn tracers",1)
263	endif
264	ENDIF !of if (co2clouds)
265
266	IF (water) THEN
267	write(,) "correction for the shape of the ice particles ?"
268	beta=0.75 ! default value
269	call getin_p("ice_shape",beta)
270	write(,) " ice_shape = ",beta
271
272	write(,) "water_param nueff Sedimentation:", nuice_sed
273	IF (activice) THEN
274	write(,) "water_param nueff Radiative:", nuice_ref
275	ENDIF
276	ENDIF
277
278	IF (rdstorm) THEN ! identifying stormdust tracers for sedimentation
279	istormdust_mass=0 ! dummy initialization
280	istormdust_number=0 ! dummy initialization
281
282	do iq=1,nq
283	if (noms(iq).eq."stormdust_mass") then
284	istormdust_mass=iq
285	write(,)"callsedim: istormdust_mass=",istormdust_mass
286	endif
287	if (noms(iq).eq."stormdust_number") then
288	istormdust_number=iq
289	write(,)"callsedim: istormdust_number=",
290	& istormdust_number
291	endif
292	enddo
293
294	! check that we did find the tracers
295	if ((istormdust_mass.eq.0).or.(istormdust_number.eq.0)) then
296	write(,) 'callsedim: error! could not identify'
297	write(,) ' tracers for stormdust mass and number mixing'
298	write(,) ' ratio and rdstorm is activated!'
299	call abort_physic(modname,"missing stormdust tracers",1)
300	endif
301	ENDIF !of if (rdstorm)
302
303	IF (slpwind) THEN ! identifying topdust tracers for sedimentation
304	itopdust_mass=0 ! dummy initialization
305	itopdust_number=0 ! dummy initialization
306
307	do iq=1,nq
308	if (noms(iq).eq."topdust_mass") then
309	itopdust_mass=iq
310	write(,)"callsedim: itopdust_mass=",itopdust_mass
311	endif
312	if (noms(iq).eq."topdust_number") then
313	itopdust_number=iq
314	write(,)"callsedim: itopdust_number=",
315	& itopdust_number
316	endif
317	enddo
318
319	! check that we did find the tracers
320	if ((itopdust_mass.eq.0).or.(itopdust_number.eq.0)) then
321	write(,) 'callsedim: error! could not identify'
322	write(,) ' tracers for topdust mass and number mixing'
323	write(,) ' ratio and slpwind is activated!'
324	call abort_physic(modname,"missing topdust tracers",1)
325	endif
326	ENDIF !of if (slpwind)
327
328	firstcall=.false.
329	ENDIF ! of IF (firstcall)
330
331	c-----------------------------------------------------------------------
332	c 1. Initialization
333	c -----------------
334
335	! zqi(1:ngrid,1:nlay,1:nqmx) = 0.
336	c Update the mass mixing ratio and temperature with the tendencies coming
337	c from other parameterizations:
338	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
339	zqi(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)
340	& +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
341	zq0(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq) !MVals: keep the input value
342	& +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
343	zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay)
344	& +pdt(1:ngrid,1:nlay)*ptimestep
345
346	c Computing the different layer properties
347	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
348	c Mass (kg.m-2), thickness(m), crossing time (s) etc.
349
350	do l=1,nlay
351	do ig=1, ngrid
352	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
353	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
354	end do
355	end do
356
357	c =================================================================
358	c Compute the geometric mean radius used for sedimentation
359
360	if (doubleq) then
361	do l=1,nlay
362	do ig=1, ngrid
363
364	call updaterdust(zqi(ig,l,igcm_dust_mass),
365	& zqi(ig,l,igcm_dust_number),r0dust(ig,l),
366	& tauscaling(ig))
367
368	end do
369	end do
370	endif
371	c rocket dust storm
372	if (rdstorm) then
373	do l=1,nlay
374	do ig=1, ngrid
375
376	call updaterdust(zqi(ig,l,igcm_stormdust_mass),
377	& zqi(ig,l,igcm_stormdust_number),r0stormdust(ig,l),
378	& tauscaling(ig))
379
380	end do
381	end do
382	endif
383	c entrainment by slope wind
384	if (slpwind) then
385	do l=1,nlay
386	do ig=1, ngrid
387
388	call updaterdust(zqi(ig,l,igcm_topdust_mass),
389	& zqi(ig,l,igcm_topdust_number),r0topdust(ig,l),
390	& tauscaling(ig))
391
392	end do
393	end do
394	endif
395	c =================================================================
396	do iq=1,nq
397	if(radius(iq).gt.1.e-9 .and.(iq.ne.ico2_ice) .and.
398	& (iq .ne. iccnco2_mass) .and. (iq .ne.
399	& iccnco2_number) .and. ! no sedim for gaz or CO2 clouds (done in microtimestep)
400	& iq .ne. igcm_hdo_ice) then !MVals: hdo is transported by h2o
401	c -----------------------------------------------------------------
402	c DOUBLEQ CASE
403	c -----------------------------------------------------------------
404	if ( doubleq.and.
405	& ((iq.eq.idust_mass).or.(iq.eq.idust_number).or.
406	& (iq.eq.istormdust_mass).or.(iq.eq.istormdust_number).or.
407	& (iq.eq.itopdust_mass).or.(iq.eq.itopdust_number)) ) then
408
409	c Computing size distribution:
410	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
411
412	if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
413	do l=1,nlay
414	do ig=1, ngrid
415	r0(ig,l)=r0dust(ig,l)
416	end do
417	end do
418	else if ((iq.eq.istormdust_mass).or.
419	& (iq.eq.istormdust_number)) then
420	do l=1,nlay
421	do ig=1, ngrid
422	r0(ig,l)=r0stormdust(ig,l)
423	end do
424	end do
425	else if ((iq.eq.itopdust_mass).or.
426	& (iq.eq.itopdust_number)) then
427	do l=1,nlay
428	do ig=1, ngrid
429	r0(ig,l)=r0topdust(ig,l)
430	end do
431	end do
432	endif
433	sigma0 = varian
434
435	c Computing mass mixing ratio for each particle size
436	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
437	IF ((iq.EQ.idust_mass).or.(iq.EQ.istormdust_mass)
438	& .or.(iq.EQ.itopdust_mass)) then
439	radpower = 2
440	ELSE ! number
441	radpower = -1
442	ENDIF
443	Sq(1:ngrid,1:nlay) = 0.
444	do ir=1,nr
445	do l=1,nlay
446	do ig=1,ngrid
447	c ****************
448	c Size distribution integration
449	c (Trapezoid Integration Method)
450	qr(ig,l,ir)=0.5(rr(2,ir)-rr(1,ir))
451	& (rr(1,ir)*radpower)
452	& exp(-(log(rr(1,ir)/r0(ig,l)))*2/(2sigma0**2))
453	do iint=2,ninter-1
454	qr(ig,l,ir)=qr(ig,l,ir) +
455	& 0.5(rr(iint+1,ir)-rr(iint-1,ir))
456	& (rr(iint,ir)*radpower)
457	& exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
458	& (2sigma0*2))
459	end do
460	qr(ig,l,ir)=qr(ig,l,ir) +
461	& 0.5(rr(ninter,ir)-rr(ninter-1,ir))
462	& (rr(ninter,ir)*radpower)
463	& exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
464	& (2sigma0*2))
465
466	c **************** old method (not recommended!)
467	c qr(ig,l,ir)=(rd(ir)*(5-3iq))*
468	c & exp( -(log(rd(ir)/r0(ig,l)))*2 / (2sigma0**2) )
469	c ******************************
470
471	Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
472	enddo
473	enddo
474	enddo
475
476	do ir=1,nr
477	do l=1,nlay
478	do ig=1,ngrid
479	qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
480	enddo
481	enddo
482	enddo
483
484	c Computing sedimentation for each tracer
485	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
486
487	zqi(1:ngrid,1:nlay,iq) = 0.
488	pdqs_sed(1:ngrid,iq) = 0.
489
490	do ir=1,nr
491	call newsedim(ngrid,nlay,1,1,ptimestep,
492	& pplev,masse,epaisseur,zt,rd(ir),(/rho_dust/),qr(1,1,ir),
493	& wq,0.5)
494
495	c Tendencies
496	c ~~~~~~~~~~
497	do ig=1,ngrid
498	pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
499	& + wq(ig,1)/ptimestep
500	end do
501	DO l = 1, nlay
502	DO ig=1,ngrid
503	zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
504	ENDDO
505	ENDDO
506	enddo ! of do ir=1,nr
507	c -----------------------------------------------------------------
508	c WATER CYCLE CASE
509	c -----------------------------------------------------------------
510	else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
511	& .or. (iq .eq. igcm_h2o_ice)) then
512	if (microphys) then
513	c water ice sedimentation
514	c ~~~~~~~~~~
515	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
516	& ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rhocloud,
517	& zqi(1,1,iq),wq,beta)
518	c Surface Tendencies
519	c ~~~~~~~~~~
520	do ig=1,ngrid
521	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
522	end do
523	else
524	c water ice sedimentation
525	c ~~~~~~~~~~
526	call newsedim(ngrid,nlay,ngrid*nlay,1,
527	& ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rho_q(iq),
528	& zqi(1,1,iq),wq,beta)
529	c Surface tendencies
530	c ~~~~~~~~~~
531	do ig=1,ngrid
532	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
533	end do
534	c Special case of isotopes
535	c ~~~~~~~~~~
536	!MVals: Loop over the sons ("fils")
537	if (nqfils(iq).gt.0) then
538	iq2=iqfils(nqfils(iq),iq) ! for now always nqfils(iq)=1 (special case of HDO only son of H2O)
539	!MVals: input paramters in vlz_fi for hdo
540	do l=1,nlay
541	do ig=1,ngrid
542	if (zq0(ig,l,iq).gt.qperemin) then
543	Ratio0(ig,l)=zq0(ig,l,iq2)/zq0(ig,l,iq)
544	else
545	Ratio0(ig,l)=0.
546	endif
547	Ratio(ig,l)=Ratio0(ig,l)
548	masseq(ig,l)=max(masse(ig,l)*zq0(ig,l,iq),masseqmin)
549	w(ig,l)=wq(ig,l) !MVals: very important: hdo is transported by h2o (see vlsplt_p.F: correction bugg 15 mai 2015)
550	enddo !ig=1,ngrid
551	enddo !l=1,nlay
552	!MVals: no need to enter newsedim as the transporting mass w has been already calculated
553	call vlz_fi(ngrid,nlay,Ratio,2.,masseq,w,wq)
554	zqi(:,nlay,iq2)=zqi(:,nlay,iq)*Ratio0(:,nlay)
555	do l=1,nlay-1
556	do ig=1,ngrid
557	newmasse=max(masseq(ig,l)+w(ig,l+1)-w(ig,l),masseqmin)
558	Ratio(ig,l)=(Ratio0(ig,l)*masseq(ig,l)
559	& +wq(ig,l+1)-wq(ig,l))/newmasse
560	zqi(ig,l,iq2)=zqi(ig,l,iq)*Ratio(ig,l)
561	enddo
562	enddo !l=1,nlay-1
563	!MVals: hdo surface tendency
564	do ig=1,ngrid
565	if (w(ig,1).gt.masseqmin) then
566	pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*(wq(ig,1)/w(ig,1))
567	else
568	pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*Ratio0(ig,1)
569	endif
570	end do
571	endif !(nqfils(iq).gt.0)
572	endif ! of if (microphys)
573
574	c -----------------------------------------------------------------
575	c GENERAL CASE
576	c -----------------------------------------------------------------
577	else
578	call newsedim(ngrid,nlay,1,1,ptimestep,
579	& pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
580	& zqi(1,1,iq),wq,1.0)
581	c Tendencies
582	c ~~~~~~~~~~
583	do ig=1,ngrid
584	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
585	end do
586	endif ! of if doubleq and if water
587	c -----------------------------------------------------------------
588
589	c Compute the final tendency:
590	c ---------------------------
591	DO l = 1, nlay
592	DO ig=1,ngrid
593	pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
594	$ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
595	!MVals: Special case of isotopes: for now only HDO
596	if (nqfils(iq).gt.0) then
597	iq2=iqfils(nqfils(iq),iq)
598	pdqsed(ig,l,iq2)=(zqi(ig,l,iq2)-
599	$ (pq(ig,l,iq2) + pdqfi(ig,l,iq2)*ptimestep))/ptimestep
600	endif
601	ENDDO
602	ENDDO
603
604	endif ! of if(radius(iq).gt.1.e-9)
605	c =================================================================
606	enddo ! of do iq=1,nq
607
608	c Update the dust particle size "rdust"
609	c -------------------------------------
610	if (doubleq) then
611	DO l = 1, nlay
612	DO ig=1,ngrid
613
614
615	call updaterdust(zqi(ig,l,igcm_dust_mass),
616	& zqi(ig,l,igcm_dust_number),rdust(ig,l),
617	& tauscaling(ig))
618
619
620	ENDDO
621	ENDDO
622	endif ! of if (doubleq)
623
624	if (rdstorm) then
625	DO l = 1, nlay
626	DO ig=1,ngrid
627	call updaterdust(zqi(ig,l,igcm_stormdust_mass),
628	& zqi(ig,l,igcm_stormdust_number),rstormdust(ig,l),
629	& tauscaling(ig))
630	ENDDO
631	ENDDO
632	endif ! of if (rdstorm)
633
634	if (slpwind) then
635	DO l = 1, nlay
636	DO ig=1,ngrid
637	call updaterdust(zqi(ig,l,igcm_topdust_mass),
638	& zqi(ig,l,igcm_topdust_number),rtopdust(ig,l),
639	& tauscaling(ig))
640	ENDDO
641	ENDDO
642	endif ! of if (slpwind)
643
644	c Update the ice particle size "rice"
645	c -------------------------------------
646	if (water) then
647	IF(microphys) THEN
648
649
650	DO l = 1, nlay
651	DO ig=1,ngrid
652
653	call updaterice_micro(zqi(ig,l,igcm_h2o_ice),
654	& zqi(ig,l,igcm_ccn_mass),zqi(ig,l,igcm_ccn_number),
655	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
656
657	ENDDO
658	ENDDO
659
660	ELSE
661
662	DO l = 1, nlay
663	DO ig=1,ngrid
664
665	call updaterice_typ(zqi(ig,l,igcm_h2o_ice),
666	& tau(ig,1),zlay(ig,l),rice(ig,l))
667
668	ENDDO
669	ENDDO
670	ENDIF ! of IF(microphys)
671	endif ! of if (water)
672
673	END SUBROUTINE callsedim
674
675	END MODULE callsedim_mod
676

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

Download in other formats:

Original Format