Context Navigation

callsedim_mod.F @ 2156

Last change on this file since 2156 was 1983, checked in by mvals, 6 years ago

Mars GCM:
Cosmetic/practical changes:

swmain and lwmain become modules swmain_mod, lwmain_mod
Addition of the intent in/out characteristics of variables in swmain_mod and lwmain_mod subroutines

Correction:

in callsedim_mod, declaration of variable tau(ngrid,nlay) corrected to tau(ngrid,naerkind)

File size: 19.3 KB

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

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

Download in other formats:

Original Format