Context Navigation

improvedCO2clouds.F @ 1619

Last change on this file since 1619 was 1617, checked in by jaudouard, 8 years ago
Added modifications for CO2 clouds scheme in physiq_mod.F and added several routines and variables for CO2 microphysics. October 2016 Joachim Audouard
Property svn:executable set to ``*
File size: 32.1 KB

Line
1	subroutine improvedCO2clouds(ngrid,nlay,ptimestep,
2	& pplay,pt,pdt,
3	& pq,pdq,pdqcloudco2,pdtcloudco2,
4	& nq,tauscaling)
5	! to use 'getin'
6	USE comcstfi_h
7	USE ioipsl_getincom
8	USE updaterad
9	use tracer_mod
10	!, only: rho_ice_co2, nuiceco2_sed, igcm_co2,
11	! & rho_ice,igcm_h2o_ice, igcm_ccn_number,
12	! & igcm_co2_ice, igcm_dust_mass,
13	! & igcm_dust_number, igcm_ccnco2_mass,
14	! & igcm_ccnco2_number
15	use conc_mod, only: mmean
16	implicit none
17
18
19	c------------------------------------------------------------------
20	c This routine is used to form CO2 clouds when a parcel of the GCM is
21	c saturated. It includes the ability to have supersaturation, a
22	c computation of the nucleation rates, growthrates and the
23	c scavenging of dust particles by clouds.
24	c It is worth noting that the amount of dust is computed using the
25	c dust optical depth computed in aeropacity.F. That's why
26	c the variable called "tauscaling" is used to convert
27	c pq(dust_mass) and pq(dust_number), which are relative
28	c quantities, to absolute and realistic quantities stored in zq.
29	c This has to be done to convert the inputs into absolute
30	c values, but also to convert the outputs back into relative
31	c values which are then used by the sedimentation and advection
32	c schemes.
33
34	c Authors of the water ice clouds microphysics
35	c J.-B. Madeleine, based on the work by Franck Montmessin
36	c (October 2011)
37	c T. Navarro, debug,correction, new scheme (October-April 2011)
38	c A. Spiga, optimization (February 2012)
39	c Adaptation for CO2 clouds by Joachim Audouard (09/16), based on the work
40	c of Constantino Listowski
41	c------------------------------------------------------------------
42	!#include "dimensions.h"
43	!#include "dimphys.h"
44	#include "callkeys.h"
45	!#include "tracer.h"
46	!#include "comgeomfi.h"
47	!#include "dimradmars.h"
48	#include "microphys.h"
49	!#include "microphysCO2.h"
50	!#include "conc.h"
51	c------------------------------------------------------------------
52	c Inputs:
53
54	INTEGER ngrid,nlay
55	integer nq ! nombre de traceurs
56	REAL ptimestep ! pas de temps physique (s)
57	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
58
59	REAL pt(ngrid,nlay) ! temperature at the middle of the
60	! layers (K)
61	REAL pdt(ngrid,nlay) ! tendance temperature des autres
62	! param.
63	REAL pq(ngrid,nlay,nq) ! traceur (kg/kg)
64	REAL pdq(ngrid,nlay,nq) ! tendance avant condensation
65	! (kg/kg.s-1)
66	REAL tauscaling(ngrid) ! Convertion factor for qdust and Ndust
67
68	real rice(ngrid,nlay) ! Water Ice mass mean radius (m)
69	! used for nucleation of CO2 on ice-coated ccns
70
71	c Outputs:
72	REAL pdqcloudco2(ngrid,nlay,nq) ! tendance de la condensation
73	! CO2 (kg/kg.s-1)
74	! condensation si igcm_co2_ice
75	REAL pdtcloudco2(ngrid,nlay) ! tendance temperature due
76	! a la chaleur latente
77
78	c------------------------------------------------------------------
79	c Local variables:
80
81	LOGICAL firstcall
82	DATA firstcall/.true./
83	SAVE firstcall
84
85	REAL*8 derf ! Error function
86	!external derf
87
88	!REAL*8 massflowrateCO2
89	!external massflowrateCO2
90
91	INTEGER ig,l,i
92
93	REAL zq(ngrid,nlay,nq) ! local value of tracers
94	REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
95	REAL zt(ngrid,nlay) ! local value of temperature
96	REAL zqsat(ngrid,nlay) ! saturation vapor pressure for CO2
97	REAL lw !Latent heat of sublimation (J.kg-1)
98	REAL l0,l1,l2,l3,l4
99	REAL cste
100	REAL dMice ! mass of condensed ice
101	DOUBLE PRECISION sumcheck
102	DOUBLE PRECISION pco2 ! Co2 vapor partial pressure (Pa)
103	DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice
104	DOUBLE PRECISION Mo,No
105	DOUBLE PRECISION Rn, Rm, dev2, n_derf, m_derf
106
107	! Radius used by the microphysical scheme (m)
108	DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin
109	DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin
110
111	DOUBLE PRECISION n_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
112	DOUBLE PRECISION m_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
113	DOUBLE PRECISION rad_h2oice(nbinco2_cld)
114
115	c REAL*8 sigco2 ! Co2-ice/air surface tension (N.m)
116	c EXTERNAL sigco2
117
118	DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM
119	DOUBLE PRECISION rate(nbinco2_cld) ! nucleation rate
120	DOUBLE PRECISION rateh2o(nbinco2_cld) ! nucleation rate
121	REAL seq
122
123	DOUBLE PRECISION riceco2(ngrid,nlay) ! CO2Ice mean radius (m)
124
125	REAL rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
126	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
127
128	c REAL res ! Resistance growth
129	DOUBLE PRECISION Ic_rice ! Mass transfer rate CO2 ice crystal
130
131
132	c Parameters of the size discretization
133	c used by the microphysical scheme
134	DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m)
135	DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-5 ! Maximum radius (m)
136	DOUBLE PRECISION, PARAMETER :: rbmin_cld =0.099e-9
137	! Minimum boundary radius (m)
138	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-2 ! Maximum boundary radius (m)
139	DOUBLE PRECISION vrat_cld ! Volume ratio
140	DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
141	SAVE rb_cldco2
142
143	DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m)
144	DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3)
145
146	DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff
147	REAL sigma_iceco2 ! Variance of the ice and CCN distributions
148	SAVE sigma_iceco2
149
150
151
152	c----------------------------------
153	c TESTS
154
155	INTEGER countcells
156
157	LOGICAL test_flag ! flag for test/debuging outputs
158	SAVE test_flag
159
160
161	REAL satubf(ngrid,nlay),satuaf(ngrid,nlay)
162	REAL res_out(ngrid,nlay)
163
164
165	c------------------------------------------------------------------
166
167	IF (firstcall) THEN
168	!=============================================================
169	! 0. Definition of the size grid
170	!=============================================================
171	c rad_cldco2 is the primary radius grid used for microphysics computation.
172	c The grid spacing is computed assuming a constant volume ratio
173	c between two consecutive bins; i.e. vrat_cld.
174	c vrat_cld is determined from the boundary values of the size grid:
175	c rmin_cld and rmax_cld.
176	c The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
177	c dr_cld is the width of each rad_cldco2 bin.
178
179	c Volume ratio between two adjacent bins
180	! vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
181	! vrat_cld = exp(vrat_cld)
182	vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
183	vrat_cld = exp(vrat_cld)
184	c write(,) "vrat_cld", vrat_cld
185
186	rb_cldco2(1) = rbmin_cld
187	rad_cldco2(1) = rmin_cld
188	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
189	! vol_cld(1) = 4./3. * pi * rmin_cldrmin_cldrmin_cld
190
191	do i=1,nbinco2_cld-1
192	rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.)
193	vol_cld(i+1) = vol_cld(i) * vrat_cld
194	enddo
195
196	do i=1,nbinco2_cld
197	rb_cldco2(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
198	& rad_cldco2(i)
199	dr_cld(i) = rb_cldco2(i+1) - rb_cldco2(i)
200	enddo
201	rb_cldco2(nbinco2_cld+1) = rbmax_cld
202	dr_cld(nbinco2_cld) = rb_cldco2(nbinco2_cld+1) -
203	& rb_cldco2(nbinco2_cld)
204
205	print*, ' '
206	print*,'Microphysics co2: size bin information:'
207	print*,'i,rb_cldco2(i), rad_cldco2(i),dr_cld(i)'
208	print*,'-----------------------------------'
209	do i=1,nbinco2_cld
210	write(*,'(i3,3x,3(e12.6,4x))') i,rb_cldco2(i), rad_cldco2(i),
211	& dr_cld(i)
212	enddo
213	write(*,'(i3,3x,e12.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1)
214	print*,'-----------------------------------'
215
216	do i=1,nbinco2_cld+1
217	rb_cldco2(i) = log(rb_cldco2(i)) !! we save that so that it is not computed
218	!! at each timestep and gridpoint
219	enddo
220
221	c Contact parameter of co2 ice on dst ( m=cos(theta) )
222	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
223	c mteta = 0.952
224	c mtetaco2 = 0.952
225	c write(,) 'co2_param contact parameter:', mtetaco2
226
227	c Volume of a co2 molecule (m3)
228	vo1 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
229	vo1co2=vo1 ! AJOUT JA
230	c Variance of the ice and CCN distributions
231	sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
232
233
234	c write(,) 'Variance of ice & CCN distribs :', sigma_iceco2
235	c write(,) 'nuice for sedimentation:', nuiceco2_sed
236	c write(,) 'Volume of a co2 molecule:', vo1co2
237
238
239	test_flag = .false.
240	firstcall=.false.
241
242	END IF
243
244
245	!=============================================================
246	! 1. Initialisation
247	!=============================================================
248	!cste = 4pirho_ice*ptimestep !not used for co2
249
250	res_out(:,:) = 0
251	rice(:,:) = 1.e-11
252	riceco2(:,:) = 1.e-11
253
254	c Initialize the tendencies
255	pdqcloudco2(1:ngrid,1:nlay,1:nq)=0.
256	pdtcloudco2(1:ngrid,1:nlay)=0.
257
258	c pt temperature layer; pdt dT.s-1
259	c pq traceur kg/kg; pdq tendance idem .s-1
260	zt(1:ngrid,1:nlay) =
261	& pt(1:ngrid,1:nlay) +
262	& pdt(1:ngrid,1:nlay) * ptimestep
263
264	zq(1:ngrid,1:nlay,1:nq) =
265	& pq(1:ngrid,1:nlay,1:nq) +
266	& pdq(1:ngrid,1:nlay,1:nq) * ptimestep
267
268
269	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
270	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
271
272	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
273
274	!=============================================================
275	! 2. Compute saturation
276	!=============================================================
277
278
279	dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
280
281	call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2)
282
283	countcells = 0
284
285	c Faire rice co2 update en n-1 puis a chaque microdt, mettre a jour riceco2
286
287	c Main loop over the GCM's grid
288	DO l=1,nlay
289	DO ig=1,ngrid
290
291
292	c Get the partial pressure of co2 vapor and its saturation ratio
293	pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/44.01) * pplay(ig,l)
294	c satu = zq(ig,l,igcm_co2) / zqsat(ig,l)
295	satu = pco2 / zqsat(ig,l)
296	!=============================================================
297	! 3. Nucleation
298	!=============================================================
299
300	c call updaterccn(zq(ig,l,igcm_dust_mass),
301	c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
302
303	IF ( satu .ge. 1d0 ) THEN ! if there is condensation
304	write(,)
305	write(,) "l, pco2, satu= ",l,pco2,satu
306	c Masse_atm=mmean(ig,l)1.e-3pplay(ig,l)/rgp/zt(ig,l) !Kg par couche
307
308
309	call updaterccn(zq(ig,l,igcm_dust_mass),
310	& zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
311	write(,) "Improved, l,Rdust = ",l,rdust(ig,l)
312
313	rdust(ig,l)= zq(ig,l,igcm_dust_mass)
314	& *0.75/pi/rho_dust
315	& / zq(ig,l,igcm_dust_number)
316	rdust(ig,l)= rdust(ig,l)**(1./3.)
317	write(,) "Improved2, l,Rdust = ",l,rdust(ig,l)
318	rdust(ig,l)=max(1.e-9,rdust(ig,l))
319	rdust(ig,l)=min(2e-6,rdust(ig,l))
320	write(,) "Improved3, l,Rdust = ",l,rdust(ig,l)
321
322	c Expand the dust moments into a binned distribution
323	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)
324	No = zq(ig,l,igcm_dust_number)* tauscaling(ig)
325	write(,) "dust number, mass = ",
326	& zq(ig,l,igcm_dust_number)* tauscaling(ig),
327	& zq(ig,l,igcm_dust_mass)* tauscaling(ig)
328	c write(,) "No, Mo = ",No, Mo
329	Rn = rdust(ig,l)
330	Rn = -log(Rn)
331	Rm = Rn - 3. * sigma_iceco2*sigma_iceco2
332	n_derf = erf( (rb_cldco2(1)+Rn) *dev2)
333	m_derf = erf( (rb_cldco2(1)+Rm) *dev2)
334
335	do i = 1, nbinco2_cld
336	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
337	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
338	n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
339	m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
340	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
341	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
342	c write(,) "i, rb_cldco2(i) = ",i, rb_cldco2(i),n_aer(i)
343
344	enddo
345
346
347	sumcheck = 0
348	do i = 1, nbinco2_cld
349	sumcheck = sumcheck + n_aer(i)
350	enddo
351	sumcheck = abs(sumcheck/No - 1)
352	if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
353	print*, "WARNING, No sumcheck PROBLEM"
354	print*, "sumcheck, No",sumcheck, No
355	print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
356	print*, "Dust binned distribution", n_aer
357	STOP
358	endif
359
360	sumcheck = 0
361	do i = 1, nbinco2_cld
362	sumcheck = sumcheck + m_aer(i)
363	enddo
364	sumcheck = abs(sumcheck/Mo - 1)
365	if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld))
366	& then
367	print*, "WARNING, Mo sumcheck PROBLEM"
368	print*, "sumcheck, Mo",sumcheck, Mo
369	print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l
370	print*, "Dust binned distribution", m_aer
371	STOP
372	endif
373
374	c Expand the water ice moments into a binned distribution
375	c For now the radius grid's bound are same as for dust
376	c min=100 nm and max=10microns
377	c might need a change if rice (water) is large (but how large?) - listo
378
379	Mo = zq(ig,l,igcm_h2o_ice)* tauscaling(ig) + 1.e-30
380	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
381	Rn = rice(ig,l)
382	Rn = -log(Rn)
383	Rm = Rn - 3. * sigma_iceco2*sigma_iceco2
384	n_derf = erf( (rb_cldco2(1)+Rn) *dev2)
385	m_derf = erf( (rb_cldco2(1)+Rm) *dev2)
386	do i = 1, nbinco2_cld
387	n_aer_h2oice(i) = -0.5 * No * n_derf !! this ith previously computed
388	m_aer_h2oice(i) = -0.5 * Mo * m_derf !! this ith previously computed
389	n_derf = derf( (rb_cldco2(i+1)+Rn) *dev2)
390	m_derf = derf( (rb_cldco2(i+1)+Rm) *dev2)
391	n_aer_h2oice(i) = n_aer(i) + 0.5 * No * n_derf ! vector not really needed - temp var - listo
392	m_aer_h2oice(i) = m_aer(i) + 0.5 * Mo * m_derf ! vector not really needed - temp var
393	rad_h2oice(i) = ((m_aer_h2oice(i)/rho_ice/
394	& n_aer_h2oice(i) + vol_cld(i)))
395	& 0.75/pi*(1./3)
396	c write(,) "before nuc, i,rad_h2o(i)= ",i,rad_h2oice(i)
397	c & ,m_aer(i),n_aer(i)
398	enddo
399
400
401
402	c Get the rates of nucleation
403	call nucleaCO2(dble(pco2),zt(ig,l),dble(satu)
404	& ,n_aer,rate,n_aer_h2oice
405	& ,rad_h2oice,rateh2o)
406	! regarder rateh20, et mettre = 0 si non nul pour le moment
407	dN = 0.
408	dM = 0.
409	dNh2o = 0.
410	dMh2o = 0.
411	do i = 1, nbinco2_cld
412	!n_aer(i) = n_aer(i)/( 1. + (rate(i)+rateh2o(i))*ptimestep)
413	!m_aer(i) = m_aer(i)/( 1. + (rate(i)+rateh2o(i))*ptimestep)
414	Proba=1.0-dexp(-1.ptimesteprate(i))
415
416
417	c dNh2o = dNh2o + n_aer_h2oice(i) * rateh2o(i)
418	c & * ptimestep
419	c dMh2o = dMh2o + m_aer_h2oice(i) * rateh2o(i)
420	c & *ptimestep
421
422	dN = dN + n_aer(i) * Proba
423	dM = dM + m_aer(i) * Proba
424	c write(,) "i, dNi, dN= ",i,n_aer(i)*Proba,dN
425	enddo
426
427	c do i=1,nbinco2_cld
428	c write(,) "i n_aer m_aer = ",i,n_aer(i),m_aer(i)
429	c enddo
430	! dM masse activée (kg) et dN nb particules par kg d'air
431
432	c Update Dust particles
433	c For CO2 ice : no subtraction from dust (neither for water ice particles)
434	! zq(ig,l,igcm_dust_mass) =
435	! & zq(ig,l,igcm_dust_mass) - dM/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
436	! zq(ig,l,igcm_dust_number) =
437	! & zq(ig,l,igcm_dust_number) - dN/ tauscaling(ig) !max(tauscaling(ig),1.e-10)
438	c write(,) " nuclea dM = ",dM/tauscaling(ig),
439	c & " nuclea dN = ", dN/tauscaling(ig)
440
441	dNN= dN/tauscaling(ig)
442	dMM= dM/tauscaling(ig)
443
444	dNN=min(dNN,abs(zq0(ig,l,igcm_dust_number)))
445	dMM=min(dMM,zq0(ig,l,igcm_dust_mass))
446
447	c Update CCNs for CO2 crystals
448	! WARNING dM dMh2o, interaction nuages eau-co2 -- h20 set to 0 for now
449	zq(ig,l,igcm_ccnco2_mass) =
450	& zq(ig,l,igcm_ccnco2_mass) + dMM
451
452	zq(ig,l,igcm_ccnco2_number) =
453	& zq(ig,l,igcm_ccnco2_number) + dNN
454
455	zq(ig,l,igcm_dust_mass) =
456	& zq(ig,l,igcm_dust_mass) - dMM
457	zq(ig,l,igcm_dust_number) =
458	& zq(ig,l,igcm_dust_number) - dNN
459
460
461	c + enlever les CCN a la distri de dust
462
463	write(,) "new dust_mass, number =",
464	& zq(ig,l,igcm_dust_mass)* tauscaling(ig),
465	& zq(ig,l,igcm_dust_number)*tauscaling(ig)
466	write(,) "new ccn mass, number =",
467	& zq(ig,l,igcm_ccnco2_mass)* tauscaling(ig)
468	& ,zq(ig,l,igcm_ccnco2_number)*tauscaling(ig)
469
470	ENDIF ! of is satu >1
471	!=============================================================
472	! 4. Ice growth: scheme for radius evolution
473	!=============================================================
474
475	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
476	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
477	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
478	c IF ( zq(ig,l,igcm_ccnco2_number)*tauscaling(ig).ge. 1.0) THEN
479
480	IF (zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) .ge. threshJA)
481	& THEN ! we trigger crystal growth
482
483	call updaterccn(zq(ig,l,igcm_dust_mass),
484	& zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
485	rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass)
486	& *0.75/pi/rho_dust
487	& / zq(ig,l,igcm_ccnco2_number)
488	rdust(ig,l)= rdust(ig,l)**(1./3.)
489
490	rdust(ig,l)=max(1.e-10,rdust(ig,l))
491	rdust(ig,l)=min(2e-6,rdust(ig,l))
492	! rdust(ig,l)=1.e-7
493
494	IF (zq(ig,l,igcm_ccnco2_mass) .lt. 0. .or.
495	& zq(ig,l,igcm_ccnco2_number) .lt. 0. .or.
496	& zq(ig,l,igcm_dust_mass) .lt. 0. .or.
497	& zq(ig,l,igcm_dust_number) .lt. 0. ) THEN
498
499	write(,) "before growth CCN N,M = "
500	& ,zq(ig,l,igcm_ccnco2_number)*tauscaling(ig)
501	& ,zq(ig,l,igcm_ccnco2_mass)*tauscaling(ig)
502
503	write(,) "before growth dust number mass = ",
504	& zq(ig,l,igcm_dust_number)*tauscaling(ig),
505	& zq(ig,l,igcm_dust_mass)*tauscaling(ig)
506	STOP
507	END IF
508
509	c write(,) "reff dN = ",reff,dN
510	c reff=reff/dble(dN)
511	c if (zq(ig,l,igcm_co2_ice) .le. 1.e-20) then
512	c riceco2(ig,l)=reff
513	c endif
514
515	c write(,) "Rdust in improved = ",rdust(ig,l)
516
517	riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/
518	& (4.0rho_ice_co2pi*zq(ig,l,igcm_ccnco2_number)
519	& tauscaling(ig)) +rdust(ig,l)rdust(ig,l)
520	& rdust(ig,l) )*(1.0/3.0)
521
522	! WATCH OUT: CO2 nuclei is supposed to be dust
523	! only when deriving rhocloud (otherwise would need to keep info on water embedded in co2) - listo
524	write(,) "Rdust before growth = ",rdust(ig,l)
525	write(,) "Riceco2 before growth = ",riceco2(ig,l)
526
527	call updaterice_microCO2(zq(ig,l,igcm_co2_ice),
528	& zq(ig,l,igcm_ccnco2_mass),zq(ig,l,igcm_ccnco2_number)
529	& ,tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
530	write(,) "Riceco2 update before growth = ",riceco2(ig,l)
531
532	No = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig)
533	& + 1.e-30
534	! No nb de particules de poussieres mis à l'échelle pour donner une opacité optique
535
536	c saturation at equilibrium
537	c rice should not be too small, otherwise seq value is not valid
538	c seq = exp(2.sigco2mco2 / (rho_ice_co2rgpzt(ig,l)*
539	c & max(riceco2(ig,l),1.e-7))) !Exponant sans unité OK
540
541	ccccccc Scheme of microphys. mass growth for CO2
542
543	call massflowrateCO2(pplay(ig,l),zt(ig,l),
544	& satu,riceco2(ig,l),mmean(ig,l),Ic_rice) ! Mass transfer rate (kg/s) for a rice particle
545	! Ic_rice mass flux kg.s-1 <0 si croissance !
546	drsurdt=-1.0/(4.0piriceco2(ig,l)*
547	& riceco2(ig,l)rho_ice_co2)Ic_rice
548	dMice = No * Ic_rice * ptimestep ! Kg par kg d'air, <0 si croissance !
549	write(,) "dMicev0 in improved = " , dMice
550
551	if (dMice .gt. 0) dMice = min(dMice,zq0(ig,l,igcm_co2_ice))
552	if (dMice .lt. 0) dMice = max(dMice,-1.*zq0(ig,l,igcm_co2))
553
554
555
556
557	riceco2(ig,l)=riceco2(ig,l)+drsurdt*ptimestep
558	c write(,) "riceco2+dr/dt = ", riceco2(ig,l)
559
560	write(,) "dMice in improved = " , dMice
561
562
563	zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice)
564	& -dMice
565	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)+dMice
566	c write(,) "zq co2 ice = ", zq(ig,l,igcm_co2_ice)
567	countcells = countcells + 1
568
569	riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/
570	& (4.0rho_ice_co2pi*zq(ig,l,igcm_ccnco2_number)
571	& tauscaling(ig)) +rdust(ig,l)rdust(ig,l)
572	& rdust(ig,l) )*(1.0/3.0)
573	write(,) "new riceco2 = ",riceco2(ig,l)
574
575	call updaterice_microCO2(zq(ig,l,igcm_co2_ice),
576	& zq(ig,l,igcm_ccnco2_mass),zq(ig,l,igcm_ccnco2_number)
577	& ,tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
578	write(,) "new riceco2 updaterad= ",riceco2(ig,l)
579
580	! latent heat release
581
582	l0=595594.
583	l1=903.111
584	l2=-11.5959
585	l3=0.0528288
586	l4=-0.000103183
587
588	lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 +
589	& l3 * zt(ig,l)*3 + l4 zt(ig,l)**4 !J.kg-1
590	c write(,) "CPP= ",cpp ! = 744.5
591
592	pdtcloudco2(ig,l)= dMicelw/cpp/ptimestep ! kg par couche J par kg /J par K / s = K par seconde
593
594	c write(,) "pdtcloudco2 = ",pdtcloudco2(ig,l)
595
596	c Peut etre -1*dMice?
597
598
599	!deltaT par condens/subli. qui remplace le dT du CO2 du newcondens pré-constantino
600	!PDT should be in K/s.
601	!=============================================================
602	! 5. Dust cores released, tendancies, latent heat, etc ...
603	!=============================================================
604
605	c If all the ice particles sublimate, all the condensation
606	c nuclei are released:
607
608	c !!! with CO2 ice nuclei: dust/H2O nuclei are not released because
609	c they were not subtracted to dust_number
610	c Their counter is just set to "0".
611	c (see end of section 3.) : On peut les enlever à dust
612
613	c interaction ho-co2 ici, dans la mise a jour des traceurs WARNING reflechir
614
615
616	ENDIF !of if Nccn>1
617	c if (riceco2(ig,l) .lt. 1.e-9) then
618
619	if (zq(ig,l,igcm_co2_ice).le.1.e-20 .or.
620	& riceco2(ig,l) .lt. 1.e-9 .or. riceco2(ig,l)
621	& .ge. 4.999e-4) then
622	! Reverser les ccn libérés dans les h2o ou dust?
623
624	c c ICI: Co2 ice devient vapeur
625	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)
626	& + zq(ig,l,igcm_co2_ice)
627
628	zq(ig,l,igcm_dust_mass) =
629	& zq(ig,l,igcm_dust_mass)
630	& + zq(ig,l,igcm_ccnco2_mass)
631	zq(ig,l,igcm_dust_number)
632	& = zq(ig,l,igcm_dust_number)
633	& + zq(ig,l,igcm_ccnco2_number)
634	c c CCNs
635	zq(ig,l,igcm_ccnco2_mass) = 0.
636	zq(ig,l,igcm_ccnco2_number) =0.
637	zq(ig,l,igcm_co2_ice) = 0.
638	riceco2(ig,l)=0.
639	endif
640
641	c write(,) "zq co2 end imp= ", zq(i g,l,igcm_co2_ice),satu
642
643
644
645	ENDDO ! of ig loop
646	ENDDO ! of nlayer loop
647
648
649	! Get cloud tendencies
650	pdqcloudco2(1:ngrid,1:nlay,igcm_co2) =
651	& (zq(1:ngrid,1:nlay,igcm_co2) -
652	& zq0(1:ngrid,1:nlay,igcm_co2))/ptimestep
653
654	pdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) =
655	& (zq(1:ngrid,1:nlay,igcm_co2_ice) -
656	& zq0(1:ngrid,1:nlay,igcm_co2_ice))/ptimestep
657	c do l=1,nlay
658	c write(,) "end imp",pdqcloudco2(1,l,igcm_co2),
659	c & pdqcloudco2(1,l,igcm_co2_ice)
660	c enddo
661	pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) =
662	& (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) -
663	& zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/ptimestep
664
665	pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) =
666	& (zq(1:ngrid,1:nlay,igcm_ccnco2_number) -
667	& zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/ptimestep
668
669
670	if (scavenging) then
671
672	pdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) =
673	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
674	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep
675	pdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) =
676	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
677	& zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
678	endif
679
680	c call WRITEdiagfi(ngrid,"Improvedb","co2 traceur","kg/kg",1,
681	c & zq0(1,:,igcm_co2_ice) )
682
683	c call WRITEdiagfi(ngrid,"Improveda","co2 traceur","kg/kg",1,
684	c & zq(1,:,igcm_co2_ice) )
685
686
687
688
689	c call WRITEdiagfi(ngrid,"satuco2","satu co2 in improved","kg/kg",1,
690	c & satu)
691
692
693	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
694	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
695	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
696	IF (test_flag) then
697
698	! error2d(:) = 0.
699	DO l=1,nlay
700	DO ig=1,ngrid
701	! error2d(ig) = max(abs(error_out(ig,l)),error2d(ig))
702	satubf(ig,l) = zq0(ig,l,igcm_co2)/zqsat(ig,l) ! att. if zqsat=mmr or psat
703	satuaf(ig,l) = zq(ig,l,igcm_co2)/zqsat(ig,l)
704	ENDDO
705	ENDDO
706
707	write(,) 'count is ',countcells, ' i.e. ',
708	& countcells100/(nlayngrid), '% for microphys computation'
709
710	! IF (ngrid.ne.1) THEN ! 3D
711	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",3,
712	! & satu_out)
713	! call WRITEDIAGFI(ngrid,"dM","ccn variation","kg/kg",3,
714	! & dM_out)
715	! call WRITEDIAGFI(ngrid,"dN","ccn variation","#",3,
716	! & dN_out)
717	! call WRITEDIAGFI(ngrid,"error","dichotomy max error","%",2,
718	! & error2d)
719	! call WRITEDIAGFI(ngrid,"zqsat","zqsat","kg",3,
720	! & zqsat)
721	! ENDIF
722
723	! IF (ngrid.eq.1) THEN ! 1D
724	! call WRITEDIAGFI(ngrid,"error","incertitude sur glace","%",1,
725	! & error_out)
726	! call WRITEdiagfi(ngrid,"resist","resistance","s/m2",1,
727	! & res_out)
728	call WRITEdiagfi(ngrid,"satu_bf","satu before","kg/kg",1,
729	& satubf)
730	call WRITEdiagfi(ngrid,"satu_af","satu after","kg/kg",1,
731	& satuaf)
732	call WRITEdiagfi(ngrid,"vapbf","CO2vap before","kg/kg",1,
733	& zq0(1,1,igcm_co2))
734	call WRITEdiagfi(ngrid,"vapaf","CO2vap after","kg/kg",1,
735	& zq(1,1,igcm_co2))
736	call WRITEdiagfi(ngrid,"icebf","CO2ice before","kg/kg",1,
737	& zq0(1,1,igcm_co2_ice))
738	call WRITEdiagfi(ngrid,"iceaf","CO2ice after","kg/kg",1,
739	& zq(1,1,igcm_co2_ice))
740	call WRITEdiagfi(ngrid,"ccnbf","ccn before","/kg",1,
741	& zq0(1,1,igcm_ccnco2_number))
742	call WRITEdiagfi(ngrid,"ccnaf","ccn after","/kg",1,
743	& zq(1,1,igcm_ccnco2_number))
744	c call WRITEDIAGFI(ngrid,"growthrate","growth rate","m^2/s",1,
745	c & gr_out)
746	c call WRITEDIAGFI(ngrid,"nuclearate","nucleation rate","",1,
747	c & rate_out)
748	c call WRITEDIAGFI(ngrid,"dM","ccn variation","kg",1,
749	c & dM_out)
750	c call WRITEDIAGFI(ngrid,"dN","ccn variation","#",1,
751	c & dN_out)
752	c call WRITEdiagfi(ngrid,"zqsat","p vap sat","kg/kg",1,
753	c & zqsat)
754	! call WRITEDIAGFI(ngrid,"satu","ratio saturation","",1,
755	! & satu_out)
756	! call WRITEDIAGFI(ngrid,"rdust_sca","rdust","m",1,
757	! & rdust)
758	! call WRITEDIAGFI(ngrid,"rsedcloud","rsedcloud","m",1,
759	! & rsedcloud)
760	! call WRITEDIAGFI(ngrid,"rhocloud","rhocloud","kg.m-3",1,
761	! & rhocloud)
762	! ENDIF
763
764	ENDIF ! endif test_flag
765	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
766	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
767	!!!!!!!!!!!!!! TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS TESTS OUTPUTS
768
769	return
770	end
771
772
773
774	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
775	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
776	c The so -called "phi" function is such as phi(r) - phi(r0) = t - t0
777	c It is an analytical solution to the ice radius growth equation,
778	c with the approximation of a constant 'reduced' cunningham correction factor
779	c (lambda in growthrate.F) taken at radius req instead of rice
780	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
781	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
782
783	c subroutine phi(rice,req,coeff1,coeff2,time)
784	c
785	c implicit none
786	c
787	c ! inputs
788	c real rice ! ice radius
789	c real req ! ice radius at equilibirum
790	c real coeff1 ! coeff for the log
791	c real coeff2 ! coeff for the arctan
792	c
793	c ! output
794	c real time
795	c
796	c !local
797	c real var
798	c
799	c ! 1.73205 is sqrt(3)
800	c
801	c var = max(
802	c & abs(rice-req) / sqrt(ricerice + ricereq + req*req),1e-30)
803	c
804	c time =
805	c & coeff1 *
806	c & log( var )
807	c & + coeff2 * 1.73205 *
808	c & atan( (2rice+req) / (1.73205req) )
809	c
810	c return
811	c end
812
813
814

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

Download in other formats:

Original Format