Context Navigation

callsedim_mod.F @ 2616

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

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

Download in other formats:

Original Format