Context Navigation

callsedim_mod.F @ 2541

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

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

Download in other formats:

Original Format