Context Navigation

improvedco2clouds_mod.F90 @ 3807

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

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

Download in other formats:

Original Format