Context Navigation

callsedim.F @ 828

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

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

Download in other formats:

Original Format