Context Navigation

improvedco2clouds_mod.F90 @ 2616

Last change on this file since 2616 was 2616, checked in by romain.vande, 3 years ago
LMDZ_MARS RV : Open_MP; Put all the "save" variables as "!$OMP THREADPRIVATE" in phymars. The code can now be tested, see README for more info
Property svn:executable set to ``*
File size: 30.3 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
136
137	character(len=23), parameter:: &
138	file_meteoritic_flux = 'Meteo_flux_Plane.dat'
139	!----------------------------------------------------------------------------------------------------------------------!
140	!----2) Saved:
141	!-------------
142	real, save :: &
143	sigma_ice ! Variance of the h2o ice and CCN distributions
144
145	!$OMP THREADPRIVATE(sigma_ice)
146
147	double precision, save :: &
148	pression_meteor(nlev_meteor), &! pressure from meteoritic flux input file
149	meteor(nlev_meteor,nbin_meteor), &! Meteoritic flux read from file uMeteor
150	dev3 ! 1. / ( sqrt(2.) * sigma_ice )
151
152	!$OMP THREADPRIVATE(meteor,dev3)
153
154	logical, save :: &
155	firstcall = .true. ! Used to compute saved variables
156
157	!$OMP THREADPRIVATE(firstcall)
158
159	!----------------------------------------------------------------------------------------------------------------------!
160	!----3) Variables:
161	!-----------------
162	integer :: &
163	ig, &! loop on ngrid
164	l, &! loop on nlay
165	i, &! loop on nbinco2
166	! ---Variables for meteoritic flux
167	ibin, &! loop on nbin_meteor
168	idx_min,&! index of min(diff_pressure)
169	read_ok ! file uMeteor iostat
170
171	real :: &
172	Nccnco2, &
173	Qccnco2, &
174	Nccnco2_h2o, &
175	Qccnco2_h2o, &
176	masse(ngrid,nlay), &! mass layer (kg.m-2)
177	rice(ngrid,nlay), &! water ice mass mean radius (m): used for nucleation of CO2 on ice-coated ccns
178	zq(ngrid,nlay,nq), &! local value of tracers (kg/kg)
179	zq0(ngrid,nlay,nq), &! local init value of tracers (kg/kg)
180	zt(ngrid,nlay), &! local value of temperature (K)
181	zqsat(ngrid,nlay), &! saturation vapor pressure for CO2 (K)
182	tcond(ngrid,nlay), &! condensation temperature of CO2 (K)
183	lw, &! Latent heat of sublimation (J.kg-1)
184	rdust(ngrid,nlay), &! Dust geometric mean radius (m)
185	rhocloud(ngrid,nlay), &! Cloud density (kg.m-3)
186	rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
187
188	real(kind=8) :: &
189	derf ! Error function
190
191	double precision :: &
192	dMice, &! mass of condensed ice
193	facteurmax, &! for energy limit on mass growth
194	pco2, &! Co2 vapor partial pressure (Pa)
195	satu, &! Co2 vapor saturation ratio over ice
196	Mo, &! mass of aerosol particles
197	No, &! number of aerosol particles
198	Rn, &! logarithm of rdust/rice
199	Rm, &! Rn * variance of ice and CCN distribution
200	n_derf, &! derf( (rb_cldco2(1)+Rn) *dev3)
201	m_derf, &! derf( (rb_cldco2(1)+Rm) *dev2)
202	n_aer(nbinco2_cld), &! Radius used by the microphysical scheme (m)
203	m_aer(nbinco2_cld), &! number concentration V-1 of particle/each size bin
204	n_aer_meteor(nbinco2_cld), &! Radius used by the microphysical scheme (m)
205	m_aer_meteor(nbinco2_cld), &! number concentration V-1 of particle/each size bin
206	n_aer_h2oice(nbinco2_cld), &! mass mixing ratio of particle/each size bin
207	m_aer_h2oice(nbinco2_cld), &! Same - for CO2 nucleation
208	Ic_rice, &! Mass transfer rate CO2 ice crystal
209	ratioh2o_ccn, &! 1./(zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
210	vo2co2, &! volume of co2 ice particle
211	dN, &! number of particle of dust used as ccn
212	dM, &! mass of dN
213	dN_meteor, &! number of particle of dust used as ccn
214	dM_meteor, &! mass of dN
215	dNh2o, &! number of particle of h2o ice used as ccn
216	dMh2o, &! mass of dNh2o
217	dNN, &! min(dN,zq(ig,l,igcm_dust_number))
218	dMM, &! min(dM,zq(ig,l,igcm_dust_mass))
219	dNN_meteor, &! number of particle of dust used as ccn
220	dMM_meteor, &! mass of dN
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
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.microtimesteprateh2o(i))
233	diff_pressure(nlev_meteor), &! abs(pression_meteor(:)-pplev(ig,l))
234	meteor_ccn(ngrid,nlay,nbinco2_cld) !
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)//'/'//file_meteoritic_flux, exist=file_ok)
256	if (.not. file_ok) then
257	call abort_physic("CO2clouds", 'file '//file_meteoritic_flux//' should be in'//trim(datadir), 1)
258	end if
259
260	! open file
261	open(unit=uMeteor,file=trim(datadir)//'/'//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	end if
512
513	call updaterice_microco2(dble(zq(ig,l,igcm_co2_ice)), dble(Qccnco2), dble(Nccnco2), dble(Qccnco2_h2o), &
514	dble(Nccnco2_h2o), zt(ig,l), tauscaling(ig), riceco2(ig,l), rhocloudco2(ig,l))
515	Ic_rice = 0.
516
517	! J.kg-1
518	lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)*2 + l3 zt(ig,l)*3 + l4 zt(ig,l)**4
519
520	facteurmax = abs(tcond(ig,l)-zt(ig,l)) * (cpp/lw)
521
522	! call scheme of microphys. mass growth for CO2 (evaporation/condensation)
523	call massflowrateCO2(pplay(ig,l), zt(ig,l), satu, riceco2(ig,l), mmean(ig,l), Ic_rice)
524
525	! Ic_rice Mass transfer rate (kg/s) for a rice particle > 0 si croissance !
526	if (isnan(Ic_rice) .or. Ic_rice == 0.) then
527	Ic_rice = 0.
528	subpdtcloudco2(ig,l) = -sum_subpdt(ig,l)
529	dMice = 0
530	else
531	! Kg par kg d'air, >0 si croissance !
532	! kg.s-1 par particule * nb particule par kg air*s = kg par kg air
533	dMice = zq(ig,l,igcm_ccnco2_number) * Ic_rice * microtimestep * tauscaling(ig)
534
535	! facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
536	! latent heat release > 0 if growth i.e. if dMice > 0
537	dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
538	dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
539
540	! kgco2/kgair* J/kgco2 * 1/(J.kgair-1.K-1)/s = K /s
541	subpdtcloudco2(ig,l) = dMice * lw / cpp / microtimestep
542
543	!Now update tracers
544	zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice) + dMice
545	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) - dMice
546	end if
547	end if ! if zq(ccnco2_number) >= 1
548	!----------------------------------------------------------------------------------------------------------------------!
549	! 4.3 Dust cores releasing if no more co2 ice
550	!----------------------------------------------------------------------------------------------------------------------!
551	! On sublime tout
552	if ((zq(ig,l,igcm_co2_ice) <= threshold).or.(zq(ig,l,igcm_ccnco2_number) * tauscaling(ig) < 1.)) then
553	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccnco2_mass)
554	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccnco2_number)
555	zq(ig,l,igcm_ccnco2_mass) = 0.
556	zq(ig,l,igcm_ccnco2_number) = 0.
557	if (meteo_flux) then
558	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - zq(ig,l,igcm_ccnco2_meteor_mass)
559	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - zq(ig,l,igcm_ccnco2_meteor_number)
560	zq(ig,l,igcm_ccnco2_meteor_mass) = 0.
561	zq(ig,l,igcm_ccnco2_meteor_number) = 0.
562	end if
563	if (co2useh2o) then
564	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - zq(ig,l,igcm_ccnco2_h2o_mass_ccn) - &
565	zq(ig,l,igcm_ccnco2_h2o_mass_ice)
566	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - zq(ig,l,igcm_ccnco2_h2o_number)
567
568	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + zq(ig,l,igcm_ccnco2_h2o_mass_ccn)
569	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccnco2_h2o_mass_ice)
570	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + zq(ig,l,igcm_ccnco2_h2o_number)
571	zq(ig,l,igcm_ccnco2_h2o_number) = 0.
572	zq(ig,l,igcm_ccnco2_h2o_mass_ccn) = 0.
573	zq(ig,l,igcm_ccnco2_h2o_mass_ice) = 0.
574	end if
575	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) + zq(ig,l,igcm_co2_ice)
576	zq(ig,l,igcm_co2_ice) = 0.
577	riceco2(ig,l) = 0.
578	end if !of if co2_ice < threshold or zq(ccnco2_number) < 1
579	end do ! of ig loop
580	end do ! of nlayer loop
581	!----------------------------------------------------------------------------------------------------------------------!
582	! 5. Get cloud tendencies
583	!----------------------------------------------------------------------------------------------------------------------!
584	subpdqcloudco2(:,:,igcm_co2) = ( zq(:,:,igcm_co2) - zq0(:,:,igcm_co2) ) / microtimestep
585
586	subpdqcloudco2(:,:,igcm_co2_ice) = ( zq(:,:,igcm_co2_ice) - zq0(:,:,igcm_co2_ice) ) / microtimestep
587
588	subpdqcloudco2(:,:,igcm_ccnco2_mass) = ( zq(:,:,igcm_ccnco2_mass) - zq0(:,:,igcm_ccnco2_mass) ) / microtimestep
589
590	subpdqcloudco2(:,:,igcm_ccnco2_number) = ( zq(:,:,igcm_ccnco2_number) - zq0(:,:,igcm_ccnco2_number))/microtimestep
591
592	subpdqcloudco2(:,:,igcm_dust_mass) = ( zq(:,:,igcm_dust_mass) - zq0(:,:,igcm_dust_mass) ) / microtimestep
593
594	subpdqcloudco2(:,:,igcm_dust_number) = ( zq(:,:,igcm_dust_number) - zq0(:,:,igcm_dust_number) ) / microtimestep
595
596	if (meteo_flux) then
597	subpdqcloudco2(:,:,igcm_ccnco2_meteor_mass) = (zq(:,:,igcm_ccnco2_meteor_mass)-zq0(:,:,igcm_ccnco2_meteor_mass) &
598	)/microtimestep
599
600	subpdqcloudco2(:,:,igcm_ccnco2_meteor_number) = ( zq(:,:,igcm_ccnco2_meteor_number) - &
601	zq0(:,:,igcm_ccnco2_meteor_number) )/microtimestep
602	end if
603
604	if (co2useh2o) then
605	subpdqcloudco2(:,:,igcm_h2o_ice) = ( zq(:,:,igcm_h2o_ice) - zq0(:,:,igcm_h2o_ice) ) / microtimestep
606
607	subpdqcloudco2(:,:,igcm_ccn_mass) = ( zq(:,:,igcm_ccn_mass) - zq0(:,:,igcm_ccn_mass) ) / microtimestep
608
609	subpdqcloudco2(:,:,igcm_ccn_number) = ( zq(:,:,igcm_ccn_number) - zq0(:,:,igcm_ccn_number) ) / microtimestep
610
611	subpdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ice) = (zq(:,:,igcm_ccnco2_h2o_mass_ice)-zq0(:,:,igcm_ccnco2_h2o_mass_ice)&
612	)/microtimestep
613
614	subpdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ccn) = (zq(:,:,igcm_ccnco2_h2o_mass_ccn)-zq0(:,:,igcm_ccnco2_h2o_mass_ccn)&
615	)/microtimestep
616
617	subpdqcloudco2(:,:,igcm_ccnco2_h2o_number) = ( zq(:,:,igcm_ccnco2_h2o_number) - zq0(:,:,igcm_ccnco2_h2o_number) &
618	)/microtimestep
619	end if
620	!======================================================================================================================!
621	! END =================================================================================================================!
622	!======================================================================================================================!
623	end subroutine improvedCO2clouds
624
625	end module improvedCO2clouds_mod
626

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

Download in other formats:

Original Format