Context Navigation

improvedclouds.F @ 524

Last change on this file since 524 was 520, checked in by tnavarro, 13 years ago

10/02/12 == TN

Major update on watercycle: a smaller integration timestep is now used

in watercloud.F, sedimentation of clouds is done in watercloud instead of
callsedim.F

Temperature-dependant contact parameter in nuclea.F
No dust lifting if CO2 ice in vdif.c
Ice integrated column opacity is written in diagfi from physiq.F, instead

of aeropacity.F. Mandatory if iradia is not 1.

New definition of permanent ice in surfini.F and possibility to have an ice

cap in it in 1d.

Update in deftank: callphys.def.outliers,run.def.1d; added traceur.def.scavenging

File size: 24.6 KB

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

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

Download in other formats:

Original Format