Context Navigation

improvedclouds.F @ 414

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

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

Download in other formats:

Original Format