Context Navigation

co2cloud.F90 @ 2800

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

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

Download in other formats:

Original Format