Context Navigation

callsedim_mod.F @ 1974

Last change on this file since 1974 was 1974, checked in by mvals, 6 years ago
Mars GCM: Integration of the detached dust layer parametrizations (rocket dust storm, slope wind lifting, CW, and dust injection scheme, DB). Still experimental, default behaviour (rdstorm=.false., dustinjection=0) identical to previous revision. NB: Updated newstart requires an updated "surface.nc" containing the "hmons" field. EM+MV
File size: 19.3 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,rdust,rstormdust,rice,
9	& rsedcloud,rhocloud,
10	& pq,pdqfi,pdqsed,pdqs_sed,nq,
11	& tau,tauscaling)
12
13	USE ioipsl_getincom, only: getin
14	USE updaterad, only: updaterdust,updaterice_micro,updaterice_typ
15	USE tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
16	& rho_dust, rho_q, radius, varian,
17	& igcm_ccn_mass, igcm_ccn_number,
18	& igcm_h2o_ice, nuice_sed, nuice_ref,
19	& igcm_ccnco2_mass,igcm_ccnco2_number,
20	& igcm_co2_ice, igcm_stormdust_mass,
21	& igcm_stormdust_number
22	USE newsedim_mod, ONLY: newsedim
23	USE comcstfi_h, ONLY: g
24	IMPLICIT NONE
25
26	c=======================================================================
27	c Sedimentation of the Martian aerosols
28	c depending on their density and radius
29	c
30	c F.Forget 1999
31	c
32	c Modified by J.-B. Madeleine 2010: Now includes the doubleq
33	c technique in order to have only one call to callsedim in
34	c physiq.F.
35	c
36	c Modified by J. Audouard 09/16: Now includes the co2clouds case
37	c If the co2 microphysics is on, then co2 theice & ccn tracers
38	c are being sedimented in the microtimestep (co2cloud.F), not
39	c in this routine.
40	c
41	c=======================================================================
42
43	c-----------------------------------------------------------------------
44	c declarations:
45	c -------------
46
47	include "callkeys.h"
48
49	c
50	c arguments:
51	c ----------
52
53	integer,intent(in) :: ngrid ! number of horizontal grid points
54	integer,intent(in) :: nlay ! number of atmospheric layers
55	real,intent(in) :: ptimestep ! physics time step (s)
56	real,intent(in) :: pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
57	real,intent(in) :: zlev(ngrid,nlay+1) ! altitude at layer boundaries
58	real,intent(in) :: zlay(ngrid,nlay) ! altitude at the middle of the layers
59	real,intent(in) :: pt(ngrid,nlay) ! temperature at mid-layer (K)
60	real,intent(in) :: pdt(ngrid,nlay) ! tendency on temperature, from
61	! previous processes (K/s)
62	c Aerosol radius provided by the water ice microphysical scheme:
63	real,intent(out) :: rdust(ngrid,nlay) ! Dust geometric mean radius (m)
64	real,intent(out) :: rstormdust(ngrid,nlay) ! Stormdust geometric mean radius (m)
65	real,intent(out) :: rice(ngrid,nlay) ! H2O Ice geometric mean radius (m)
66	c Sedimentation radius of water ice
67	real,intent(in) :: rsedcloud(ngrid,nlay)
68	c Cloud density (kg.m-3)
69	real,intent(inout) :: rhocloud(ngrid,nlay)
70	c Traceurs :
71	real,intent(in) :: pq(ngrid,nlay,nq) ! tracers (kg/kg)
72	real,intent(in) :: pdqfi(ngrid,nlay,nq) ! tendency before sedimentation (kg/kg.s-1)
73	real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
74	real,intent(out) :: pdqs_sed(ngrid,nq) ! flux at surface (kg.m-2.s-1)
75	integer,intent(in) :: nq ! number of tracers
76	real,intent(in) :: tau(ngrid,nlay) ! dust opacity
77	real,intent(in) :: tauscaling(ngrid)
78
79	c local:
80	c ------
81
82	INTEGER l,ig, iq
83
84	real zqi(ngrid,nlay,nq) ! to locally store tracers
85	real zt(ngrid,nlay) ! to locally store temperature
86	real masse (ngrid,nlay) ! Layer mass (kg.m-2)
87	real epaisseur (ngrid,nlay) ! Layer thickness (m)
88	real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2)
89	real r0(ngrid,nlay) ! geometric mean radius used for
90	! sedimentation (m)
91	real r0dust(ngrid,nlay) ! geometric mean radius used for
92	! dust (m)
93	real r0stormdust(ngrid,nlay) ! Geometric mean radius used for stormdust (m)
94	! ! CCNs (m)
95	real,save :: beta ! correction for the shape of the ice particles (cf. newsedim)
96	c for ice radius computation
97	REAL Mo,No
98	REAl ccntyp
99
100
101
102	c Discrete size distributions (doubleq)
103	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
104	c 1) Parameters used to represent the changes in fall
105	c velocity as a function of particle size;
106	integer ir
107	integer,parameter :: nr=12 !(nr=7) ! number of bins
108	real,save :: rd(nr)
109	real qr(ngrid,nlay,nr)
110	real,save :: rdi(nr+1) ! extreme and intermediate radii
111	real Sq(ngrid,nlay)
112	real,parameter :: rdmin=1.e-8
113	real,parameter :: rdmax=30.e-6
114	real,parameter :: rdimin=1.e-8 ! 1.e-7
115	real,parameter :: rdimax=1.e-4
116
117	c 2) Second size distribution for the log-normal integration
118	c (the mass mixing ratio is computed for each radius)
119
120	integer iint
121	integer,parameter :: ninter=4 ! number of points between each rdi radii
122	real,save :: rr(ninter,nr)
123	integer radpower
124	real sigma0
125
126	c 3) Other local variables used in doubleq
127
128	INTEGER,SAVE :: idust_mass ! index of tracer containing dust mass
129	! mix. ratio
130	INTEGER,SAVE :: idust_number ! index of tracer containing dust number
131	! mix. ratio
132	INTEGER,SAVE :: iccn_mass ! index of tracer containing CCN mass
133	! mix. ratio
134	INTEGER,SAVE :: iccn_number ! index of tracer containing CCN number
135	! mix. ratio
136	INTEGER,SAVE :: istormdust_mass ! index of tracer containing
137	!stormdust mass mix. ratio
138	INTEGER,SAVE :: istormdust_number ! index of tracer containing
139	!stormdust number mix. ratio
140	INTEGER,SAVE :: iccnco2_number ! index of tracer containing CCN number
141	INTEGER,SAVE :: iccnco2_mass ! index of tracer containing CCN number
142	INTEGER,SAVE :: ico2_ice ! index of tracer containing CCN number
143
144
145	LOGICAL,SAVE :: firstcall=.true.
146
147
148
149	c ** un petit test de coherence
150	c --------------------------
151	! AS: firstcall OK absolute
152	IF (firstcall) THEN
153
154	c Doubleq: initialization
155	IF (doubleq) THEN
156	do ir=1,nr
157	rd(ir)= rdmin(rdmax/rdmin)*(float(ir-1)/float(nr-1))
158	end do
159	rdi(1)=rdimin
160	do ir=2,nr
161	rdi(ir)= sqrt(rd(ir-1)*rd(ir))
162	end do
163	rdi(nr+1)=rdimax
164
165	do ir=1,nr
166	do iint=1,ninter
167	rr(iint,ir)=
168	& rdi(ir)*
169	& (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
170	c write(,) rr(iint,ir)
171	end do
172	end do
173
174	! identify tracers corresponding to mass mixing ratio and
175	! number mixing ratio
176	idust_mass=0 ! dummy initialization
177	idust_number=0 ! dummy initialization
178
179	do iq=1,nq
180	if (noms(iq).eq."dust_mass") then
181	idust_mass=iq
182	write(,)"callsedim: idust_mass=",idust_mass
183	endif
184	if (noms(iq).eq."dust_number") then
185	idust_number=iq
186	write(,)"callsedim: idust_number=",idust_number
187	endif
188	enddo
189
190	! check that we did find the tracers
191	if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
192	write(,) 'callsedim: error! could not identify'
193	write(,) ' tracers for dust mass and number mixing'
194	write(,) ' ratio and doubleq is activated!'
195	stop
196	endif
197	ENDIF !of if (doubleq)
198
199	IF (microphys) THEN
200	iccn_mass=0
201	iccn_number=0
202	do iq=1,nq
203	if (noms(iq).eq."ccn_mass") then
204	iccn_mass=iq
205	write(,)"callsedim: iccn_mass=",iccn_mass
206	endif
207	if (noms(iq).eq."ccn_number") then
208	iccn_number=iq
209	write(,)"callsedim: iccn_number=",iccn_number
210	endif
211	enddo
212	! check that we did find the tracers
213	if ((iccn_mass.eq.0).or.(iccn_number.eq.0)) then
214	write(,) 'callsedim: error! could not identify'
215	write(,) ' tracers for ccn mass and number mixing'
216	write(,) ' ratio and microphys is activated!'
217	stop
218	endif
219	ENDIF !of if (microphys)
220
221	IF (co2clouds) THEN
222	iccnco2_mass=0
223	iccnco2_number=0
224	ico2_ice=0
225	do iq=1,nq
226	if (noms(iq).eq."ccnco2_mass") then
227	iccnco2_mass=iq
228	write(,)"callsedim: iccnco2_mass=",iccnco2_mass
229	endif
230	if (noms(iq).eq."co2_ice") then
231	ico2_ice=iq
232	write(,)"callsedim: ico2_ice=",ico2_ice
233	endif
234	if (noms(iq).eq."ccnco2_number") then
235	iccnco2_number=iq
236	write(,)"callsedim: iccnco2_number=",iccnco2_number
237	endif
238	enddo
239	! check that we did find the tracers
240	if ((iccnco2_mass.eq.0).or.(iccnco2_number.eq.0)) then
241	write(,) 'callsedim: error! could not identify'
242	write(,) ' tracers for ccn co2 mass and number mixing'
243	write(,) ' ratio and co2clouds are activated!'
244	stop
245	endif
246	ENDIF !of if (co2clouds)
247
248	IF (water) THEN
249	write(,) "correction for the shape of the ice particles ?"
250	beta=0.75 ! default value
251	call getin("ice_shape",beta)
252	write(,) " ice_shape = ",beta
253
254	write(,) "water_param nueff Sedimentation:", nuice_sed
255	IF (activice) THEN
256	write(,) "water_param nueff Radiative:", nuice_ref
257	ENDIF
258	ENDIF
259
260	IF (rdstorm) THEN ! identifying stormdust tracers for sedimentation
261	istormdust_mass=0 ! dummy initialization
262	istormdust_number=0 ! dummy initialization
263
264	do iq=1,nq
265	if (noms(iq).eq."stormdust_mass") then
266	istormdust_mass=iq
267	write(,)"callsedim: istormdust_mass=",istormdust_mass
268	endif
269	if (noms(iq).eq."stormdust_number") then
270	istormdust_number=iq
271	write(,)"callsedim: istormdust_number=",
272	& istormdust_number
273	endif
274	enddo
275
276	! check that we did find the tracers
277	if ((istormdust_mass.eq.0).or.(istormdust_number.eq.0)) then
278	write(,) 'callsedim: error! could not identify'
279	write(,) ' tracers for stormdust mass and number mixing'
280	write(,) ' ratio and rdstorm is activated!'
281	stop
282	endif
283	ENDIF !of if (rdstorm)
284
285	firstcall=.false.
286	ENDIF ! of IF (firstcall)
287
288	c-----------------------------------------------------------------------
289	c 1. Initialization
290	c -----------------
291
292	! zqi(1:ngrid,1:nlay,1:nqmx) = 0.
293	c Update the mass mixing ratio and temperature with the tendencies coming
294	c from other parameterizations:
295	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
296	zqi(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)
297	& +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
298	zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay)
299	& +pdt(1:ngrid,1:nlay)*ptimestep
300
301	c Computing the different layer properties
302	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
303	c Mass (kg.m-2), thickness(m), crossing time (s) etc.
304
305	do l=1,nlay
306	do ig=1, ngrid
307	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
308	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
309	end do
310	end do
311
312	c =================================================================
313	c Compute the geometric mean radius used for sedimentation
314
315	if (doubleq) then
316	do l=1,nlay
317	do ig=1, ngrid
318
319	call updaterdust(zqi(ig,l,igcm_dust_mass),
320	& zqi(ig,l,igcm_dust_number),r0dust(ig,l),
321	& tauscaling(ig))
322
323	end do
324	end do
325	endif
326	! rocket dust storm
327	if (rdstorm) then
328	do l=1,nlay
329	do ig=1, ngrid
330
331	call updaterdust(zqi(ig,l,igcm_stormdust_mass),
332	& zqi(ig,l,igcm_stormdust_number),r0stormdust(ig,l),
333	& tauscaling(ig))
334
335	end do
336	end do
337	endif
338	c =================================================================
339	do iq=1,nq
340	if(radius(iq).gt.1.e-9 .and.(iq.ne.ico2_ice) .and.
341	& (iq .ne. iccnco2_mass) .and. (iq .ne.
342	& iccnco2_number)) then ! no sedim for gaz or CO2 clouds (done in microtimestep)
343
344	c -----------------------------------------------------------------
345	c DOUBLEQ CASE
346	c -----------------------------------------------------------------
347	if ((doubleq.and.
348	& ((iq.eq.idust_mass).or.(iq.eq.idust_number).or.
349	& (iq.eq.istormdust_mass).or.(iq.eq.istormdust_number)))) then
350
351	c Computing size distribution:
352	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
353
354	if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
355	do l=1,nlay
356	do ig=1, ngrid
357	r0(ig,l)=r0dust(ig,l)
358	end do
359	end do
360	else if ((iq.eq.istormdust_mass).or.
361	& (iq.eq.istormdust_number)) then
362	do l=1,nlay
363	do ig=1, ngrid
364	r0(ig,l)=r0stormdust(ig,l)
365	end do
366	end do
367	endif
368	sigma0 = varian
369
370	c Computing mass mixing ratio for each particle size
371	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
372	IF ((iq.EQ.idust_mass).or.(iq.EQ.istormdust_mass)) then
373	radpower = 2
374	ELSE ! number
375	radpower = -1
376	ENDIF
377	Sq(1:ngrid,1:nlay) = 0.
378	do ir=1,nr
379	do l=1,nlay
380	do ig=1,ngrid
381	c ****************
382	c Size distribution integration
383	c (Trapezoid Integration Method)
384	qr(ig,l,ir)=0.5(rr(2,ir)-rr(1,ir))
385	& (rr(1,ir)*radpower)
386	& exp(-(log(rr(1,ir)/r0(ig,l)))*2/(2sigma0**2))
387	do iint=2,ninter-1
388	qr(ig,l,ir)=qr(ig,l,ir) +
389	& 0.5(rr(iint+1,ir)-rr(iint-1,ir))
390	& (rr(iint,ir)*radpower)
391	& exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
392	& (2sigma0*2))
393	end do
394	qr(ig,l,ir)=qr(ig,l,ir) +
395	& 0.5(rr(ninter,ir)-rr(ninter-1,ir))
396	& (rr(ninter,ir)*radpower)
397	& exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
398	& (2sigma0*2))
399
400	c **************** old method (not recommended!)
401	c qr(ig,l,ir)=(rd(ir)*(5-3iq))*
402	c & exp( -(log(rd(ir)/r0(ig,l)))*2 / (2sigma0**2) )
403	c ******************************
404
405	Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
406	enddo
407	enddo
408	enddo
409
410	do ir=1,nr
411	do l=1,nlay
412	do ig=1,ngrid
413	qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
414	enddo
415	enddo
416	enddo
417
418	c Computing sedimentation for each tracer
419	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
420
421	zqi(1:ngrid,1:nlay,iq) = 0.
422	pdqs_sed(1:ngrid,iq) = 0.
423
424	do ir=1,nr
425	call newsedim(ngrid,nlay,1,1,ptimestep,
426	& pplev,masse,epaisseur,zt,rd(ir),(/rho_dust/),qr(1,1,ir),
427	& wq,0.5)
428
429	c Tendencies
430	c ~~~~~~~~~~
431	do ig=1,ngrid
432	pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
433	& + wq(ig,1)/ptimestep
434	end do
435	DO l = 1, nlay
436	DO ig=1,ngrid
437	zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
438	ENDDO
439	ENDDO
440	enddo ! of do ir=1,nr
441	c -----------------------------------------------------------------
442	c WATER CYCLE CASE
443	c -----------------------------------------------------------------
444	else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
445	& .or. (iq .eq. igcm_h2o_ice)) then
446	if (microphys) then
447	! water ice sedimentation
448	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
449	& ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rhocloud,
450	& zqi(1,1,iq),wq,beta)
451	else
452	! water ice sedimentation
453	call newsedim(ngrid,nlay,ngrid*nlay,1,
454	& ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rho_q(iq),
455	& zqi(1,1,iq),wq,beta)
456	endif ! of if (microphys)
457	c Tendencies
458	c ~~~~~~~~~~
459	do ig=1,ngrid
460	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
461	end do
462	c -----------------------------------------------------------------
463	c GENERAL CASE
464	c -----------------------------------------------------------------
465	else
466	call newsedim(ngrid,nlay,1,1,ptimestep,
467	& pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
468	& zqi(1,1,iq),wq,1.0)
469	c Tendencies
470	c ~~~~~~~~~~
471	do ig=1,ngrid
472	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
473	end do
474	endif ! of if doubleq and if water
475	c -----------------------------------------------------------------
476
477	c Compute the final tendency:
478	c ---------------------------
479	DO l = 1, nlay
480	DO ig=1,ngrid
481	pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
482	$ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
483	ENDDO
484	ENDDO
485
486	endif ! of if(radius(iq).gt.1.e-9)
487	c =================================================================
488	enddo ! of do iq=1,nq
489
490	c Update the dust particle size "rdust"
491	c -------------------------------------
492	if (doubleq) then
493	DO l = 1, nlay
494	DO ig=1,ngrid
495
496
497	call updaterdust(zqi(ig,l,igcm_dust_mass),
498	& zqi(ig,l,igcm_dust_number),rdust(ig,l),
499	& tauscaling(ig))
500
501
502	ENDDO
503	ENDDO
504	endif ! of if (doubleq)
505
506	if (rdstorm) then
507	DO l = 1, nlay
508	DO ig=1,ngrid
509	call updaterdust(zqi(ig,l,igcm_stormdust_mass),
510	& zqi(ig,l,igcm_stormdust_number),rstormdust(ig,l),
511	& tauscaling(ig))
512	ENDDO
513	ENDDO
514	endif ! of if (rdstorm)
515
516	c Update the ice particle size "rice"
517	c -------------------------------------
518	if (water) then
519	IF(microphys) THEN
520
521
522	DO l = 1, nlay
523	DO ig=1,ngrid
524
525	call updaterice_micro(zqi(ig,l,igcm_h2o_ice),
526	& zqi(ig,l,igcm_ccn_mass),zqi(ig,l,igcm_ccn_number),
527	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
528
529	ENDDO
530	ENDDO
531
532	ELSE
533
534	DO l = 1, nlay
535	DO ig=1,ngrid
536
537	call updaterice_typ(zqi(ig,l,igcm_h2o_ice),
538	& tau(ig,1),zlay(ig,l),rice(ig,l))
539
540	ENDDO
541	ENDDO
542	ENDIF ! of IF(microphys)
543	endif ! of if (water)
544
545	END SUBROUTINE callsedim
546
547	END MODULE callsedim_mod
548

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

Download in other formats:

Original Format