Context Navigation

improvedclouds.F @ 631

Last change on this file since 631 was 626, checked in by tnavarro, 13 years ago
New scheme for the clouds, no more sub-timestep. Clouds sedimentation is done with the dust one in callsedim.F like it was before. Added latent heat for sublimating ground ice. Bugs corrected. THIS VERSION OF THE WATER CYCLE SHOULD NOT BE USED WITH THERMALS DUE TO NEGATIVE TRACERS ISSUES.
File size: 31.6 KB

Line
1	subroutine improvedclouds(ngrid,nlay,ptimestep,
2	& pplev,pplay,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	c A word about the radius growth ...
24	c A word about nucleation and ice growth strategies ...
25
26	c Authors: J.-B. Madeleine, based on the work by Franck Montmessin
27	c (October 2011)
28	c T. Navarro, debug,correction, new scheme (October-April 2011)
29	c A. Spiga, optimization (February 2012)
30	c------------------------------------------------------------------
31	#include "dimensions.h"
32	#include "dimphys.h"
33	#include "comcstfi.h"
34	#include "callkeys.h"
35	#include "tracer.h"
36	#include "comgeomfi.h"
37	#include "dimradmars.h"
38	#include "microphys.h"
39	c------------------------------------------------------------------
40	c Inputs:
41
42	INTEGER ngrid,nlay
43	integer nq ! nombre de traceurs
44	REAL ptimestep ! pas de temps physique (s)
45	REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
46	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
47
48	REAL pt(ngrid,nlay) ! temperature at the middle of the
49	! layers (K)
50	REAL pdt(ngrid,nlay) ! tendance temperature des autres
51	! param.
52	REAL pq(ngrid,nlay,nq) ! traceur (kg/kg)
53	REAL pdq(ngrid,nlay,nq) ! tendance avant condensation
54	! (kg/kg.s-1)
55	REAL tauscaling(ngridmx) ! Convertion factor for qdust and Ndust
56	REAL rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
57
58	c Outputs:
59	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
60	! (r_c in montmessin_2004)
61	REAL nuice(ngrid,nlay) ! Estimated effective variance
62	! of the size distribution
63	REAL rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
64	REAL rhocloud(ngridmx,nlayermx) ! Cloud density (kg.m-3)
65	REAL pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
66	! H2O(kg/kg.s-1)
67	REAL pdtcloud(ngrid,nlay) ! tendance temperature due
68	! a la chaleur latente
69
70	c------------------------------------------------------------------
71	c Local variables:
72
73	LOGICAL firstcall
74	DATA firstcall/.true./
75	SAVE firstcall
76
77	REAL*8 derf ! Error function
78	!external derf
79
80	REAL CBRT
81	EXTERNAL CBRT
82
83	INTEGER ig,l,i
84
85
86	REAL zq(ngridmx,nlayermx,nqmx) ! local value of tracers
87	REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
88	REAL zt(ngridmx,nlayermx) ! local value of temperature
89	REAL zqsat(ngridmx,nlayermx) ! saturation
90	REAL lw !Latent heat of sublimation (J.kg-1)
91	REAL Cste
92	REAL dMice ! mass of condensated ice
93	REAL sumcheck
94	REAL*8 ph2o ! Water vapor partial pressure (Pa)
95	REAL*8 satu ! Water vapor saturation ratio over ice
96	REAL*8 Mo,No
97	REAL*8 dN,dM,newvap
98	REAL*8 Rn, Rm, dev2, yeah, n_derf, m_derf
99	REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
100	REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
101	REAL*8 rate(nbin_cld) ! nucleation rate
102	!REAL*8 up,dwn,Ctot,gr
103	REAl*8 seq
104	REAL*8 sig ! Water-ice/air surface tension (N.m)
105	EXTERNAL sig
106
107	c Parameters of the size discretization
108	c used by the microphysical scheme
109	DOUBLE PRECISION, PARAMETER :: rmin_cld = 0.1e-6 ! Minimum radius (m)
110	DOUBLE PRECISION, PARAMETER :: rmax_cld = 10.e-6 ! Maximum radius (m)
111	DOUBLE PRECISION, PARAMETER :: rbmin_cld = 0.0001e-6
112	! Minimum boundary radius (m)
113	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
114	DOUBLE PRECISION vrat_cld ! Volume ratio
115	DOUBLE PRECISION rb_cld(nbin_cld+1)! boundary values of each rad_cld bin (m)
116	SAVE rb_cld
117	DOUBLE PRECISION dr_cld(nbin_cld) ! width of each rad_cld bin (m)
118	DOUBLE PRECISION vol_cld(nbin_cld) ! particle volume for each bin (m3)
119
120	REAL sigma_ice ! Variance of the ice and CCN distributions
121	SAVE sigma_ice
122
123	REAL tdicho, tmax, rmin, rmax, req, rdicho
124	REAL coeff0, coeff1, coeff2
125	REAL error_out(ngridmx,nlayermx)
126	REAL error2d(ngridmx)
127
128	REAL var1,var2,var3
129
130	REAL rccn, epsilon
131
132	c----------------------------------
133	c some outputs for 1D -- TESTS
134	REAL satu_out(ngridmx,nlayermx) ! satu ratio for output
135	REAL dN_out(ngridmx,nlayermx) ! mass variation for output
136	REAL dM_out(ngridmx,nlayermx) ! number variation for output
137	REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)
138	REAL gr_out(ngridmx,nlayermx) ! for 1d output
139	REAL rice_out(ngridmx,nlayermx) ! ice radius change
140	REAL rate_out(ngridmx,nlayermx) ! nucleation rate
141	REAL satubf(ngridmx,nlayermx),satuaf(ngridmx,nlayermx)
142	REAL ccnbf(ngridmx,nlayermx),ccnaf(ngridmx,nlayermx)
143	INTEGER count
144
145	LOGICAL test_flag ! flag for test/debuging outputs
146
147	test_flag = .false.
148	c----------------------------------
149	c----------------------------------
150
151	c------------------------------------------------------------------
152
153	IF (firstcall) THEN
154	!=============================================================
155	! 0. Definition of the size grid
156	!=============================================================
157	c rad_cld is the primary radius grid used for microphysics computation.
158	c The grid spacing is computed assuming a constant volume ratio
159	c between two consecutive bins; i.e. vrat_cld.
160	c vrat_cld is determined from the boundary values of the size grid:
161	c rmin_cld and rmax_cld.
162	c The rb_cld array contains the boundary values of each rad_cld bin.
163	c dr_cld is the width of each rad_cld bin.
164
165	c Volume ratio between two adjacent bins
166	vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbin_cld-1) *3.
167	vrat_cld = dexp(vrat_cld)
168	write(,) "vrat_cld", vrat_cld
169
170	rb_cld(1) = rbmin_cld
171	rad_cld(1) = rmin_cld
172	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
173
174	do i=1,nbin_cld-1
175	rad_cld(i+1) = rad_cld(i) * vrat_cld**(1./3.)
176	vol_cld(i+1) = vol_cld(i) * vrat_cld
177	enddo
178
179	do i=1,nbin_cld
180	rb_cld(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
181	& rad_cld(i)
182	dr_cld(i) = rb_cld(i+1) - rb_cld(i)
183	enddo
184	rb_cld(nbin_cld+1) = rbmax_cld
185	dr_cld(nbin_cld) = rb_cld(nbin_cld+1) - rb_cld(nbin_cld)
186
187	print*, ' '
188	print*,'Microphysics: size bin information:'
189	print*,'i,rb_cld(i), rad_cld(i),dr_cld(i)'
190	print*,'-----------------------------------'
191	do i=1,nbin_cld
192	write(*,'(i2,3x,3(e12.6,4x))') i,rb_cld(i), rad_cld(i),
193	& dr_cld(i)
194	enddo
195	write(*,'(i2,3x,e12.6)') nbin_cld+1,rb_cld(nbin_cld+1)
196	print*,'-----------------------------------'
197
198	do i=1,nbin_cld+1
199	rb_cld(i) = dlog(rb_cld(i)) !! we save that so that it is not computed
200	!! at each timestep and gridpoint
201	enddo
202
203	c Contact parameter of water ice on dust ( m=cos(theta) )
204	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
205	! mteta = 0.95
206	write(,) 'water_param contact parameter:', mteta
207
208	c Volume of a water molecule (m3)
209	vo1 = mh2o / dble(rho_ice)
210	c Variance of the ice and CCN distributions
211	sigma_ice = sqrt(log(1.+nuice_sed))
212
213	write(,) 'Variance of ice & CCN distribs :', sigma_ice
214	write(,) 'nuice for sedimentation:', nuice_sed
215	write(,) 'Volume of a water molecule:', vo1
216
217	firstcall=.false.
218	END IF
219
220	!=============================================================
221	! 1. Initialisation
222	!=============================================================
223
224	c Initialize the tendencies
225	pdqcloud(1:ngrid,1:nlay,1:nq)=0
226	pdtcloud(1:ngrid,1:nlay)=0
227
228	c Update the needed variables
229	do l=1,nlay
230	do ig=1,ngrid
231	c Temperature
232	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
233	c Dust mass mixing ratio
234	zq(ig,l,igcm_dust_mass) =
235	& pq(ig,l,igcm_dust_mass) +
236	& pdq(ig,l,igcm_dust_mass) * ptimestep
237	zq0(ig,l,igcm_dust_mass)=zq(ig,l,igcm_dust_mass)
238	c Dust particle number
239	zq(ig,l,igcm_dust_number) =
240	& pq(ig,l,igcm_dust_number) +
241	& pdq(ig,l,igcm_dust_number) * ptimestep
242	zq0(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)
243	c Update rdust from last tendencies
244	rdust(ig,l)=
245	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
246	& max(zq(ig,l,igcm_dust_number),0.01))
247	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
248	c CCN mass mixing ratio
249	zq(ig,l,igcm_ccn_mass)=
250	& pq(ig,l,igcm_ccn_mass)+pdq(ig,l,igcm_ccn_mass)*ptimestep
251	zq(ig,l,igcm_ccn_mass)=max(zq(ig,l,igcm_ccn_mass),1.E-30)
252	zq0(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
253	c CCN particle number
254	zq(ig,l,igcm_ccn_number)=
255	& pq(ig,l,igcm_ccn_number)+pdq(ig,l,igcm_ccn_number)*ptimestep
256	zq(ig,l,igcm_ccn_number)=max(zq(ig,l,igcm_ccn_number),1.E-30)
257	zq0(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
258	c Water vapor
259	zq(ig,l,igcm_h2o_vap)=
260	& pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
261	zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004
262	zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
263	c Water ice
264	zq(ig,l,igcm_h2o_ice)=
265	& pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
266	zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),1e-30) ! FF 12/2004
267	zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
268	enddo
269	enddo
270
271	!=============================================================
272	! 2. Compute saturation
273	!=============================================================
274
275	dev2 = 1. / ( sqrt(2.) * sigma_ice )
276	error_out(:,:) = 0.
277
278
279	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
280
281	count = 0
282
283	c Main loop over the GCM's grid
284	DO l=1,nlay
285	DO ig=1,ngrid
286
287	c Get the partial pressure of water vapor and its saturation ratio
288	ph2o = zq(ig,l,igcm_h2o_vap) * (44./18.) * pplay(ig,l)
289	satu = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
290
291	IF (( satu .ge. 1. ) .or. ! if there is condensation
292	& ( zq(ig,l,igcm_ccn_number)*tauscaling(ig).ge. 1.)) THEN ! or sublimation
293
294
295	!=============================================================
296	! 3. Nucleation
297	!=============================================================
298
299	c Expand the dust moments into a binned distribution
300	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig) !+ 1.e-30
301	No = zq(ig,l,igcm_dust_number)* tauscaling(ig) + 1.e-30
302	Rn = rdust(ig,l)
303	Rn = -dlog(Rn)
304	Rm = Rn - 3. * sigma_ice*sigma_ice
305	n_derf = derf( (rb_cld(1)+Rn) *dev2)
306	m_derf = derf( (rb_cld(1)+Rm) *dev2)
307	do i = 1, nbin_cld
308	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
309	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
310	n_derf = derf( (rb_cld(i+1)+Rn) *dev2)
311	m_derf = derf( (rb_cld(i+1)+Rm) *dev2)
312	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
313	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
314	enddo
315	!!! MORE EFFICIENT COMPUTATIONALLY THAN
316	! Rn = rdust(ig,l)
317	! Rm = Rn * exp( 3. * sigma_ice**2. )
318	! Rn = 1. / Rn
319	! Rm = 1. / Rm
320	! do i = 1, nbin_cld
321	! n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)Rn) dev2 )
322	! & -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
323	! m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)Rm) dev2 )
324	! & -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
325	! enddo
326	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
327
328	! sumcheck = 0
329	! do i = 1, nbin_cld
330	! sumcheck = sumcheck + n_aer(i)
331	! enddo
332	! sumcheck = abs(sumcheck/No - 1)
333	! if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
334	! print*, "WARNING, No sumcheck PROBLEM"
335	! print*, "sumcheck, No",sumcheck, No
336	! print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
337	! print*, "Dust binned distribution", n_aer
338	! endif
339	!
340	! sumcheck = 0
341	! do i = 1, nbin_cld
342	! sumcheck = sumcheck + m_aer(i)
343	! enddo
344	! sumcheck = abs(sumcheck/Mo - 1)
345	! if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld)) then
346	! print*, "WARNING, Mo sumcheck PROBLEM"
347	! print*, "sumcheck, Mo",sumcheck, Mo
348	! print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig, l
349	! print*, "Dust binned distribution", m_aer
350	! endif
351
352	c Get the rates of nucleation
353	call nuclea(ph2o,zt(ig,l),satu,n_aer,rate)
354
355
356	dN = 0.
357	dM = 0.
358	rate_out(ig,l) = 0
359	do i = 1, nbin_cld
360	n_aer(i) = n_aer(i) / ( 1. + rate(i)*ptimestep )
361	m_aer(i) = m_aer(i) / ( 1. + rate(i)*ptimestep )
362	dN = dN + n_aer(i) * rate(i) * ptimestep
363	dM = dM + m_aer(i) * rate(i) * ptimestep
364	!rate_out(ig,l)=rate_out(ig,l)+rate(i)
365	enddo
366
367	c Update CCNs, can also be done after the radius growth
368	c Dust particles
369	zq(ig,l,igcm_dust_mass) =
370	& zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig)
371	zq(ig,l,igcm_dust_number) =
372	& zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig)
373	c CCNs
374	zq(ig,l,igcm_ccn_mass) =
375	& zq(ig,l,igcm_ccn_mass) + dM/ tauscaling(ig)
376	zq(ig,l,igcm_ccn_number) =
377	& zq(ig,l,igcm_ccn_number) + dN/ tauscaling(ig)
378
379
380	!=============================================================
381	! 4. Ice growth: scheme for radius evolution
382	!=============================================================
383
384	Mo = zq(ig,l,igcm_h2o_ice) +
385	& zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
386	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
387
388	rhocloud(ig,l) = zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
389	& + zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
390	& / Mo * rho_dust
391	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
392
393	c nuclei radius
394	rccn = CBRT(zq(ig,l,igcm_ccn_mass)/
395	& (pirho_dustzq(ig,l,igcm_ccn_number)*4./3.))
396
397	c ice crystal radius
398	rice (ig,l) =
399	& CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) )
400
401	c enforce physical value : crystal cannot be smaller than its nuclei !
402	rice(ig,l) = max(rice(ig,l), rccn)
403
404	c saturation at equilibrium
405	seq = exp( 2.sig(zt(ig,l))mh2o /
406	& (rho_icergpzt(ig,l)*rice(ig,l)) )
407
408
409	c If there is more than on nuclei, we peform ice growth
410	var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN
411	IF (var1 .ge. -1) THEN
412
413
414	if (test_flag) then
415	print*, ' '
416	print*, ptimestep
417	print*, 'satu,seq', real(satu), real(seq), ig,l
418	print*, 'dN,dM', real(dN),real(dM)
419	print*,'rccn', rccn
420	print, 'Nccn, Mccn', zq(ig,l,igcm_ccn_number)tauscaling(ig),
421	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
422	endif
423
424	c crystal radius to reach saturation at equilibrium (i.e. satu=seq)
425	req = ( zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap)
426	& + zq(ig,l,igcm_ccn_mass)tauscaling(ig)
427	& rho_ice/rho_dust - seq * zqsat(ig,l))
428	& / ( zq(ig,l,igcm_ccn_number)tauscaling(ig)
429	& pirho_ice4./3. )
430	req = CBRT(req)
431
432	c compute ice radius growth resistances (diffusive and latent heat resistancea)
433	call growthrate(zt(ig,l),pplay(ig,l),
434	& ph2o/satu,seq,req,coeff1,coeff2)
435
436	coeff0 = -zqsat(ig,l) / (4.pireq*rho_ice
437	& * zq(ig,l,igcm_ccn_number)*tauscaling(ig))
438
439	c compute tmax, the time needed to reach req
440	call phi(rice(ig,l),req,coeff1,coeff2,tmax)
441
442	if (test_flag) then
443	print*, 'coeffs', coeff0,coeff1,coeff2
444	print, 'req,tmax', req,tmaxcoeff0
445	print*, 'i,rmin,rdicho,rmax,tdicho'
446	endif
447
448	c rmin is rice if r increases (satu >1) or req if it decreases (satu<1)
449	c if req is lower than rccn (ie not enough ice to reach saturation), rmin is forced to rccn
450	if (satu .ge. seq) then
451	! crystal size is increasing
452	rmin = max(min(rice(ig,l),req),rccn)
453	rmax = max(rice(ig,l),req)
454	else
455	! crystal size is decreasing
456	rmax = max(min(rice(ig,l),req),rccn)
457	rmin = max(rice(ig,l),req)
458	endif
459	!rmax = min(rmax,1.e-3) ! au max on a des rayons de 1 mm pour la dicho ...
460	rdicho = 0.5*(rmin+rmax)
461
462	! for output
463	var1 = rice(ig,l)
464	var2 = rmin
465	var3 = rmax
466
467	c Given the phi function is monotonous, we perform a simple dichotomy to find the raidus at t+1
468	do i = 1,10 ! dichotomy loop
469
470	c compute tdicho, the time needed to reach rdicho
471	call phi(rdicho,req,coeff1,coeff2,tdicho)
472	!print*, tdicho,tmax
473	tdicho = coeff0*(tdicho - tmax)
474
475	if (test_flag) print*, i,rmin,rdicho,rmax,tdicho
476
477	if (tdicho .ge. ptimestep) then
478	rmax = rdicho
479	else
480	rmin = rdicho
481	endif
482
483	rdicho = 0.5*(rmin+rmax)
484
485	enddo ! of dichotomy loop
486
487	if (test_flag) then
488	if (abs(rdicho - rccn) .ge. 1e-15) then ! to avoid infinite values
489	epsilon = (rmax - rmin)/(2**10)
490	error_out(ig,l) =
491	& 100(epsilon3 +3epsilon*2rdicho +3epsilonrdicho**2)
492	& / (rdicho3-rccn3)
493	endif
494	print*, 'error masse glace %', error_out(ig,l)
495	print*, 'rice,ice,vap bf',
496	& rice(ig,l),zq0(ig,l,igcm_h2o_ice),zq0(ig,l,igcm_h2o_vap)
497	endif
498
499	c If the initial condition is subsaturated and there is not enough ice available for sublimation
500	c to reach equilibrium, req is neagtive. Therefore, enforce physical value.
501	rice(ig,l) = max(rdicho,rccn)
502
503	!!!!! Water ice mass is computed with radius at t+1 and their current number
504	!!!!! Nccn is at t or t+1, depending on what has been done before
505	! zq(ig,l,igcm_h2o_ice) =
506	! & (pirho_icezq(ig,l,igcm_ccn_number)*4./3.
507	! & * rdichordichordicho -
508	! & zq(ig,l,igcm_ccn_mass)*rho_ice/rho_dust)
509	! & * tauscaling(ig)
510
511	!!!!! Water ice mass is computed with radius at t+1 and their number at t+1
512	!!!!! that is dirty, but this enforces the use of Nccn at t+1, whatever is done before
513	!!!!! TO BE CLEANED ONE DAY
514	var1 = zq0(ig,l,igcm_ccn_number)*tauscaling(ig) + dN
515	var2 = zq0(ig,l,igcm_ccn_mass)*tauscaling(ig) + dM
516	zq(ig,l,igcm_h2o_ice) =
517	& (pirho_icevar1*4./3.
518	& * rdichordichordicho -
519	& var2*rho_ice/rho_dust)
520
521
522	!!!! enforce realistic values (if the case of growth on Nccn(t) and condensation on Nccn(t+1))
523	zq(ig,l,igcm_h2o_ice) =
524	& min(zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap),
525	& zq(ig,l,igcm_h2o_ice))
526	zq(ig,l,igcm_h2o_ice) =
527	& max(1e-30,zq(ig,l,igcm_h2o_ice))
528
529	zq(ig,l,igcm_h2o_vap) =
530	& zq0(ig,l,igcm_h2o_ice) + zq0(ig,l,igcm_h2o_vap)
531	& - zq(ig,l,igcm_h2o_ice)
532
533
534	if (test_flag) then
535	print*, 'rice,ice,vap af',
536	& rice(ig,l),zq(ig,l,igcm_h2o_ice),zq(ig,l,igcm_h2o_vap)
537	print*, 'satu bf, af',
538	& zq0(ig,l,igcm_h2o_vap)/zqsat(ig,l),
539	& zq(ig,l,igcm_h2o_vap)/zqsat(ig,l)
540	endif
541
542
543	!!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
544	if ((zq(ig,l,igcm_h2o_ice).le. -1e-8)
545	& .or. (zq(ig,l,igcm_h2o_vap).le. -1e-8)) then
546	print*, 'NEGATIVE WATER'
547	print*, 'ig,l', ig,l
548	print*, 'satu', satu
549	print*, 'vap, ice bf',
550	& zq0(ig,l,igcm_h2o_vap), zq0(ig,l,igcm_h2o_ice)
551	print*, 'vap, ice af',
552	& zq(ig,l,igcm_h2o_vap), zq(ig,l,igcm_h2o_ice)
553	print*, 'ccn N,M bf',
554	& zq0(ig,l,igcm_ccn_number), zq0(ig,l,igcm_ccn_mass)
555	print*, 'ccn N,M af',
556	& zq(ig,l,igcm_ccn_number), zq(ig,l,igcm_ccn_mass)
557	print*, 'tauscaling',
558	& tauscaling(ig)
559	print*, 'req,rccn,rice bf,rice af',
560	& req,rccn,var1,rice(ig,l)
561	print*, 'rmin, rmax', var2,var3
562	print*, 'error_out,tdicho,ptimestep',
563	& error_out(ig,l),tdicho,ptimestep
564	print*, ' '
565	endif
566	!!!!!!!!!!!! TEST TEST TEST TEST TEST TEST TEST TEST TEST TEST
567
568
569	ENDIF !of if Nccn >1
570
571
572	!=============================================================
573	! 5. Dust cores released, tendancies, latent heat, etc ...
574	!=============================================================
575
576	c If all the ice particles sublimate, all the condensation
577	c nuclei are released:
578	if (zq(ig,l,igcm_h2o_ice).le.1e-10) then
579	! for coherence
580	! dM = 0
581	! dN = 0
582	! dN_out(ig,l) = - zq(ig,l,igcm_ccn_number)*tauscaling(ig)
583	! dM_out(ig,l) = - zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
584	c Water
585	zq(ig,l,igcm_h2o_vap) = zq(ig,l,igcm_h2o_vap)
586	& + zq(ig,l,igcm_h2o_ice)
587	zq(ig,l,igcm_h2o_ice) = 0.
588	c Dust particles
589	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass)
590	& + zq(ig,l,igcm_ccn_mass)
591	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number)
592	& + zq(ig,l,igcm_ccn_number)
593	c CCNs
594	zq(ig,l,igcm_ccn_mass) = 0.
595	zq(ig,l,igcm_ccn_number) = 0.
596	endif
597
598
599	! dN = dN/ tauscaling(ig)
600	! dM = dM/ tauscaling(ig)
601	!c Dust particles
602	! zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dM
603	! zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dN
604	!c CCNs
605	! zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + dM
606	! zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + dN
607
608
609	pdqcloud(ig,l,igcm_dust_mass)=(zq(ig,l,igcm_dust_mass)
610	& -zq0(ig,l,igcm_dust_mass))/ptimestep
611	pdqcloud(ig,l,igcm_dust_number)=(zq(ig,l,igcm_dust_number)
612	& -zq0(ig,l,igcm_dust_number))/ptimestep
613	pdqcloud(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
614	& -zq0(ig,l,igcm_ccn_mass))/ptimestep
615	pdqcloud(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
616	& -zq0(ig,l,igcm_ccn_number))/ptimestep
617	pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
618	& -zq0(ig,l,igcm_h2o_vap))/ptimestep
619	pdqcloud(ig,l,igcm_h2o_ice)=(zq(ig,l,igcm_h2o_ice)
620	& -zq0(ig,l,igcm_h2o_ice))/ptimestep
621
622
623	count = count +1
624
625
626	ELSE ! for coherence (rdust, rice computations etc ..)
627	zq(ig,l,igcm_dust_mass) = zq0(ig,l,igcm_dust_mass)
628	zq(ig,l,igcm_dust_number) = zq0(ig,l,igcm_dust_number)
629	zq(ig,l,igcm_ccn_mass) = zq0(ig,l,igcm_ccn_mass)
630	zq(ig,l,igcm_ccn_number) = zq0(ig,l,igcm_ccn_number)
631	zq(ig,l,igcm_h2o_ice) = zq0(ig,l,igcm_h2o_ice)
632	zq(ig,l,igcm_h2o_vap) = zq0(ig,l,igcm_h2o_vap)
633
634	! pour les sorties de test
635	! satu_out(ig,l) = satu
636	! gr_out(ig,l) = 0
637	! dN_out(ig,l) = 0
638	! dM_out(ig,l) = 0
639
640	ENDIF ! end if (saturation ratio > 1) or (there is h2o_ice)
641
642	! ccnbf(ig,l) = zq0(ig,l,igcm_ccn_number) * tauscaling(ig)
643	! ccnaf(ig,l) = zq(ig,l,igcm_ccn_number) * tauscaling(ig)
644	!
645	! satubf(ig,l) = zq0(ig,l,igcm_h2o_vap) / zqsat(ig,l)
646	! satuaf(ig,l) = zq(ig,l,igcm_h2o_vap) / zqsat(ig,l)
647
648
649	c-----update temperature (latent heat relase)
650	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
651	pdtcloud(ig,l)= -pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp
652
653	c----- update rice & rhocloud AFTER microphysic
654	Mo = zq(ig,l,igcm_h2o_ice) +
655	& zq(ig,l,igcm_ccn_mass)* tauscaling(ig) + 1.e-30
656	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig)+ 1e-30
657	rhocloud(ig,l) = zq(ig,l,igcm_h2o_ice) / Mo * rho_ice
658	& +zq(ig,l,igcm_ccn_mass)* tauscaling(ig)
659	& / Mo * rho_dust
660	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
661
662	if ((Mo.lt.1.e-20) .or. (No.le.1)) then
663	rice(ig,l) = 1.e-8
664	else
665	rice(ig,l) =
666	& CBRT( real(Mo)/real(No) * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
667	endif
668
669	nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
670
671	c----- update rdust and sedimentation radius
672	rdust(ig,l)=
673	& CBRT(r3n_q*zq(ig,l,igcm_dust_mass)/
674	& max(zq(ig,l,igcm_dust_number),0.01))
675	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
676
677	rsedcloud(ig,l)=max( rice(ig,l)(1.+nuice_sed)
678	& (1.+nuice_sed)*(1.+nuice_sed),
679	& rdust(ig,l) )
680	rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
681
682	ENDDO ! of ig loop
683	ENDDO ! of nlayer loop
684
685
686	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
687	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
688	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
689	IF (test_flag) then
690
691	error2d(:) = 0.
692	DO l=1,nlay
693	DO ig=1,ngrid
694	error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
695	ENDDO
696	ENDDO
697
698	print*, 'count is ',count, ' i.e. ',
699	& count100/(nlayngrid), '% for microphys computation'
700
701	IF (ngrid.ne.1) THEN ! 3D
702	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
703	! & satu_out)
704	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
705	! & dM_out)
706	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
707	! & dN_out)
708	call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
709	& error2d)
710	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
711	! & zqsat)
712	ENDIF
713
714	IF (ngrid.eq.1) THEN ! 1D
715	call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
716	& error_out)
717	call WRITEDIAGFI(ngrid,"satu_bf","satu before","kg/kg",1,
718	& satubf)
719	call WRITEDIAGFI(ngrid,"satu_af","satu after","kg/kg",1,
720	& satuaf)
721	call WRITEDIAGFI(ngrid,"vapbf","h2ovap before","kg/kg",1,
722	& zq0(1,:,igcm_h2o_vap))
723	call WRITEDIAGFI(ngrid,"vapaf","h2ovap after","kg/kg",1,
724	& zq(1,:,igcm_h2o_vap))
725	call WRITEDIAGFI(ngrid,"icebf","h2oice before","kg/kg",1,
726	& zq0(1,:,igcm_h2o_ice))
727	call WRITEDIAGFI(ngrid,"iceaf","h2oice after","kg/kg",1,
728	& zq(1,:,igcm_h2o_ice))
729	call WRITEDIAGFI(ngrid,"ccnbf","ccn before","/kg",1,
730	& ccnbf)
731	call WRITEDIAGFI(ngrid,"ccnaf","ccn after","/kg",1,
732	& ccnaf)
733	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
734	c & gr_out)
735	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
736	c & rate_out)
737	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
738	c & dM_out)
739	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
740	c & dN_out)
741	call WRITEDIAGFI(ngrid,"zqsat","p vap sat","kg/kg",1,
742	& zqsat)
743	call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
744	& satu_out)
745	call WRITEDIAGFI(ngrid,"rice_sca","ice radius","m",1,
746	& rice)
747	call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
748	& rdust)
749	call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
750	& rsedcloud)
751	call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
752	& rhocloud)
753	ENDIF
754
755	ENDIF ! endif test_flag
756	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
757	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
758	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
759
760	return
761	end
762
763
764
765	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
766	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
767	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
768	c It is an analytical solution to the ice radius growth equation,
769	c with the approximation of a constant 'reduced' cunningham correction factor
770	c (lambda in growthrate.F) taken at radius req instead of rice
771	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
772	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
773
774	subroutine phi(rice,req,coeff1,coeff2,time)
775
776	implicit none
777
778	! inputs
779	real rice ! ice radius
780	real req ! ice radius at equilibirum
781	real coeff1 ! coeff for the log
782	real coeff2 ! coeff for the arctan
783
784	! output
785	real time
786
787	!local
788	real var
789
790	! 1.73205 is sqrt(3)
791
792	var = max(
793	& abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
794
795	time =
796	& coeff1 *
797	& log( var )
798	& + coeff2 * 1.73205 *
799	& atan( (2rice+req) / (1.73205req) )
800
801	return
802	end
803
804
805

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

Download in other formats:

Original Format