Context Navigation

callsedim.F @ 1009

Last change on this file since 1009 was 1005, checked in by emillour, 11 years ago

Mars GCM:

Bug fix: when running with photochemistry, ccns did not sediment! Fixed initracer.F. Also added that callsedim/newsedim use updated temperatures.
Converted run0 and run_mcd scripts to bash.

File size: 15.1 KB

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

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

Download in other formats:

Original Format