Context Navigation

callsedim_mod.F @ 2562

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

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

Download in other formats:

Original Format