Context Navigation

improvedco2clouds_mod.F90 @ 2997

Last change on this file since 2997 was 2660, checked in by cmathe, 3 years ago
correction on rsedcloudco2/opacities/riceco2 computations, and more comments
Property svn:executable set to ``*
File size: 30.6 KB

Line
1	!======================================================================================================================!
2	! Module: Scheme of co2 cloud formation ===============================================================================!
3	!----------------------------------------------------------------------------------------------------------------------!
4	! Authors: Joaquim Audouard, Constantino Listowski, Anni Määttänen
5	! Date: 09/2016
6	!----------------------------------------------------------------------------------------------------------------------!
7	! Contains subroutines:
8	! - improvedco2clouds_mod: nucleation
9	!======================================================================================================================!
10	module improvedco2clouds_mod
11
12	implicit none
13
14	contains
15	!======================================================================================================================!
16	! SUBROUTINE: improvedco2clouds =======================================================================================!
17	!======================================================================================================================!
18	! Subject:
19	!---------
20	! This routine is used to form CO2 clouds when a parcel of the GCM is saturated.
21	!----------------------------------------------------------------------------------------------------------------------!
22	! Comments:
23	!----------
24	! Adaptation for CO2 clouds based on improvedclouds_mod.F
25	!
26	! It includes the ability to have supersaturation, a computation of the nucleation rates, growthrates and the
27	! scavenging of dust particles by clouds. It is worth noting that the amount of dust is computed using the dust
28	! optical depth computed in aeropacity.F.
29	! That's why the variable called "tauscaling" is used to convert pq(dust_mass) and pq(dust_number), which are relative
30	! quantities, to absolute and realistic quantities stored in zq. This has to be done to convert the inputs into
31	! absolute values, but also to convert the outputs back into relative values which are then used by the sedimentation
32	! and advection schemes.
33	! CO2 ice particles can nucleate on both dust and water ice particles. When CO2 ice is deposited onto a water ice
34	! particles, the particle is removed from the water tracers. Memory of the origin of the co2 particles is kept and
35	! thus the water cycle shouldn't be modified by this.
36	! There is an energy limit to how much co2 can sublimate/condensate. It is defined by the difference of the GCM
37	! temperature with the co2 condensation temperature.
38	!
39	! /!\ WARNING /!\
40	! No sedimentation of the water ice origin is performed in the microphysical timestep in co2cloud.F.
41	!
42	! If meteoritic particles are activated and turn into co2 ice particles, then they will be reversed in the dust
43	! tracers if the cloud sublimates.
44	!----------------------------------------------------------------------------------------------------------------------!
45	! Paper:
46	!-------
47	! see co2cloud.F90
48	!----------------------------------------------------------------------------------------------------------------------!
49	! Algorithm:
50	!-----------
51	! 0. Firstcall
52	! 0.1. Bonus: meteoritic component, extract data
53	! 1. Initialization
54	! 2. Compute saturation
55	! 3. Bonus: additional meteoritic particles for nucleation
56	! 4. Actual microphysics: Main loop over the GCM's grid
57	! 4.1 Nucleation
58	! 4.2. Ice growth: scheme for radius evolution
59	! 4.3 Dust cores releasing if no more co2 ice
60	! 5. Get cloud tendencies
61	!======================================================================================================================!
62	subroutine improvedCO2clouds(ngrid, nlay, microtimestep, pplay, pplev, pteff, sum_subpdt, pqeff, sum_subpdq, &
63	subpdqcloudco2, subpdtcloudco2, nq, tauscaling, rb_cldco2, sigma_iceco2, dev2)
64
65	use comcstfi_h, only: pi, g, cpp
66
67	use updaterad, only: updaterice_micro, updaterice_microco2, updaterccnCO2
68
69	use tracer_mod, only: igcm_dust_mass, igcm_dust_number, rho_dust, igcm_h2o_ice, igcm_ccn_mass, igcm_ccn_number, &
70	nuice_sed, igcm_co2, igcm_co2_ice, igcm_ccnco2_mass, igcm_ccnco2_number, &
71	igcm_ccnco2_h2o_mass_ice, igcm_ccnco2_h2o_mass_ccn, igcm_ccnco2_h2o_number, nuiceco2_sed, &
72	nuiceco2_ref, igcm_ccnco2_meteor_mass, igcm_ccnco2_meteor_number
73
74	use conc_mod, only: mmean
75	use time_phylmdz_mod, only: daysec
76	use nucleaco2_mod, only: nucleaco2
77	use datafile_mod, only: datadir
78
79	use density_co2_ice_mod, only: density_co2_ice
80
81	implicit none
82
83	include "callkeys.h"
84	include "microphys.h"
85	!----------------------------------------------------------------------------------------------------------------------!
86	! VARIABLES DECLARATION
87	!----------------------------------------------------------------------------------------------------------------------!
88	! Input arguments:
89	!-----------------
90	integer, intent(in) :: &
91	nq, &! number of tracers
92	ngrid, &! number of point grid
93	nlay ! number of layer
94
95	real, intent(in) :: &
96	microtimestep, &! physics time step (s)
97	pplay(ngrid,nlay), &! mid-layer pressure (Pa)
98	pplev(ngrid,nlay+1), &! inter-layer pressure (Pa)
99	pteff(ngrid,nlay), &! temperature at the middle of the layers (K)
100	sum_subpdt(ngrid,nlay), &! tendency on temperature from previous physical parametrizations
101	pqeff(ngrid,nlay,nq), &! tracers (kg/kg)
102	tauscaling(ngrid), &! convertion factor for qdust and Ndust
103	sum_subpdq(ngrid,nlay,nq) ! tendencies on tracers before condensation (kg/kg.s-1)
104
105	real, intent(in) :: &
106	sigma_iceco2 ! Variance of the co2 ice and CCN distributions
107
108	double precision, intent(in) :: &
109	rb_cldco2(nbinco2_cld+1), & ! boundary values of each rad_cldco2
110	dev2
111	!----------------------------------------------------------------------------------------------------------------------!
112	! Output arguments:
113	!------------------
114	real, intent(out) :: &
115	subpdtcloudco2(ngrid,nlay), &! tendency on tracers due to CO2 condensation (K/s)
116	subpdqcloudco2(ngrid,nlay,nq) ! tendency on tracers due to CO2 condensation (kg/kg.s-1)
117	!----------------------------------------------------------------------------------------------------------------------!
118	! Local:
119	!-------
120	!----1) Parameters:
121	!------------------
122	integer, parameter :: &
123	! ---Meteoritic flux input file
124	nbin_meteor = 100, &! Dimension 1
125	nlev_meteor = 130, &! Dimension 2
126	uMeteor = 666, &! File unit
127	! ---Latent heat computation
128	l0 = 595594d0, &! coeff from: Azreg-Aïnou (2005)
129	l1 = 903.111d0, &! Title: "Low-temperature data for carbon dioxide"
130	l2 = -11.5959d0, &! Pulication: eprint arXiv:1403.4403
131	l3 = 0.0528288d0, &!
132	l4 = -0.000103183d0 !
133
134	real, parameter :: &
135	threshold = 1e-30, & ! limit value considering as zero
136	threshold_2 = 1e-13 ! limit value considering the value is physical (below this value => computer noise for dble)
137
138	character(len=50), parameter:: &
139	file_meteoritic_flux = 'Meteo_flux_Plane.dat'
140	!----------------------------------------------------------------------------------------------------------------------!
141	!----2) Saved:
142	!-------------
143	real, save :: &
144	sigma_ice ! Variance of the h2o ice and CCN distributions
145
146	!$OMP THREADPRIVATE(sigma_ice)
147
148	double precision, save :: &
149	pression_meteor(nlev_meteor), &! pressure from meteoritic flux input file
150	meteor(nlev_meteor,nbin_meteor), &! Meteoritic flux read from file uMeteor
151	dev3 ! 1. / ( sqrt(2.) * sigma_ice )
152
153	!$OMP THREADPRIVATE(meteor,dev3)
154
155	logical, save :: &
156	firstcall = .true. ! Used to compute saved variables
157
158	!$OMP THREADPRIVATE(firstcall)
159
160	!----------------------------------------------------------------------------------------------------------------------!
161	!----3) Variables:
162	!-----------------
163	integer :: &
164	ig, &! loop on ngrid
165	l, &! loop on nlay
166	i, &! loop on nbinco2
167	! ---Variables for meteoritic flux
168	ibin, &! loop on nbin_meteor
169	idx_min,&! index of min(diff_pressure)
170	read_ok ! file uMeteor iostat
171
172	real :: &
173	Nccnco2, &
174	Qccnco2, &
175	Nccnco2_h2o, &
176	Qccnco2_h2o, &
177	masse(ngrid,nlay), &! mass layer (kg.m-2)
178	rice(ngrid,nlay), &! water ice mass mean radius (m): used for nucleation of CO2 on ice-coated ccns
179	zq(ngrid,nlay,nq), &! local value of tracers (kg/kg)
180	zq0(ngrid,nlay,nq), &! local init value of tracers (kg/kg)
181	zt(ngrid,nlay), &! local value of temperature (K)
182	zqsat(ngrid,nlay), &! saturation vapor pressure for CO2 (K)
183	tcond(ngrid,nlay), &! condensation temperature of CO2 (K)
184	lw, &! Latent heat of sublimation (J.kg-1)
185	rdust(ngrid,nlay), &! Dust geometric mean radius (m)
186	rhocloud(ngrid,nlay), &! Cloud density (kg.m-3)
187	rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
188
189	real(kind=8) :: &
190	derf ! Error function
191
192	double precision :: &
193	dMice, &! mass of condensed ice
194	facteurmax, &! for energy limit on mass growth
195	pco2, &! Co2 vapor partial pressure (Pa)
196	satu, &! Co2 vapor saturation ratio over ice
197	Mo, &! mass of dust particles
198	No, &! number of dust particles
199	Rn, &! logarithm of rdust/rice
200	Rm, &! Rn * variance of ice and CCN distribution
201	n_derf, &! derf( (rb_cldco2(1)+Rn) *dev3)
202	m_derf, &! derf( (rb_cldco2(1)+Rm) *dev2)
203	n_aer(nbinco2_cld), &! Radius used by the microphysical scheme (m) for dust
204	m_aer(nbinco2_cld), &! number concentration V-1 of particle/each size bin for dust
205	n_aer_meteor(nbinco2_cld), &! Radius used by the microphysical scheme (m) for meteor
206	m_aer_meteor(nbinco2_cld), &! number concentration V-1 of particle/each size bin for meteor
207	n_aer_h2oice(nbinco2_cld), &! mass mixing ratio of particle/each size bin for h2o ice
208	m_aer_h2oice(nbinco2_cld), &! Same - for CO2 nucleation for h2o ice
209	Ic_rice, &! Mass transfer rate CO2 ice crystal
210	ratioh2o_ccn, &! 1./(zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
211	vo2co2, &! volume of co2 ice particle
212	dN, &! number of particle of dust used as ccn
213	dM, &! mass of dN
214	dN_meteor, &! number of particle of meteoritic particles used as ccn
215	dM_meteor, &! mass of dN_meteor
216	dNh2o, &! number of particle of h2o ice used as ccn
217	dMh2o, &! mass of dNh2o
218	dNN, &! min(dN,zq(ig,l,igcm_dust_number))
219	dMM, &! min(dM,zq(ig,l,igcm_dust_mass))
220	dNNh2o, &! min(dNNh2o,zq(ig,l,igcm_ccn_number))
221	dMh2o_ice, &! min(dMh2ozq(ig,l,igcm_h2o_ice)ratioh2o_ccn, zq(ig,l,igcm_h2o_ice))
222	dMh2o_ccn, &! min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
223	rate(nbinco2_cld), &! nucleation rate
224	rate_meteor(nbinco2_cld), &! nucleation rate for meteor
225	rateh2o(nbinco2_cld), &! nucleation rate for h2o
226	rho_ice_co2T, &! density of co2 ice Temperature-dependent
227	riceco2(ngrid,nlay), &! CO2 ice mean radius (m)
228	vrat_cld, &! Volume ratio
229	Proba, &! 1.0 - exp(-1.microtimesteprate(i))
230	Probah2o, &! 1.0 - exp(-1.microtimesteprateh2o(i))
231	Proba_meteor, &! 1.0 - exp(-1.microtimesteprate_meteor(i))
232	diff_pressure(nlev_meteor), &! abs(pression_meteor(:)-pplev(ig,l))
233	meteor_ccn(ngrid,nlay,nbinco2_cld) ! input flux of meteoritc particles
234
235	logical :: &
236	file_ok ! test if meteoritic input file exists
237	!======================================================================================================================!
238	! BEGIN ===============================================================================================================!
239	!======================================================================================================================!
240	! 0. Firstcall
241	!----------------------------------------------------------------------------------------------------------------------!
242	if (firstcall) then
243	firstcall = .false.
244
245	! Variance of the ice and CCN distributions
246	sigma_ice = sqrt(log(1.+nuice_sed))
247	dev3 = 1. / ( sqrt(2.) * sigma_ice )
248	!----------------------------------------------------------------------------------------------------------------------!
249	! 0.1. Bonus: meteoritic component, extract data
250	!----------------------------------------------------------------------------------------------------------------------!
251
252	if (meteo_flux) then
253	! Check if file exists
254	inquire(file=trim(datadir)//'/'//trim(file_meteoritic_flux), exist=file_ok)
255	if (.not. file_ok) then
256	call abort_physic("CO2clouds", 'file '//trim(file_meteoritic_flux)//' should be in'//trim(datadir), 1)
257	end if
258	write(,)'Meteoritic flux file used: ', trim(file_meteoritic_flux)
259
260	! open file
261	open(unit=uMeteor,file=trim(datadir)//'/'//trim(file_meteoritic_flux), FORM='formatted')
262
263	!skip 1 line
264	read(uMeteor,*)
265
266	! extract pressure_meteor
267	do i = 1, nlev_meteor
268	read(uMeteor,*)pression_meteor(i)
269	end do
270
271	!skip 1 line
272	read(uMeteor,*)
273
274	! extract meteor flux
275	do i = 1, nlev_meteor
276	! les mêmes 100 bins size que la distri nuclea : on touche pas
277	do ibin = 1, nbin_meteor
278	read(uMeteor,'(F12.6)')meteor(i,ibin)
279	end do
280	end do
281
282	! close file
283	close(uMeteor)
284	end if ! of if meteo_flux
285	end if ! firstcall
286	!----------------------------------------------------------------------------------------------------------------------!
287	! 1. Initialization
288	!----------------------------------------------------------------------------------------------------------------------!
289	rdust(:,:) = 0.
290	meteor_ccn(:,:,:) = 0d0
291	rice(:,:) = 1.e-8
292	riceco2(:,:) = 1.e-11
293
294	! Initialize the tendencies
295	subpdqcloudco2(:,:,:) = 0.
296	subpdtcloudco2(:,:) = 0.
297
298	! pteff temperature layer; sum_subpdt dT.s-1
299	zt(1:ngrid,1:nlay) = 0.
300	zt(:,:) = pteff(:,:) + sum_subpdt(:,:) * microtimestep
301
302	! pqeff traceur kg/kg; sum_subpdq tendance idem .s-1
303	zq(:,:,:) = pqeff(:,:,:) + sum_subpdq(:,:,:) * microtimestep
304	where( zq(:,:,:) < threshold ) zq(:,:,:) = threshold
305
306	zq0(:,:,:) = zq(:,:,:)
307	zqsat(:,:) = 0.
308	!----------------------------------------------------------------------------------------------------------------------!
309	! 2. Compute saturation
310	!----------------------------------------------------------------------------------------------------------------------!
311	call co2sat(ngrid*nlay,zt,zqsat)
312	call tcondco2(ngrid,nlay,pplay, zq(:,:,igcm_co2) + zq(:,:,igcm_co2_ice),tcond)
313	!----------------------------------------------------------------------------------------------------------------------!
314	! 3. Bonus: additional meteoritic particles for nucleation
315	!----------------------------------------------------------------------------------------------------------------------!
316	if (meteo_flux) then
317	do l = 1, nlay
318	do ig = 1, ngrid
319	masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1)) / g
320
321	diff_pressure(:) = abs(pression_meteor(:)-pplev(ig,l))
322
323	idx_min = minloc(diff_pressure, DIM=1)
324	! Par kg air csi par m carre, x epaisseur/masse pour par kg/air. Check original unit with J. Plane
325	meteor_ccn(ig,l,:) = meteor(idx_min,:)/masse(ig,l)
326	end do
327	end do
328	end if
329	!----------------------------------------------------------------------------------------------------------------------!
330	! 4. Actual microphysics: Main loop over the GCM's grid
331	!----------------------------------------------------------------------------------------------------------------------!
332	do l = 1, nlay
333	do ig = 1, ngrid
334
335	! Get the partial pressure of co2 vapor and its saturation ratio
336	pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/(mco21e3)) pplay(ig,l) ! mco2 (kg/mol) => mmean (g/mol)
337	satu = pco2 / zqsat(ig,l)
338
339	!T-dependant CO2 ice density
340	call density_co2_ice(dble(zt(ig,l)), rho_ice_co2T)
341
342	vo2co2 = m0co2 / rho_ice_co2T
343	!----------------------------------------------------------------------------------------------------------------------!
344	! 4.1 Nucleation
345	!----------------------------------------------------------------------------------------------------------------------!
346	! if there is condensation
347	if ( satu >= 1 ) then
348	call updaterccnCO2(zq(ig,l,igcm_dust_mass), zq(ig,l,igcm_dust_number), rdust(ig,l), tauscaling(ig))
349
350	! Expand the dust moments into a binned distribution
351	n_aer(:) = 0d0 ! number of aerosol particles
352	m_aer(:) = 0d0 ! mass of aerosol particles
353
354	No = zq(ig,l,igcm_dust_number) * tauscaling(ig)
355
356	Mo = (4./3.) * pi * rho_dust * No * rdust(ig,l)*3 exp(9.*nuiceco2_ref/2.)
357
358	if (No > threshold) then
359	Rn = -log(rdust(ig,l))
360
361	Rm = Rn - 3. * sigma_iceco2 * sigma_iceco2
362
363	n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
364	m_derf = derf( (rb_cldco2(1)+Rm) *dev2)
365
366	do i = 1, nbinco2_cld
367	n_aer(i) = -0.5 * No * n_derf
368	m_aer(i) = -0.5 * Mo * m_derf
369
370	n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
371	m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
372
373	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
374	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
375	end do
376
377	! Call to nucleation routine
378	call nucleaco2(dble(pco2), zt(ig,l), dble(satu), n_aer, rate, vo2co2, mtetaco2)
379	dN = 0.
380	dM = 0.
381
382	do i = 1, nbinco2_cld
383	Proba = 1.0 - exp(-1.microtimesteprate(i))
384	dN = dN + n_aer(i) * Proba
385	dM = dM + m_aer(i) * Proba
386	end do
387
388	! Now increment CCN tracers and update dust tracers
389	dNN = min(dN, zq(ig,l,igcm_dust_number)) ! dNN est devenu DN
390	dMM = min(dM, zq(ig,l,igcm_dust_mass)) ! idem dans le min
391
392	zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dMM /tauscaling(ig)
393	zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dNN /tauscaling(ig)
394
395	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dMM /tauscaling(ig)
396	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dNN /tauscaling(ig)
397	end if ! of if No > 1e-30
398
399	! Ajout meteor_ccn particles aux particules de poussière background
400	if (meteo_flux) then
401	n_aer_meteor(1:nbinco2_cld) = 0d0
402	m_aer_meteor(1:nbinco2_cld) = 0d0
403	do i = 1, nbinco2_cld
404	n_aer_meteor(i) = meteor_ccn(ig,l,i)
405	m_aer_meteor(i) = (4./3.) * pi * rho_dust * meteor_ccn(ig,l,i) * rad_cldco2(i)**3
406	end do
407	! Call to nucleation routine
408	rate_meteor(1:nbinco2_cld) = 0d0
409	call nucleaco2(dble(pco2), zt(ig,l), dble(satu), n_aer_meteor, rate_meteor, vo2co2, mtetaco2)
410
411	dN_meteor = 0.
412	dM_meteor = 0.
413	do i = 1, nbinco2_cld
414	Proba_meteor = 1.0 - exp(-1.microtimesteprate_meteor(i))
415	dN_meteor = dN_meteor + n_aer_meteor(i) * Proba_meteor
416	dM_meteor = dM_meteor + m_aer_meteor(i) * Proba_meteor
417	end do
418	! Now increment CCN tracers and update dust tracers
419	zq(ig,l,igcm_ccnco2_meteor_mass) = zq(ig,l,igcm_ccnco2_meteor_mass) + dM_meteor
420	zq(ig,l,igcm_ccnco2_meteor_number) = zq(ig,l,igcm_ccnco2_meteor_number) + dN_meteor
421
422	zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dM_meteor
423	zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dN_meteor
424	end if
425
426	! Same but with h2o particles as CCN only if co2useh2o = .true.
427	if (co2useh2o) then
428	call updaterice_micro(zq(ig,l,igcm_h2o_ice), zq(ig,l,igcm_ccn_mass), zq(ig,l,igcm_ccn_number), &
429	tauscaling(ig), rice(ig,l), rhocloud(ig,l))
430
431	! Total mass of H20 crystals,CCN included
432	Mo = zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccn_mass) * tauscaling(ig) + threshold
433
434	No = zq(ig,l,igcm_ccn_number) * tauscaling(ig) + threshold
435	if (No > threshold) then
436	Rn = -log(rice(ig,l))
437
438	Rm = Rn - 3. * sigma_ice * sigma_ice
439
440	n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
441	m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
442
443	n_aer_h2oice(:)=0.
444	m_aer_h2oice(:)=0.
445	do i = 1, nbinco2_cld
446	n_aer_h2oice(i) = -0.5 * No * n_derf
447	m_aer_h2oice(i) = -0.5 * Mo * m_derf
448
449	n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
450	m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
451
452	n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf
453	m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf
454	end do
455
456	call nucleaco2(dble(pco2), zt(ig,l), dble(satu), n_aer_h2oice, rateh2o, vo2co2, mteta)
457	dNh2o = 0.
458	dMh2o = 0.
459
460	do i = 1, nbinco2_cld
461	Probah2o = 1.0 - exp(-1.microtimesteprateh2o(i))
462	dNh2o = dNh2o + n_aer_h2oice(i) * Probah2o
463	dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o
464	end do
465
466	! Update CCN for CO2 nucleating on H2O CCN : Warning: must keep memory of it
467	dNNh2o = dNh2o/tauscaling(ig)
468	dNNh2o = min(dNNh2o,zq(ig,l,igcm_ccn_number))
469
470	ratioh2o_ccn = 1./(zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
471
472	dMh2o_ccn = dMh2o * zq(ig,l,igcm_ccn_mass) * tauscaling(ig) * ratioh2o_ccn
473	dMh2o_ccn = dMh2o_ccn/tauscaling(ig)
474	dMh2o_ccn = min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
475
476	dMh2o_ice = dMh2o * zq(ig,l,igcm_h2o_ice) * ratioh2o_ccn
477	dMh2o_ice = min(dMh2o_ice,zq(ig,l,igcm_h2o_ice))
478
479	zq(ig,l,igcm_ccnco2_h2o_mass_ice) = zq(ig,l,igcm_ccnco2_h2o_mass_ice) + dMh2o_ice
480	zq(ig,l,igcm_ccnco2_h2o_mass_ccn) = zq(ig,l,igcm_ccnco2_h2o_mass_ccn) + dMh2o_ccn
481	zq(ig,l,igcm_ccnco2_h2o_number) = zq(ig,l,igcm_ccnco2_h2o_number) + dNNh2o
482
483	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) - dNNh2o
484	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) - dMh2o_ice
485	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) - dMh2o_ccn
486
487	zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dMh2o_ice + dMh2o_ccn
488	zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dNNh2o
489	end if
490	end if ! of if co2useh2o
491	end if ! of is satu > 1
492	!----------------------------------------------------------------------------------------------------------------------!
493	! 4.2. Ice growth: scheme for radius evolution
494	!----------------------------------------------------------------------------------------------------------------------!
495	! We trigger crystal growth if and only if there is at least one nuclei (N>1). Indeed, if we are supersaturated
496	! and still don't have at least one nuclei, we should better wait to avoid unrealistic value for nuclei radius
497	! and so on for cases that remain negligible.
498	!----------------------------------------------------------------------------------------------------------------------!
499	! we trigger crystal growth
500	if ((zq(ig,l,igcm_ccnco2_number)) * tauscaling(ig) + threshold >= 1) then
501	Nccnco2 = dble(zq(ig,l,igcm_ccnco2_number))
502	Qccnco2 = dble(zq(ig,l,igcm_ccnco2_mass))
503	Nccnco2_h2o = 0.
504	Qccnco2_h2o = 0.
505
506	if (co2useh2o) then
507	Nccnco2_h2o = zq(ig,l,igcm_ccnco2_h2o_number)
508	Qccnco2_h2o = zq(ig,l,igcm_ccnco2_h2o_mass_ice) + zq(ig,l,igcm_ccnco2_h2o_mass_ccn)
509	Nccnco2 = Nccnco2 - Nccnco2_h2o
510	Qccnco2 = Qccnco2 - Qccnco2_h2o
511	if (Nccnco2 <= 0.) then
512	Nccnco2 = threshold
513	Qccnco2 = threshold
514	end if
515	end if
516
517	call updaterice_microco2(dble(zq(ig,l,igcm_co2_ice)), dble(Qccnco2), dble(Nccnco2), dble(Qccnco2_h2o), &
518	dble(Nccnco2_h2o), zt(ig,l), tauscaling(ig), riceco2(ig,l), rhocloudco2(ig,l))
519	Ic_rice = 0.
520
521	! J.kg-1
522	lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)*2 + l3 zt(ig,l)*3 + l4 zt(ig,l)**4
523
524	facteurmax = abs(tcond(ig,l)-zt(ig,l)) * (cpp/lw)
525
526	! call scheme of microphys. mass growth for CO2 (evaporation/condensation)
527	call massflowrateCO2(pplay(ig,l), zt(ig,l), satu, riceco2(ig,l), mmean(ig,l), Ic_rice)
528
529	! Ic_rice Mass transfer rate (kg/s) for a rice particle > 0 si croissance !
530	if (isnan(Ic_rice) .or. Ic_rice == 0.) then
531	Ic_rice = 0.
532	subpdtcloudco2(ig,l) = -sum_subpdt(ig,l)
533	dMice = 0
534	else
535	! Kg par kg d'air, >0 si croissance !
536	! kg.s-1 par particule * nb particule par kg air*s = kg par kg air
537	dMice = zq(ig,l,igcm_ccnco2_number) * Ic_rice * microtimestep * tauscaling(ig)
538
539	! facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
540	! latent heat release > 0 if growth i.e. if dMice > 0
541	dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
542	dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
543
544	! kgco2/kgair* J/kgco2 * 1/(J.kgair-1.K-1)/s = K /s
545	subpdtcloudco2(ig,l) = dMice * lw / cpp / microtimestep
546
547	!Now update tracers
548	zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice) + dMice
549	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) - dMice
550	end if
551	end if ! if zq(ccnco2_number) >= 1
552	!----------------------------------------------------------------------------------------------------------------------!
553	! 4.3 Dust cores releasing if no more co2 ice
554	!----------------------------------------------------------------------------------------------------------------------!
555	! On sublime tout
556	if ((zq(ig,l,igcm_co2_ice) < threshold_2).or.(zq(ig,l,igcm_ccnco2_number) * tauscaling(ig) < 1.)) then
557	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccnco2_mass)
558	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccnco2_number)
559	zq(ig,l,igcm_ccnco2_mass) = 0.
560	zq(ig,l,igcm_ccnco2_number) = 0.
561	if (meteo_flux) then
562	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - zq(ig,l,igcm_ccnco2_meteor_mass)
563	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - zq(ig,l,igcm_ccnco2_meteor_number)
564	zq(ig,l,igcm_ccnco2_meteor_mass) = 0.
565	zq(ig,l,igcm_ccnco2_meteor_number) = 0.
566	end if
567	if (co2useh2o) then
568	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - zq(ig,l,igcm_ccnco2_h2o_mass_ccn) - &
569	zq(ig,l,igcm_ccnco2_h2o_mass_ice)
570	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - zq(ig,l,igcm_ccnco2_h2o_number)
571
572	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + zq(ig,l,igcm_ccnco2_h2o_mass_ccn)
573	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccnco2_h2o_mass_ice)
574	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + zq(ig,l,igcm_ccnco2_h2o_number)
575	zq(ig,l,igcm_ccnco2_h2o_number) = 0.
576	zq(ig,l,igcm_ccnco2_h2o_mass_ccn) = 0.
577	zq(ig,l,igcm_ccnco2_h2o_mass_ice) = 0.
578	end if
579	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) + zq(ig,l,igcm_co2_ice)
580	zq(ig,l,igcm_co2_ice) = 0.
581	riceco2(ig,l) = 0.
582	end if !of if co2_ice < threshold or zq(ccnco2_number) < 1
583	end do ! of ig loop
584	end do ! of nlayer loop
585	!----------------------------------------------------------------------------------------------------------------------!
586	! 5. Get cloud tendencies
587	!----------------------------------------------------------------------------------------------------------------------!
588	subpdqcloudco2(:,:,igcm_co2) = ( zq(:,:,igcm_co2) - zq0(:,:,igcm_co2) ) / microtimestep
589
590	subpdqcloudco2(:,:,igcm_co2_ice) = ( zq(:,:,igcm_co2_ice) - zq0(:,:,igcm_co2_ice) ) / microtimestep
591
592	subpdqcloudco2(:,:,igcm_ccnco2_mass) = ( zq(:,:,igcm_ccnco2_mass) - zq0(:,:,igcm_ccnco2_mass) ) / microtimestep
593
594	subpdqcloudco2(:,:,igcm_ccnco2_number) = ( zq(:,:,igcm_ccnco2_number) - zq0(:,:,igcm_ccnco2_number))/microtimestep
595
596	subpdqcloudco2(:,:,igcm_dust_mass) = ( zq(:,:,igcm_dust_mass) - zq0(:,:,igcm_dust_mass) ) / microtimestep
597
598	subpdqcloudco2(:,:,igcm_dust_number) = ( zq(:,:,igcm_dust_number) - zq0(:,:,igcm_dust_number) ) / microtimestep
599
600	if (meteo_flux) then
601	subpdqcloudco2(:,:,igcm_ccnco2_meteor_mass) = (zq(:,:,igcm_ccnco2_meteor_mass)-zq0(:,:,igcm_ccnco2_meteor_mass) &
602	)/microtimestep
603
604	subpdqcloudco2(:,:,igcm_ccnco2_meteor_number) = ( zq(:,:,igcm_ccnco2_meteor_number) - &
605	zq0(:,:,igcm_ccnco2_meteor_number) )/microtimestep
606	end if
607
608	if (co2useh2o) then
609	subpdqcloudco2(:,:,igcm_h2o_ice) = ( zq(:,:,igcm_h2o_ice) - zq0(:,:,igcm_h2o_ice) ) / microtimestep
610
611	subpdqcloudco2(:,:,igcm_ccn_mass) = ( zq(:,:,igcm_ccn_mass) - zq0(:,:,igcm_ccn_mass) ) / microtimestep
612
613	subpdqcloudco2(:,:,igcm_ccn_number) = ( zq(:,:,igcm_ccn_number) - zq0(:,:,igcm_ccn_number) ) / microtimestep
614
615	subpdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ice) = (zq(:,:,igcm_ccnco2_h2o_mass_ice)-zq0(:,:,igcm_ccnco2_h2o_mass_ice)&
616	)/microtimestep
617
618	subpdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ccn) = (zq(:,:,igcm_ccnco2_h2o_mass_ccn)-zq0(:,:,igcm_ccnco2_h2o_mass_ccn)&
619	)/microtimestep
620
621	subpdqcloudco2(:,:,igcm_ccnco2_h2o_number) = ( zq(:,:,igcm_ccnco2_h2o_number) - zq0(:,:,igcm_ccnco2_h2o_number) &
622	)/microtimestep
623	end if
624	!======================================================================================================================!
625	! END =================================================================================================================!
626	!======================================================================================================================!
627	end subroutine improvedCO2clouds
628
629	end module improvedCO2clouds_mod
630

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

Download in other formats:

Original Format