Context Navigation

callsedim.F @ 1292

Last change on this file since 1292 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

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

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

Download in other formats:

Original Format