Context Navigation

improvedclouds.F @ 532

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

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

Download in other formats:

Original Format