Context Navigation

callsedim.F @ 357

Last change on this file since 357 was 82, checked in by aslmd, 14 years ago

mars + LMD_MM_MARS : modifications pour cycle de l'eau, valeurs tunees JBM

used settings reached by JBM to obtain his PhD results

alb_surfice = 0.45 --- in physiq.F and meso_physiq.F
ccn_factor = 4.5 --- in watercloud.F
nuice_sed = 0.45 --- in callsedim.F

NB: this is supposed to be further refined in the future

File size: 11.7 KB

Line
1	SUBROUTINE callsedim(ngrid,nlay, ptimestep,
2	$ pplev,zlev, pt, rdust, rice,
3	& pq, pdqfi, pdqsed,pdqs_sed,nq)
4	IMPLICIT NONE
5
6	c=======================================================================
7	c Sedimentation of the Martian aerosols
8	c depending on their density and radius
9	c
10	c F.Forget 1999
11	c
12	c Modified by J.-B. Madeleine 2010: Now includes the doubleq
13	c technique in order to have only one call to callsedim in
14	c physiq.F.
15	c
16	c=======================================================================
17
18	c-----------------------------------------------------------------------
19	c declarations:
20	c -------------
21
22	#include "dimensions.h"
23	#include "dimphys.h"
24	#include "comcstfi.h"
25	#include "tracer.h"
26	#include "callkeys.h"
27
28	c
29	c arguments:
30	c ----------
31
32	INTEGER ngrid ! number of horizontal grid points
33	INTEGER nlay ! number of atmospheric layers
34	REAL ptimestep ! physics time step (s)
35	REAL pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
36	REAL pt(ngrid,nlay) ! temperature at mid-layer (K)
37	REAL zlev(ngrid,nlay+1) ! altitude at layer boundaries
38	c Aerosol radius provided by the water ice microphysical scheme:
39	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
40	REAL rice(ngrid,nlay) ! Ice geometric mean radius (m)
41
42	c Traceurs :
43	integer nq ! number of tracers
44	real pq(ngrid,nlay,nq) ! tracers (kg/kg)
45	real pdqfi(ngrid,nlay,nq) ! tendency before sedimentation (kg/kg.s-1)
46	real pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
47	real pdqs_sed(ngrid,nq) ! flux at surface (kg.m-2.s-1)
48
49	c local:
50	c ------
51
52	REAL CBRT
53	EXTERNAL CBRT
54
55	INTEGER l,ig, iq
56
57	real zqi(ngridmx,nlayermx,nqmx) ! to locally store tracers
58	real masse (ngridmx,nlayermx) ! Layer mass (kg.m-2)
59	real epaisseur (ngridmx,nlayermx) ! Layer thickness (m)
60	real wq(ngridmx,nlayermx+1) ! displaced tracer mass (kg.m-2)
61	real r0(ngridmx,nlayermx) ! geometric mean doubleq radius (m)
62	c Sedimentation radius of water ice
63	real rfall(ngridmx,nlayermx)
64	c Sedimentation effective variance of water ice
65	REAL, PARAMETER :: nuice_sed = 0.45 !! TESTS_JB !! 0.1 avant
66
67	c Discrete size distributions (doubleq)
68	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
69	c 1) Parameters used to represent the changes in fall
70	c velocity as a function of particle size;
71	integer nr,ir
72	parameter (nr=12) !(nr=7) ! number of bins
73	real rd(nr),qr(ngridmx,nlayermx,nr)
74	real rdi(nr+1) ! extreme and intermediate radii
75	real Sq(ngridmx,nlayermx)
76	real rdmin,rdmax,rdimin,rdimax
77	data rdmin/1.e-8/ !/1.e-7/
78	data rdmax/30.e-6/
79	data rdimin/1.e-10/
80	data rdimax/1e-4/
81	save rd, rdi
82
83	c 2) Second size distribution for the log-normal integration
84	c (the mass mixing ratio is computed for each radius)
85
86	integer ninter, iint
87	parameter (ninter=4) ! nombre de point entre chaque rayon rdi
88	real rr(ninter,nr)
89	save rr
90	integer radpower
91
92	c 3) Other local variables used in doubleq
93
94	real reff(ngridmx,nlayermx,2) ! for diagnostic only
95	INTEGER idust_mass ! index of tracer containing dust mass
96	! mix. ratio
97	INTEGER idust_number ! index of tracer containing dust number
98	! mix. ratio
99	SAVE idust_mass,idust_number
100
101	c Firstcall:
102
103	LOGICAL firstcall
104	SAVE firstcall
105	DATA firstcall/.true./
106
107	c ** un petit test de coherence
108	c --------------------------
109
110	IF (firstcall) THEN
111	IF(ngrid.NE.ngridmx) THEN
112	PRINT*,'STOP dans callsedim'
113	PRINT*,'probleme de dimensions :'
114	PRINT*,'ngrid =',ngrid
115	PRINT*,'ngridmx =',ngridmx
116	STOP
117	ENDIF
118
119	c Doubleq: initialization
120	IF (doubleq) THEN
121	do ir=1,nr
122	rd(ir)= rdmin(rdmax/rdmin)*(float(ir-1)/float(nr-1))
123	end do
124	rdi(1)=rdimin
125	do ir=2,nr
126	rdi(ir)= sqrt(rd(ir-1)*rd(ir))
127	end do
128	rdi(nr+1)=rdimax
129
130	do ir=1,nr
131	do iint=1,ninter
132	rr(iint,ir)=
133	& rdi(ir)*
134	& (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
135	c write(,) rr(iint,ir)
136	end do
137	end do
138
139	! identify tracers corresponding to mass mixing ratio and
140	! number mixing ratio
141	idust_mass=0 ! dummy initialization
142	idust_number=0 ! dummy initialization
143
144	do iq=1,nq
145	if (noms(iq).eq."dust_mass") then
146	idust_mass=iq
147	endif
148	if (noms(iq).eq."dust_number") then
149	idust_number=iq
150	endif
151	enddo
152
153	! check that we did find the tracers
154	if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
155	write(,) 'callsedim: error! could not identify'
156	write(,) ' tracers for dust mass and number mixing'
157	write(,) ' ratio and doubleq is activated!'
158	stop
159	endif
160	ENDIF !of if (doubleq)
161
162	IF (water) THEN
163	write(,) "water_param nueff Sedimentation:", nuice_sed
164	IF (activice) THEN
165	write(,) "water_param nueff Radiative:", nuice_ref
166	ENDIF
167	ENDIF
168
169	firstcall=.false.
170	ENDIF ! of IF (firstcall)
171
172	c-----------------------------------------------------------------------
173	c 1. Initialization
174	c -----------------
175
176	zqi(1:ngrid,1:nlay,1:nqmx) = 0.
177	c Updating the mass mixing ratio with the tendencies coming
178	c from other parameterizations:
179	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
180
181	do iq=1,nq
182	do l=1,nlay
183	do ig=1,ngrid
184	zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
185	enddo
186	enddo
187	enddo
188
189	c Computing the different layer properties
190	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
191	c Mass (kg.m-2), thickness(m), crossing time (s) etc.
192
193	do l=1,nlay
194	do ig=1, ngrid
195	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
196	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
197	end do
198	end do
199
200	c =================================================================
201	do iq=1,nq
202	if(radius(iq).gt.1.e-9) then ! no sedim for gaz
203
204	c -----------------------------------------------------------------
205	c DOUBLEQ CASE
206	c -----------------------------------------------------------------
207	if (doubleq.and.
208	& ((iq.eq.idust_mass).or.
209	& (iq.eq.idust_number))) then
210
211	c Computing size distribution:
212	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
213
214	do l=1,nlay
215	do ig=1, ngrid
216	r0(ig,l)=
217	& CBRT(r3n_q*zqi(ig,l,idust_mass)/
218	& max(zqi(ig,l,idust_number),0.01))
219	r0(ig,l)=min(max(r0(ig,l),1.e-10),500.e-6)
220	end do
221	end do
222
223	c Computing mass mixing ratio for each particle size
224	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
225	IF (iq.EQ.idust_mass) then
226	radpower = 2
227	ELSE
228	radpower = -1
229	ENDIF
230	Sq(1:ngrid,1:nlay) = 0.
231	do ir=1,nr
232	do l=1,nlay
233	do ig=1,ngrid
234	c ****************
235	c Size distribution integration
236	c (Trapezoid Integration Method)
237	qr(ig,l,ir)=0.5(rr(2,ir)-rr(1,ir))
238	& (rr(1,ir)*radpower)
239	& exp(-(log(rr(1,ir)/r0(ig,l)))*2/(2varian**2))
240	do iint=2,ninter-1
241	qr(ig,l,ir)=qr(ig,l,ir) +
242	& 0.5(rr(iint+1,ir)-rr(iint-1,ir))
243	& (rr(iint,ir)*radpower)
244	& exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
245	& (2varian*2))
246	end do
247	qr(ig,l,ir)=qr(ig,l,ir) +
248	& 0.5(rr(ninter,ir)-rr(ninter-1,ir))
249	& (rr(ninter,ir)*radpower)
250	& exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
251	& (2varian*2))
252
253	c **************** old method (not recommended!)
254	c qr(ig,l,ir)=(rd(ir)*(5-3iq))*
255	c & exp( -(log(rd(ir)/r0(ig,l)))*2 / (2varian**2) )
256	c ******************************
257
258	Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
259	enddo
260	enddo
261	enddo
262
263	do ir=1,nr
264	do l=1,nlay
265	do ig=1,ngrid
266	qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
267	enddo
268	enddo
269	enddo
270
271	c Computing sedimentation for each tracer
272	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
273
274	c call zerophys(ngridmx*nlayermx,zqi(1,1,iq))
275	zqi(1:ngrid,1:nlay,iq) = 0.
276	c call zerophys(ngridmx,pdqs_sed(1,iq))
277	pdqs_sed(1:ngrid,iq) = 0.
278
279	do ir=1,nr
280	call newsedim(ngrid,nlay,1,ptimestep,
281	& pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),
282	& wq,0.5)
283
284	c Tendencies
285	c ~~~~~~~~~~
286	do ig=1,ngrid
287	pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
288	& + wq(ig,1)/ptimestep
289	end do
290	DO l = 1, nlay
291	DO ig=1,ngrid
292	zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
293	ENDDO
294	ENDDO
295	enddo ! of do ir=1,nr
296	c -----------------------------------------------------------------
297	c WATER CYCLE CASE
298	c -----------------------------------------------------------------
299	else if (water.and.(iq.eq.igcm_h2o_ice)) then
300	c if (doubleq) then
301	c do l=1,nlay
302	c do ig=1,ngrid
303	c rfall(ig,l)=max( rice(ig,l),rdust(ig,l) )
304	c rfall(ig,l)=min(rfall(ig,l),1.e-4)
305	c enddo
306	c enddo ! of do l=1,nlay
307	c else
308	do l=1,nlay
309	do ig=1,ngrid
310	c For water cycle, a typical dust size is assumed:
311	c r(z)=r0*exp(-z/H) with r0=0.8 micron and H=18 km.
312	c rfall(ig,l)=max( rice(ig,l)*1.5,rdust(ig,l) )
313	rfall(ig,l)=max( rice(ig,l)(1.+nuice_sed)*3.,
314	& rdust(ig,l) )
315	c modif FranckMM pour ameliorer cycle H2O: rfall= 20 microns
316	c mars commente pour l'instant rfall(ig,l)=20.e-6
317	rfall(ig,l)=min(rfall(ig,l),1.e-4)
318	enddo
319	enddo ! of do l=1,nlay
320	c endif
321	call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
322	& pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi(1,1,iq),
323	& wq,1.0)
324	c Tendencies
325	c ~~~~~~~~~~
326	do ig=1,ngrid
327	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
328	end do
329	c -----------------------------------------------------------------
330	c GENERAL CASE
331	c -----------------------------------------------------------------
332	else
333	call newsedim(ngrid,nlay,1,ptimestep,
334	& pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
335	& zqi(1,1,iq),wq,1.0)
336	c Tendencies
337	c ~~~~~~~~~~
338	do ig=1,ngrid
339	pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
340	end do
341	endif ! of if doubleq and if water
342	c -----------------------------------------------------------------
343
344	DO l = 1, nlay
345	DO ig=1,ngrid
346	pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
347	$ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
348	ENDDO
349	ENDDO
350
351	endif ! of if(radius(iq).gt.1.e-9)
352	c =================================================================
353	enddo ! of do iq=1,nq
354
355	RETURN
356	END
357

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

Download in other formats:

Original Format