Context Navigation

improvedclouds.F @ 426

Last change on this file since 426 was 420, checked in by tnavarro, 14 years ago

24/11/11 == TN

corrected minor bug in updatereffrad.F : reffdust was not saved

ccn_factor as not correctly used in sedimentation.

It is now initialized in inifis.F, declared in tracer.h and
used in both simpleclouds.F & callsedim.F to update ice radius.

Commented diagfi outputs in aeropacity.F & improvedclouds.F for non scavenging users.

File size: 23.2 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	INTEGER count
129
130	LOGICAL output_sca ! scavenging outputs flag for tests
131	output_sca = .false.
132	c----------------------------------
133	c----------------------------------
134
135	c------------------------------------------------------------------
136
137	IF (firstcall) THEN
138
139	c Definition of the size grid
140	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~
141	c rad_cld is the primary radius grid used for microphysics computation.
142	c The grid spacing is computed assuming a constant volume ratio
143	c between two consecutive bins; i.e. vrat_cld.
144	c vrat_cld is determined from the boundary values of the size grid:
145	c rmin_cld and rmax_cld.
146	c The rb_cld array contains the boundary values of each rad_cld bin.
147	c dr_cld is the width of each rad_cld bin.
148
149	c Volume ratio between two adjacent bins
150	vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
151	vrat_cld = dexp(vrat_cld)
152	write(,) "vrat_cld", vrat_cld
153
154	rb_cld(1) = rbmin_cld
155	rad_cld(1) = rmin_cld
156	vol_cld(1) = 4./3. * dble(pi) * rmin_cld**3.
157
158	do i=1,nbin_cld-1
159	rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
160	vol_cld(i+1) = vol_cld(i) * vrat_cld
161	enddo
162
163	do i=1,nbin_cld
164	rb_cld(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
165	& rad_cld(i)
166	dr_cld(i) = rb_cld(i+1) - rb_cld(i)
167	enddo
168	rb_cld(nbin_cld+1) = rbmax_cld
169	dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
170
171	print*, ' '
172	print*,'Microphysics: size bin information:'
173	print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
174	print*,'-----------------------------------'
175	do i=1,nbin_cld
176	write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
177	& dr_cld(i)
178	enddo
179	write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
180	print*,'-----------------------------------'
181
182	c Contact parameter of water ice on dust ( m=cos(theta) )
183	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
184	! mteta = 0.95
185	write(,) 'water_param contact parameter:', mteta
186
187	c Volume of a water molecule (m3)
188	vo1 = mh2o / dble(rho_ice)
189	c Variance of the ice and CCN distributions
190	sigma_ice = sqrt(log(1.+nuice_sed))
191
192	write(,) 'Variance of ice & CCN distribs :', sigma_ice
193	write(,) 'Volume of a water molecule:', vo1
194
195	firstcall=.false.
196	END IF
197
198	c-----------------------------------------------------------------------
199	c 1. Initialization
200	c-----------------------------------------------------------------------
201	c Initialize the tendencies
202	pdqscloud(1:ngrid,1:nq)=0
203	pdqcloud(1:ngrid,1:nlay,1:nq)=0
204	pdtcloud(1:ngrid,1:nlay)=0
205
206	c Update the needed variables
207	do l=1,nlay
208	do ig=1,ngrid
209	c Temperature
210	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
211	c Dust mass mixing ratio
212	zq(ig,l,igcm_dust_mass) =
213	& pq(ig,l,igcm_dust_mass) +
214	& pdq(ig,l,igcm_dust_mass) * ptimestep
215	zq0(ig,l,igcm_dust_mass)=zq(ig,l,igcm_dust_mass)
216	c Dust particle number
217	zq(ig,l,igcm_dust_number) =
218	& pq(ig,l,igcm_dust_number) +
219	& pdq(ig,l,igcm_dust_number) * ptimestep
220	zq0(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)
221	c Update rdust from last tendencies
222	rdust(ig,l)=
223	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
224	& max(zq(ig,l,igcm_dust_number),0.01))
225	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
226	c CCN mass mixing ratio
227	zq(ig,l,igcm_ccn_mass)=
228	& pq(ig,l,igcm_ccn_mass)+pdq(ig,l,igcm_ccn_mass)*ptimestep
229	zq(ig,l,igcm_ccn_mass)=max(zq(ig,l,igcm_ccn_mass),1.E-30)
230	zq0(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
231	c CCN particle number
232	zq(ig,l,igcm_ccn_number)=
233	& pq(ig,l,igcm_ccn_number)+pdq(ig,l,igcm_ccn_number)*ptimestep
234	zq(ig,l,igcm_ccn_number)=max(zq(ig,l,igcm_ccn_number),1.E-30)
235	zq0(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
236	c Water vapor
237	zq(ig,l,igcm_h2o_vap)=
238	& pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
239	zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004
240	zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
241	c Water ice
242	zq(ig,l,igcm_h2o_ice)=
243	& pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
244	zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004
245	zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
246	enddo
247	enddo
248
249
250	c------------------------------------------------------------------
251	c Cloud microphysical scheme
252	c------------------------------------------------------------------
253
254	Cste = ptimestep * 4. * pi * rho_ice
255
256	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
257
258	count = 0
259
260	c Main loop over the GCM's grid
261	DO l=1,nlay
262	DO ig=1,ngrid
263
264	c Get the partial pressure of water vapor and its saturation ratio
265	ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l)
266	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
267
268	IF ((satu .ge. 1)! ) THEN ! if we have condensation
269	& .or. ( zq(ig,l,igcm_ccn_number)
270	& .ge. 10) ) THEN ! or sublimation
271
272
273	c Expand the dust moments into a binned distribution
274	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)
275	No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+ 1.e-30
276	Rn = rdust(ig,l)
277	Rm = Rn * exp( 3. * sigma_ice**2. )
278	Rn = 1. / Rn
279	Rm = 1. / Rm
280	dev2 = 1. / ( sqrt(2.) * sigma_ice )
281	do i = 1, nbin_cld
282	n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)Rn) dev2 )
283	& -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
284	m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)Rm) dev2 )
285	& -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
286	enddo
287
288	! sumcheck = 0
289	! do i = 1, nbin_cld
290	! sumcheck = sumcheck + n_aer(i)
291	! enddo
292	! sumcheck = abs(sumcheck/No - 1)
293	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
294	! print*, "WARNING, No sumcheck PROBLEM"
295	! print*, "sumcheck, No",sumcheck, No
296	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
297	! print*, "Dust binned distribution", n_aer
298	! endif
299	!
300	! sumcheck = 0
301	! do i = 1, nbin_cld
302	! sumcheck = sumcheck + m_aer(i)
303	! enddo
304	! sumcheck = abs(sumcheck/Mo - 1)
305	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
306	! print*, "WARNING, Mo sumcheck PROBLEM"
307	! print*, "sumcheck, Mo",sumcheck, Mo
308	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
309	! print*, "Dust binned distribution", m_aer
310	! endif
311
312	c Get the rates of nucleation
313	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
314
315
316	dN = 0.
317	dM = 0.
318	do i = 1, nbin_cld
319	n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
320	m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
321	dN = dN + n_aer(i) * rate(i) * ptimestep
322	dM = dM + m_aer(i) * rate(i) * ptimestep
323	enddo
324
325	! dN = min( max(dN,-zq(ig,l,igcm_ccn_number) ),
326	! & zq(ig,l,igcm_dust_number) )
327	!
328	! dM = min( max(dM,-zq(ig,l,igcm_ccn_mass) ),
329	! & zq(ig,l,igcm_dust_mass) )
330
331	Mo = zq0(ig,l,igcm_h2o_ice) +
332	& zq0(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
333	No = zq0(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
334	!write(,) "l,cloud particles,cloud mass",l, No, Mo
335	rhocloud(ig,l) = zq0(ig,l,igcm_h2o_ice) / Mo * rho_ice
336	& +zq0(ig,l,igcm_ccn_mass)* tauscaling(ig)
337	& / Mo * rho_dust
338	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
339	rice(ig,l) =
340	& ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
341	c nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
342	if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
343	seq = exp( 2.sig(zt(ig,l))mh2o /
344	& (rho_icergpzt(ig,l)*rice(ig,l)) )
345
346	call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
347	& ph2o,ph2o/satu,seq,rice(ig,l),gr)
348
349	up = Cste * gr * rice(ig,l) * No * seq +
350	& zq(ig,l,igcm_h2o_vap)
351	dwn = Cste * gr * rice(ig,l) * No / zqsat(ig,l)+ 1.
352
353	Ctot = zq0(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_h2o_vap)
354
355	newvap = min(up/dwn,Ctot)
356
357	gr = gr * ( newvap/zqsat(ig,l) - seq )
358
359
360	dMice = min( max(Cste * No * rice(ig,l) * gr,
361	& -zq(ig,l,igcm_h2o_ice) ),
362	& zq(ig,l,igcm_h2o_vap) )
363
364	c----------- TESTS 1D output ---------
365	satu_out(ig,l) = satu
366	Mcon_out(ig,l) = dMice
367	newvap_out(ig,l) = newvap
368	gr_out(ig,l) = gr
369	dN_out(ig,l) = dN
370	dM_out(ig,l) = dM
371	c-------------------------------------
372
373	c Water ice
374	zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice) +
375	& dMice
376
377	c Water vapor
378	zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) -
379	& dMice
380
381	c If all the ice particles sublimate, all the condensation
382	c nuclei are released:
383	if (zq(ig,l,igcm_h2o_ice).le.1e-30) then
384	c for coherence
385	dM = 0
386	dN = 0
387	dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig)
388	dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
389	c Water ice particles
390	zq(ig,l,igcm_h2o_ice) = 0.
391	c Dust particles
392	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) +
393	& zq(ig,l,igcm_ccn_mass)
394	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) +
395	& zq(ig,l,igcm_ccn_number)
396	c CCNs
397	zq(ig,l,igcm_ccn_mass) = 0.
398	zq(ig,l,igcm_ccn_number) = 0.
399	endif
400
401	dN = dN/ tauscaling(ig)
402	dM = dM/ tauscaling(ig)
403	c Dust particles
404	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM
405	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN
406	c CCNs
407	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
408	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
409
410	pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass)
411	& -zq0(ig,l,igcm_dust_mass))/ptimestep
412	pdqcloud(ig,l,igcm_dust_number)=(zq(ig,l,igcm_dust_number)
413	& -zq0(ig,l,igcm_dust_number))/ptimestep
414	pdqcloud(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
415	& -zq0(ig,l,igcm_ccn_mass))/ptimestep
416	pdqcloud(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
417	& -zq0(ig,l,igcm_ccn_number))/ptimestep
418	pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
419	& -zq0(ig,l,igcm_h2o_vap))/ptimestep
420	pdqcloud(ig,l,igcm_h2o_ice)=(zq(ig,l,igcm_h2o_ice)
421	& -zq0(ig,l,igcm_h2o_ice))/ptimestep
422
423	count = count +1
424	ELSE ! for coherence (rdust, rice computations etc ..)
425	zq(ig,l,igcm_dust_mass) = zq0(ig,l,igcm_dust_mass)
426	zq(ig,l,igcm_dust_number) = zq0(ig,l,igcm_dust_number)
427	zq(ig,l,igcm_ccn_mass) = zq0(ig,l,igcm_ccn_mass)
428	zq(ig,l,igcm_ccn_number) = zq0(ig,l,igcm_ccn_number)
429	zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice)
430	zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap)
431
432	! pour les sorties de test
433	satu_out(ig,l) = satu
434	Mcon_out(ig,l) = 0
435	newvap_out(ig,l) = zq(ig,l,igcm_h2o_vap)
436	gr_out(ig,l) = 0
437	dN_out(ig,l) = 0
438	dM_out(ig,l) = 0
439
440	ENDIF ! end if (saturation ratio > 1) or (there is h2o_ice)
441
442	c-----update temperature
443	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
444	pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
445
446	c----- update rice & rhocloud AFTER scavenging
447	Mo = zq(ig,l,igcm_h2o_ice) +
448	& zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
449	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
450	rhocloud(ig,l) = zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
451	& +zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
452	& / Mo * rho_dust
453	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
454	rice(ig,l) =
455	& ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
456	if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
457
458	nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
459
460	c----- update rdust and sedimentation radius
461	rdust(ig,l)=
462	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
463	& max(zq(ig,l,igcm_dust_number),0.01))
464	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
465
466	rsedcloud(ig,l)=max( rice(ig,l)(1.+nuice_sed)*3.,
467	& rdust(ig,l) )
468	rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
469
470	ENDDO
471	ENDDO
472
473
474	c------------------------------------------------------------------
475	c------------------------------------------------------------------
476	c------------------------------------------------------------------
477	c------------------------------------------------------------------
478	c------------------------------------------------------------------
479	c TESTS
480
481	IF (output_sca) then
482
483	print*, 'count is ',count, ' i.e. ',
484	& count100/(nlayngrid), '% for microphys computation'
485
486	dM_col(:) = 0
487	dN_col(:) = 0
488	Mdust_col(:) = 0
489	Ndust_col(:) = 0
490	Mccn_col(:) = 0
491	Nccn_col(:) = 0
492	do l=1, nlay
493	do ig=1,ngrid
494	dM_col(ig) = dM_col(ig) +
495	& dM_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g
496	dN_col(ig) = dN_col(ig) +
497	& dN_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g
498	Mdust_col(ig) = Mdust_col(ig) +
499	& zq(ig,l,igcm_dust_mass)*tauscaling(ig)
500	& *(pplev(ig,l) - pplev(ig,l+1)) / g
501	Ndust_col(ig) = Ndust_col(ig) +
502	& zq(ig,l,igcm_dust_number)*tauscaling(ig)
503	& *(pplev(ig,l) - pplev(ig,l+1)) / g
504	Mccn_col(ig) = Mccn_col(ig) +
505	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
506	& *(pplev(ig,l) - pplev(ig,l+1)) / g
507	Nccn_col(ig) = Nccn_col(ig) +
508	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
509	& *(pplev(ig,l) - pplev(ig,l+1)) / g
510	enddo ! of do ig=1,ngrid
511	enddo ! of do l=1,nlay
512
513
514	IF (ngrid.ne.1) THEN ! 3D
515	call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
516	& satu_out)
517	call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",2,
518	& dM_col)
519	call WRITEDIAGFI(ngrid,"dN_col","dN column","N",2,
520	& dN_col)
521	call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",2,
522	& Ndust_col)
523	call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",2,
524	& Mdust_col)
525	call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",2,
526	& Nccn_col)
527	call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",2,
528	& Mccn_col)
529	call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
530	& dM_out)
531	call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
532	& dN_out)
533	ENDIF
534
535	IF (ngrid.eq.1) THEN ! 1D
536
537	call WRITEDIAGFI(ngrid,"newvap","h2o newvap","kg",1,
538	& newvap_out)
539	call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
540	& gr_out)
541	call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
542	& dM_out)
543	call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
544	& dN_out)
545	call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
546	& Mcon_out)
547	call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
548	& zqsat)
549	call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
550	& satu_out)
551	call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1,
552	& rice)
553	call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
554	& rdust)
555	call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
556	& rsedcloud)
557	call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
558	& rhocloud)
559	call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",0,
560	& dM_col)
561	call WRITEDIAGFI(ngrid,"dN_col","dN column","N",0,
562	& dN_col)
563	call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",0,
564	& Ndust_col)
565	call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",0,
566	& Mdust_col)
567	call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",0,
568	& Nccn_col)
569	call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",0,
570	& Mccn_col)
571	ENDIF
572	ENDIF ! endif output_sca
573	c------------------------------------------------------------------
574	return
575	end

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

Download in other formats:

Original Format