Context Navigation

callsedim_mod.F @ 3804

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

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

Download in other formats:

Original Format