Context Navigation

callsedim.F @ 414

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

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

Download in other formats:

Original Format