Context Navigation

improvedclouds.F @ 623

Last change on this file since 623 was 541, checked in by aslmd, 13 years ago
LMDZ.MARS. Added small optimizations in watercycle. 10% improvement.
File size: 25.7 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	do i=1,nbin_cld+1
196	rb_cld(i) = dlog(rb_cld(i)) !! we save that so that it is not computed
197	!! at each timestep and gridpoint
198	enddo
199
200	c Contact parameter of water ice on dust ( m=cos(theta) )
201	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
202	! mteta = 0.95
203	write(,) 'water_param contact parameter:', mteta
204
205	c Volume of a water molecule (m3)
206	vo1 = mh2o / dble(rho_ice)
207	c Variance of the ice and CCN distributions
208	sigma_ice = sqrt(log(1.+nuice_sed))
209
210	write(,) 'Variance of ice & CCN distribs :', sigma_ice
211	write(,) 'nuice for sedimentation:', nuice_sed
212	write(,) 'Volume of a water molecule:', vo1
213
214	firstcall=.false.
215	END IF
216
217	c-----------------------------------------------------------------------
218	c 1. Initialization
219	c-----------------------------------------------------------------------
220	c Initialize the tendencies
221	pdqcloud(1:ngrid,1:nlay,1:nq)=0
222	pdtcloud(1:ngrid,1:nlay)=0
223
224	c Update the needed variables
225	do l=1,nlay
226	do ig=1,ngrid
227	c Temperature
228	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
229	c Dust mass mixing ratio
230	zq(ig,l,igcm_dust_mass) =
231	& pq(ig,l,igcm_dust_mass) +
232	& pdq(ig,l,igcm_dust_mass) * ptimestep
233	zq0(ig,l,igcm_dust_mass)=zq(ig,l,igcm_dust_mass)
234	c Dust particle number
235	zq(ig,l,igcm_dust_number) =
236	& pq(ig,l,igcm_dust_number) +
237	& pdq(ig,l,igcm_dust_number) * ptimestep
238	zq0(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)
239	c Update rdust from last tendencies
240	rdust(ig,l)=
241	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
242	& max(zq(ig,l,igcm_dust_number),0.01))
243	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
244	c CCN mass mixing ratio
245	zq(ig,l,igcm_ccn_mass)=
246	& pq(ig,l,igcm_ccn_mass)+pdq(ig,l,igcm_ccn_mass)*ptimestep
247	zq(ig,l,igcm_ccn_mass)=max(zq(ig,l,igcm_ccn_mass),1.E-30)
248	zq0(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
249	c CCN particle number
250	zq(ig,l,igcm_ccn_number)=
251	& pq(ig,l,igcm_ccn_number)+pdq(ig,l,igcm_ccn_number)*ptimestep
252	zq(ig,l,igcm_ccn_number)=max(zq(ig,l,igcm_ccn_number),1.E-30)
253	zq0(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
254	c Water vapor
255	zq(ig,l,igcm_h2o_vap)=
256	& pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
257	zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004
258	zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
259	c Water ice
260	zq(ig,l,igcm_h2o_ice)=
261	& pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
262	zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004
263	zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
264	enddo
265	enddo
266
267	c------------------------------------------------------------------
268	c Cloud microphysical scheme
269	c------------------------------------------------------------------
270
271	Cste = ptimestep * 4. * pi * rho_ice
272	dev2 = 1. / ( sqrt(2.) * sigma_ice )
273
274	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
275
276	count = 0
277
278	c Main loop over the GCM's grid
279	DO l=1,nlay
280	DO ig=1,ngrid
281
282	c Get the partial pressure of water vapor and its saturation ratio
283	ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l)
284	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
285
286	IF ((satu .ge. 1)! ) THEN ! if we have condensation
287	& .or. ( zq(ig,l,igcm_ccn_number)*tauscaling(ig)
288	& .ge. 2) ) THEN ! or sublimation
289
290
291	c Expand the dust moments into a binned distribution
292	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)
293	No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+ 1.e-30
294	Rn = rdust(ig,l)
295	Rn = -dlog(Rn)
296	Rm = Rn - 3. * sigma_ice*sigma_ice
297	yeahn = derf( (rb_cld(1)+Rn) *dev2)
298	yeahm = derf( (rb_cld(1)+Rm) *dev2)
299	do i = 1, nbin_cld
300	n_aer(i) = -0.5 * No * yeahn !! this ith previously computed
301	m_aer(i) = -0.5 * Mo * yeahm !! this ith previously computed
302	yeahn = derf( (rb_cld(i+1)+Rn) *dev2)
303	yeahm = derf( (rb_cld(i+1)+Rm) *dev2)
304	n_aer(i) = n_aer(i) + 0.5 * No * yeahn
305	m_aer(i) = m_aer(i) + 0.5 * Mo * yeahm
306	enddo
307	!!! MORE EFFICIENT COMPUTATIONALLY THAN
308	! Rn = rdust(ig,l)
309	! Rm = Rn * exp( 3. * sigma_ice**2. )
310	! Rn = 1. / Rn
311	! Rm = 1. / Rm
312	! do i = 1, nbin_cld
313	! n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)Rn) dev2 )
314	! & -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
315	! m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)Rm) dev2 )
316	! & -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
317	! enddo
318	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
319
320	! sumcheck = 0
321	! do i = 1, nbin_cld
322	! sumcheck = sumcheck + n_aer(i)
323	! enddo
324	! sumcheck = abs(sumcheck/No - 1)
325	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
326	! print*, "WARNING, No sumcheck PROBLEM"
327	! print*, "sumcheck, No",sumcheck, No
328	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
329	! print*, "Dust binned distribution", n_aer
330	! endif
331	!
332	! sumcheck = 0
333	! do i = 1, nbin_cld
334	! sumcheck = sumcheck + m_aer(i)
335	! enddo
336	! sumcheck = abs(sumcheck/Mo - 1)
337	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
338	! print*, "WARNING, Mo sumcheck PROBLEM"
339	! print*, "sumcheck, Mo",sumcheck, Mo
340	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
341	! print*, "Dust binned distribution", m_aer
342	! endif
343
344	c Get the rates of nucleation
345	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
346
347
348	dN = 0.
349	dM = 0.
350	rate_out(ig,l) = 0
351	do i = 1, nbin_cld
352	! schema simple
353	n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
354	m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
355	dN = dN + n_aer(i) * rate(i) * ptimestep
356	dM = dM + m_aer(i) * rate(i) * ptimestep
357	!rate_out(ig,l)=rate_out(ig,l)+rate(i)
358	enddo
359
360	Mo = zq0(ig,l,igcm_h2o_ice) +
361	& zq0(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
362	No = zq0(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
363	!write(,) "l,cloud particles,cloud mass",l, No, Mo
364	rhocloud(ig,l) = zq0(ig,l,igcm_h2o_ice) / Mo * rho_ice
365	& +zq0(ig,l,igcm_ccn_mass)* tauscaling(ig)
366	& / Mo * rho_dust
367	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
368	if ((Mo.lt.1.e-20) .or. (No.le.1)) then
369	rice(ig,l) = 1.e-8
370	else
371	rice(ig,l) =
372	& CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
373	endif
374	c nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
375	seq = exp( 2.sig(zt(ig,l))mh2o /
376	& (rho_icergpzt(ig,l)*rice(ig,l)) )
377
378	call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
379	& ph2o,ph2o/satu,seq,rice(ig,l),gr)
380
381	up = Cste * gr * rice(ig,l) * No * seq +
382	& zq(ig,l,igcm_h2o_vap)
383	dwn = Cste * gr * rice(ig,l) * No / zqsat(ig,l)+ 1.
384
385	Ctot = zq0(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_h2o_vap)
386
387	newvap = min(up/dwn,Ctot)
388
389	gr = gr * ( newvap/zqsat(ig,l) - seq )
390
391
392	dMice = min( max(Cste * No * rice(ig,l) * gr,
393	& -zq(ig,l,igcm_h2o_ice) ),
394	& zq(ig,l,igcm_h2o_vap) )
395
396	c----------- TESTS 1D output ---------
397	satu_out(ig,l) = satu
398	Mcon_out(ig,l) = dMice
399	newvap_out(ig,l) = newvap
400	gr_out(ig,l) = gr
401	dN_out(ig,l) = dN
402	dM_out(ig,l) = dM
403	c-------------------------------------
404
405	c Water ice
406	zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice) +
407	& dMice
408
409	c Water vapor
410	zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap) -
411	& dMice
412
413	c If all the ice particles sublimate, all the condensation
414	c nuclei are released:
415	if (zq(ig,l,igcm_h2o_ice).le.1e-30) then
416	c for coherence
417	dM = 0
418	dN = 0
419	dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig)
420	dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
421	c Water ice particles
422	zq(ig,l,igcm_h2o_ice) = 0.
423	c Dust particles
424	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) +
425	& zq(ig,l,igcm_ccn_mass)
426	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) +
427	& zq(ig,l,igcm_ccn_number)
428	c CCNs
429	zq(ig,l,igcm_ccn_mass) = 0.
430	zq(ig,l,igcm_ccn_number) = 0.
431	endif
432
433	dN = dN/ tauscaling(ig)
434	dM = dM/ tauscaling(ig)
435	c Dust particles
436	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM
437	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN
438	c CCNs
439	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
440	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
441
442
443	pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass)
444	& -zq0(ig,l,igcm_dust_mass))/ptimestep
445	pdqcloud(ig,l,igcm_dust_number)=(zq(ig,l,igcm_dust_number)
446	& -zq0(ig,l,igcm_dust_number))/ptimestep
447	pdqcloud(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
448	& -zq0(ig,l,igcm_ccn_mass))/ptimestep
449	pdqcloud(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
450	& -zq0(ig,l,igcm_ccn_number))/ptimestep
451	pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
452	& -zq0(ig,l,igcm_h2o_vap))/ptimestep
453	pdqcloud(ig,l,igcm_h2o_ice)=(zq(ig,l,igcm_h2o_ice)
454	& -zq0(ig,l,igcm_h2o_ice))/ptimestep
455
456
457	count = count +1
458	ELSE ! for coherence (rdust, rice computations etc ..)
459	zq(ig,l,igcm_dust_mass) = zq0(ig,l,igcm_dust_mass)
460	zq(ig,l,igcm_dust_number) = zq0(ig,l,igcm_dust_number)
461	zq(ig,l,igcm_ccn_mass) = zq0(ig,l,igcm_ccn_mass)
462	zq(ig,l,igcm_ccn_number) = zq0(ig,l,igcm_ccn_number)
463	zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice)
464	zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap)
465
466	! pour les sorties de test
467	satu_out(ig,l) = satu
468	Mcon_out(ig,l) = 0
469	newvap_out(ig,l) = zq(ig,l,igcm_h2o_vap)
470	gr_out(ig,l) = 0
471	dN_out(ig,l) = 0
472	dM_out(ig,l) = 0
473
474	ENDIF ! end if (saturation ratio > 1) or (there is h2o_ice)
475
476	c-----update temperature
477	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
478	pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
479	c pdtcloud(ig,l)=pdt(ig,l)-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
480
481	c----- update rice & rhocloud AFTER microphysic
482	Mo = zq(ig,l,igcm_h2o_ice) +
483	& zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
484	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
485	rhocloud(ig,l) = zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
486	& +zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
487	& / Mo * rho_dust
488	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
489
490	rice_out(ig,l)=rice(ig,l)
491	if ((Mo.lt.1.e-20) .or. (No.le.1)) then
492	rice(ig,l) = 1.e-8
493	else
494	rice(ig,l) =
495	& CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
496	endif
497	rice_out(ig,l)=rice(ig,l)-rice_out(ig,l)
498
499	nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
500
501	c----- update rdust and sedimentation radius
502	rdust(ig,l)=
503	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
504	& max(zq(ig,l,igcm_dust_number),0.01))
505	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
506
507	rsedcloud(ig,l)=max( rice(ig,l)(1.+nuice_sed)
508	& (1.+nuice_sed)*(1.+nuice_sed),
509	& rdust(ig,l) )
510	rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
511
512	ENDDO
513	ENDDO
514
515	! print*, 'SATU'
516	! print*, satu_out(1,:)
517
518	c------------------------------------------------------------------
519	c------------------------------------------------------------------
520	c------------------------------------------------------------------
521	c------------------------------------------------------------------
522	c------------------------------------------------------------------
523	c TESTS
524
525	IF (output_sca) then
526
527	print*, 'count is ',count, ' i.e. ',
528	& count100/(nlayngrid), '% for microphys computation'
529
530	dM_col(:) = 0
531	dN_col(:) = 0
532	Mdust_col(:) = 0
533	Ndust_col(:) = 0
534	Mccn_col(:) = 0
535	Nccn_col(:) = 0
536	dMice_col(:) = 0
537	drice_col(:) = 0
538	icetot(:) = 0
539	do l=1, nlay
540	do ig=1,ngrid
541	dM_col(ig) = dM_col(ig) +
542	& dM_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g
543	dN_col(ig) = dN_col(ig) +
544	& dN_out(ig,l)*(pplev(ig,l) - pplev(ig,l+1)) / g
545	Mdust_col(ig) = Mdust_col(ig) +
546	& zq(ig,l,igcm_dust_mass)*tauscaling(ig)
547	& *(pplev(ig,l) - pplev(ig,l+1)) / g
548	Ndust_col(ig) = Ndust_col(ig) +
549	& zq(ig,l,igcm_dust_number)*tauscaling(ig)
550	& *(pplev(ig,l) - pplev(ig,l+1)) / g
551	Mccn_col(ig) = Mccn_col(ig) +
552	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
553	& *(pplev(ig,l) - pplev(ig,l+1)) / g
554	Nccn_col(ig) = Nccn_col(ig) +
555	& zq(ig,l,igcm_ccn_number)*tauscaling(ig)
556	& *(pplev(ig,l) - pplev(ig,l+1)) / g
557	dMice_col(ig) = dMice_col(ig) +
558	& Mcon_out(ig,l)
559	& *(pplev(ig,l) - pplev(ig,l+1)) / g
560	drice_col(ig) = drice_col(ig) +
561	& rice_out(ig,l)*zq(ig,l,igcm_h2o_ice)
562	& *(pplev(ig,l) - pplev(ig,l+1)) / g
563	icetot(ig) = icetot(ig) +
564	& zq(ig,l,igcm_h2o_ice)
565	& *(pplev(ig,l) - pplev(ig,l+1)) / g
566	enddo ! of do ig=1,ngrid
567	enddo ! of do l=1,nlay
568
569	drice_col(ig) = drice_col(ig)/icetot(ig)
570
571	IF (ngrid.ne.1) THEN ! 3D
572	call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
573	& satu_out)
574	call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",2,
575	& dM_col)
576	call WRITEDIAGFI(ngrid,"dN_col","dN column","N",2,
577	& dN_col)
578	call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",2,
579	& Ndust_col)
580	call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",2,
581	& Mdust_col)
582	call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",2,
583	& Nccn_col)
584	call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",2,
585	& Mccn_col)
586	call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
587	& dM_out)
588	call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
589	& dN_out)
590	ENDIF
591
592	IF (ngrid.eq.1) THEN ! 1D
593
594	call WRITEDIAGFI(ngrid,"newvap","h2o newvap","kg",1,
595	& newvap_out)
596	call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
597	& gr_out)
598	call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
599	& rate_out)
600	call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
601	& dM_out)
602	call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
603	& dN_out)
604	call WRITEDIAGFI(ngrid,"dicetot","ice col var","kg/m2",0,
605	& dMice_col)
606	call WRITEDIAGFI(ngrid,"dricetot","ice col var","m",0,
607	& drice_col)
608	call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
609	& Mcon_out)
610	call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
611	& zqsat)
612	call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
613	& satu_out)
614	call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1,
615	& rice)
616	call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
617	& rdust)
618	call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
619	& rsedcloud)
620	call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
621	& rhocloud)
622	call WRITEDIAGFI(ngrid,"dM_col","dM column","kg",0,
623	& dM_col)
624	call WRITEDIAGFI(ngrid,"dN_col","dN column","N",0,
625	& dN_col)
626	call WRITEDIAGFI(ngrid,"Ndust_col","M column","N",0,
627	& Ndust_col)
628	call WRITEDIAGFI(ngrid,"Mdust_col","N column","kg",0,
629	& Mdust_col)
630	call WRITEDIAGFI(ngrid,"Nccn_col","M column","N",0,
631	& Nccn_col)
632	call WRITEDIAGFI(ngrid,"Mccn_col","N column","kg",0,
633	& Mccn_col)
634	ENDIF
635	ENDIF ! endif output_sca
636	c------------------------------------------------------------------
637	return
638	end

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

Download in other formats:

Original Format