Context Navigation

improvedco2clouds_mod.F90 @ 3118

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

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

Download in other formats:

Original Format