Context Navigation

improvedclouds.F @ 517

Last change on this file since 517 was 455, checked in by tnavarro, 14 years ago
tiny change : nuice_sed initialisation is now done in inifis. Also changed initracer and improvedclouds.
File size: 23.5 KB

Line
1	subroutine improvedclouds(ngrid,nlay,ptimestep,
2	& pplev,pplay,pt,pdt,
3	& pq,pdq,pdqcloud,pdqscloud,pdtcloud,
4	& nq,tauscaling,rdust,rice,nuice,
5	& rsedcloud,rhocloud)
6	implicit none
7	c------------------------------------------------------------------
8	c This routine is used to form clouds when a parcel of the GCM is
9	c saturated. It includes the ability to have supersaturation, a
10	c computation of the nucleation rates, growthrates and the
11	c scavenging of dust particles by clouds.
12	c It is worth noting that the amount of dust is computed using the
13	c dust optical depth computed in aeropacity.F. That's why
14	c the variable called "tauscaling" is used to convert
15	c pq(dust_mass) and pq(dust_number), which are relative
16	c quantities, to absolute and realistic quantities stored in zq.
17	c This has to be done to convert the inputs into absolute
18	c values, but also to convert the outputs back into relative
19	c values which are then used by the sedimentation and advection
20	c schemes.
21
22	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
23	c (October 2011)
24	c------------------------------------------------------------------
25	#include "dimensions.h"
26	#include "dimphys.h"
27	#include "comcstfi.h"
28	#include "callkeys.h"
29	#include "tracer.h"
30	#include "comgeomfi.h"
31	#include "dimradmars.h"
32	#include "microphys.h"
33	c------------------------------------------------------------------
34	c Inputs:
35
36	INTEGER ngrid,nlay
37	integer nq ! nombre de traceurs
38	REAL ptimestep ! pas de temps physique (s)
39	REAL pplev(ngrid,nlay+1),pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
40	REAL pt(ngrid,nlay) ! temperature at the middle of the
41	! layers (K)
42	REAL pdt(ngrid,nlay) ! tendance temperature des autres
43	! param.
44	REAL pq(ngrid,nlay,nq) ! traceur (kg/kg)
45	REAL pdq(ngrid,nlay,nq) ! tendance avant condensation
46	! (kg/kg.s-1)
47	REAL tauscaling(ngridmx) ! Convertion factor for qdust and Ndust
48	REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
49
50	c Outputs:
51	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
52	! (r_c in montmessin_2004)
53	REAL nuice(ngrid,nlay) ! Estimated effective variance
54	! of the size distribution
55	REAL rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
56	REAL rhocloud(ngridmx,nlayermx) ! Cloud density (kg.m-3)
57	REAL pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
58	! H2O(kg/kg.s-1)
59	REAL pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1)
60	REAL pdtcloud(ngrid,nlay) ! tendance temperature due
61	! a la chaleur latente
62
63	c------------------------------------------------------------------
64	c Local variables:
65
66	LOGICAL firstcall
67	DATA firstcall/.true./
68	SAVE firstcall
69
70	REAL*8 derf ! Error function
71	!external derf
72
73	REAL CBRT
74	EXTERNAL CBRT
75
76	INTEGER ig,l,i
77
78	REAL zq(ngridmx,nlayermx,nqmx) ! local value of tracers
79	REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
80	REAL zt(ngridmx,nlayermx) ! local value of temperature
81	REAL zqsat(ngridmx,nlayermx) ! saturation
82	REAL lw !Latent heat of sublimation (J.kg-1)
83	REAL Cste
84	REAL dMice ! mass of condensated ice
85	REAL sumcheck
86	REAL*8 ph2o ! Water vapor partial pressure (Pa)
87	REAL*8 satu ! Water vapor saturation ratio over ice
88	REAL*8 Mo,No
89	REAL*8 dN,dM,newvap
90	REAL*8 Rn, Rm, dev2
91	REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
92	REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
93	REAL*8 rate(nbin_cld) ! nucleation rate
94	REAL*8 up,dwn,Ctot,gr,seq
95	REAL*8 sig ! Water-ice/air surface tension (N.m)
96	EXTERNAL sig
97
98	c Parameters of the size discretization
99	c used by the microphysical scheme
100	DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
101	DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
102	DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
103	! Minimum boundary radius (m)
104	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
105	DOUBLE PRECISION vrat_cld ! Volume ratio
106	DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
107	SAVE rb_cld
108	DOUBLE PRECISION dr_cld(nbin_cld)! width of each rad_cld bin (m)
109	DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3)
110
111	REAL sigma_ice ! Variance of the ice and CCN distributions
112	SAVE sigma_ice
113
114	c----------------------------------
115	c some outputs for 1D -- TEMPORARY
116	REAL satu_out(ngridmx,nlayermx) ! satu ratio for output
117	REAL dN_out(ngridmx,nlayermx) ! mass variation for output
118	REAL dM_out(ngridmx,nlayermx) ! number variation for output
119	REAL dM_col(ngridmx) ! total mass condensed in column
120	REAL dN_col(ngridmx) ! total mass condensed in column
121	REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
122	REAL gr_out(ngridmx,nlayermx) ! for 1d output
123	REAL newvap_out(ngridmx,nlayermx) ! for 1d output
124	REAL Mdust_col(ngridmx) ! total column dust mass
125	REAL Ndust_col(ngridmx) ! total column dust mass
126	REAL Mccn_col(ngridmx) ! total column ccn mass
127	REAL Nccn_col(ngridmx) ! total column ccn mass
128	REAL rate_out(ngridmx,nlayermx) ! nucleation rate
129	INTEGER count
130
131	LOGICAL output_sca ! scavenging outputs flag for tests
132	output_sca = .false.
133	c----------------------------------
134	c----------------------------------
135
136	c------------------------------------------------------------------
137
138	IF (firstcall) THEN
139
140	c Definition of the size grid
141	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~
142	c rad_cld is the primary radius grid used for microphysics computation.
143	c The grid spacing is computed assuming a constant volume ratio
144	c between two consecutive bins; i.e. vrat_cld.
145	c vrat_cld is determined from the boundary values of the size grid:
146	c rmin_cld and rmax_cld.
147	c The rb_cld array contains the boundary values of each rad_cld bin.
148	c dr_cld is the width of each rad_cld bin.
149
150	c Volume ratio between two adjacent bins
151	vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
152	vrat_cld = dexp(vrat_cld)
153	write(,) "vrat_cld", vrat_cld
154
155	rb_cld(1) = rbmin_cld
156	rad_cld(1) = rmin_cld
157	vol_cld(1) = 4./3. * dble(pi) * rmin_cld**3.
158
159	do i=1,nbin_cld-1
160	rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
161	vol_cld(i+1) = vol_cld(i) * vrat_cld
162	enddo
163
164	do i=1,nbin_cld
165	rb_cld(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
166	& rad_cld(i)
167	dr_cld(i) = rb_cld(i+1) - rb_cld(i)
168	enddo
169	rb_cld(nbin_cld+1) = rbmax_cld
170	dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
171
172	print*, ' '
173	print*,'Microphysics: size bin information:'
174	print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
175	print*,'-----------------------------------'
176	do i=1,nbin_cld
177	write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
178	& dr_cld(i)
179	enddo
180	write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
181	print*,'-----------------------------------'
182
183	c Contact parameter of water ice on dust ( m=cos(theta) )
184	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
185	! mteta = 0.95
186	write(,) 'water_param contact parameter:', mteta
187
188	c Volume of a water molecule (m3)
189	vo1 = mh2o / dble(rho_ice)
190	c Variance of the ice and CCN distributions
191	sigma_ice = sqrt(log(1.+nuice_sed))
192
193	write(,) 'Variance of ice & CCN distribs :', sigma_ice
194	write(,) 'nuice for sedimentation:', nuice_sed
195	write(,) 'Volume of a water molecule:', vo1
196
197	firstcall=.false.
198	END IF
199
200	c-----------------------------------------------------------------------
201	c 1. Initialization
202	c-----------------------------------------------------------------------
203	c Initialize the tendencies
204	pdqscloud(1:ngrid,1:nq)=0
205	pdqcloud(1:ngrid,1:nlay,1:nq)=0
206	pdtcloud(1:ngrid,1:nlay)=0
207
208	c Update the needed variables
209	do l=1,nlay
210	do ig=1,ngrid
211	c Temperature
212	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
213	c Dust mass mixing ratio
214	zq(ig,l,igcm_dust_mass) =
215	& pq(ig,l,igcm_dust_mass) +
216	& pdq(ig,l,igcm_dust_mass) * ptimestep
217	zq0(ig,l,igcm_dust_mass)=zq(ig,l,igcm_dust_mass)
218	c Dust particle number
219	zq(ig,l,igcm_dust_number) =
220	& pq(ig,l,igcm_dust_number) +
221	& pdq(ig,l,igcm_dust_number) * ptimestep
222	zq0(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)
223	c Update rdust from last tendencies
224	rdust(ig,l)=
225	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
226	& max(zq(ig,l,igcm_dust_number),0.01))
227	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
228	c CCN mass mixing ratio
229	zq(ig,l,igcm_ccn_mass)=
230	& pq(ig,l,igcm_ccn_mass)+pdq(ig,l,igcm_ccn_mass)*ptimestep
231	zq(ig,l,igcm_ccn_mass)=max(zq(ig,l,igcm_ccn_mass),1.E-30)
232	zq0(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
233	c CCN particle number
234	zq(ig,l,igcm_ccn_number)=
235	& pq(ig,l,igcm_ccn_number)+pdq(ig,l,igcm_ccn_number)*ptimestep
236	zq(ig,l,igcm_ccn_number)=max(zq(ig,l,igcm_ccn_number),1.E-30)
237	zq0(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
238	c Water vapor
239	zq(ig,l,igcm_h2o_vap)=
240	& pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
241	zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004
242	zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
243	c Water ice
244	zq(ig,l,igcm_h2o_ice)=
245	& pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
246	zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004
247	zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
248	enddo
249	enddo
250
251
252	c------------------------------------------------------------------
253	c Cloud microphysical scheme
254	c------------------------------------------------------------------
255
256	Cste = ptimestep * 4. * pi * rho_ice
257
258	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
259
260	count = 0
261
262	c Main loop over the GCM's grid
263	DO l=1,nlay
264	DO ig=1,ngrid
265
266	c Get the partial pressure of water vapor and its saturation ratio
267	ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l)
268	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
269
270	IF ((satu .ge. 1)! ) THEN ! if we have condensation
271	& .or. ( zq(ig,l,igcm_ccn_number)
272	& .ge. 10) ) THEN ! or sublimation
273
274
275	c Expand the dust moments into a binned distribution
276	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)
277	No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+ 1.e-30
278	Rn = rdust(ig,l)
279	Rm = Rn * exp( 3. * sigma_ice**2. )
280	Rn = 1. / Rn
281	Rm = 1. / Rm
282	dev2 = 1. / ( sqrt(2.) * sigma_ice )
283	do i = 1, nbin_cld
284	n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)Rn) dev2 )
285	& -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
286	m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)Rm) dev2 )
287	& -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
288	enddo
289
290	! sumcheck = 0
291	! do i = 1, nbin_cld
292	! sumcheck = sumcheck + n_aer(i)
293	! enddo
294	! sumcheck = abs(sumcheck/No - 1)
295	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
296	! print*, "WARNING, No sumcheck PROBLEM"
297	! print*, "sumcheck, No",sumcheck, No
298	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
299	! print*, "Dust binned distribution", n_aer
300	! endif
301	!
302	! sumcheck = 0
303	! do i = 1, nbin_cld
304	! sumcheck = sumcheck + m_aer(i)
305	! enddo
306	! sumcheck = abs(sumcheck/Mo - 1)
307	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
308	! print*, "WARNING, Mo sumcheck PROBLEM"
309	! print*, "sumcheck, Mo",sumcheck, Mo
310	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
311	! print*, "Dust binned distribution", m_aer
312	! endif
313
314	c Get the rates of nucleation
315	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
316
317
318	dN = 0.
319	dM = 0.
320	rate_out(ig,l) = 0
321	do i = 1, nbin_cld
322	n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
323	m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
324	dN = dN + n_aer(i) * rate(i) * ptimestep
325	dM = dM + m_aer(i) * rate(i) * ptimestep
326	rate_out(ig,l)=rate_out(ig,l)+rate(i)
327	enddo
328
329	! dN = min( max(dN,-zq(ig,l,igcm_ccn_number) ),
330	! & zq(ig,l,igcm_dust_number) )
331	!
332	! dM = min( max(dM,-zq(ig,l,igcm_ccn_mass) ),
333	! & zq(ig,l,igcm_dust_mass) )
334
335	Mo = zq0(ig,l,igcm_h2o_ice) +
336	& zq0(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
337	No = zq0(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
338	!write(,) "l,cloud particles,cloud mass",l, No, Mo
339	rhocloud(ig,l) = zq0(ig,l,igcm_h2o_ice) / Mo * rho_ice
340	& +zq0(ig,l,igcm_ccn_mass)* tauscaling(ig)
341	& / Mo * rho_dust
342	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
343	rice(ig,l) =
344	& ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
345	c nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
346	if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
347	seq = exp( 2.sig(zt(ig,l))mh2o /
348	& (rho_icergpzt(ig,l)*rice(ig,l)) )
349
350	call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
351	& ph2o,ph2o/satu,seq,rice(ig,l),gr)
352
353	up = Cste * gr * rice(ig,l) * No * seq +
354	& zq(ig,l,igcm_h2o_vap)
355	dwn = Cste * gr * rice(ig,l) * No / zqsat(ig,l)+ 1.
356
357	Ctot = zq0(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_h2o_vap)
358
359	newvap = min(up/dwn,Ctot)
360
361	gr = gr * ( newvap/zqsat(ig,l) - seq )
362
363
364	dMice = min( max(Cste * No * rice(ig,l) * gr,
365	& -zq(ig,l,igcm_h2o_ice) ),
366	& zq(ig,l,igcm_h2o_vap) )
367
368	c----------- TESTS 1D output ---------
369	satu_out(ig,l) = satu
370	Mcon_out(ig,l) = dMice
371	newvap_out(ig,l) = newvap
372	gr_out(ig,l) = gr
373	dN_out(ig,l) = dN
374	dM_out(ig,l) = dM
375	c-------------------------------------
376
377	c Water ice
378	zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice) +
379	& dMice
380
381	c Water vapor
382	zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) -
383	& dMice
384
385	c If all the ice particles sublimate, all the condensation
386	c nuclei are released:
387	if (zq(ig,l,igcm_h2o_ice).le.1e-30) then
388	c for coherence
389	dM = 0
390	dN = 0
391	dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig)
392	dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
393	c Water ice particles
394	zq(ig,l,igcm_h2o_ice) = 0.
395	c Dust particles
396	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) +
397	& zq(ig,l,igcm_ccn_mass)
398	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) +
399	& zq(ig,l,igcm_ccn_number)
400	c CCNs
401	zq(ig,l,igcm_ccn_mass) = 0.
402	zq(ig,l,igcm_ccn_number) = 0.
403	endif
404
405	dN = dN/ tauscaling(ig)
406	dM = dM/ tauscaling(ig)
407	c Dust particles
408	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM
409	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN
410	c CCNs
411	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
412	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
413
414	pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass)
415	& -zq0(ig,l,igcm_dust_mass))/ptimestep
416	pdqcloud(ig,l,igcm_dust_number)=(zq(ig,l,igcm_dust_number)
417	& -zq0(ig,l,igcm_dust_number))/ptimestep
418	pdqcloud(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
419	& -zq0(ig,l,igcm_ccn_mass))/ptimestep
420	pdqcloud(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
421	& -zq0(ig,l,igcm_ccn_number))/ptimestep
422	pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
423	& -zq0(ig,l,igcm_h2o_vap))/ptimestep
424	pdqcloud(ig,l,igcm_h2o_ice)=(zq(ig,l,igcm_h2o_ice)
425	& -zq0(ig,l,igcm_h2o_ice))/ptimestep
426
427	count = count +1
428	ELSE ! for coherence (rdust, rice computations etc ..)
429	zq(ig,l,igcm_dust_mass) = zq0(ig,l,igcm_dust_mass)
430	zq(ig,l,igcm_dust_number) = zq0(ig,l,igcm_dust_number)
431	zq(ig,l,igcm_ccn_mass) = zq0(ig,l,igcm_ccn_mass)
432	zq(ig,l,igcm_ccn_number) = zq0(ig,l,igcm_ccn_number)
433	zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice)
434	zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap)
435
436	! pour les sorties de test
437	satu_out(ig,l) = satu
438	Mcon_out(ig,l) = 0
439	newvap_out(ig,l) = zq(ig,l,igcm_h2o_vap)
440	gr_out(ig,l) = 0
441	dN_out(ig,l) = 0
442	dM_out(ig,l) = 0
443
444	ENDIF ! end if (saturation ratio > 1) or (there is h2o_ice)
445
446	c-----update temperature
447	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
448	pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
449
450	c----- update rice & rhocloud AFTER scavenging
451	Mo = zq(ig,l,igcm_h2o_ice) +
452	& zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
453	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
454	rhocloud(ig,l) = zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
455	& +zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
456	& / Mo * rho_dust
457	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
458	rice(ig,l) =
459	& ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
460	if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
461
462	nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
463
464	c----- update rdust and sedimentation radius
465	rdust(ig,l)=
466	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
467	& max(zq(ig,l,igcm_dust_number),0.01))
468	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
469
470	rsedcloud(ig,l)=max( rice(ig,l)(1.+nuice_sed)*3.,
471	& rdust(ig,l) )
472	rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
473
474	ENDDO
475	ENDDO
476
477
478	c------------------------------------------------------------------
479	c------------------------------------------------------------------
480	c------------------------------------------------------------------
481	c------------------------------------------------------------------
482	c------------------------------------------------------------------
483	c TESTS
484
485	IF (output_sca) then
486
487	print*, 'count is ',count, ' i.e. ',
488	& count100/(nlayngrid), '% for microphys computation'
489
490	dM_col(:) = 0
491	dN_col(:) = 0
492	Mdust_col(:) = 0
493	Ndust_col(:) = 0
494	Mccn_col(:) = 0
495	Nccn_col(:) = 0
496	do l=1, nlay
497	do ig=1,ngrid
498	dM_col(ig) = dM_col(ig) +
499	& dM_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g
500	dN_col(ig) = dN_col(ig) +
501	& dN_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g
502	Mdust_col(ig) = Mdust_col(ig) +
503	& zq(ig,l,igcm_dust_mass)*tauscaling(ig)
504	& *(pplev(ig,l) - pplev(ig,l+1)) / g
505	Ndust_col(ig) = Ndust_col(ig) +
506	& zq(ig,l,igcm_dust_number)*tauscaling(ig)
507	& *(pplev(ig,l) - pplev(ig,l+1)) / g
508	Mccn_col(ig) = Mccn_col(ig) +
509	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
510	& *(pplev(ig,l) - pplev(ig,l+1)) / g
511	Nccn_col(ig) = Nccn_col(ig) +
512	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
513	& *(pplev(ig,l) - pplev(ig,l+1)) / g
514	enddo ! of do ig=1,ngrid
515	enddo ! of do l=1,nlay
516
517
518	IF (ngrid.ne.1) THEN ! 3D
519	call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
520	& satu_out)
521	call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",2,
522	& dM_col)
523	call WRITEDIAGFI(ngrid,"dN_col","dN column","N",2,
524	& dN_col)
525	call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",2,
526	& Ndust_col)
527	call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",2,
528	& Mdust_col)
529	call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",2,
530	& Nccn_col)
531	call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",2,
532	& Mccn_col)
533	call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
534	& dM_out)
535	call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
536	& dN_out)
537	ENDIF
538
539	IF (ngrid.eq.1) THEN ! 1D
540
541	call WRITEDIAGFI(ngrid,"newvap","h2o newvap","kg",1,
542	& newvap_out)
543	call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
544	& gr_out)
545	call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
546	& rate_out)
547	call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
548	& dM_out)
549	call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
550	& dN_out)
551	call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
552	& Mcon_out)
553	call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
554	& zqsat)
555	call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
556	& satu_out)
557	call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1,
558	& rice)
559	call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
560	& rdust)
561	call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
562	& rsedcloud)
563	call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
564	& rhocloud)
565	call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",0,
566	& dM_col)
567	call WRITEDIAGFI(ngrid,"dN_col","dN column","N",0,
568	& dN_col)
569	call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",0,
570	& Ndust_col)
571	call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",0,
572	& Mdust_col)
573	call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",0,
574	& Nccn_col)
575	call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",0,
576	& Mccn_col)
577	ENDIF
578	ENDIF ! endif output_sca
579	c------------------------------------------------------------------
580	return
581	end

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

Download in other formats:

Original Format