Context Navigation

co2cloud.F90 @ 3325

Last change on this file since 3325 was 3008, checked in by emillour, 16 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: 61.3 KB

Line
1	!======================================================================================================================!
2	! Module: CO2 clouds formation ========================================================================================!
3	!----------------------------------------------------------------------------------------------------------------------!
4	! Authors: Joachim Audouard, Constantino Listowski, Anni Määttänen
5	! Date: 09/2016
6	!----------------------------------------------------------------------------------------------------------------------!
7	! Contains subroutines:
8	! - co2cloud: of co2 cloud microphysics
9	!======================================================================================================================!
10	module co2cloud_mod
11
12	implicit none
13
14	contains
15	!======================================================================================================================!
16	! SUBROUTINE: co2cloud ================================================================================================!
17	!======================================================================================================================!
18	! Subject:
19	!---------
20	! Main subroutine of co2 cloud microphysics
21	!----------------------------------------------------------------------------------------------------------------------!
22	! Comments:
23	!----------
24	! - Adaptation of the water ice clouds scheme (with specific microphysics) of Montmessin, Navarro et al.
25	!
26	! - Microphysics subroutine is improvedCO2clouds.F
27	!
28	! - There is a time loop specific to cloud formation due to timescales smaller than the GCM integration timestep
29	!
30	! - The microphysics time step is a fraction of the physical one
31	!
32	! - The CO2 clouds tracers (co2_ice, ccn mass and concentration) are sedimented at each microtimestep. pdqs_sedco2
33	! keeps track of the CO2 flux at the surface
34	!
35	! - The subgrid Temperature distribution is modulated (0 or 1) by Spiga et al. (GRL 2012)
36	!
37	! - Saturation Index to account for GW propagation or dissipation upwards
38	!
39	! - 4D and column opacities are computed using Qext values at 1µm
40	!----------------------------------------------------------------------------------------------------------------------!
41	! Papers:
42	!--------
43	! "Near-pure vapor condensation in the Martian atmosphere: CO2 ice crystal growth", Listowski et al. (2013), JGRE
44	! "Modeling the microphysics of CO2 ice clouds within wave-induced cold pockets in the martian mesosphere", Listowski
45	! et al. (2014), Icarus
46	! "Global climate modeling of the Martian water cycle with improved microphysics and radiatively active water ice
47	! clouds", Navarro et al. (2014), JGRE
48	! "Martian GCM with complete CO2 clouds microphysics", Audouard et al. (2017), EPSC abstract
49	!----------------------------------------------------------------------------------------------------------------------!
50	! Algorithm:
51	!-----------
52	! 0. Firstcall
53	! 0.1. Initialization of microtimestep from imicroco2
54	! 0.2. Compute the radius grid of CO2 ice particles (rb_cldco2)
55	! 0.3. Read file 'optprop_co2ice_1mic.dat' to extract optical properties of CO2 ice at 1 micron (Qext)
56	! 0.4. Interpole the radius grid (rb_cldco2) to get the corresponding exctinction coefficients (Qext1bins)
57	! 0.5. Save the radius grid of CO2 particles (rb_cldco2)
58	! 1. Initialization
59	! 2. Compute mass and thickness layers
60	! 3. Define the sub-grid cloud (CLFvaryingCO2)
61	! 3.1. Representation of sub-grid CO2 ice clouds (CLFvaryingCO2 = True)
62	! 3.1.a. Saturation index CO2
63	! 3.1.b. Compute tcond
64	! 3.1.c. Compute cloud fraction in cells
65	! 3.2. No sub-grid cloud representation (CLFvaryingCO2 = False)
66	! 4. Microphysics of CO2 cloud formation
67	! 4.1. Effective tracers quantities in the cloud
68	! 4.2. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry
69	! 4.3. Main call to the cloud scheme
70	! 4.4. Updating tendencies after cloud scheme
71	! 4.5. Gravitational sedimentation
72	! 4.5.a. Compute cloud density
73	! 4.5.b. Save actualized tracer values to compute sedimentation tendancies
74	! 4.5.c. Sedimentation of co2 ice
75	! 4.5.d. Sedimentation for other tracers
76	! 4.5.e. Compute tendencies due to the sedimation process
77	! 5. Compute final tendencies after time loop
78	! 6. Update clouds physical values in the cloud (for output)
79	! 6.1. Update density of co2 ice, riceco2 and opacity
80	! 6.2. Update rice and rdust
81	! 7. Correction if a lot of subliming CO2 fills the 1st layer
82	! 8. CO2 saturated quantities
83	! 8.1 Compute CO2 saturation in layers
84	! 8.2 Compute CO2 saturated quantities in layers
85	! 9. Everything modified by CO2 microphysics must be wrt co2cloudfrac
86	! 10. Compute opacity at 1 micron
87	! 11. Write outputs in diagfi.nc
88	!======================================================================================================================!
89	subroutine co2cloud(ngrid, nlay, ptimestep, pplev, pplay, pdpsrf, pzlay, pt, pdt, pq, pdq, pdqcloudco2, pdtcloudco2, &
90	nq, tau, tauscaling, rdust, rice, riceco2, nuice, rhocloud, rsedcloudco2, rhocloudco2, pzlev,&
91	pdqs_sedco2, pdqs_sedccn, pdu, pu, pcondicea)
92
93	use ioipsl_getincom, only: getin
94	use dimradmars_mod, only: naerkind
95	use comcstfi_h, only: pi, g, cpp
96	use updaterad, only: updaterice_microco2, updaterice_micro, updaterdust
97	use conc_mod, only: mmean, rnew
98	use tracer_mod, only: igcm_co2, igcm_co2_ice, igcm_dust_mass, igcm_dust_number, igcm_h2o_ice, &
99	igcm_ccn_mass, igcm_ccn_number, igcm_ccnco2_mass, igcm_ccnco2_number, &
100	igcm_ccnco2_h2o_number, igcm_ccnco2_h2o_mass_ice, igcm_ccnco2_h2o_mass_ccn, rho_dust, &
101	nuiceco2_sed, nuiceco2_ref, r3n_q, rho_ice, nuice_sed, igcm_ccnco2_meteor_mass, &
102	igcm_ccnco2_meteor_number
103
104	use newsedim_mod, only: newsedim
105
106	use datafile_mod, only: datadir
107	use density_co2_ice_mod, only: density_co2_ice
108
109	use improvedCO2clouds_mod, only: improvedCO2clouds
110	use microphys_h, only: nbinco2_cld, rad_cldco2, mco2
111	use write_output_mod, only: write_output
112	#ifndef MESOSCALE
113	use vertical_layers_mod, only: ap, bp
114	#endif
115
116	implicit none
117
118	include "callkeys.h"
119	!----------------------------------------------------------------------------------------------------------------------!
120	! VARIABLES DECLARATION
121	!----------------------------------------------------------------------------------------------------------------------!
122	! Input arguments:
123	!-----------------
124	integer, intent(in) ::&
125	ngrid, &! Number of grid points
126	nlay, &! Number of layers
127	nq ! Number of tracers
128
129	real, intent(in) :: &
130	ptimestep, &! Physical time step (s)
131	pplev(ngrid,nlay+1), &! Inter-layer pressures (Pa)
132	pplay(ngrid,nlay), &! Mid-layer pressures (Pa)
133	pdpsrf(ngrid), &! Tendency on surface pressure
134	pzlay(ngrid,nlay), &! Altitude at the middle of the layers
135	pt(ngrid,nlay), &! Temperature at the middle of the layers (K)
136	pdt(ngrid,nlay), &! Tendency on temperature from other parametrizations
137	pq(ngrid,nlay,nq), &! Tracers (kg/kg)
138	pdq(ngrid,nlay,nq), &! Tendencies before condensation (kg/kg.s-1)
139	tau(ngrid,naerkind), &! Column dust optical depth at each point
140	tauscaling(ngrid), &! Convertion factor for dust amount
141	pu(ngrid,nlay), &! Zonal Wind: zu = pu + (pdu * ptimestep)
142	pdu(ngrid,nlay), &! Tendency of zonal wind before condensation
143	pzlev(ngrid,nlay+1), &! Altitude at the boundaries of the layers
144	nuice(ngrid,nlay) ! Estimated effective variance of the size distribution
145	!----------------------------------------------------------------------------------------------------------------------!
146	! Output arguments:
147	!------------------
148	real, intent(out) :: &
149	rice(ngrid,nlay), & ! Water Ice mass mean radius (m)
150	! rsedcloud(ngrid,nlay), & ! Water Cloud sedimentation radius
151	rhocloud(ngrid,nlay), & ! Water Cloud density (kg.m-3)
152	pdqs_sedco2(ngrid), & ! CO2 flux at the surface
153	pdqs_sedccn(ngrid,nq), & ! CCN flux at the surface
154	pdqcloudco2(ngrid,nlay,nq),& ! Tendency due to CO2 condensation (kg/kg.s-1)
155	pcondicea(ngrid,nlay), & ! Rate of condensation/sublimation of co2 ice in layers
156	pdtcloudco2(ngrid,nlay), & ! Tendency on temperature due to latent heat
157	rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius
158
159	real, intent(inout) :: &
160	rdust(ngrid,nlay) ! Dust geometric mean radius (m)
161
162	double precision, intent(out) :: &
163	riceco2(ngrid,nlay) ! Ice mass mean radius (m) r_c in Montmessin et al. (2004)
164	!----------------------------------------------------------------------------------------------------------------------!
165	! Local:
166	!-------
167	!-----1) Parameters:
168	!-------------------
169	integer, parameter :: &
170	uQext = 555, &! file_qext unit ID
171	var_dim_qext = 10000 ! Exact dimension of radv and qextv1mic from file_qext
172
173	real, parameter :: &
174	mincloud = 0.1, &! Minimum cloud fraction
175	beta = 0.85, &! correction for the shape of the particles (see Murphy et al. JGR 1990 vol.95):
176	! beta = 1 for spheres
177	! beta = 0.85 for irregular particles
178	! beta = 0.5 for disk shaped particles
179	threshold = 1e-30, & ! limit value considering as zero
180	threshold_2 = 1e-13 ! limit value considering the value is physical (below this value => computer noise for dble)
181
182	double precision, parameter :: &
183	rmin_cld = 1.e-9, &! Minimum radius (m)
184	rmax_cld = 5.e-6, &! Maximum radius (m)
185	rbmin_cld = 1.e-10,&! Minimum boundary radius (m)
186	rbmax_cld = 2.e-4, &! Maximum boundary radius (m)
187	Fo = 7.5e-7, &! for sat index (J.m-3)
188	lambdaH = 150.e3 ! for sat index (km)
189
190	character(len=23), parameter :: &
191	file_qext = 'optprop_co2ice_1mic.dat' ! File extinction coefficients of CO2 particles
192	!----------------------------------------------------------------------------------------------------------------------!
193	!-----2) Saved:
194	!--------------
195	integer, save :: &
196	imicroco2 ! Time subsampling for coupled water microphysics sedimentation microtimestep timeloop for microphysics:
197	! if imicroco2 = 1, subpdt is the same as pdt
198	real, save :: &
199	sigma_iceco2, &! Variance of the ice and CCN distributions
200	microtimestep ! Integration timestep for coupled water microphysics & sedimentation
201
202	double precision, save :: &
203	dev2, &! 1. / ( sqrt(2.) * sigma_iceco2 )
204	Qext1bins(nbinco2_cld), &! Extinction coefficients for rb_cldco2 radius of CO2 ice particles
205	Qextv1mic(var_dim_qext), &
206	radv(var_dim_qext), & ! radius of CO2 ice at 1 µm (read from file_qext)
207	rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
208	logical, save :: &
209	firstcall = .true. ! Used to compute saved variables
210
211	!$OMP THREADPRIVATE(imicroco2,sigma_iceco2,microtimestep)
212	!$OMP THREADPRIVATE(dev2,Qext1bins,Qextv1mic,radv,rb_cldco2)
213	!$OMP THREADPRIVATE(firstcall)
214
215	!----------------------------------------------------------------------------------------------------------------------!
216	!-----3) Variables:
217	!------------------
218	integer :: &
219	iq, &! loop on tracers
220	ig, &! loop on grid points
221	l, &! loop on layers
222	i, &! loop on nbinco2_cld
223	nelem, &! number of point between lebon1 and lebon2 => interpolation
224	lebon1, &! bound limit for the interpolation
225	lebon2, &! bound limit for the interpolation
226	microstep ! Time subsampling step variable
227
228	real :: &
229	! ---Tendency given by clouds inside the micro loop
230	subpdqcloudco2(ngrid,nlay,nq), &! On tracers, cf. pdqcloud
231	subpdtcloudco2(ngrid,nlay), &! On temperature, cf. pdtcloud
232	! ---Global tendency (clouds+physics)
233	sum_subpdq(ngrid,nlay,nq), &! On tracers, cf. pdqcloud
234	sum_subpdt(ngrid,nlay), &! On temperature, cf. pdtcloud
235	! ---Sedimentation
236	ztsed(ngrid,nlay), &! Temperature with real-time value in microtimeloop
237	zqsed(ngrid,nlay,nq), &! Tracers with real-time value in µloop
238	zqsed0(ngrid,nlay,nq), &! For sedimentation tendancy
239	subpdqsed(ngrid,nlay,nq), &! Tendancy due to sedimentation
240	sum_subpdqs_sedco2(ngrid), &! CO2 flux at the surface
241	sum_subpdqs_sedccn(ngrid,nq), &!
242	! ---For sub grid T distribution
243	zt(ngrid,nlay), &! Local value of temperature
244	zq_co2vap(ngrid, nlay), &! Local value of CO2 vap
245	rhocloudco2t(ngrid, nlay), &! Cloud density (kg.m-3)
246	! ---For Saturation Index computation
247	zdelt, &! Delta T for the temperature distribution
248	co2cloudfrac(ngrid,nlay), &! Cloud fraction used only with CLFvarying is true
249	! ---Misc
250	rhocloudco2(ngrid, nlay), &! Cloud density (kg.m-3)
251	Nccnco2, &! buffer: number of ccn used for co2 condensation
252	Qccnco2, &! buffer: mass of ccn used for co2 condensation
253	Nccnco2_h2o, &
254	Qccnco2_h2o, &
255	Niceco2, &! buffer: mmr co2 ice
256	epaisseur(ngrid,nlay), &! Layer thickness (m)
257	masse(ngrid,nlay), &! Layer mass (kg.m-2)
258	pteff(ngrid, nlay), &! Effective temperature in the cloud
259	pqeff(ngrid, nlay, nq), &! Effective tracers quantities in the cloud
260	wq(ngrid,nlay+1), &! Displaced tracer mass (kg.m-2) during microtimestep
261	satuco2(ngrid,nlay), &! CO2 satu ratio for output diagfi
262	zqsatco2(ngrid,nlay), &! Saturation co2
263	availco2,&
264	masslayer, &
265	tmp, a,b, &
266	new_pdq(ngrid,nlay)
267
268	double precision :: &
269	! ---Extinction coefficients at 1 micron of CO2 particles
270	vrat_cld, &! Volume ratio
271	n_derf, &! derf( (rb_cldco2(1)-log(riceco2(ig,l))) *dev2)
272	Qtemp, &! mean value in the interval during the interpolation
273	ltemp1(var_dim_qext), &! abs(radv(:)-rb_cldco2(i))
274	ltemp2(var_dim_qext), &! abs(radv(:)-rb_cldco2(i+1))
275	n_aer(nbinco2_cld), &! -0.5 * Nccnco2tauscaling(ig) n_derf
276	tau1mic(ngrid), &! CO2 ice column opacity at 1µm
277	Qext1bins2(ngrid,nlay), &! CO2 ice opacities
278	! ---For Saturation Index computation
279	rho, &! background density
280	zu, &! absolute value of zonal wind field
281	NN, &! N^2 static stability
282	gradT, &! thermal gradient
283	SatIndex(ngrid,nlay), &! Saturation index S in Spiga 2012 paper, assuming like in the paper GW phase speed
284	! (stationary waves): c = 0 m.s-1, lambdaH = 150 km, Fo = 7.5e-7 J.m-3
285	SatIndexmap(ngrid), &! maxval(SatIndex(ig,12:26))
286	! ---Misc
287	myT, &! Temperature scalar for co2 density computation
288	tcond(ngrid,nlay) ! CO2 condensation temperature
289
290	logical :: &
291	file_qext_ok ! Check if file_qext exists
292	!======================================================================================================================!
293	! BEGIN ===============================================================================================================!
294	!======================================================================================================================!
295	! 0. Firstcall
296	!----------------------------------------------------------------------------------------------------------------------!
297	if (firstcall) then
298	firstcall=.false.
299	!----------------------------------------------------------------------------------------------------------------------!
300	! 0.1. Initialization of microtimestep from imicroco2
301	!----------------------------------------------------------------------------------------------------------------------!
302	#ifdef MESOSCALE
303	imicroco2 = 2
304	#else
305	imicroco2 = 30
306	#endif
307	call getin("imicroco2", imicroco2)
308	microtimestep = ptimestep/real(imicroco2)
309	sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
310	dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
311	!----------------------------------------------------------------------------------------------------------------------!
312	! 0.2. Compute the radius grid of CO2 ice particles (rb_cldco2)
313	! > the grid spacing is computed assuming a constant volume ratio between two consecutive bins; i.e. vrat_cld.
314	! - rad_cldco2 is the primary radius grid used for microphysics computation.
315	! - The grid spacing is computed assuming a constant volume ratio between two consecutive bins; i.e. vrat_cld.
316	! - vrat_cld is determined from the boundary values of the size grid: rmin_cld and rmax_cld.
317	! - The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
318	!----------------------------------------------------------------------------------------------------------------------!
319	! vrat_cld is determined from the boundary values of the size grid: rmin_cld and rmax_cld.
320	vrat_cld = exp(log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) * 3.)
321
322	! rad_cldco2 is the primary radius grid used for microphysics computation.
323	rad_cldco2(1) = rmin_cld
324	do i = 1, nbinco2_cld-1
325	rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.)
326	end do
327
328	! rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
329	rb_cldco2(1) = rbmin_cld
330	do i = 1, nbinco2_cld
331	rb_cldco2(i+1) = ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.) rad_cldco2(i)
332	end do
333	rb_cldco2(nbinco2_cld+1) = rbmax_cld
334	!----------------------------------------------------------------------------------------------------------------------!
335	! 0.3. Read file 'optprop_co2ice_1mic.dat' to extract optical properties of CO2 ice at 1 micron (Qext)
336	!----------------------------------------------------------------------------------------------------------------------!
337	! get information about file_qext
338	inquire(file=trim(datadir)//'/'//file_qext, exist=file_qext_ok)
339
340	! if file_qext is missing then stop
341	if (.not. file_qext_ok) then
342	write(,)'file'//file_qext//'should be in ', trim(datadir)
343	call abort_physic('co2cloud', 'file missing', 1)
344	end if
345
346	! read file_qext
347	open(unit=uQext,file=trim(datadir)//'/'//file_qext, form='formatted')
348
349	! skip 1 line
350	read(uQext,*)
351
352	! extract radv
353	do i = 1, var_dim_qext
354	read(uQext,'(E12.5)')radv(i)
355	end do
356
357	! skip 1 line
358	read(uQext,*)
359
360	! Qextv1mic
361	do i = 1 , var_dim_qext
362	read(uQext,'(E12.5)')Qextv1mic(i)
363	end do
364
365	! close file_qext
366	close(uQext)
367	!----------------------------------------------------------------------------------------------------------------------!
368	! 0.4. Interpole the radius grid (rb_cldco2) to get the corresponding exctinction coefficients (Qext1bins), using
369	! file_qext values (radv, Qextv1mic)
370	!----------------------------------------------------------------------------------------------------------------------!
371	do i = 1, nbinco2_cld
372	ltemp1 = abs(radv(:)-rb_cldco2(i))
373	ltemp2 = abs(radv(:)-rb_cldco2(i+1))
374	lebon1 = minloc(ltemp1,DIM=1)
375	lebon2 = min(minloc(ltemp2,DIM=1), var_dim_qext)
376	nelem = lebon2 - lebon1 + 1.
377
378	! mean value in the interval
379	Qtemp = 0d0
380	do l = 0, nelem
381	Qtemp = Qtemp + Qextv1mic(min(lebon1+l, var_dim_qext))
382	end do
383
384	Qext1bins(i) = Qtemp / nelem
385	end do
386
387	Qext1bins(:) = Qext1bins(:) * pi * (rad_cldco2(:)**2)
388
389	! print result of the interpolation
390	write(,)'--------------------------------------------'
391	write(,)'Microphysics co2: size bin-Qext information:'
392	write(,)' i, rad_cldco2(i), Qext1bins(i)'
393	do i = 1, nbinco2_cld
394	write(*,'(i3,3x,3(e13.6,4x))')i, rad_cldco2(i), Qext1bins(i)
395	end do
396	write(,)'--------------------------------------------'
397	!----------------------------------------------------------------------------------------------------------------------!
398	! 0.5. Save the radius grid of CO2 particles (rb_cldco2)
399	!----------------------------------------------------------------------------------------------------------------------!
400	do i = 1, nbinco2_cld+1
401	rb_cldco2(i) = log(rb_cldco2(i))
402	end do
403	end if ! of IF (firstcall)
404	!----------------------------------------------------------------------------------------------------------------------!
405	! 1. Initialization
406	!----------------------------------------------------------------------------------------------------------------------!
407	sum_subpdq(1:ngrid,1:nlay,1:nq) = 0.
408	sum_subpdt(1:ngrid,1:nlay) = 0.
409
410	subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0.
411	subpdtcloudco2(1:ngrid,1:nlay) = 0.
412
413	pdqcloudco2(1:ngrid,1:nlay,1:nq) = 0.
414	pdtcloudco2(1:ngrid,1:nlay) = 0.
415
416	! default value if no ice
417	rhocloudco2(1:ngrid,1:nlay) = rho_dust
418	rhocloudco2t(1:ngrid,1:nlay) = rho_dust
419	epaisseur(1:ngrid,1:nlay) = 0.
420	masse(1:ngrid,1:nlay) = 0.
421	riceco2(1:ngrid, 1:nlay) = 0.
422	zqsed0(1:ngrid,1:nlay,1:nq) = 0.
423	sum_subpdqs_sedco2(1:ngrid) = 0.
424	sum_subpdqs_sedccn(1:ngrid,1:nq) = 0.
425	subpdqsed(1:ngrid,1:nlay,1:nq) = 0.
426	!----------------------------------------------------------------------------------------------------------------------!
427	! 2. Compute mass and thickness layers
428	!----------------------------------------------------------------------------------------------------------------------!
429	do l = 1, nlay
430	do ig = 1, ngrid
431	#ifndef MESOSCALE
432	masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1) + (bp(l)-bp(l+1)) ) / g
433	#else
434	masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1)) / g
435	#endif
436	epaisseur(ig,l) = pzlev(ig,l+1) - pzlev(ig,l)
437	end do
438	end do
439	!----------------------------------------------------------------------------------------------------------------------!
440	! 3. Define the sub-grid cloud (CLFvaryingCO2)
441	!----------------------------------------------------------------------------------------------------------------------!
442	! 3.1. Representation of sub-grid CO2 ice clouds (CLFvaryingCO2 = True)
443	!----------------------------------------------------------------------------------------------------------------------!
444	if (CLFvaryingCO2) then
445	! effective temperature
446	pteff(:,:) = pt(:,:)
447
448	! min co2cloudfrac when there is ice
449	co2cloudfrac(:,:) = mincloud
450
451	! temperature
452	do l=1,nlay
453	do ig=1,ngrid
454	zt(ig,l) = pt(ig,l) + pdt(ig,l)*ptimestep
455	end do
456	end do
457
458	! Quantities of traceurs
459	if (igcm_co2 /= 0) then
460	do l = 1, nlay
461	do ig = 1, ngrid
462	zq_co2vap(ig,l) = pq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2)*ptimestep
463	end do
464	end do
465	end if
466	!----------------------------------------------------------------------------------------------------------------------!
467	! 3.1.a. Saturation index CO2
468	!----------------------------------------------------------------------------------------------------------------------!
469	! if saturation index co2 is true
470	if (satindexco2) then
471	! layers 12 --> 26 ~ 12->85 km
472	do l = 12, 26
473	do ig = 1, ngrid
474	! compute N^2 static stability
475	gradT = (zt(ig,l+1)-zt(ig,l))/(pzlev(ig,l+1)-pzlev(ig,l))
476	NN = sqrt(g/zt(iq,l) * (max(gradT,-g/cpp) + g/cpp))
477
478	! compute absolute value of zonal wind field
479	zu = abs(pu(ig,l) + pdu(ig,l)*ptimestep)
480
481	! compute background density
482	rho = pplay(ig,l) / (rnew(ig,l)*zt(ig,l))
483
484	! saturation index: Modulate the DeltaT by GW propagation index:
485	! --------------------------------------------------------------
486	SatIndex(ig,l) = sqrt(FolambdaH/(2.pi)NN / (rhozu**3) )
487	end do
488	end do
489
490	! Then compute Satindex map in layers 12 --> 26 ~ 12->85 km
491	do ig = 1, ngrid
492	SatIndexmap(ig) = maxval(SatIndex(ig,12:26))
493	end do
494
495	! Write outputs in diagfi.nc
496	call write_output( "SatIndex", "SatIndex", " ", SatIndex(:,:))
497
498	call write_output( "SatIndexmap", "SatIndexmap", "km", SatIndexmap(:))
499	!------------------------------------------------------------------------------------------------------------------!
500	! if saturation index co2 is false, set saturation index to 0.05
501	!------------------------------------------------------------------------------------------------------------------!
502	else
503	do ig = 1, ngrid
504	SatIndexmap(ig)=0.05
505	end do
506	end if ! of if (satindexco2)
507	!----------------------------------------------------------------------------------------------------------------------!
508	! 3.1.b. Compute tcond
509	!----------------------------------------------------------------------------------------------------------------------!
510	call tcondco2(ngrid,nlay,pplay,zq_co2vap,tcond)
511	!----------------------------------------------------------------------------------------------------------------------!
512	! 3.1.c. Compute cloud fraction in cells
513	!----------------------------------------------------------------------------------------------------------------------!
514	do ig = 1, ngrid
515	if (SatIndexmap(ig) <= 0.1) then
516	do l = 1, nlay-1
517
518	! The entire fraction is saturated
519	if (tcond(ig,l) >= (zt(ig,l)+zdelt) .or. tcond(ig,l) <= 0.) then
520	pteff(ig,l) = zt(ig,l)
521	co2cloudfrac(ig,l) = 1.
522
523	! No saturation at all
524	else if (tcond(ig,l) <= (zt(ig,l)-zdelt)) then
525	pteff(ig,l) = zt(ig,l) - zdelt
526	co2cloudfrac(ig,l) = mincloud
527
528	! Mean temperature of the cloud fraction
529	else
530	pteff(ig,l) = (tcond(ig,l)+zt(ig,l)-zdelt) / 2.
531	co2cloudfrac(ig,l) = (tcond(ig,l)-zt(ig,l)+zdelt) / (2.0*zdelt)
532	end if
533
534	pteff(ig,l) = pteff(ig,l) - pdt(ig,l)*ptimestep
535
536	! check boundary values of co2cloudfrac
537	if (co2cloudfrac(ig,l) <= mincloud) then
538	co2cloudfrac(ig,l) = mincloud
539	else if (co2cloudfrac(ig,l)> 1) then
540	co2cloudfrac(ig,l) = 1.
541	end if
542	end do
543
544	! SatIndex not favorable for GW: leave pt untouched
545	else
546	pteff(ig,l) = pt(ig,l)
547	co2cloudfrac(ig,l) = mincloud
548	end if ! of if (SatIndexmap <= 0.1)
549	end do ! of ngrid
550	! TODO: Totalcloud frac of the column missing here
551	!----------------------------------------------------------------------------------------------------------------------!
552	! 3.2. No sub-grid cloud representation (CLFvaryingCO2 = False)
553	!----------------------------------------------------------------------------------------------------------------------!
554	else
555	do l = 1, nlay
556	do ig = 1, ngrid
557	pteff(ig,l) = pt(ig,l)
558	end do
559	end do
560	end if ! end if (CLFvaryingco2)
561	!----------------------------------------------------------------------------------------------------------------------!
562	! 4. Microphysics of CO2 cloud formation
563	!----------------------------------------------------------------------------------------------------------------------!
564	pqeff(:,:,:) = pq(:,:,:)
565	pteff(:,:) = pt(:,:)
566	!----------------------------------------------------------------------------------------------------------------------!
567	! 4.1. Effective tracers quantities in the cloud
568	!----------------------------------------------------------------------------------------------------------------------!
569	if (CLFvaryingCO2) then
570	pqeff(:,:,igcm_ccnco2_mass) = pq(:,:,igcm_ccnco2_mass) / co2cloudfrac(:,:)
571
572	pqeff(:,:,igcm_ccnco2_number) = pq(:,:,igcm_ccnco2_number) / co2cloudfrac(:,:)
573
574	pqeff(:,:,igcm_co2_ice) = pq(:,:,igcm_co2_ice) / co2cloudfrac(:,:)
575	end if
576	!----------------------------------------------------------------------------------------------------------------------!
577	do microstep = 1, imicroco2
578	!----------------------------------------------------------------------------------------------------------------------!
579	! 4.2. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry
580	!----------------------------------------------------------------------------------------------------------------------!
581	do l = 1, nlay
582	do ig = 1, ngrid
583	! on temperature
584	sum_subpdt(ig,l) = sum_subpdt(ig,l) + pdt(ig,l)
585
586	! on tracers
587	sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + pdq(ig,l,igcm_dust_mass)
588
589	sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + pdq(ig,l,igcm_dust_number)
590
591	sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + pdq(ig,l,igcm_ccnco2_mass)
592
593	sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + pdq(ig,l,igcm_ccnco2_number)
594
595	sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + pdq(ig,l,igcm_co2_ice)
596
597	sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + pdq(ig,l,igcm_co2)
598
599	if (meteo_flux) then
600	sum_subpdq(ig,l,igcm_ccnco2_meteor_number) = sum_subpdq(ig,l,igcm_ccnco2_meteor_number) + &
601	pdq(ig,l,igcm_ccnco2_meteor_number)
602
603	sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) = sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) + &
604	pdq(ig,l,igcm_ccnco2_meteor_mass)
605	end if
606	if (co2useh2o) then
607	sum_subpdq(ig,l,igcm_h2o_ice) = sum_subpdq(ig,l,igcm_h2o_ice) + pdq(ig,l,igcm_h2o_ice)
608
609	sum_subpdq(ig,l,igcm_ccn_mass) = sum_subpdq(ig,l,igcm_ccn_mass) + pdq(ig,l,igcm_ccn_mass)
610
611	sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + pdq(ig,l,igcm_ccn_number)
612
613	sum_subpdq(ig,l,igcm_ccnco2_h2o_number) = sum_subpdq(ig,l,igcm_ccnco2_h2o_number) + &
614	pdq(ig,l,igcm_ccnco2_h2o_number)
615
616	sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) + &
617	pdq(ig,l,igcm_ccnco2_h2o_mass_ice)
618
619	sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) + &
620	pdq(ig,l,igcm_ccnco2_h2o_mass_ccn)
621	end if
622	end do ! ngrid
623	end do ! nlay
624	!----------------------------------------------------------------------------------------------------------------------!
625	! 4.3. Main call to the cloud scheme
626	!----------------------------------------------------------------------------------------------------------------------!
627	call improvedco2clouds(ngrid, nlay, microtimestep, pplay, pplev, pteff, sum_subpdt, pqeff, sum_subpdq, &
628	subpdqcloudco2, subpdtcloudco2, nq, tauscaling, rb_cldco2, sigma_iceco2, dev2)
629
630	do l = 1, nlay
631	do ig = 1, ngrid
632	if(pq(ig,l,igcm_co2_ice) + microtimestep*(sum_subpdq(ig,l,igcm_co2_ice)+subpdqcloudco2(ig,l,igcm_co2_ice)) &
633	<= 1.e-12) then
634	subpdqcloudco2(ig,l,igcm_co2_ice) = -pq(ig,l,igcm_co2_ice)/microtimestep - sum_subpdq(ig,l,igcm_co2_ice)
635	subpdqcloudco2(ig,l,igcm_co2) = -subpdqcloudco2(ig,l,igcm_co2_ice)
636	end if
637
638	if(pq(ig,l,igcm_co2) + microtimestep*(sum_subpdq(ig,l,igcm_co2)+subpdqcloudco2(ig,l,igcm_co2)) <= 1.e-12) then
639	subpdqcloudco2(ig,l,igcm_co2) = - pq(ig,l,igcm_co2)/microtimestep - sum_subpdq(ig,l,igcm_co2)
640	subpdqcloudco2(ig,l,igcm_co2_ice) = -subpdqcloudco2(ig,l,igcm_co2)
641	end if
642
643	! ccnco2_number and ccnco2_mass
644	if (((pq(ig,l,igcm_ccnco2_number)+(sum_subpdq(ig,l,igcm_ccnco2_number)+subpdqcloudco2(ig,l,igcm_ccnco2_number)) &
645	*microtimestep)<=1.) .or. &
646	(pq(ig,l,igcm_ccnco2_mass)+(sum_subpdq(ig,l,igcm_ccnco2_mass)+subpdqcloudco2(ig,l,igcm_ccnco2_mass)) &
647	*microtimestep<=1e-20)) then
648	subpdqcloudco2(ig,l,igcm_ccnco2_number) = - pq(ig,l,igcm_ccnco2_number)/microtimestep + 1. &
649	- sum_subpdq(ig,l,igcm_ccnco2_number)
650	subpdqcloudco2(ig,l,igcm_dust_number) = - subpdqcloudco2(ig,l,igcm_ccnco2_number)
651
652	subpdqcloudco2(ig,l,igcm_ccnco2_mass) = - pq(ig,l,igcm_ccnco2_mass)/microtimestep + 1e-20 &
653	- sum_subpdq(ig,l,igcm_ccnco2_mass)
654	subpdqcloudco2(ig,l,igcm_dust_mass) = - subpdqcloudco2(ig,l,igcm_ccnco2_mass)
655	end if
656
657	! ccnco2_meteor_number and ccnco2_meteor_mass
658	if (meteo_flux) then
659	if (((pq(ig,l,igcm_ccnco2_meteor_number)+(sum_subpdq(ig,l,igcm_ccnco2_meteor_number)+ &
660	subpdqcloudco2(ig,l,igcm_ccnco2_meteor_number))*microtimestep)<=1.) .or. &
661	(pq(ig,l,igcm_ccnco2_meteor_mass)+(sum_subpdq(ig,l,igcm_ccnco2_meteor_mass)+ &
662	subpdqcloudco2(ig,l,igcm_ccnco2_meteor_mass))*microtimestep<=1e-20)) then
663	subpdqcloudco2(ig,l,igcm_ccnco2_meteor_number) = - pq(ig,l,igcm_ccnco2_meteor_number)/microtimestep + 1. &
664	- sum_subpdq(ig,l,igcm_ccnco2_meteor_number)
665	subpdqcloudco2(ig,l,igcm_ccnco2_meteor_mass) = - pq(ig,l,igcm_ccnco2_meteor_mass)/microtimestep + 1e-20 &
666	- sum_subpdq(ig,l,igcm_ccnco2_meteor_mass)
667	end if
668	end if
669
670	! ccnco2_h2o_number and masses
671	if (co2useh2o) then
672	if (((pq(ig,l,igcm_ccnco2_h2o_number) + (sum_subpdq(ig,l,igcm_ccnco2_h2o_number) + &
673	subpdqcloudco2(ig,l,igcm_ccnco2_h2o_number)) &
674	*microtimestep)<=1.) .or. &
675	(pq(ig,l,igcm_ccnco2_h2o_mass_ice)+pq(ig,l,igcm_ccnco2_h2o_mass_ccn) +&
676	(sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice)+subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) +&
677	sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn)+subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn)) &
678	*microtimestep<=1e-20)) then
679
680	subpdqcloudco2(ig,l,igcm_ccnco2_h2o_number) = - pq(ig,l,igcm_ccnco2_h2o_number)/microtimestep + 1. &
681	- sum_subpdq(ig,l,igcm_ccnco2_h2o_number)
682	subpdqcloudco2(ig,l,igcm_ccn_number) = - subpdqcloudco2(ig,l,igcm_ccnco2_h2o_number)
683
684	subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) = - pq(ig,l,igcm_ccnco2_h2o_mass_ice)/microtimestep + 1e-20 &
685	- sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice)
686
687	subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) = - pq(ig,l,igcm_ccnco2_h2o_mass_ccn)/microtimestep + 1e-20 &
688	- sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn)
689
690	subpdqcloudco2(ig,l,igcm_ccn_mass) = - subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) &
691	- subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn)
692	end if
693	end if
694	end do
695	end do
696	!----------------------------------------------------------------------------------------------------------------------!
697	! 4.4. Updating tendencies after cloud scheme
698	!----------------------------------------------------------------------------------------------------------------------!
699	do l = 1, nlay
700	do ig = 1, ngrid
701	sum_subpdt(ig,l) = sum_subpdt(ig,l) + subpdtcloudco2(ig,l)
702
703	sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + subpdqcloudco2(ig,l,igcm_dust_mass)
704
705	sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + subpdqcloudco2(ig,l,igcm_dust_number)
706
707	sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
708
709	sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + &
710	subpdqcloudco2(ig,l,igcm_ccnco2_number)
711
712	sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + subpdqcloudco2(ig,l,igcm_co2_ice)
713
714	sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + subpdqcloudco2(ig,l,igcm_co2)
715
716	if (meteo_flux) then
717	sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) = sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) + &
718	subpdqcloudco2(ig,l,igcm_ccnco2_meteor_mass)
719
720	sum_subpdq(ig,l,igcm_ccnco2_meteor_number) = sum_subpdq(ig,l,igcm_ccnco2_meteor_number) + &
721	subpdqcloudco2(ig,l,igcm_ccnco2_meteor_number)
722	end if
723	if (co2useh2o) then
724	sum_subpdq(ig,l,igcm_h2o_ice) = sum_subpdq(ig,l,igcm_h2o_ice) + subpdqcloudco2(ig,l,igcm_h2o_ice)
725
726	sum_subpdq(ig,l,igcm_ccn_mass) = sum_subpdq(ig,l,igcm_ccn_mass) + subpdqcloudco2(ig,l,igcm_ccn_mass)
727
728	sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + subpdqcloudco2(ig,l,igcm_ccn_number)
729
730	sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) + &
731	subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice)
732
733	sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) + &
734	subpdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn)
735
736	sum_subpdq(ig,l,igcm_ccnco2_h2o_number) = sum_subpdq(ig,l,igcm_ccnco2_h2o_number) + &
737	subpdqcloudco2(ig,l,igcm_ccnco2_number)
738
739	end if
740	end do ! ngrid
741	end do ! nlay
742	!----------------------------------------------------------------------------------------------------------------------!
743	! 4.5. Gravitational sedimentation
744	!----------------------------------------------------------------------------------------------------------------------!
745	if (sedimentation) then
746	!----------------------------------------------------------------------------------------------------------------------!
747	! 4.5.a. Compute cloud density
748	!----------------------------------------------------------------------------------------------------------------------!
749	do l = 1, nlay
750	do ig = 1, ngrid
751	! temperature during the sedimentation process
752	ztsed(ig,l) = pteff(ig,l) + sum_subpdt(ig,l) * microtimestep
753
754	! quantities tracers during the sedimentation process
755	zqsed(ig,l,:) = pqeff(ig,l,:) + sum_subpdq(ig,l,:) * microtimestep
756
757	! assure positive value of co2_ice mmr, ccnco2 number, ccnco2 mass
758	! meteoritic particle are considered like dust, rho_dust
759	Niceco2 = max(zqsed(ig,l,igcm_co2_ice), threshold)
760	Nccnco2 = max(zqsed(ig,l,igcm_ccnco2_number), threshold)
761	Qccnco2 = max(zqsed(ig,l,igcm_ccnco2_mass), threshold)
762
763	Nccnco2_h2o = 0.
764	Qccnco2_h2o = 0.
765	if (co2useh2o) then
766	Nccnco2_h2o = max(zqsed(ig,l,igcm_ccnco2_h2o_number), threshold)
767	Qccnco2_h2o = max(zqsed(ig,l,igcm_ccnco2_h2o_mass_ice)+zqsed(ig,l,igcm_ccnco2_h2o_mass_ccn), threshold)
768	Nccnco2 = Nccnco2 - Nccnco2_h2o
769	Qccnco2 = Qccnco2 - Qccnco2_h2o
770	if (Nccnco2 < 0.) then
771	Nccnco2 = threshold
772	Qccnco2 = threshold
773	end if
774	end if
775
776	! Get density cloud and co2 ice particle radius
777	if (Niceco2/=0d0) then
778	call updaterice_microco2(dble(Niceco2), dble(Qccnco2), dble(Nccnco2), dble(Qccnco2_h2o), dble(Nccnco2_h2o),&
779	ztsed(ig,l), tauscaling(ig), riceco2(ig,l), rhocloudco2t(ig,l))
780	else
781	riceco2(ig,l) = 0.
782	rhocloudco2t(ig,l) = 0.
783	end if
784	end do ! ngrid
785	end do ! nlay
786	!----------------------------------------------------------------------------------------------------------------------!
787	! 4.5.b. Save actualized tracer values to compute sedimentation tendancies
788	!----------------------------------------------------------------------------------------------------------------------!
789	zqsed0(:,:,igcm_co2_ice) = zqsed(:,:,igcm_co2_ice)
790	zqsed0(:,:,igcm_ccnco2_mass) = zqsed(:,:,igcm_ccnco2_mass)
791	zqsed0(:,:,igcm_ccnco2_number) = zqsed(:,:,igcm_ccnco2_number)
792
793	if (meteo_flux) then
794	zqsed0(:,:,igcm_ccnco2_meteor_mass) = zqsed(:,:,igcm_ccnco2_meteor_mass)
795	zqsed0(:,:,igcm_ccnco2_meteor_number) = zqsed(:,:,igcm_ccnco2_meteor_number)
796	end if
797
798	if (co2useh2o) then
799	zqsed0(:,:,igcm_ccnco2_h2o_number) = zqsed(:,:,igcm_ccnco2_h2o_number)
800	zqsed0(:,:,igcm_ccnco2_h2o_mass_ice) = zqsed(:,:,igcm_ccnco2_h2o_mass_ice)
801	zqsed0(:,:,igcm_ccnco2_h2o_mass_ccn) = zqsed(:,:,igcm_ccnco2_h2o_mass_ccn)
802	end if
803	!----------------------------------------------------------------------------------------------------------------------!
804	! 4.5.c. Sedimentation of co2 ice
805	!----------------------------------------------------------------------------------------------------------------------!
806	do ig = 1, ngrid
807	do l = 1, nlay
808	rsedcloudco2(ig,l) = max( riceco2(ig,l)(1.+nuiceco2_sed)(1.+nuiceco2_sed)*(1.+nuiceco2_sed), &
809	rdust(ig,l) )
810	end do
811	end do
812
813	wq(:,:) = 0.
814	call newsedim(ngrid, nlay, ngridnlay, ngridnlay, microtimestep, pplev, masse, epaisseur, ztsed, &
815	rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_co2_ice), wq, beta)
816
817	do ig = 1, ngrid
818	sum_subpdqs_sedco2(ig) = sum_subpdqs_sedco2(ig) + wq(ig,1) / microtimestep !wq est en kg.m-2
819	end do
820	!----------------------------------------------------------------------------------------------------------------------!
821	! 4.5.d. Sedimentation for other tracers
822	!----------------------------------------------------------------------------------------------------------------------!
823	wq(:,:) = 0.
824	! for ccnco2_mass
825	call newsedim(ngrid, nlay, ngridnlay, ngridnlay, microtimestep, pplev, masse, epaisseur, ztsed, &
826	rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_mass), wq, beta)
827	do ig = 1, ngrid
828	sum_subpdqs_sedccn(ig,igcm_ccnco2_mass) = sum_subpdqs_sedccn(ig,igcm_ccnco2_mass) + wq(ig,1) / microtimestep
829	end do
830
831	wq(:,:) = 0.
832	! for ccnco2_number
833	call newsedim(ngrid, nlay, ngridnlay, ngridnlay,microtimestep, pplev, masse, epaisseur, ztsed, &
834	rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_number), wq, beta)
835	do ig = 1, ngrid
836	sum_subpdqs_sedccn(ig,igcm_ccnco2_number) = sum_subpdqs_sedccn(ig,igcm_ccnco2_number) + wq(ig,1) / microtimestep
837	end do
838
839	if (meteo_flux) then
840	wq(:,:) = 0.
841	! for ccnco2_meteor_mass
842	call newsedim(ngrid, nlay, ngridnlay, ngridnlay, microtimestep, pplev, masse, epaisseur, ztsed, &
843	rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_meteor_mass), wq, beta)
844	do ig = 1, ngrid
845	sum_subpdqs_sedccn(ig,igcm_ccnco2_meteor_mass) = sum_subpdqs_sedccn(ig,igcm_ccnco2_meteor_mass) + &
846	wq(ig,1)/ microtimestep
847	end do
848
849	wq(:,:) = 0.
850	! for ccnco2_meteor_number
851	call newsedim(ngrid, nlay, ngridnlay, ngridnlay,microtimestep, pplev, masse, epaisseur, ztsed, &
852	rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_meteor_number), wq, beta)
853	do ig = 1, ngrid
854	sum_subpdqs_sedccn(ig,igcm_ccnco2_meteor_number) = sum_subpdqs_sedccn(ig,igcm_ccnco2_meteor_number) + &
855	wq(ig,1) / microtimestep
856	end do
857	end if
858	! for ccnco2_h2o_mass_ice
859	if (co2useh2o) then
860	wq(:,:) = 0.
861	call newsedim(ngrid, nlay, ngridnlay, ngridnlay,microtimestep, pplev, masse, epaisseur, ztsed, &
862	rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_h2o_mass_ice), wq, beta)
863	do ig = 1, ngrid
864	sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_mass_ice) = sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_mass_ice) + &
865	wq(ig,1) / microtimestep
866	end do
867
868	wq(:,:) = 0.
869	! for ccnco2_h2o_mass_ccn
870	call newsedim(ngrid, nlay, ngridnlay, ngridnlay,microtimestep, pplev, masse, epaisseur, ztsed, &
871	rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_h2o_mass_ccn), wq, beta)
872	do ig = 1, ngrid
873	sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_mass_ccn) = sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_mass_ccn) + &
874	wq(ig,1) / microtimestep
875	end do
876
877	wq(:,:) = 0.
878	! for ccnco2_h2o_number
879	call newsedim(ngrid, nlay, ngridnlay, ngridnlay,microtimestep, pplev, masse, epaisseur, ztsed, &
880	rsedcloudco2, rhocloudco2t, zqsed(:,:,igcm_ccnco2_h2o_number), wq, beta)
881	do ig = 1, ngrid
882	sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_number) = sum_subpdqs_sedccn(ig,igcm_ccnco2_h2o_number) + &
883	wq(ig,1) / microtimestep
884	end do
885	end if
886	!----------------------------------------------------------------------------------------------------------------------!
887	! 4.5.e. Compute tendencies due to the sedimation process
888	!----------------------------------------------------------------------------------------------------------------------!
889	do l = 1, nlay
890	do ig = 1, ngrid
891	subpdqsed(ig,l,igcm_ccnco2_mass) = ( zqsed(ig,l,igcm_ccnco2_mass) - zqsed0(ig,l,igcm_ccnco2_mass) ) &
892	/ microtimestep
893
894	subpdqsed(ig,l,igcm_ccnco2_number) = ( zqsed(ig,l,igcm_ccnco2_number) - zqsed0(ig,l,igcm_ccnco2_number) )&
895	/ microtimestep
896
897	subpdqsed(ig,l,igcm_co2_ice) = ( zqsed(ig,l,igcm_co2_ice) - zqsed0(ig,l,igcm_co2_ice) ) / microtimestep
898
899	if (meteo_flux) then
900	subpdqsed(ig,l,igcm_ccnco2_meteor_mass) = ( zqsed(ig,l,igcm_ccnco2_meteor_mass) - &
901	zqsed0(ig,l,igcm_ccnco2_meteor_mass) ) / microtimestep
902
903	subpdqsed(ig,l,igcm_ccnco2_meteor_number) = ( zqsed(ig,l,igcm_ccnco2_meteor_number) - &
904	zqsed0(ig,l,igcm_ccnco2_meteor_number) ) / microtimestep
905	end if
906	if (co2useh2o) then
907	subpdqsed(ig,l,igcm_ccnco2_h2o_number) = ( zqsed(ig,l,igcm_ccnco2_h2o_number) - &
908	zqsed0(ig,l,igcm_ccnco2_h2o_number) ) / microtimestep
909
910	subpdqsed(ig,l,igcm_ccnco2_h2o_mass_ice) = ( zqsed(ig,l,igcm_ccnco2_h2o_mass_ice) - &
911	zqsed0(ig,l,igcm_ccnco2_h2o_mass_ice) ) / microtimestep
912
913	subpdqsed(ig,l,igcm_ccnco2_h2o_mass_ccn) = ( zqsed(ig,l,igcm_ccnco2_h2o_mass_ccn) - &
914	zqsed0(ig,l,igcm_ccnco2_h2o_mass_ccn) ) / microtimestep
915	end if
916	end do
917	end do
918	! update subtimestep tendencies with sedimentation input
919	do l = 1, nlay
920	do ig = 1, ngrid
921	sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + subpdqsed(ig,l,igcm_ccnco2_mass)
922
923	sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + subpdqsed(ig,l,igcm_ccnco2_number)
924
925	sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + subpdqsed(ig,l,igcm_co2_ice)
926	if (meteo_flux) then
927	sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) = sum_subpdq(ig,l,igcm_ccnco2_meteor_mass) + &
928	subpdqsed(ig,l,igcm_ccnco2_meteor_mass)
929
930	sum_subpdq(ig,l,igcm_ccnco2_meteor_number) = sum_subpdq(ig,l,igcm_ccnco2_meteor_number) + &
931	subpdqsed(ig,l,igcm_ccnco2_meteor_number)
932	end if
933	if (co2useh2o) then
934	sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) + &
935	subpdqsed(ig,l,igcm_ccnco2_h2o_mass_ice)
936
937	sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) = sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) + &
938	subpdqsed(ig,l,igcm_ccnco2_h2o_mass_ccn)
939
940	sum_subpdq(ig,l,igcm_ccnco2_h2o_number) = sum_subpdq(ig,l,igcm_ccnco2_h2o_number) + &
941	subpdqsed(ig,l,igcm_ccnco2_h2o_number)
942	end if
943	end do
944	end do
945	end if !(end if sedimentation)
946	end do ! of do microstep = 1, imicroco2
947	!----------------------------------------------------------------------------------------------------------------------!
948	! 5. Compute final tendencies after time loop
949	!----------------------------------------------------------------------------------------------------------------------!
950	! condensation/sublimation rate in the atmosphere (kg.kg-1.s-1)
951	do l = nlay, 1, -1
952	do ig = 1, ngrid
953	pcondicea(ig,l) = sum_subpdq(ig,l,igcm_co2_ice) / real(imicroco2)
954	end do
955	end do
956
957	! CO2 flux at surface (kg.m-2.s-1)
958	do ig = 1, ngrid
959	pdqs_sedco2(ig) = sum_subpdqs_sedco2(ig) / real(imicroco2)
960	pdqs_sedccn(ig,:) = sum_subpdqs_sedccn(ig,:) / real(imicroco2)
961	end do
962
963	! temperature tendency (T.s-1)
964	do l = 1, nlay
965	do ig = 1, ngrid
966	pdtcloudco2(ig,l) = ( sum_subpdt(ig,l)/real(imicroco2) ) - pdt(ig,l)
967	if (pdtcloudco2(ig,l) /= pdtcloudco2(ig,l)) then
968	write(,)ig,l,pdtcloudco2(ig,l), pdt(ig,l), sum_subpdt(ig,l)
969	call abort_physic('co2clouds', 'ptdcloudco2 is NaN', 1)
970	end if
971	end do
972	end do
973
974	! tracers tendencies
975	do l = 1, nlay
976	do ig = 1, ngrid
977	pdqcloudco2(ig,l,igcm_co2) = 0. ! here is the trick, this tendency is computed in co2condens
978
979	pdqcloudco2(ig,l,igcm_co2_ice) = ( sum_subpdq(ig,l,igcm_co2_ice) / real(imicroco2) ) - pdq(ig,l,igcm_co2_ice)
980
981	pdqcloudco2(ig,l,igcm_ccnco2_mass) = ( sum_subpdq(ig,l,igcm_ccnco2_mass)/real(imicroco2) ) - &
982	pdq(ig,l,igcm_ccnco2_mass)
983
984	pdqcloudco2(ig,l,igcm_ccnco2_number) = ( sum_subpdq(ig,l,igcm_ccnco2_number) / real(imicroco2) ) - &
985	pdq(ig,l,igcm_ccnco2_number)
986
987	pdqcloudco2(ig,l,igcm_dust_mass) = ( sum_subpdq(ig,l,igcm_dust_mass) / real(imicroco2) ) - &
988	pdq(ig,l,igcm_dust_mass)
989
990	pdqcloudco2(ig,l,igcm_dust_number) = ( sum_subpdq(ig,l,igcm_dust_number)/real(imicroco2) ) - &
991	pdq(ig,l,igcm_dust_number)
992
993	if (meteo_flux) then
994	pdqcloudco2(ig,l,igcm_ccnco2_meteor_mass) = ( sum_subpdq(ig,l,igcm_ccnco2_meteor_mass)/real(imicroco2) ) &
995	- pdq(ig,l,igcm_ccnco2_meteor_mass)
996
997	pdqcloudco2(ig,l,igcm_ccnco2_meteor_number) = ( sum_subpdq(ig,l,igcm_ccnco2_meteor_number) / real(imicroco2) ) &
998	- pdq(ig,l,igcm_ccnco2_meteor_number)
999	end if
1000	if (co2useh2o) then
1001	pdqcloudco2(ig,l,igcm_h2o_ice) = ( sum_subpdq(ig,l,igcm_h2o_ice) / real(imicroco2) ) - &
1002	pdq(ig,l,igcm_h2o_ice)
1003
1004	pdqcloudco2(ig,l,igcm_ccn_mass) = ( sum_subpdq(ig,l,igcm_ccn_mass) / real(imicroco2) ) - &
1005	pdq(ig,l,igcm_ccn_mass)
1006
1007	pdqcloudco2(ig,l,igcm_ccn_number) = ( sum_subpdq(ig,l,igcm_ccn_number) / real(imicroco2) ) - &
1008	pdq(ig,l,igcm_ccn_number)
1009
1010	pdqcloudco2(ig,l,igcm_ccnco2_h2o_number) = ( sum_subpdq(ig,l,igcm_ccnco2_h2o_number) / real(imicroco2) ) - &
1011	pdq(ig,l,igcm_ccnco2_h2o_number)
1012
1013	pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) = ( sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ice) / real(imicroco2) ) - &
1014	pdq(ig,l,igcm_ccnco2_h2o_mass_ice)
1015
1016	pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) = ( sum_subpdq(ig,l,igcm_ccnco2_h2o_mass_ccn) / real(imicroco2) )- &
1017	pdq(ig,l,igcm_ccnco2_h2o_mass_ccn)
1018
1019	end if
1020	end do ! ngrid
1021	end do ! nlay
1022	!----------------------------------------------------------------------------------------------------------------------!
1023	! 6. Update clouds physical values in the cloud (for output)
1024	!----------------------------------------------------------------------------------------------------------------------!
1025	! 6.1. Update density of co2 ice, riceco2 and opacity
1026	!----------------------------------------------------------------------------------------------------------------------!
1027	do l = 1, nlay
1028	do ig = 1, ngrid
1029	Niceco2 = pqeff(ig,l,igcm_co2_ice) + ( pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice) ) * ptimestep
1030	Niceco2 = max(Niceco2, threshold)
1031
1032	! meteoritic particles are considered like dust, => rho_dust
1033	Nccnco2 = max( (pqeff(ig,l,igcm_ccnco2_number) + (pdq(ig,l,igcm_ccnco2_number) + &
1034	pdqcloudco2(ig,l, igcm_ccnco2_number)) * ptimestep) &
1035	, threshold)
1036
1037	Qccnco2 = max( (pqeff(ig,l,igcm_ccnco2_mass) + (pdq(ig,l,igcm_ccnco2_mass) + &
1038	pdqcloudco2(ig,l, igcm_ccnco2_mass)) * ptimestep)&
1039	, threshold)
1040
1041	myT = pteff(ig,l) + (pdt(ig,l)+pdtcloudco2(ig,l))*ptimestep
1042
1043	Nccnco2_h2o = 0.
1044	Qccnco2_h2o = 0.
1045	if (co2useh2o) then
1046	Nccnco2_h2o = max( (pqeff(ig,l,igcm_ccnco2_h2o_number) + (pdq(ig,l,igcm_ccnco2_h2o_number) + &
1047	pdqcloudco2(ig,l, igcm_ccnco2_h2o_number)) * ptimestep) &
1048	, threshold)
1049
1050	Qccnco2_h2o = max( (pqeff(ig,l,igcm_ccnco2_h2o_mass_ice) + pqeff(ig,l,igcm_ccnco2_h2o_mass_ccn) + &
1051	(pdq(ig,l,igcm_ccnco2_h2o_mass_ice) + pdq(ig,l,igcm_ccnco2_h2o_mass_ccn) + &
1052	pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) + pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn)) * &
1053	ptimestep)&
1054	, threshold)
1055	Nccnco2 = Nccnco2 - Nccnco2_h2o
1056	Qccnco2 = Qccnco2 - Qccnco2_h2o
1057	if (Nccnco2 <= 0) then
1058	Nccnco2 = threshold
1059	Qccnco2 = threshold
1060	end if
1061	end if
1062
1063	! Compute particle size
1064	call updaterice_microco2(dble(Niceco2), dble(Qccnco2), dble(Nccnco2), dble(Qccnco2_h2o), dble(Nccnco2_h2o), myT, &
1065	tauscaling(ig), riceco2(ig,l), rhocloudco2(ig,l))
1066
1067	! Compute opacities
1068	if ( (Niceco2 >= threshold_2 .and. (Nccnco2+Nccnco2_h2o)*tauscaling(ig) >= 1.) ) then
1069	n_derf = derf( (rb_cldco2(1)-log(riceco2(ig,l))) *dev2)
1070	Qext1bins2(ig,l) = 0.
1071	do i = 1, nbinco2_cld
1072	n_aer(i) = -0.5 * (Nccnco2+Nccnco2_h2o)tauscaling(ig) n_derf
1073
1074	n_derf = derf((rb_cldco2(i+1)-log(riceco2(ig,l))) *dev2)
1075	n_aer(i) = n_aer(i) + (0.5 * (Nccnco2+Nccnco2_h2o)tauscaling(ig) n_derf)
1076
1077	Qext1bins2(ig,l) = Qext1bins2(ig,l) + Qext1bins(i) * n_aer(i)
1078	end do
1079	else
1080	riceco2(ig,l) = 0.
1081	Qext1bins2(ig,l) = 0.
1082	end if
1083	!----------------------------------------------------------------------------------------------------------------------!
1084	! 6.2. Update rice and rdust
1085	!----------------------------------------------------------------------------------------------------------------------!
1086	! update rice water only if co2 use h2o ice as CCN
1087	if (co2useh2o) then
1088	call updaterice_micro( &
1089	pqeff(ig,l,igcm_h2o_ice) + (pdq(ig,l,igcm_h2o_ice) + pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep, &
1090	pqeff(ig,l,igcm_ccn_mass) + (pdq(ig,l,igcm_ccn_mass) + pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep, &
1091	pqeff(ig,l,igcm_ccn_number) + (pdq(ig,l,igcm_ccn_number) + pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, &
1092	tauscaling(ig),rice(ig,l),rhocloud(ig,l))
1093	end if
1094
1095	! update rdust
1096	call updaterdust( &
1097	pqeff(ig,l,igcm_dust_mass) + (pdq(ig,l,igcm_dust_mass) + pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, &
1098	pqeff(ig,l,igcm_dust_number) + (pdq(ig,l,igcm_dust_number) + pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, &
1099	rdust(ig,l))
1100	end do ! ngrid
1101	end do ! nlay
1102	!----------------------------------------------------------------------------------------------------------------------!
1103	! 7. A correction if a lot of subliming CO2 fills the 1st layer FF (04/2005). Then that should not affect the ice
1104	! particle radius
1105	!----------------------------------------------------------------------------------------------------------------------!
1106	do ig = 1, ngrid
1107	if ( pdpsrf(ig)ptimestep > 0.9(pplev(ig,1)-pplev(ig,2))) then
1108
1109	if ( pdpsrf(ig)ptimestep > 0.9(pplev(ig,1)-pplev(ig,3)) ) then
1110	riceco2(ig,2) = riceco2(ig,3)
1111	end if
1112
1113	riceco2(ig,1) = riceco2(ig,2)
1114	end if
1115	end do
1116	!----------------------------------------------------------------------------------------------------------------------!
1117	! 8. CO2 saturated quantities
1118	!----------------------------------------------------------------------------------------------------------------------!
1119	! 8.1 Compute CO2 saturation in layers
1120	!----------------------------------------------------------------------------------------------------------------------!
1121	call co2sat(ngridnlay, pteff+(pdt+pdtcloudco2)ptimestep, zqsatco2)
1122	!----------------------------------------------------------------------------------------------------------------------!
1123	! 8.2 Compute CO2 saturated quantities in layers
1124	!----------------------------------------------------------------------------------------------------------------------!
1125	do l = 1, nlay
1126	do ig = 1, ngrid
1127	satuco2(ig,l) = ( pqeff(ig,l,igcm_co2) + (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2))ptimestep ) &
1128	(mmean(ig,l)/(mco21e3)) pplay(ig,l) / zqsatco2(ig,l)
1129	end do
1130	end do
1131	!----------------------------------------------------------------------------------------------------------------------!
1132	! 9. Everything modified by CO2 microphysics must be wrt co2cloudfrac
1133	!----------------------------------------------------------------------------------------------------------------------!
1134	if (CLFvaryingCO2) then
1135	do l = 1, nlay
1136	do ig = 1, ngrid
1137	pdqcloudco2(ig,l,igcm_ccnco2_mass) = pdqcloudco2(ig,l,igcm_ccnco2_mass) * co2cloudfrac(ig,l)
1138
1139	pdqcloudco2(ig,l,igcm_ccnco2_number) = pdqcloudco2(ig,l,igcm_ccnco2_number) * co2cloudfrac(ig,l)
1140
1141	pdqcloudco2(ig,l,igcm_dust_mass) = pdqcloudco2(ig,l,igcm_dust_mass) * co2cloudfrac(ig,l)
1142
1143	pdqcloudco2(ig,l,igcm_dust_number) = pdqcloudco2(ig,l,igcm_dust_number) * co2cloudfrac(ig,l)
1144
1145	pdqcloudco2(ig,l,igcm_co2_ice) = pdqcloudco2(ig,l,igcm_co2_ice) * co2cloudfrac(ig,l)
1146
1147	pdqcloudco2(ig,l,igcm_co2) = pdqcloudco2(ig,l,igcm_co2) * co2cloudfrac(ig,l)
1148
1149	pdtcloudco2(ig,l) = pdtcloudco2(ig,l) * co2cloudfrac(ig,l)
1150
1151	Qext1bins2(ig,l) = Qext1bins2(ig,l) * co2cloudfrac(ig,l)
1152
1153	if (meteo_flux) then
1154	pdqcloudco2(ig,l,igcm_ccnco2_meteor_mass) = pdqcloudco2(ig,l,igcm_ccnco2_meteor_mass) * co2cloudfrac(ig,l)
1155
1156	pdqcloudco2(ig,l,igcm_ccnco2_meteor_number) = pdqcloudco2(ig,l,igcm_ccnco2_meteor_number) * co2cloudfrac(ig,l)
1157	end if
1158	if (co2useh2o) then
1159	pdqcloudco2(ig,l,igcm_h2o_ice) = pdqcloudco2(ig,l,igcm_h2o_ice) * co2cloudfrac(ig,l)
1160
1161	pdqcloudco2(ig,l,igcm_ccn_mass) = pdqcloudco2(ig,l,igcm_ccn_mass) * co2cloudfrac(ig,l)
1162
1163	pdqcloudco2(ig,l,igcm_ccn_number) = pdqcloudco2(ig,l,igcm_ccn_number) * co2cloudfrac(ig,l)
1164
1165	pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) = pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ice) * co2cloudfrac(ig,l)
1166
1167	pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) = pdqcloudco2(ig,l,igcm_ccnco2_h2o_mass_ccn) * co2cloudfrac(ig,l)
1168
1169	pdqcloudco2(ig,l,igcm_ccnco2_h2o_number) = pdqcloudco2(ig,l,igcm_ccnco2_h2o_number) * co2cloudfrac(ig,l)
1170	end if
1171	end do ! ngrid
1172	end do ! nlay
1173	end if ! if CLFvaryingCO2 is true
1174	!----------------------------------------------------------------------------------------------------------------------!
1175	! 10. Compute opacity at 1 micron: Opacity in mesh ig is the sum over l of Qext1bins2. Is this true ?
1176	!----------------------------------------------------------------------------------------------------------------------!
1177	tau1mic(:)=0.
1178	do l = 1, nlay
1179	do ig = 1, ngrid
1180	tau1mic(ig) = tau1mic(ig) + Qext1bins2(ig,l)
1181	end do
1182	end do
1183	!----------------------------------------------------------------------------------------------------------------------!
1184	! 11. Write outputs in diagfi.nc
1185	!----------------------------------------------------------------------------------------------------------------------!
1186	call write_output( "satuco2", "vap in satu", " ", satuco2(:,:))
1187
1188	call write_output( "precip_co2_ice_rate", "surface deposition rate of co2 ice", "kg.m-2.s-1", pdqs_sedco2(:))
1189
1190	call write_output( "co2ice_cond_rate", "CO2 condensation rate in atm layers", "kg.kg-1.s-1", pcondicea(:,:))
1191
1192	call write_output( "pdtcloudco2", "temperature variation of CO2 latent heat", "K.s-1", pdtcloudco2(:,:))
1193
1194	call write_output( "riceco2", "ice radius", "m", riceco2(:,:))
1195
1196	call write_output( "Tau3D1mic", " co2 ice opacities", " ", Qext1bins2(:,:))
1197
1198	call write_output( "tau1mic", "co2 ice opacity 1 micron", " ", tau1mic(:))
1199
1200	if (CLFvaryingCO2) then
1201	call write_output( "co2cloudfrac", "co2 cloud fraction", " ", co2cloudfrac(:,:))
1202	end if
1203	!======================================================================================================================!
1204	! END =================================================================================================================!
1205	!======================================================================================================================!
1206	end subroutine co2cloud
1207
1208	end module co2cloud_mod

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

Download in other formats:

Original Format