Context Navigation

callsedim_mod.F @ 2596

Last change on this file since 2596 was 2589, checked in by cmathe, 4 years ago
GCM Mars: add meteoritic particles as CN for CO2 microphysics
File size: 27.8 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_co2_ice, igcm_stormdust_mass,
22	& igcm_stormdust_number,igcm_topdust_mass,
23	& igcm_topdust_number,
24	& nqfils,qperemin,masseqmin ! MVals: variables isotopes
25	USE newsedim_mod, ONLY: newsedim
26	USE comcstfi_h, ONLY: g
27	USE dimradmars_mod, only: naerkind
28	USE dust_param_mod, ONLY: doubleq
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 :: iccnco2_h2o_number ! index of tracer containing CCN number
161	INTEGER,SAVE :: iccnco2_h2o_mass_ice ! index of tracer containing CCN number
162	INTEGER,SAVE :: iccnco2_h2o_mass_ccn ! index of tracer containing CCN number
163	INTEGER,SAVE :: iccnco2_meteor_number ! index of tracer containing CCN number
164	INTEGER,SAVE :: iccnco2_meteor_mass ! index of tracer containing CCN number
165	INTEGER,SAVE :: ico2_ice ! index of tracer containing CCN number
166
167
168	LOGICAL,SAVE :: firstcall=.true.
169
170
171
172	c ** un petit test de coherence
173	c --------------------------
174	! AS: firstcall OK absolute
175	IF (firstcall) THEN
176
177	c Doubleq: initialization
178	IF (doubleq) THEN
179	do ir=1,nr
180	rd(ir)= rdmin(rdmax/rdmin)*(float(ir-1)/float(nr-1))
181	end do
182	rdi(1)=rdimin
183	do ir=2,nr
184	rdi(ir)= sqrt(rd(ir-1)*rd(ir))
185	end do
186	rdi(nr+1)=rdimax
187
188	do ir=1,nr
189	do iint=1,ninter
190	rr(iint,ir)=
191	& rdi(ir)*
192	& (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
193	c write(,) rr(iint,ir)
194	end do
195	end do
196
197	! identify tracers corresponding to mass mixing ratio and
198	! number mixing ratio
199	idust_mass=0 ! dummy initialization
200	idust_number=0 ! dummy initialization
201
202	do iq=1,nq
203	if (noms(iq).eq."dust_mass") then
204	idust_mass=iq
205	write(,)"callsedim: idust_mass=",idust_mass
206	endif
207	if (noms(iq).eq."dust_number") then
208	idust_number=iq
209	write(,)"callsedim: idust_number=",idust_number
210	endif
211	enddo
212
213	! check that we did find the tracers
214	if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
215	write(,) 'callsedim: error! could not identify'
216	write(,) ' tracers for dust mass and number mixing'
217	write(,) ' ratio and doubleq is activated!'
218	call abort_physic(modname,"missing dust tracers",1)
219	endif
220	ENDIF !of if (doubleq)
221
222	IF (microphys) THEN
223	iccn_mass=0
224	iccn_number=0
225	do iq=1,nq
226	if (noms(iq).eq."ccn_mass") then
227	iccn_mass=iq
228	write(,)"callsedim: iccn_mass=",iccn_mass
229	endif
230	if (noms(iq).eq."ccn_number") then
231	iccn_number=iq
232	write(,)"callsedim: iccn_number=",iccn_number
233	endif
234	enddo
235	! check that we did find the tracers
236	if ((iccn_mass.eq.0).or.(iccn_number.eq.0)) then
237	write(,) 'callsedim: error! could not identify'
238	write(,) ' tracers for ccn mass and number mixing'
239	write(,) ' ratio and microphys is activated!'
240	call abort_physic(modname,"missing ccn tracers",1)
241	endif
242	ENDIF !of if (microphys)
243
244	IF (co2clouds) THEN
245	iccnco2_mass=0
246	iccnco2_number=0
247	iccnco2_h2o_mass_ice=0
248	iccnco2_h2o_mass_ccn=0
249	iccnco2_h2o_number=0
250	iccnco2_meteor_mass=0
251	iccnco2_meteor_number=0
252	ico2_ice=0
253	do iq=1,nq
254	if (noms(iq).eq."ccnco2_mass") then
255	iccnco2_mass=iq
256	write(,)"callsedim: iccnco2_mass=",iccnco2_mass
257	endif
258	if (noms(iq).eq."co2_ice") then
259	ico2_ice=iq
260	write(,)"callsedim: ico2_ice=",ico2_ice
261	endif
262	if (noms(iq).eq."ccnco2_number") then
263	iccnco2_number=iq
264	write(,)"callsedim: iccnco2_number=",iccnco2_number
265	endif
266	if (noms(iq).eq."ccnco2_h2o_number") then
267	iccnco2_h2o_number=iq
268	write(,)"callsedim: iccnco2_h2o_number=",
269	& iccnco2_h2o_number
270	endif
271	if (noms(iq).eq."ccnco2_h2o_mass_ice") then
272	iccnco2_h2o_mass_ice=iq
273	write(,)"callsedim: iccnco2_h2o_mass_ice=",
274	& iccnco2_h2o_mass_ice
275	endif
276	if (noms(iq).eq."ccnco2_h2o_mass_ccn") then
277	iccnco2_h2o_mass_ccn=iq
278	write(,)"callsedim: iccnco2_h2o_mass_ccn=",
279	& iccnco2_h2o_mass_ccn
280	endif
281	if (noms(iq).eq."ccnco2_meteor_number") then
282	iccnco2_meteor_number=iq
283	write(,)"callsedim: iccnco2_meteor_number=",
284	& iccnco2_meteor_number
285	endif
286	if (noms(iq).eq."ccnco2_meteor_mass") then
287	iccnco2_meteor_mass=iq
288	write(,)"callsedim: iccnco2_meteor_mass=",
289	& iccnco2_meteor_mass
290	endif
291	enddo
292	! check that we did find the tracers
293	if ((iccnco2_mass.eq.0).or.(iccnco2_number.eq.0)) then
294	write(,) 'callsedim: error! could not identify'
295	write(,) ' tracers for ccn co2 mass and number mixing'
296	write(,) ' ratio and co2clouds are activated!'
297	call abort_physic(modname,"missing co2 ccn tracers",1)
298	endif
299	if (co2useh2o) then
300	if((iccnco2_h2o_mass_ccn.eq.0).or.
301	& (iccnco2_h2o_mass_ice.eq.0).or.
302	& (iccnco2_h2o_number.eq.0)) then
303	write(,) 'callsedim: error! could not identify'
304	write(,) ' tracers for ccn co2 mass and number mixing'
305	write(,) ' ratio and co2clouds are activated!'
306	call abort_physic(modname,"missing co2 ccn tracers",1)
307	end if
308	end if
309	if (meteo_flux) then
310	if((iccnco2_meteor_mass.eq.0).or.
311	& (iccnco2_meteor_number.eq.0)) then
312	write(,) 'callsedim: error! could not identify'
313	write(,) ' tracers for ccn co2 mass and number mixing'
314	write(,) ' ratio and co2clouds are activated!'
315	call abort_physic(modname,"missing co2 ccn tracers",1)
316	end if
317	end if
318
319	ENDIF !of if (co2clouds)
320
321	IF (water) THEN
322	write(,) "correction for the shape of the ice particles ?"
323	beta=0.75 ! default value
324	call getin_p("ice_shape",beta)
325	write(,) " ice_shape = ",beta
326
327	write(,) "water_param nueff Sedimentation:", nuice_sed
328	IF (activice) THEN
329	write(,) "water_param nueff Radiative:", nuice_ref
330	ENDIF
331	ENDIF
332
333	IF (rdstorm) THEN ! identifying stormdust tracers for sedimentation
334	istormdust_mass=0 ! dummy initialization
335	istormdust_number=0 ! dummy initialization
336
337	do iq=1,nq
338	if (noms(iq).eq."stormdust_mass") then
339	istormdust_mass=iq
340	write(,)"callsedim: istormdust_mass=",istormdust_mass
341	endif
342	if (noms(iq).eq."stormdust_number") then
343	istormdust_number=iq
344	write(,)"callsedim: istormdust_number=",
345	& istormdust_number
346	endif
347	enddo
348
349	! check that we did find the tracers
350	if ((istormdust_mass.eq.0).or.(istormdust_number.eq.0)) then
351	write(,) 'callsedim: error! could not identify'
352	write(,) ' tracers for stormdust mass and number mixing'
353	write(,) ' ratio and rdstorm is activated!'
354	call abort_physic(modname,"missing stormdust tracers",1)
355	endif
356	ENDIF !of if (rdstorm)
357
358	IF (slpwind) THEN ! identifying topdust tracers for sedimentation
359	itopdust_mass=0 ! dummy initialization
360	itopdust_number=0 ! dummy initialization
361
362	do iq=1,nq
363	if (noms(iq).eq."topdust_mass") then
364	itopdust_mass=iq
365	write(,)"callsedim: itopdust_mass=",itopdust_mass
366	endif
367	if (noms(iq).eq."topdust_number") then
368	itopdust_number=iq
369	write(,)"callsedim: itopdust_number=",
370	& itopdust_number
371	endif
372	enddo
373
374	! check that we did find the tracers
375	if ((itopdust_mass.eq.0).or.(itopdust_number.eq.0)) then
376	write(,) 'callsedim: error! could not identify'
377	write(,) ' tracers for topdust mass and number mixing'
378	write(,) ' ratio and slpwind is activated!'
379	call abort_physic(modname,"missing topdust tracers",1)
380	endif
381	ENDIF !of if (slpwind)
382
383	firstcall=.false.
384	ENDIF ! of IF (firstcall)
385
386	c-----------------------------------------------------------------------
387	c 1. Initialization
388	c -----------------
389
390	! zqi(1:ngrid,1:nlay,1:nqmx) = 0.
391	c Update the mass mixing ratio and temperature with the tendencies coming
392	c from other parameterizations:
393	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
394	zqi(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)
395	& +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
396	zq0(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq) !MVals: keep the input value
397	& +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
398	zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay)
399	& +pdt(1:ngrid,1:nlay)*ptimestep
400
401	c Computing the different layer properties
402	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
403	c Mass (kg.m-2), thickness(m), crossing time (s) etc.
404
405	do l=1,nlay
406	do ig=1, ngrid
407	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
408	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
409	end do
410	end do
411
412	c =================================================================
413	c Compute the geometric mean radius used for sedimentation
414
415	if (doubleq) then
416	do l=1,nlay
417	do ig=1, ngrid
418
419	call updaterdust(zqi(ig,l,igcm_dust_mass),
420	& zqi(ig,l,igcm_dust_number),r0dust(ig,l),
421	& tauscaling(ig))
422
423	end do
424	end do
425	endif
426	c rocket dust storm
427	if (rdstorm) then
428	do l=1,nlay
429	do ig=1, ngrid
430
431	call updaterdust(zqi(ig,l,igcm_stormdust_mass),
432	& zqi(ig,l,igcm_stormdust_number),r0stormdust(ig,l),
433	& tauscaling(ig))
434
435	end do
436	end do
437	endif
438	c entrainment by slope wind
439	if (slpwind) then
440	do l=1,nlay
441	do ig=1, ngrid
442
443	call updaterdust(zqi(ig,l,igcm_topdust_mass),
444	& zqi(ig,l,igcm_topdust_number),r0topdust(ig,l),
445	& tauscaling(ig))
446
447	end do
448	end do
449	endif
450	c =================================================================
451	do iq=1,nq
452	if(radius(iq).gt.1.e-9 .and.(iq.ne.ico2_ice) .and.
453	& (iq .ne. iccnco2_mass) .and. (iq .ne.
454	& iccnco2_number) .and.
455	& (iq.ne.iccnco2_h2o_number).and.
456	& (iq.ne.iccnco2_h2o_mass_ice).and.
457	& (iq.ne.iccnco2_h2o_mass_ccn).and.
458	& (iq.ne.iccnco2_meteor_mass).and.
459	& (iq.ne.iccnco2_meteor_number).and. ! no sedim for gaz or CO2 clouds (done in microtimestep)
460	& iq .ne. igcm_hdo_ice) then !MVals: hdo is transported by h2o
461	c -----------------------------------------------------------------
462	c DOUBLEQ CASE
463	c -----------------------------------------------------------------
464	if ( doubleq.and.
465	& ((iq.eq.idust_mass).or.(iq.eq.idust_number).or.
466	& (iq.eq.istormdust_mass).or.(iq.eq.istormdust_number).or.
467	& (iq.eq.itopdust_mass).or.(iq.eq.itopdust_number)) ) then
468
469	c Computing size distribution:
470	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
471
472	if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
473	do l=1,nlay
474	do ig=1, ngrid
475	r0(ig,l)=r0dust(ig,l)
476	end do
477	end do
478	else if ((iq.eq.istormdust_mass).or.
479	& (iq.eq.istormdust_number)) then
480	do l=1,nlay
481	do ig=1, ngrid
482	r0(ig,l)=r0stormdust(ig,l)
483	end do
484	end do
485	else if ((iq.eq.itopdust_mass).or.
486	& (iq.eq.itopdust_number)) then
487	do l=1,nlay
488	do ig=1, ngrid
489	r0(ig,l)=r0topdust(ig,l)
490	end do
491	end do
492	endif
493	sigma0 = varian
494
495	c Computing mass mixing ratio for each particle size
496	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
497	IF ((iq.EQ.idust_mass).or.(iq.EQ.istormdust_mass)
498	& .or.(iq.EQ.itopdust_mass)) then
499	radpower = 2
500	ELSE ! number
501	radpower = -1
502	ENDIF
503	Sq(1:ngrid,1:nlay) = 0.
504	do ir=1,nr
505	do l=1,nlay
506	do ig=1,ngrid
507	c ****************
508	c Size distribution integration
509	c (Trapezoid Integration Method)
510	qr(ig,l,ir)=0.5(rr(2,ir)-rr(1,ir))
511	& (rr(1,ir)*radpower)
512	& exp(-(log(rr(1,ir)/r0(ig,l)))*2/(2sigma0**2))
513	do iint=2,ninter-1
514	qr(ig,l,ir)=qr(ig,l,ir) +
515	& 0.5(rr(iint+1,ir)-rr(iint-1,ir))
516	& (rr(iint,ir)*radpower)
517	& exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
518	& (2sigma0*2))
519	end do
520	qr(ig,l,ir)=qr(ig,l,ir) +
521	& 0.5(rr(ninter,ir)-rr(ninter-1,ir))
522	& (rr(ninter,ir)*radpower)
523	& exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
524	& (2sigma0*2))
525
526	c **************** old method (not recommended!)
527	c qr(ig,l,ir)=(rd(ir)*(5-3iq))*
528	c & exp( -(log(rd(ir)/r0(ig,l)))*2 / (2sigma0**2) )
529	c ******************************
530
531	Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
532	enddo
533	enddo
534	enddo
535
536	do ir=1,nr
537	do l=1,nlay
538	do ig=1,ngrid
539	qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
540	enddo
541	enddo
542	enddo
543
544	c Computing sedimentation for each tracer
545	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
546
547	zqi(1:ngrid,1:nlay,iq) = 0.
548	pdqs_sed(1:ngrid,iq) = 0.
549
550	do ir=1,nr
551	call newsedim(ngrid,nlay,1,1,ptimestep,
552	& pplev,masse,epaisseur,zt,rd(ir),(/rho_dust/),qr(1,1,ir),
553	& wq,0.5)
554
555	c Tendencies
556	c ~~~~~~~~~~
557	do ig=1,ngrid
558	pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
559	& + wq(ig,1)/ptimestep
560	end do
561	DO l = 1, nlay
562	DO ig=1,ngrid
563	zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
564	ENDDO
565	ENDDO
566	enddo ! of do ir=1,nr
567	c -----------------------------------------------------------------
568	c WATER CYCLE CASE
569	c -----------------------------------------------------------------
570	else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
571	& .or. (iq .eq. igcm_h2o_ice)) then
572	if (microphys) then
573	c water ice sedimentation
574	c ~~~~~~~~~~
575	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
576	& ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rhocloud,
577	& zqi(1,1,iq),wq,beta)
578	else
579	c water ice sedimentation
580	c ~~~~~~~~~~
581	call newsedim(ngrid,nlay,ngrid*nlay,1,
582	& ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rho_q(iq),
583	& zqi(1,1,iq),wq,beta)
584	endif ! of if (microphys)
585	c Surface tendencies
586	c ~~~~~~~~~~
587	do ig=1,ngrid
588	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
589	end do
590	c Special case of isotopes
591	c ~~~~~~~~~~
592	!MVals: for now only water can have an isotope, a son ("fils"), and it has to be hdo
593	if (nqfils(iq).gt.0) then
594	if (iq.eq.igcm_h2o_ice) then
595	iq2=igcm_hdo_ice
596	else
597	call abort_physic("callsedim_mod","invalid isotope",1)
598	endif
599	!MVals: input parameters in vlz_fi for hdo
600	do l=1,nlay
601	do ig=1,ngrid
602	if (zq0(ig,l,iq).gt.qperemin) then
603	Ratio0(ig,l)=zq0(ig,l,iq2)/zq0(ig,l,iq)
604	else
605	Ratio0(ig,l)=0.
606	endif
607	Ratio(ig,l)=Ratio0(ig,l)
608	masseq(ig,l)=max(masse(ig,l)*zq0(ig,l,iq),masseqmin)
609	w(ig,l)=wq(ig,l) !MVals: very important: hdo is transported by h2o (see vlsplt_p.F: correction bugg 15 mai 2015)
610	enddo !ig=1,ngrid
611	enddo !l=1,nlay
612	!MVals: no need to enter newsedim as the transporting mass w has been already calculated
613	call vlz_fi(ngrid,nlay,Ratio,2.,masseq,w,wq)
614	zqi(:,nlay,iq2)=zqi(:,nlay,iq)*Ratio0(:,nlay)
615	do l=1,nlay-1
616	do ig=1,ngrid
617	newmasse=max(masseq(ig,l)+w(ig,l+1)-w(ig,l),masseqmin)
618	Ratio(ig,l)=(Ratio0(ig,l)*masseq(ig,l)
619	& +wq(ig,l+1)-wq(ig,l))/newmasse
620	zqi(ig,l,iq2)=zqi(ig,l,iq)*Ratio(ig,l)
621	enddo
622	enddo !l=1,nlay-1
623	!MVals: hdo surface tendency
624	do ig=1,ngrid
625	if (w(ig,1).gt.masseqmin) then
626	pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*(wq(ig,1)/w(ig,1))
627	else
628	pdqs_sed(ig,iq2)=pdqs_sed(ig,iq)*Ratio0(ig,1)
629	endif
630	end do
631	endif !(nqfils(iq).gt.0)
632	c -----------------------------------------------------------------
633	c GENERAL CASE
634	c -----------------------------------------------------------------
635	else
636	call newsedim(ngrid,nlay,1,1,ptimestep,
637	& pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
638	& zqi(1,1,iq),wq,1.0)
639	c Tendencies
640	c ~~~~~~~~~~
641	do ig=1,ngrid
642	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
643	end do
644	endif ! of if doubleq and if water
645	c -----------------------------------------------------------------
646
647	c Compute the final tendency:
648	c ---------------------------
649	DO l = 1, nlay
650	DO ig=1,ngrid
651	pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
652	$ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
653	!MVals: Special case of isotopes: for now only HDO
654	if (nqfils(iq).gt.0) then
655	if (iq.eq.igcm_h2o_ice) then
656	iq2=igcm_hdo_ice
657	else
658	call abort_physic("callsedim_mod","invalid isotope",1)
659	endif
660	pdqsed(ig,l,iq2)=(zqi(ig,l,iq2)-
661	$ (pq(ig,l,iq2) + pdqfi(ig,l,iq2)*ptimestep))/ptimestep
662	endif
663	ENDDO
664	ENDDO
665
666	endif ! of if(radius(iq).gt.1.e-9)
667	c =================================================================
668	enddo ! of do iq=1,nq
669
670	c Update the dust particle size "rdust"
671	c -------------------------------------
672	if (doubleq) then
673	DO l = 1, nlay
674	DO ig=1,ngrid
675
676
677	call updaterdust(zqi(ig,l,igcm_dust_mass),
678	& zqi(ig,l,igcm_dust_number),rdust(ig,l),
679	& tauscaling(ig))
680
681
682	ENDDO
683	ENDDO
684	endif ! of if (doubleq)
685
686	if (rdstorm) then
687	DO l = 1, nlay
688	DO ig=1,ngrid
689	call updaterdust(zqi(ig,l,igcm_stormdust_mass),
690	& zqi(ig,l,igcm_stormdust_number),rstormdust(ig,l),
691	& tauscaling(ig))
692	ENDDO
693	ENDDO
694	endif ! of if (rdstorm)
695
696	if (slpwind) then
697	DO l = 1, nlay
698	DO ig=1,ngrid
699	call updaterdust(zqi(ig,l,igcm_topdust_mass),
700	& zqi(ig,l,igcm_topdust_number),rtopdust(ig,l),
701	& tauscaling(ig))
702	ENDDO
703	ENDDO
704	endif ! of if (slpwind)
705
706	c Update the ice particle size "rice"
707	c -------------------------------------
708	if (water) then
709	IF(microphys) THEN
710
711
712	DO l = 1, nlay
713	DO ig=1,ngrid
714
715	call updaterice_micro(zqi(ig,l,igcm_h2o_ice),
716	& zqi(ig,l,igcm_ccn_mass),zqi(ig,l,igcm_ccn_number),
717	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
718
719	ENDDO
720	ENDDO
721
722	ELSE
723
724	DO l = 1, nlay
725	DO ig=1,ngrid
726
727	call updaterice_typ(zqi(ig,l,igcm_h2o_ice),
728	& tau(ig,1),zlay(ig,l),rice(ig,l))
729
730	ENDDO
731	ENDDO
732	ENDIF ! of IF(microphys)
733	endif ! of if (water)
734	END SUBROUTINE callsedim
735
736	END MODULE callsedim_mod
737

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

Download in other formats:

Original Format