Context Navigation

callsedim.F @ 463

Last change on this file since 463 was 420, checked in by tnavarro, 14 years ago

24/11/11 == TN

corrected minor bug in updatereffrad.F : reffdust was not saved

ccn_factor as not correctly used in sedimentation.

It is now initialized in inifis.F, declared in tracer.h and
used in both simpleclouds.F & callsedim.F to update ice radius.

Commented diagfi outputs in aeropacity.F & improvedclouds.F for non scavenging users.

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

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

Download in other formats:

Original Format