Context Navigation

callsedim.F @ 706

Last change on this file since 706 was 635, checked in by emillour, 13 years ago

Mars GCM: Update of the chemistry package, including:

93 reactions are accounted for (instead of 22); tracking 28 species (instead of 11)
computation of photoabsorption using raytracing
improved time stepping in the photochemistry
updated parameters (cross-sections); with this new version input files

are in 'EUV/param_v5' of "datafile" directory.

transition between lower and upper atmosphere chemistry set to 0.1 Pa (calchim.F90)
Lots of code clean-up: removed obsolete files column.F, param_v3.h, flujo.F, phdisrate.F, ch.F, interpfast.F, paramfoto.F, getch.F Converted chemtermos.F -> chemthermos.F90 and euvheat.F -> euvheat.F90. Added paramfoto_compact.F , param_v4.h and iono.h
Upadted surfacearea.F
Cleaned initracer.F and callkeys.h (removed use of obsolete "nqchem" and "oldnames" case when initializing tracers).
Minor correction in "callsedim": compute "rdust" and/or "rice" only when it makes sense.

FGG+FL+EM

File size: 16.0 KB

Line
1	SUBROUTINE callsedim(ngrid,nlay, ptimestep,
2	& pplev,zlev, zlay, pt, rdust, rice,
3	& rsedcloud,rhocloud,
4	& pq, pdqfi, pdqsed,pdqs_sed,nq,
5	& tau, tauscaling)
6	! to use 'getin'
7	USE ioipsl_getincom
8	IMPLICIT NONE
9
10	c=======================================================================
11	c Sedimentation of the Martian aerosols
12	c depending on their density and radius
13	c
14	c F.Forget 1999
15	c
16	c Modified by J.-B. Madeleine 2010: Now includes the doubleq
17	c technique in order to have only one call to callsedim in
18	c physiq.F.
19	c
20	c=======================================================================
21
22	c-----------------------------------------------------------------------
23	c declarations:
24	c -------------
25
26	#include "dimensions.h"
27	#include "dimphys.h"
28	#include "comcstfi.h"
29	#include "tracer.h"
30	#include "callkeys.h"
31
32	c
33	c arguments:
34	c ----------
35
36	INTEGER ngrid ! number of horizontal grid points
37	INTEGER nlay ! number of atmospheric layers
38	REAL ptimestep ! physics time step (s)
39	REAL pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
40	REAL pt(ngrid,nlay) ! temperature at mid-layer (K)
41	REAL zlev(ngrid,nlay+1) ! altitude at layer boundaries
42	c Aerosol radius provided by the water ice microphysical scheme:
43	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
44	REAL rice(ngrid,nlay) ! Ice geometric mean radius (m)
45
46	c Traceurs :
47	integer nq ! number of tracers
48	real pq(ngrid,nlay,nq) ! tracers (kg/kg)
49	real pdqfi(ngrid,nlay,nq) ! tendency before sedimentation (kg/kg.s-1)
50	real pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
51	real pdqs_sed(ngrid,nq) ! flux at surface (kg.m-2.s-1)
52
53	c local:
54	c ------
55
56	REAL CBRT
57	EXTERNAL CBRT
58
59	INTEGER l,ig, iq
60
61	real zqi(ngridmx,nlayermx,nqmx) ! to locally store tracers
62	real masse (ngridmx,nlayermx) ! Layer mass (kg.m-2)
63	real epaisseur (ngridmx,nlayermx) ! Layer thickness (m)
64	real wq(ngridmx,nlayermx+1) ! displaced tracer mass (kg.m-2)
65	real r0(ngridmx,nlayermx) ! geometric mean radius used for
66	! sedimentation (m)
67	real r0dust(ngridmx,nlayermx) ! geometric mean radius used for
68	! dust (m)
69	real r0ccn(ngridmx,nlayermx) ! geometric mean radius used for
70	! CCNs (m)
71	c Sedimentation radius of water ice
72	real rsedcloud(ngridmx,nlayermx)
73	real beta ! correction for the shape of the ice particles (cf. newsedim)
74	save beta
75	c Cloud density (kg.m-3)
76	real rhocloud(ngridmx,nlayermx)
77
78	c for ice radius computation
79	REAL Mo,No
80	REAL tau(ngrid,nlay), tauscaling(ngrid)
81	REAL zlay(ngrid,nlay) ! altitude at the middle of the layers
82	REAl ccntyp
83
84
85
86	c Discrete size distributions (doubleq)
87	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
88	c 1) Parameters used to represent the changes in fall
89	c velocity as a function of particle size;
90	integer nr,ir
91	parameter (nr=12) !(nr=7) ! number of bins
92	real rd(nr),qr(ngridmx,nlayermx,nr)
93	real rdi(nr+1) ! extreme and intermediate radii
94	real Sq(ngridmx,nlayermx)
95	real rdmin,rdmax,rdimin,rdimax
96	data rdmin/1.e-8/ !/1.e-7/
97	data rdmax/30.e-6/
98	data rdimin/1.e-10/
99	data rdimax/1e-4/
100	save rd, rdi
101
102	c 2) Second size distribution for the log-normal integration
103	c (the mass mixing ratio is computed for each radius)
104
105	integer ninter, iint
106	parameter (ninter=4) ! nombre de point entre chaque rayon rdi
107	real rr(ninter,nr)
108	save rr
109	integer radpower
110	real sigma0
111
112	c 3) Other local variables used in doubleq
113
114	INTEGER idust_mass ! index of tracer containing dust mass
115	! mix. ratio
116	INTEGER idust_number ! index of tracer containing dust number
117	! mix. ratio
118	INTEGER iccn_mass ! index of tracer containing CCN mass
119	! mix. ratio
120	INTEGER iccn_number ! index of tracer containing CCN number
121	! mix. ratio
122	SAVE idust_mass,idust_number
123	SAVE iccn_mass,iccn_number
124
125
126	c Firstcall:
127
128	LOGICAL firstcall
129	SAVE firstcall
130	DATA firstcall/.true./
131
132	c ** un petit test de coherence
133	c --------------------------
134
135	IF (firstcall) THEN
136
137	IF(ngrid.NE.ngridmx) THEN
138	PRINT*,'STOP dans callsedim'
139	PRINT*,'probleme de dimensions :'
140	PRINT*,'ngrid =',ngrid
141	PRINT*,'ngridmx =',ngridmx
142	STOP
143	ENDIF
144
145	c Doubleq: initialization
146	IF (doubleq) THEN
147	do ir=1,nr
148	rd(ir)= rdmin(rdmax/rdmin)*(float(ir-1)/float(nr-1))
149	end do
150	rdi(1)=rdimin
151	do ir=2,nr
152	rdi(ir)= sqrt(rd(ir-1)*rd(ir))
153	end do
154	rdi(nr+1)=rdimax
155
156	do ir=1,nr
157	do iint=1,ninter
158	rr(iint,ir)=
159	& rdi(ir)*
160	& (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
161	c write(,) rr(iint,ir)
162	end do
163	end do
164
165	! identify tracers corresponding to mass mixing ratio and
166	! number mixing ratio
167	idust_mass=0 ! dummy initialization
168	idust_number=0 ! dummy initialization
169
170	do iq=1,nq
171	if (noms(iq).eq."dust_mass") then
172	idust_mass=iq
173	endif
174	if (noms(iq).eq."dust_number") then
175	idust_number=iq
176	endif
177	enddo
178
179	! check that we did find the tracers
180	if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
181	write(,) 'callsedim: error! could not identify'
182	write(,) ' tracers for dust mass and number mixing'
183	write(,) ' ratio and doubleq is activated!'
184	stop
185	endif
186	ENDIF !of if (doubleq)
187
188	IF (scavenging) THEN
189	iccn_mass=0
190	iccn_number=0
191	do iq=1,nq
192	if (noms(iq).eq."ccn_mass") then
193	iccn_mass=iq
194	endif
195	if (noms(iq).eq."ccn_number") then
196	iccn_number=iq
197	endif
198	enddo
199	! check that we did find the tracers
200	if ((iccn_mass.eq.0).or.(iccn_number.eq.0)) then
201	write(,) 'callsedim: error! could not identify'
202	write(,) ' tracers for ccn mass and number mixing'
203	write(,) ' ratio and scavenging is activated!'
204	stop
205	endif
206	ENDIF !of if (scavenging)
207
208	IF (water) THEN
209	write(,) "correction for the shape of the ice particles ?"
210	beta=0.75 ! default value
211	call getin("ice_shape",beta)
212	write(,) " ice_shape = ",beta
213
214	write(,) "water_param nueff Sedimentation:", nuice_sed
215	IF (activice) THEN
216	write(,) "water_param nueff Radiative:", nuice_ref
217	ENDIF
218	ENDIF
219
220	firstcall=.false.
221	ENDIF ! of IF (firstcall)
222
223	c-----------------------------------------------------------------------
224	c 1. Initialization
225	c -----------------
226
227	zqi(1:ngrid,1:nlay,1:nqmx) = 0.
228	c Updating the mass mixing ratio with the tendencies coming
229	c from other parameterizations:
230	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
231
232	do iq=1,nq
233	do l=1,nlay
234	do ig=1,ngrid
235	zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
236	enddo
237	enddo
238	enddo
239
240	c Computing the different layer properties
241	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
242	c Mass (kg.m-2), thickness(m), crossing time (s) etc.
243
244	do l=1,nlay
245	do ig=1, ngrid
246	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
247	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
248	end do
249	end do
250
251	c =================================================================
252	c Compute the geometric mean radius used for sedimentation
253
254	if (doubleq) then
255	do l=1,nlay
256	do ig=1, ngrid
257	r0dust(ig,l) =
258	& CBRT(r3n_q*zqi(ig,l,idust_mass)/
259	& max(zqi(ig,l,idust_number),0.01))
260	r0dust(ig,l)=min(max(r0dust(ig,l),1.e-10),500.e-6)
261	end do
262	end do
263	endif
264	if (scavenging) then
265	do l=1,nlay
266	do ig=1, ngrid
267	r0ccn(ig,l) = rsedcloud(ig,l)/(1.+nuice_sed)**4.5
268	end do
269	end do
270	endif
271
272	c =================================================================
273	do iq=1,nq
274	if(radius(iq).gt.1.e-9) then ! no sedim for gaz
275
276	c -----------------------------------------------------------------
277	c DOUBLEQ CASE
278	c -----------------------------------------------------------------
279	if ((doubleq.and.
280	& ((iq.eq.idust_mass).or.
281	& (iq.eq.idust_number)))) then !.or.
282	! & (scavenging.and.
283	! & ((iq.eq.iccn_mass).or.
284	! & (iq.eq.iccn_number)))) then
285
286	c Computing size distribution:
287	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
288
289	c if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
290	do l=1,nlay
291	do ig=1, ngrid
292	r0(ig,l)=r0dust(ig,l)
293	end do
294	end do
295	sigma0 = varian
296	c else ! ccn
297	c do l=1,nlay
298	c do ig=1, ngrid
299	c r0(ig,l)=r0ccn(ig,l)
300	c end do
301	c end do
302	c sigma0 = sqrt(log(1.+nuice_sed))
303	c endif
304
305	c Computing mass mixing ratio for each particle size
306	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
307	IF ((iq.EQ.idust_mass).or.(iq.EQ.iccn_mass)) then
308	radpower = 2
309	ELSE ! number
310	radpower = -1
311	ENDIF
312	Sq(1:ngrid,1:nlay) = 0.
313	do ir=1,nr
314	do l=1,nlay
315	do ig=1,ngrid
316	c ****************
317	c Size distribution integration
318	c (Trapezoid Integration Method)
319	qr(ig,l,ir)=0.5(rr(2,ir)-rr(1,ir))
320	& (rr(1,ir)*radpower)
321	& exp(-(log(rr(1,ir)/r0(ig,l)))*2/(2sigma0**2))
322	do iint=2,ninter-1
323	qr(ig,l,ir)=qr(ig,l,ir) +
324	& 0.5(rr(iint+1,ir)-rr(iint-1,ir))
325	& (rr(iint,ir)*radpower)
326	& exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
327	& (2sigma0*2))
328	end do
329	qr(ig,l,ir)=qr(ig,l,ir) +
330	& 0.5(rr(ninter,ir)-rr(ninter-1,ir))
331	& (rr(ninter,ir)*radpower)
332	& exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
333	& (2sigma0*2))
334
335	c **************** old method (not recommended!)
336	c qr(ig,l,ir)=(rd(ir)*(5-3iq))*
337	c & exp( -(log(rd(ir)/r0(ig,l)))*2 / (2sigma0**2) )
338	c ******************************
339
340	Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
341	enddo
342	enddo
343	enddo
344
345	do ir=1,nr
346	do l=1,nlay
347	do ig=1,ngrid
348	qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
349	enddo
350	enddo
351	enddo
352
353	c Computing sedimentation for each tracer
354	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
355
356	zqi(1:ngrid,1:nlay,iq) = 0.
357	pdqs_sed(1:ngrid,iq) = 0.
358
359	do ir=1,nr
360	! IF ((iq.EQ.idust_mass).or.(iq.EQ.idust_number)) then
361	! Dust sedimentation
362	call newsedim(ngrid,nlay,1,1,ptimestep,
363	& pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),
364	& wq,0.5)
365	! ELSE
366	! ! CCN sedimentation
367	! call newsedim(ngrid,nlay,1,ngrid*nlay,ptimestep,
368	! & pplev,masse,epaisseur,pt,rd(ir),rhocloud,qr(1,1,ir),
369	! & wq,beta)
370	! ENDIF
371
372	c Tendencies
373	c ~~~~~~~~~~
374	do ig=1,ngrid
375	pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
376	& + wq(ig,1)/ptimestep
377	end do
378	DO l = 1, nlay
379	DO ig=1,ngrid
380	zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
381	ENDDO
382	ENDDO
383	enddo ! of do ir=1,nr
384	c -----------------------------------------------------------------
385	c WATER CYCLE CASE
386	c -----------------------------------------------------------------
387	c else if (water.and.(iq.eq.igcm_h2o_ice)) then
388	else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
389	& .or. (iq .eq. igcm_h2o_ice)) then
390	if (microphys) then
391	! water ice sedimentation
392	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
393	& ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud,
394	& zqi(1,1,iq),wq,beta)
395	else
396	! water ice sedimentation
397	call newsedim(ngrid,nlay,ngrid*nlay,1,
398	& ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq),
399	& zqi(1,1,iq),wq,beta)
400	endif ! of if (microphys)
401	c Tendencies
402	c ~~~~~~~~~~
403	do ig=1,ngrid
404	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
405	end do
406	c -----------------------------------------------------------------
407	c GENERAL CASE
408	c -----------------------------------------------------------------
409	else
410	call newsedim(ngrid,nlay,1,1,ptimestep,
411	& pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
412	& zqi(1,1,iq),wq,1.0)
413	c Tendencies
414	c ~~~~~~~~~~
415	do ig=1,ngrid
416	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
417	end do
418	endif ! of if doubleq and if water
419	c -----------------------------------------------------------------
420
421	c Compute the final tendency:
422	c ---------------------------
423	DO l = 1, nlay
424	DO ig=1,ngrid
425	pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
426	$ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
427	ENDDO
428	ENDDO
429
430	endif ! of if(radius(iq).gt.1.e-9)
431	c =================================================================
432	enddo ! of do iq=1,nq
433
434	c Update the dust particle size "rdust"
435	c -------------------------------------
436	if (doubleq) then
437	DO l = 1, nlay
438	DO ig=1,ngrid
439	rdust(ig,l)=
440	& CBRT(r3n_q*zqi(ig,l,idust_mass)/
441	& max(zqi(ig,l,idust_number),0.01))
442	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
443	ENDDO
444	ENDDO
445	endif ! of if (doubleq)
446
447	c Update the ice particle size "rice"
448	c -------------------------------------
449	if (water) then
450	IF(scavenging) THEN
451	DO l = 1, nlay
452	DO ig=1,ngrid
453	Mo = zqi(ig,l,igcm_h2o_ice) +
454	& zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30
455	No = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30
456	rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice
457	& +zqi(ig,l,iccn_mass)* tauscaling(ig)
458	& / Mo * rho_dust
459	rhocloud(ig,l) =
460	& min(max(rhocloud(ig,l),rho_ice),rho_dust)
461	rice(ig,l) =
462	& ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
463	if ((Mo.lt.1.e-15) .or. (No.le.1)) rice(ig,l) = 1.e-8
464	! print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No
465	! print*, "rice, rho apres", rice(ig,l), rhocloud(ig,l)
466
467	ENDDO
468	ENDDO
469	ELSE
470	DO l = 1, nlay
471	DO ig=1,ngrid
472	ccntyp =
473	& 1.3e+8max(tau(ig,1),0.001)/0.1exp(-zlay(ig,l)/10000.)
474	ccntyp = ccntyp /ccn_factor
475	rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice
476	& +ccntyp(4./3.)pirdust(ig,l)*3.)
477	& /(ccntyp4./3.pi) ), rdust(ig,l))
478	ENDDO
479	ENDDO
480	ENDIF ! of IF(scavenging)
481	endif ! of if (water)
482
483	RETURN
484	END
485

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

Download in other formats:

Original Format