Context Navigation

callsedim.F @ 503

Last change on this file since 503 was 472, checked in by emillour, 14 years ago

Mars GCM:

fixed bug in callsedim.F: recomputation of rdust and rice should only be done

if the appropriate tracers and flags are set.

commented out call to uertst in ch.F (uertst crashes with gfortran and is

moreover puerly diagnostic messages).

In nltecool.F: Enforced using igcm_* indexes to identify tracers and not

hard coded fixed numbers (which would lead to a disaster in cases where the
order of tracers was changed).

File size: 15.3 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	if (doubleq) then
421	DO l = 1, nlay
422	DO ig=1,ngrid
423	rdust(ig,l)=
424	& CBRT(r3n_q*zqi(ig,l,idust_mass)/
425	& max(zqi(ig,l,idust_number),0.01))
426	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
427	ENDDO
428	ENDDO
429	endif !of if (doubleq)
430
431	c Update the ice particle size "rice"
432	c -------------------------------------
433	IF(scavenging) THEN
434	DO l = 1, nlay
435	DO ig=1,ngrid
436	Mo = zqi(ig,l,igcm_h2o_ice) +
437	& zqi(ig,l,iccn_mass)* tauscaling(ig) + 1.e-30
438	No = zqi(ig,l,iccn_number)* tauscaling(ig)+ 1.e-30
439	rhocloud(ig,l) = zqi(ig,l,igcm_h2o_ice) / Mo * rho_ice
440	& +zqi(ig,l,iccn_mass)* tauscaling(ig)
441	& / Mo * rho_dust
442	rhocloud(ig,l) =
443	& min(max(rhocloud(ig,l),rho_ice),rho_dust)
444	rice(ig,l) =
445	& ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
446	if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
447	c print*, "Mice,Mo, No",zqi(ig,l,igcm_h2o_ice),Mo, No
448	c print*, "rice, rho apres", rice(ig,l), rhocloud(ig,l)
449
450	ENDDO
451	ENDDO
452	ELSE
453	if (water) then
454	DO l = 1, nlay
455	DO ig=1,ngrid
456	ccntyp =
457	& 1.3e+8max(tau(ig,1),0.001)/0.1exp(-zlay(ig,l)/10000.)
458	ccntyp = ccntyp /ccn_factor
459	rice(ig,l)=max( CBRT ( (zqi(ig,l,igcm_h2o_ice)/rho_ice
460	& +ccntyp(4./3.)pirdust(ig,l)*3.)
461	& /(ccntyp4./3.pi) ), rdust(ig,l))
462	ENDDO
463	ENDDO
464	endif ! of if (water)
465	ENDIF ! of IF (scavenging)
466
467	RETURN
468	END
469

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

Download in other formats:

Original Format