Context Navigation

improvedCO2clouds_mod.F @ 1944

Last change on this file since 1944 was 1939, checked in by aslmd, 6 years ago
rename improvedCO2clouds in improvedCO2clouds_mod since it is a module
Property svn:executable set to ``*
File size: 29.6 KB

Line
1	MODULE improvedCO2clouds_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	subroutine improvedCO2clouds(ngrid,nlay,microtimestep,
8	& pplay,pplev,pteff,sum_subpdt,
9	& pqeff,sum_subpdq,subpdqcloudco2,subpdtcloudco2,
10	& nq,tauscaling,
11	& mem_Mccn_co2,mem_Mh2o_co2,mem_Nccn_co2)
12	USE comcstfi_h, only: pi, g, cpp
13	USE updaterad, only: updaterice_micro, updaterice_microco2
14	use tracer_mod, only: igcm_dust_mass, igcm_dust_number, rho_dust,
15	& igcm_h2o_ice, igcm_ccn_mass,
16	& igcm_ccn_number, nuice_sed,
17	& igcm_co2, igcm_co2_ice, igcm_ccnco2_mass,
18	& igcm_ccnco2_number, nuiceco2_sed,
19	& rho_ice_co2
20	use conc_mod, only: mmean
21	use datafile_mod, only: datadir
22
23	implicit none
24
25	c------------------------------------------------------------------
26	c This routine is used to form CO2 clouds when a parcel of the GCM is
27	c saturated. It includes the ability to have supersaturation, a
28	c computation of the nucleation rates, growthrates and the
29	c scavenging of dust particles by clouds.
30	c It is worth noting that the amount of dust is computed using the
31	c dust optical depth computed in aeropacity.F. That's why
32	c the variable called "tauscaling" is used to convert
33	c pq(dust_mass) and pq(dust_number), which are relative
34	c quantities, to absolute and realistic quantities stored in zq.
35	c This has to be done to convert the inputs into absolute
36	c values, but also to convert the outputs back into relative
37	c values which are then used by the sedimentation and advection
38	c schemes.
39	c CO2 ice particles can nucleate on both dust and on water ice particles
40	c When CO2 ice is deposited onto a water ice particles, the particle is
41	c removed from the water tracers.
42	c Memory of the origin of the co2 particles is kept and thus the
43	c water cycle shouldn't be modified by this.
44	c WARNING: no sedimentation of the water ice origin is performed
45	c in the microphysical timestep in co2cloud.F.
46
47	c Authors of the water ice clouds microphysics
48	c J.-B. Madeleine, based on the work by Franck Montmessin
49	c (October 2011)
50	c T. Navarro, debug,correction, new scheme (October-April 2011)
51	c A. Spiga, optimization (February 2012)
52	c Adaptation for CO2 clouds by Joachim Audouard (09/16), based on the work
53	c of Constantino Listowski
54	c There is an energy limit to how much co2 can sublimate/condensate. It is
55	c defined by the difference of the GCM temperature with the co2 condensation
56	c temperature.
57	c Warning:
58	c If meteoritic particles are activated and turn into co2 ice particles,
59	c then they will be reversed in the dust tracers if the cloud sublimates
60
61	c------------------------------------------------------------------
62	include "callkeys.h"
63	include "microphys.h"
64	c------------------------------------------------------------------
65	c Arguments:
66
67	INTEGER,INTENT(in) :: ngrid,nlay
68	integer,intent(in) :: nq ! number of tracers
69	REAL,INTENT(in) :: microtimestep ! physics time step (s)
70	REAL,INTENT(in) :: pplay(ngrid,nlay) ! mid-layer pressure (Pa)
71	REAL,INTENT(in) :: pplev(ngrid,nlay+1) ! inter-layer pressure (Pa)
72	REAL,INTENT(in) :: pteff(ngrid,nlay) ! temperature at the middle of the
73	! layers (K)
74	REAL,INTENT(in) :: sum_subpdt(ngrid,nlay) ! tendency on temperature from
75	! previous physical parametrizations
76	REAL,INTENT(in) :: pqeff(ngrid,nlay,nq) ! tracers (kg/kg)
77	REAL,INTENT(in) :: sum_subpdq(ngrid,nlay,nq) ! tendencies on tracers
78	! before condensation (kg/kg.s-1)
79	REAL,INTENT(in) :: tauscaling(ngrid) ! Convertion factor for qdust and Ndust
80	c Outputs:
81	REAL,INTENT(out) :: subpdqcloudco2(ngrid,nlay,nq) ! tendency on tracers
82	! due to CO2 condensation (kg/kg.s-1)
83	! condensation si igcm_co2_ice
84	REAL,INTENT(out) :: subpdtcloudco2(ngrid,nlay) ! tendency on temperature due
85	! to latent heat
86
87	c------------------------------------------------------------------
88	c Local variables:
89	LOGICAL,SAVE :: firstcall=.true.
90	REAL*8 derf ! Error function
91	INTEGER ig,l,i
92
93	real masse (ngrid,nlay) ! Layer mass (kg.m-2)
94	REAL rice(ngrid,nlay) ! Water Ice mass mean radius (m)
95	! used for nucleation of CO2 on ice-coated ccns
96	REAL rccnh2o(ngrid,nlay) ! Water Ice mass mean radius (m)
97	REAL zq(ngrid,nlay,nq) ! local value of tracers
98	REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
99	REAL zt(ngrid,nlay) ! local value of temperature
100	REAL zqsat(ngrid,nlay) ! saturation vapor pressure for CO2
101	real tcond(ngrid,nlay)
102	real zqco2(ngrid,nlay)
103	REAL lw !Latent heat of sublimation (J.kg-1)
104	REAL,save :: l0,l1,l2,l3,l4
105	DOUBLE PRECISION dMice ! mass of condensed ice
106	DOUBLE PRECISION sumcheck
107	DOUBLE PRECISION facteurmax!for energy limit on mass growth
108	DOUBLE PRECISION pco2,psat ! Co2 vapor partial pressure (Pa)
109	DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice
110	DOUBLE PRECISION Mo,No,No_dust,Mo_dust
111	DOUBLE PRECISION Rn, Rm, dev2,dev3, n_derf, m_derf
112	DOUBLE PRECISION mem_Mccn_co2(ngrid,nlay) ! Memory of CCN mass of H2O and dust used by CO2
113	DOUBLE PRECISION mem_Mh2o_co2(ngrid,nlay) ! Memory of H2O mass integred into CO2 crystal
114	DOUBLE PRECISION mem_Nccn_co2(ngrid,nlay) ! Memory of CCN number of H2O and dust used by CO2
115	DOUBLE PRECISION interm1,interm2,interm3
116
117	! Radius used by the microphysical scheme (m)
118	DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin
119	DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin
120	DOUBLE PRECISION m_aer_h2oice2(nbinco2_cld) ! mass mixing ratio of particle/each size bin
121
122	DOUBLE PRECISION n_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
123	DOUBLE PRECISION m_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
124	DOUBLE PRECISION rad_h2oice(nbinco2_cld)
125
126	c REAL*8 sigco2 ! Co2-ice/air surface tension (N.m)
127	c EXTERNAL sigco2
128
129	DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM,dNNh2o,dMMh2o
130	DOUBLE PRECISION dMh2o_ice,dMh2o_ccn
131
132	DOUBLE PRECISION rate(nbinco2_cld) ! nucleation rate
133	DOUBLE PRECISION rateh2o(nbinco2_cld) ! nucleation rate
134	REAL seq
135	DOUBLE PRECISION rho_ice_co2T(ngrid,nlay)
136	DOUBLE PRECISION riceco2(ngrid,nlay) ! CO2Ice mean radius (m)
137	REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
138
139	REAL rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
140	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
141
142	c REAL res ! Resistance growth
143	DOUBLE PRECISION Ic_rice ! Mass transfer rate CO2 ice crystal
144	DOUBLE PRECISION ratioh2o_ccn
145	DOUBLE PRECISION vo2co2
146
147	c Parameters of the size discretization used by the microphysical scheme
148	DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m)
149	DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6 ! Maximum radius (m)
150	DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-10 ! Minimum bounary radius (m)
151	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4 ! Maximum boundary radius (m)
152	DOUBLE PRECISION vrat_cld ! Volume ratio
153	DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
154	SAVE rb_cldco2
155	DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m)
156	DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3)
157
158	DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o
159	REAL sigma_iceco2 ! Variance of the co2 ice and CCN distributions
160	SAVE sigma_iceco2
161	REAL sigma_ice ! Variance of the h2o ice and CCN distributions
162	SAVE sigma_ice
163	DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2
164
165	! Variables for the meteoritic flux:
166	integer,parameter :: nbin_meteor=100
167	integer,parameter :: nlev_meteor=130
168	double precision meteor_ccn(ngrid,nlay,100) !100=nbinco2_cld !!!
169	double precision,save :: meteor(130,100)
170	double precision mtemp(100),pression_meteor(130)
171	logical file_ok
172	integer read_ok
173	integer nelem,lebon1,lebon2
174	double precision :: ltemp1(130),ltemp2(130)
175	integer ibin,j
176	integer,parameter :: uMeteor=666
177
178	IF (firstcall) THEN
179	!=============================================================
180	! 0. Definition of the size grid
181	!=============================================================
182	c rad_cldco2 is the primary radius grid used for microphysics computation.
183	c The grid spacing is computed assuming a constant volume ratio
184	c between two consecutive bins; i.e. vrat_cld.
185	c vrat_cld is determined from the boundary values of the size grid:
186	c rmin_cld and rmax_cld.
187	c The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
188	c dr_cld is the width of each rad_cldco2 bin.
189
190	vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
191	vrat_cld = dexp(vrat_cld)
192	rb_cldco2(1) = rbmin_cld
193	rad_cldco2(1) = rmin_cld
194	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
195	do i=1,nbinco2_cld-1
196	rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.)
197	vol_cld(i+1) = vol_cld(i) * vrat_cld
198	enddo
199	do i=1,nbinco2_cld
200	rb_cldco2(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
201	& rad_cldco2(i)
202	dr_cld(i) = rb_cldco2(i+1) - rb_cldco2(i)
203	enddo
204	rb_cldco2(nbinco2_cld+1) = rbmax_cld
205	dr_cld(nbinco2_cld) = rb_cldco2(nbinco2_cld+1) -
206	& rb_cldco2(nbinco2_cld)
207	print*, ' '
208	print*,'Microphysics co2: size bin information:'
209	print*,'i,rb_cldco2(i), rad_cldco2(i),dr_cld(i)'
210	print*,'-----------------------------------'
211	do i=1,nbinco2_cld
212	write(*,'(i3,3x,3(e13.6,4x))') i,rb_cldco2(i), rad_cldco2(i),
213	& dr_cld(i)
214	enddo
215	write(*,'(i3,3x,e13.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1)
216	print*,'-----------------------------------'
217	do i=1,nbinco2_cld+1
218	rb_cldco2(i) = dlog(rb_cldco2(i)) !! we save that so that it is not computed
219	!! at each timestep and gridpoint
220	enddo
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	vo1co2 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
229
230	c Variance of the ice and CCN distributions
231	sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
232	sigma_ice = sqrt(log(1.+nuice_sed))
233
234
235	write(,) 'Variance of ice & CCN CO2 distribs :', sigma_iceco2
236	write(,) 'nuice for co2 ice sedimentation:', nuiceco2_sed
237	write(,) 'Volume of a co2 molecule:', vo1co2
238
239	if (co2useh2o) then
240	write(,)
241	write(,) "co2useh2o=.true. in callphys.def"
242	write(,) "This means water ice particles can "
243	write(,) "serve as CCN for CO2 microphysics"
244	endif
245
246	if (meteo_flux) then
247	write(,)
248	write(,) "meteo_flux=.true. in callphys.def"
249	write(,) "meteoritic dust particles are available"
250	write(,) "for co2 ice nucleation! "
251	write(,) "Flux given by J. Plane (pressions,size bins)"
252	! Initialisation of the flux: it is constant and is it saved
253	!We must interpolate the table to the GCM pressures
254	INQUIRE(FILE=TRIM(datadir)//
255	& '/Meteo_flux_Plane.dat', EXIST=file_ok)
256	IF (.not. file_ok) THEN
257	write(,) 'file Meteo_flux_Plane.dat should be in '
258	& ,trim(datadir)
259	STOP
260	endif
261	!used Variables
262	open(unit=uMeteor,file=trim(datadir)//
263	& '/Meteo_flux_Plane.dat'
264	& ,FORM='formatted')
265	!13000 records (130 pressions x 100 bin sizes)
266	read(uMeteor,*) !skip 1 line
267	do i=1,130
268	read(uMeteor,*,iostat=read_ok) pression_meteor(i)
269	if (read_ok==0) then
270	write(,) pression_meteor(i)
271	else
272	write(,) 'Error reading Meteo_flux_Plane.dat'
273	call abort_physic("CO2clouds",
274	& "Error reading Meteo_flux_Plane.dat",1)
275	endif
276	enddo
277	read(uMeteor,*) !skip 1 line
278	do i=1,130
279	do j=1,100 ! les mêmes 100 bins size que la distri nuclea: on touche pas
280	read(uMeteor,'(F12.6)',iostat=read_ok) meteor(i,j)
281	if (read_ok/=0) then
282	write(,) 'Error reading Meteo_flux_Plane.dat'
283	call abort_physic("CO2clouds",
284	& "Error reading Meteo_flux_Plane.dat",1)
285	endif
286	enddo
287	!On doit maintenant réinterpoler le tableau(130,100) sur les pressions du GCM (nlay,100)
288	enddo
289	close(uMeteor)
290	write(,) "File meteo_flux read, end of firstcall in co2 micro"
291	endif !of if meteo_flux
292
293	! Parameter values for Latent heat computation
294	l0=595594d0
295	l1=903.111d0
296	l2=-11.5959d0
297	l3=0.0528288d0
298	l4=-0.000103183d0
299	firstcall=.false.
300	END IF
301	!=============================================================
302	! 1. Initialisation
303	!=============================================================
304
305	meteor_ccn(:,:,:)=0.
306	rice(:,:) = 1.e-8
307	riceco2(:,:) = 1.e-11
308
309	c Initialize the tendencies
310	subpdqcloudco2(1:ngrid,1:nlay,1:nq)=0.
311	subpdtcloudco2(1:ngrid,1:nlay)=0.
312
313	c pteff temperature layer; sum_subpdt dT.s-1
314	c pqeff traceur kg/kg; sum_subpdq tendance idem .s-1
315	zt(1:ngrid,1:nlay) =
316	& pteff(1:ngrid,1:nlay) +
317	& sum_subpdt(1:ngrid,1:nlay) * microtimestep
318	zq(1:ngrid,1:nlay,1:nq) =
319	& pqeff(1:ngrid,1:nlay,1:nq) +
320	& sum_subpdq(1:ngrid,1:nlay,1:nq) * microtimestep
321	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
322	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
323
324	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
325	!=============================================================
326	! 2. Compute saturation
327	!=============================================================
328	dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
329	dev3 = 1. / ( sqrt(2.) * sigma_ice )
330	call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2)
331	zqco2=zq(:,:,igcm_co2)+zq(:,:,igcm_co2_ice)
332	CALL tcondco2(ngrid,nlay,pplay,zqco2,tcond)
333	!=============================================================
334	! 3. Bonus: additional meteoritic particles for nucleation
335	!=============================================================
336	if (meteo_flux) then
337	!pression_meteo(130)
338	!pplev(ngrid,nlay+1)
339	!meteo(130,100)
340	!resultat: meteo_ccn(ngrid,nlay,100)
341	do l=1,nlay
342	do ig=1,ngrid
343	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
344	ltemp1=abs(pression_meteor(:)-pplev(ig,l))
345	ltemp2=abs(pression_meteor(:)-pplev(ig,l+1))
346	lebon1=minloc(ltemp1,DIM=1)
347	lebon2=minloc(ltemp2,DIM=1)
348	nelem=lebon2-lebon1+1.
349	mtemp(:)=0d0 !mtemp(100) : valeurs pour les 100bins
350	do ibin=1,100
351	mtemp(ibin)=sum(meteor(lebon1:lebon2,ibin))
352	enddo
353	meteor_ccn(ig,l,:)=mtemp(:)/nelem/masse(ig,l) !Par kg air
354	csi par m carre, x epaisseur/masse pour par kg/air
355	!write(,) "masse air ig l=",masse(ig,l)
356	!check original unit with J. Plane
357	enddo
358	enddo
359	endif
360	c ------------------------------------------------------------------------
361	c --------- Actual microphysics : Main loop over the GCM's grid ---------
362	c ------------------------------------------------------------------------
363	DO l=1,nlay
364	DO ig=1,ngrid
365	c Get the partial pressure of co2 vapor and its saturation ratio
366	pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/44.01) * pplay(ig,l)
367	satu = pco2 / zqsat(ig,l)
368
369	rho_ice_co2T(ig,l)=1000.(1.72391-2.53e-4zt(ig,l)
370	& -2.87e-6zt(ig,l)zt(ig,l)) !T-dependant CO2 ice density
371	vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l))
372	rho_ice_co2=rho_ice_co2T(ig,l)
373
374	!=============================================================
375	!4. Nucleation
376	!=============================================================
377	IF ( satu .ge. 1 ) THEN ! if there is condensation
378
379	c We do rdust computation "by hand" because we don't want to
380	c change the mininumum rdust value in updaterad... It would
381	c mess up various part of the GCM !
382
383	rdust(ig,l)= zq(ig,l,igcm_dust_mass)
384	& 0.75/(pirho_dust
385	& * zq(ig,l,igcm_dust_number))
386	rdust(ig,l)= rdust(ig,l)**(1./3.)
387	if (zq(ig,l,igcm_dust_mass)*tauscaling(ig) .le. 1.e-20
388	& .or. zq(ig,l,igcm_dust_number)*tauscaling(ig) .le.1
389	& .or. rdust(ig,l) .le. 1.e-9) then
390	rdust(ig,l)=1.e-9
391	endif
392	rdust(ig,l)=min(5.e-4,rdust(ig,l))
393
394	c Expand the dust moments into a binned distribution
395
396	n_aer(:)=0.
397	m_aer(:)=0.
398
399	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)+1.e-30
400	No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30
401
402	No_dust=No
403	Mo_dust=Mo
404
405	Rn = rdust(ig,l)
406	Rn = -dlog(Rn)
407	Rm = Rn - 3. * sigma_iceco2*sigma_iceco2
408	n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
409	m_derf = derf( (rb_cldco2(1)+Rm) *dev2)
410
411	do i = 1, nbinco2_cld
412	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
413	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
414	n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
415	m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
416	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
417	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
418	enddo
419
420	c Ajout meteor_ccn particles aux particules de poussière background
421	if (meteo_flux) then
422	do i = 1, nbinco2_cld
423	n_aer(i) = n_aer(i) + meteor_ccn(ig,l,i)
424	m_aer(i) = m_aer(i) + 4./3.pirho_dust
425	& meteor_ccn(ig,l,i)rad_cldco2(i)*rad_cldco2(i)
426	& *rad_cldco2(i)
427	enddo
428	endif
429
430	c Same but with h2o particles as CCN only if co2useh2o=.true.
431
432	n_aer_h2oice(:)=0.
433	m_aer_h2oice(:)=0.
434
435	if (co2useh2o) then
436	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
437	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
438	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
439	Mo = zq(ig,l,igcm_h2o_ice) +
440	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)+1.e-30
441	! Total mass of H20 crystals,CCN included
442	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
443	Rn = rice(ig,l)
444	Rn = -dlog(Rn)
445	Rm = Rn - 3. * sigma_ice*sigma_ice
446	n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
447	m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
448	do i = 1, nbinco2_cld
449	n_aer_h2oice(i) = -0.5 * No * n_derf
450	m_aer_h2oice(i) = -0.5 * Mo * m_derf
451	n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
452	m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
453	n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf
454	m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf
455	rad_h2oice(i) = rad_cldco2(i)
456	enddo
457	endif
458
459
460	! Call to nucleation routine
461	call nucleaCO2(dble(pco2),zt(ig,l),dble(satu)
462	& ,n_aer,rate,n_aer_h2oice
463	& ,rad_h2oice,rateh2o,vo2co2)
464	dN = 0.
465	dM = 0.
466	dNh2o = 0.
467	dMh2o = 0.
468	do i = 1, nbinco2_cld
469	Proba =1.0-dexp(-1.microtimesteprate(i))
470	dN = dN + n_aer(i) * Proba
471	dM = dM + m_aer(i) * Proba
472	enddo
473	if (co2useh2o) then
474	do i = 1, nbinco2_cld
475	Probah2o = 1.0-dexp(-1.microtimesteprateh2o(i))
476	dNh2o = dNh2o + n_aer_h2oice(i) * Probah2o
477	dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o
478	enddo
479	endif
480
481	! dM masse activée (kg) et dN nb particules par kg d'air
482	! Now increment CCN tracers and update dust tracers
483	dNN= dN/tauscaling(ig)
484	dMM= dM/tauscaling(ig)
485	dNN=min(dNN,zq(ig,l,igcm_dust_number))
486	dMM=min(dMM,zq(ig,l,igcm_dust_mass))
487	zq(ig,l,igcm_ccnco2_mass) =
488	& zq(ig,l,igcm_ccnco2_mass) + dMM
489	zq(ig,l,igcm_ccnco2_number) =
490	& zq(ig,l,igcm_ccnco2_number) + dNN
491	zq(ig,l,igcm_dust_mass)= zq(ig,l,igcm_dust_mass)-dMM
492	zq(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)-dNN
493
494
495	c Update CCN for CO2 nucleating on H2O CCN :
496	! Warning: must keep memory of it
497	if (co2useh2o) then
498	dNNh2o=dNh2o/tauscaling(ig)
499	dNNh2o=min(dNNh2o,zq(ig,l,igcm_ccn_number))
500	ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice)
501	& +zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
502	dMh2o_ice=dMh2ozq(ig,l,igcm_h2o_ice)ratioh2o_ccn
503	dMh2o_ccn=dMh2ozq(ig,l,igcm_ccn_mass)
504	& tauscaling(ig)*ratioh2o_ccn
505	dMh2o_ccn=dMh2o_ccn/tauscaling(ig)
506	dMh2o_ccn=min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
507	dMh2o_ice=min(dMh2o_ice,zq(ig,l,igcm_h2o_ice))
508	zq(ig,l,igcm_ccnco2_mass) =
509	& zq(ig,l,igcm_ccnco2_mass) + dMh2o_ice+dMh2o_ccn
510	zq(ig,l,igcm_ccnco2_number) =
511	& zq(ig,l,igcm_ccnco2_number) + dNNh2o
512	zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)-dNNh2o
513	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)-dMh2o_ice
514	zq(ig,l,igcm_ccn_mass)= zq(ig,l,igcm_ccn_mass)-dMh2o_ccn
515	mem_Mh2o_co2(ig,l)=mem_Mh2o_co2(ig,l)+dMh2o_ice
516	mem_Mccn_co2(ig,l)=mem_Mccn_co2(ig,l)+dMh2o_ccn
517	mem_Nccn_co2(ig,l)=mem_Nccn_co2(ig,l)+dNNh2o
518	endif ! of if co2useh2o
519	ENDIF ! of is satu >1
520
521	!=============================================================
522	! 5. Ice growth: scheme for radius evolution
523	!=============================================================
524
525	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
526	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
527	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
528	No = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig)+1.e-30
529	IF (No .ge. 1)THEN ! we trigger crystal growth
530	interm1 = DBLE(zq(ig,l,igcm_co2_ice))
531	interm2 = DBLE(zq(ig,l,igcm_ccnco2_mass))
532	interm3 = DBLE(zq(ig,l,igcm_ccnco2_number))
533	call updaterice_microco2(interm1,
534	& interm2,interm3,
535	& tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
536
537	Ic_rice=0.
538	lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 +
539	& l3 * zt(ig,l)*3 + l4 zt(ig,l)**4 !J.kg-1
540	facteurmax=abs(Tcond(ig,l)-zt(ig,l))*cpp/lw
541	!specific heat of co2 ice = 1000 J.kg-1.K-1
542	!specific heat of atm cpp = 744.5 J.kg-1.K-1
543
544	c call scheme of microphys. mass growth for CO2
545	call massflowrateCO2(pplay(ig,l),zt(ig,l),
546	& satu,riceco2(ig,l),mmean(ig,l),Ic_rice)
547	c Ic_rice Mass transfer rate (kg/s) for a rice particle >0 si croissance !
548
549	if (isnan(Ic_rice) .or. Ic_rice .eq. 0.) then
550	Ic_rice=0.
551	subpdtcloudco2(ig,l)=-sum_subpdt(ig,l)
552	dMice=0
553
554	else
555	dMice=zq(ig,l,igcm_ccnco2_number)Ic_ricemicrotimestep
556	& *tauscaling(ig) ! Kg par kg d'air, >0 si croissance !
557	!kg.s-1 par particule * nb particule par kg air*s
558	! = kg par kg air
559
560	dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
561	dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
562	! facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
563	! latent heat release >0 if growth i.e. if dMice >0
564	subpdtcloudco2(ig,l)=dMice*lw/cpp/microtimestep
565	! kgco2/kgair* J/kgco2 * 1/(J.kgair-1.K-1)/s= K par seconde
566	!Now update tracers
567	zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice)+dMice
568	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)-dMice
569	endif
570
571	!=============================================================
572	! 6. Dust cores releasing if no more co2 ice :
573	!=============================================================
574
575	if (zq(ig,l,igcm_co2_ice).le. 1.e-25)THEN
576	! On sublime tout
577	if (co2useh2o) then
578	if (mem_Mccn_co2(ig,l) .gt. 0) then
579	zq(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
580	& +mem_Mccn_co2(ig,l)
581	endif
582	if (mem_Mh2o_co2(ig,l) .gt. 0) then
583	zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
584	& +mem_Mh2o_co2(ig,l)
585	endif
586
587	if (mem_Nccn_co2(ig,l) .gt. 0) then
588	zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
589	& +mem_Nccn_co2(ig,l)
590	endif
591	endif
592	zq(ig,l,igcm_dust_mass) =
593	& zq(ig,l,igcm_dust_mass)
594	& + zq(ig,l,igcm_ccnco2_mass)-
595	& (mem_Mh2o_co2(ig,l)+mem_Mccn_co2(ig,l))
596	zq(ig,l,igcm_dust_number) =
597	& zq(ig,l,igcm_dust_number)
598	& + zq(ig,l,igcm_ccnco2_number)
599	& -mem_Nccn_co2(ig,l)
600
601	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)
602	& + zq(ig,l,igcm_co2_ice)
603
604	zq(ig,l,igcm_ccnco2_mass)=0.
605	zq(ig,l,igcm_co2_ice)=0.
606	zq(ig,l,igcm_ccnco2_number)=0.
607	mem_Nccn_co2(ig,l)=0.
608	mem_Mh2o_co2(ig,l)=0.
609	mem_Mccn_co2(ig,l)=0.
610	riceco2(ig,l)=0.
611
612	endif !of if co2_ice <1e-25
613
614	ENDIF ! of if NCCN > 1
615	ENDDO ! of ig loop
616	ENDDO ! of nlayer loop
617
618	!=============================================================
619	! 7. END: get cloud tendencies
620	!=============================================================
621
622	! Get cloud tendencies
623	subpdqcloudco2(1:ngrid,1:nlay,igcm_co2) =
624	& (zq(1:ngrid,1:nlay,igcm_co2) -
625	& zq0(1:ngrid,1:nlay,igcm_co2))/microtimestep
626	subpdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) =
627	& (zq(1:ngrid,1:nlay,igcm_co2_ice) -
628	& zq0(1:ngrid,1:nlay,igcm_co2_ice))/microtimestep
629
630	if (co2useh2o) then
631	subpdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
632	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
633	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/microtimestep
634	subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
635	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
636	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/microtimestep
637	subpdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
638	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
639	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/microtimestep
640	endif
641
642	subpdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) =
643	& (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) -
644	& zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/microtimestep
645	subpdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) =
646	& (zq(1:ngrid,1:nlay,igcm_ccnco2_number) -
647	& zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/microtimestep
648	subpdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) =
649	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
650	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/microtimestep
651	subpdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) =
652	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
653	& zq0(1:ngrid,1:nlay,igcm_dust_number))/microtimestep
654
655
656	c TEST D.BARDET
657	call WRITEDIAGFI(ngrid,"No_dust","Nombre particules de poussière"
658	& ,"part/kg",3,No_dust)
659	call WRITEDIAGFI(ngrid,"Mo_dust","Masse particules de poussière"
660	& ,"kg/kg ",3,Mo_dust)
661
662	END SUBROUTINE improvedCO2clouds
663
664	END MODULE improvedCO2clouds_mod
665

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

Download in other formats:

Original Format