Context Navigation

co2cloud.F90 @ 3026

Last change on this file since 3026 was 3008, checked in by emillour, 18 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

Rev	Line
[2362]	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
[2386]	65	! 3.2. No sub-grid cloud representation (CLFvaryingCO2 = False)
[2362]	66	! 4. Microphysics of CO2 cloud formation
[2660]	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
[2362]	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
[2660]	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
[2362]	88	!======================================================================================================================!
	89	subroutine co2cloud(ngrid, nlay, ptimestep, pplev, pplay, pdpsrf, pzlay, pt, pdt, pq, pdq, pdqcloudco2, pdtcloudco2, &
[2562]	90	nq, tau, tauscaling, rdust, rice, riceco2, nuice, rhocloud, rsedcloudco2, rhocloudco2, pzlev,&
[2826]	91	pdqs_sedco2, pdqs_sedccn, pdu, pu, pcondicea)
[2362]	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, &
[2562]	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, &
[2589]	101	nuiceco2_sed, nuiceco2_ref, r3n_q, rho_ice, nuice_sed, igcm_ccnco2_meteor_mass, &
	102	igcm_ccnco2_meteor_number
[2362]	103
	104	use newsedim_mod, only: newsedim
	105
	106	use datafile_mod, only: datadir
[2660]	107	use density_co2_ice_mod, only: density_co2_ice
[2362]	108
[2494]	109	use improvedCO2clouds_mod, only: improvedCO2clouds
[3008]	110	use microphys_h, only: nbinco2_cld, rad_cldco2, mco2
[2932]	111	use write_output_mod, only: write_output
[2447]	112	#ifndef MESOSCALE
	113	use vertical_layers_mod, only: ap, bp
	114	#endif
	115
[2362]	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
[2826]	144	nuice(ngrid,nlay) ! Estimated effective variance of the size distribution
[2362]	145	!----------------------------------------------------------------------------------------------------------------------!
	146	! Output arguments:
	147	!------------------
	148	real, intent(out) :: &
	149	rice(ngrid,nlay), & ! Water Ice mass mean radius (m)
[2447]	150	! rsedcloud(ngrid,nlay), & ! Water Cloud sedimentation radius
[2362]	151	rhocloud(ngrid,nlay), & ! Water Cloud density (kg.m-3)
	152	pdqs_sedco2(ngrid), & ! CO2 flux at the surface
[2562]	153	pdqs_sedccn(ngrid,nq), & ! CCN flux at the surface
[2362]	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
[2447]	157	rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius
	158
	159	real, intent(inout) :: &
	160	rdust(ngrid,nlay) ! Dust geometric mean radius (m)
	161
[2456]	162	double precision, intent(out) :: &
	163	riceco2(ngrid,nlay) ! Ice mass mean radius (m) r_c in Montmessin et al. (2004)
[2362]	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 :: &
[2660]	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)
[2362]	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
[2494]	205	Qextv1mic(var_dim_qext), &
	206	radv(var_dim_qext), & ! radius of CO2 ice at 1 µm (read from file_qext)
[2362]	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
[2616]	210
	211	!$OMP THREADPRIVATE(imicroco2,sigma_iceco2,microtimestep)
	212	!$OMP THREADPRIVATE(dev2,Qext1bins,Qextv1mic,radv,rb_cldco2)
	213	!$OMP THREADPRIVATE(firstcall)
	214
[2362]	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
[2562]	241	sum_subpdqs_sedccn(ngrid,nq), &!
	242	! ---For sub grid T distribution
[2362]	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
[2562]	253	Nccnco2_h2o, &
	254	Qccnco2_h2o, &
[2362]	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
[2447]	262	zqsatco2(ngrid,nlay), &! Saturation co2
	263	availco2,&
	264	masslayer, &
	265	tmp, a,b, &
	266	new_pdq(ngrid,nlay)
	267
[2362]	268	double precision :: &
	269	! ---Extinction coefficients at 1 micron of CO2 particles
	270	vrat_cld, &! Volume ratio
[2456]	271	n_derf, &! derf( (rb_cldco2(1)-log(riceco2(ig,l))) *dev2)
[2362]	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
[2494]	288	tcond(ngrid,nlay) ! CO2 condensation temperature
[2362]	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.
[2494]	421	riceco2(1:ngrid, 1:nlay) = 0.
[2362]	422	zqsed0(1:ngrid,1:nlay,1:nq) = 0.
	423	sum_subpdqs_sedco2(1:ngrid) = 0.
[2562]	424	sum_subpdqs_sedccn(1:ngrid,1:nq) = 0.
[2362]	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
[2460]	431	#ifndef MESOSCALE
[2447]	432	masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1) + (bp(l)-bp(l+1)) ) / g
[2460]	433	#else
	434	masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1)) / g
	435	#endif
[2362]	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
[2932]	496	call write_output( "SatIndex", "SatIndex", " ", SatIndex(:,:))
[2362]	497
[2932]	498	call write_output( "SatIndexmap", "SatIndexmap", "km", SatIndexmap(:))
[2362]	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	!----------------------------------------------------------------------------------------------------------------------!
[2660]	552	! 3.2. No sub-grid cloud representation (CLFvaryingCO2 = False)
[2362]	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	!----------------------------------------------------------------------------------------------------------------------!
[2447]	564	pqeff(:,:,:) = pq(:,:,:)
	565	pteff(:,:) = pt(:,:)
	566	!----------------------------------------------------------------------------------------------------------------------!
[2660]	567	! 4.1. Effective tracers quantities in the cloud
[2447]	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	!----------------------------------------------------------------------------------------------------------------------!
[2362]	577	do microstep = 1, imicroco2
	578	!----------------------------------------------------------------------------------------------------------------------!
[2660]	579	! 4.2. Stepped entry for tendancies: At each micro timestep we add pdt in order to have a stepped entry
[2362]	580	!----------------------------------------------------------------------------------------------------------------------!
[2447]	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
[2589]	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
[2447]	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)
[2562]	612
[2589]	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)
[2562]	615
[2589]	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)
[2562]	618
[2589]	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)
[2447]	621	end if
	622	end do ! ngrid
	623	end do ! nlay
	624	!----------------------------------------------------------------------------------------------------------------------!
[2660]	625	! 4.3. Main call to the cloud scheme
[2447]	626	!----------------------------------------------------------------------------------------------------------------------!
	627	call improvedco2clouds(ngrid, nlay, microtimestep, pplay, pplev, pteff, sum_subpdt, pqeff, sum_subpdq, &
[2562]	628	subpdqcloudco2, subpdtcloudco2, nq, tauscaling, rb_cldco2, sigma_iceco2, dev2)
[2494]	629
[2447]	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)) &
[2562]	633	<= 1.e-12) then
[2447]	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
[2562]	638	if(pq(ig,l,igcm_co2) + microtimestep*(sum_subpdq(ig,l,igcm_co2)+subpdqcloudco2(ig,l,igcm_co2)) <= 1.e-12) then
[2447]	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
[2562]	642
[2447]	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)) &
[2562]	645	*microtimestep)<=1.) .or. &
[2447]	646	(pq(ig,l,igcm_ccnco2_mass)+(sum_subpdq(ig,l,igcm_ccnco2_mass)+subpdqcloudco2(ig,l,igcm_ccnco2_mass)) &
[2562]	647	*microtimestep<=1e-20)) then
[2447]	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
[2562]	656
[2589]	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
[2562]	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
[2447]	694	end do
	695	end do
	696	!----------------------------------------------------------------------------------------------------------------------!
[2660]	697	! 4.4. Updating tendencies after cloud scheme
[2447]	698	!----------------------------------------------------------------------------------------------------------------------!
[2362]	699	do l = 1, nlay
	700	do ig = 1, ngrid
[2447]	701	sum_subpdt(ig,l) = sum_subpdt(ig,l) + subpdtcloudco2(ig,l)
[2362]	702
[2447]	703	sum_subpdq(ig,l,igcm_dust_mass) = sum_subpdq(ig,l,igcm_dust_mass) + subpdqcloudco2(ig,l,igcm_dust_mass)
[2362]	704
[2447]	705	sum_subpdq(ig,l,igcm_dust_number) = sum_subpdq(ig,l,igcm_dust_number) + subpdqcloudco2(ig,l,igcm_dust_number)
[2362]	706
[2447]	707	sum_subpdq(ig,l,igcm_ccnco2_mass) = sum_subpdq(ig,l,igcm_ccnco2_mass) + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
[2362]	708
[2447]	709	sum_subpdq(ig,l,igcm_ccnco2_number) = sum_subpdq(ig,l,igcm_ccnco2_number) + &
	710	subpdqcloudco2(ig,l,igcm_ccnco2_number)
[2362]	711
[2447]	712	sum_subpdq(ig,l,igcm_co2_ice) = sum_subpdq(ig,l,igcm_co2_ice) + subpdqcloudco2(ig,l,igcm_co2_ice)
[2362]	713
[2447]	714	sum_subpdq(ig,l,igcm_co2) = sum_subpdq(ig,l,igcm_co2) + subpdqcloudco2(ig,l,igcm_co2)
[2362]	715
[2589]	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
[2362]	723	if (co2useh2o) then
[2447]	724	sum_subpdq(ig,l,igcm_h2o_ice) = sum_subpdq(ig,l,igcm_h2o_ice) + subpdqcloudco2(ig,l,igcm_h2o_ice)
[2362]	725
[2447]	726	sum_subpdq(ig,l,igcm_ccn_mass) = sum_subpdq(ig,l,igcm_ccn_mass) + subpdqcloudco2(ig,l,igcm_ccn_mass)
[2362]	727
[2447]	728	sum_subpdq(ig,l,igcm_ccn_number) = sum_subpdq(ig,l,igcm_ccn_number) + subpdqcloudco2(ig,l,igcm_ccn_number)
[2562]	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
[2362]	739	end if
	740	end do ! ngrid
	741	end do ! nlay
	742	!----------------------------------------------------------------------------------------------------------------------!
[2660]	743	! 4.5. Gravitational sedimentation
[2362]	744	!----------------------------------------------------------------------------------------------------------------------!
	745	if (sedimentation) then
	746	!----------------------------------------------------------------------------------------------------------------------!
[2660]	747	! 4.5.a. Compute cloud density
[2362]	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
[2589]	758	! meteoritic particle are considered like dust, rho_dust
[2362]	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
[2562]	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
[2660]	770	if (Nccnco2 < 0.) then
	771	Nccnco2 = threshold
	772	Qccnco2 = threshold
	773	end if
[2562]	774	end if
	775
[2362]	776	! Get density cloud and co2 ice particle radius
[2562]	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))
[2447]	780	else
	781	riceco2(ig,l) = 0.
	782	rhocloudco2t(ig,l) = 0.
[2362]	783	end if
	784	end do ! ngrid
	785	end do ! nlay
	786	!----------------------------------------------------------------------------------------------------------------------!
[2660]	787	! 4.5.b. Save actualized tracer values to compute sedimentation tendancies
[2362]	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)
[2562]	792
[2589]	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
[2562]	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
[2362]	803	!----------------------------------------------------------------------------------------------------------------------!
[2660]	804	! 4.5.c. Sedimentation of co2 ice
[2362]	805	!----------------------------------------------------------------------------------------------------------------------!
[2447]	806	do ig = 1, ngrid
	807	do l = 1, nlay
[2660]	808	rsedcloudco2(ig,l) = max( riceco2(ig,l)(1.+nuiceco2_sed)(1.+nuiceco2_sed)*(1.+nuiceco2_sed), &
[2447]	809	rdust(ig,l) )
	810	end do
	811	end do
	812
[2362]	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
[2447]	818	sum_subpdqs_sedco2(ig) = sum_subpdqs_sedco2(ig) + wq(ig,1) / microtimestep !wq est en kg.m-2
[2362]	819	end do
	820	!----------------------------------------------------------------------------------------------------------------------!
[2660]	821	! 4.5.d. Sedimentation for other tracers
[2362]	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)
[2562]	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
[2362]	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)
[2562]	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
[2589]	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
[2562]	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
[2362]	886	!----------------------------------------------------------------------------------------------------------------------!
[2660]	887	! 4.5.e. Compute tendencies due to the sedimation process
[2362]	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
[2562]	898
[2589]	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
[2562]	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
[2362]	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)
[2589]	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
[2562]	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
[2362]	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	!----------------------------------------------------------------------------------------------------------------------!
[2447]	950	! condensation/sublimation rate in the atmosphere (kg.kg-1.s-1)
[2362]	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)
[2562]	960	pdqs_sedccn(ig,:) = sum_subpdqs_sedccn(ig,:) / real(imicroco2)
[2362]	961	end do
[2562]	962
[2447]	963	! temperature tendency (T.s-1)
[2362]	964	do l = 1, nlay
	965	do ig = 1, ngrid
	966	pdtcloudco2(ig,l) = ( sum_subpdt(ig,l)/real(imicroco2) ) - pdt(ig,l)
[2660]	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
[2362]	971	end do
	972	end do
	973
	974	! tracers tendencies
	975	do l = 1, nlay
	976	do ig = 1, ngrid
[2660]	977	pdqcloudco2(ig,l,igcm_co2) = 0. ! here is the trick, this tendency is computed in co2condens
[2362]	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) ) - &
[2447]	988	pdq(ig,l,igcm_dust_mass)
[2362]	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
[2589]	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
[2362]	1000	if (co2useh2o) then
[2447]	1001	pdqcloudco2(ig,l,igcm_h2o_ice) = ( sum_subpdq(ig,l,igcm_h2o_ice) / real(imicroco2) ) - &
	1002	pdq(ig,l,igcm_h2o_ice)
[2362]	1003
	1004	pdqcloudco2(ig,l,igcm_ccn_mass) = ( sum_subpdq(ig,l,igcm_ccn_mass) / real(imicroco2) ) - &
[2447]	1005	pdq(ig,l,igcm_ccn_mass)
[2362]	1006
	1007	pdqcloudco2(ig,l,igcm_ccn_number) = ( sum_subpdq(ig,l,igcm_ccn_number) / real(imicroco2) ) - &
	1008	pdq(ig,l,igcm_ccn_number)
[2562]	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
[2362]	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
[2660]	1030	Niceco2 = max(Niceco2, threshold)
[2362]	1031
[2589]	1032	! meteoritic particles are considered like dust, => rho_dust
[2362]	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
[2562]	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
[2660]	1057	if (Nccnco2 <= 0) then
	1058	Nccnco2 = threshold
	1059	Qccnco2 = threshold
	1060	end if
[2562]	1061	end if
[2660]	1062
[2362]	1063	! Compute particle size
[2562]	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))
[2362]	1066
	1067	! Compute opacities
[2660]	1068	if ( (Niceco2 >= threshold_2 .and. (Nccnco2+Nccnco2_h2o)*tauscaling(ig) >= 1.) ) then
[2456]	1069	n_derf = derf( (rb_cldco2(1)-log(riceco2(ig,l))) *dev2)
[2362]	1070	Qext1bins2(ig,l) = 0.
	1071	do i = 1, nbinco2_cld
[2660]	1072	n_aer(i) = -0.5 * (Nccnco2+Nccnco2_h2o)tauscaling(ig) n_derf
[2362]	1073
[2456]	1074	n_derf = derf((rb_cldco2(i+1)-log(riceco2(ig,l))) *dev2)
[2660]	1075	n_aer(i) = n_aer(i) + (0.5 * (Nccnco2+Nccnco2_h2o)tauscaling(ig) n_derf)
[2362]	1076
	1077	Qext1bins2(ig,l) = Qext1bins2(ig,l) + Qext1bins(i) * n_aer(i)
	1078	end do
[2660]	1079	else
	1080	riceco2(ig,l) = 0.
	1081	Qext1bins2(ig,l) = 0.
[2362]	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	!----------------------------------------------------------------------------------------------------------------------!
[2660]	1117	! 8. CO2 saturated quantities
[2362]	1118	!----------------------------------------------------------------------------------------------------------------------!
[2660]	1119	! 8.1 Compute CO2 saturation in layers
[2362]	1120	!----------------------------------------------------------------------------------------------------------------------!
	1121	call co2sat(ngridnlay, pteff+(pdt+pdtcloudco2)ptimestep, zqsatco2)
	1122	!----------------------------------------------------------------------------------------------------------------------!
[2660]	1123	! 8.2 Compute CO2 saturated quantities in layers
[2362]	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	!----------------------------------------------------------------------------------------------------------------------!
[2660]	1132	! 9. Everything modified by CO2 microphysics must be wrt co2cloudfrac
[2362]	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
[2589]	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
[2362]	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)
[2562]	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)
[2362]	1170	end if
	1171	end do ! ngrid
	1172	end do ! nlay
	1173	end if ! if CLFvaryingCO2 is true
	1174	!----------------------------------------------------------------------------------------------------------------------!
[2660]	1175	! 10. Compute opacity at 1 micron: Opacity in mesh ig is the sum over l of Qext1bins2. Is this true ?
[2362]	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	!----------------------------------------------------------------------------------------------------------------------!
[2660]	1184	! 11. Write outputs in diagfi.nc
[2362]	1185	!----------------------------------------------------------------------------------------------------------------------!
[2932]	1186	call write_output( "satuco2", "vap in satu", " ", satuco2(:,:))
[2362]	1187
[2932]	1188	call write_output( "precip_co2_ice_rate", "surface deposition rate of co2 ice", "kg.m-2.s-1", pdqs_sedco2(:))
[2524]	1189
[2932]	1190	call write_output( "co2ice_cond_rate", "CO2 condensation rate in atm layers", "kg.kg-1.s-1", pcondicea(:,:))
[2524]	1191
[2932]	1192	call write_output( "pdtcloudco2", "temperature variation of CO2 latent heat", "K.s-1", pdtcloudco2(:,:))
[2524]	1193
[2932]	1194	call write_output( "riceco2", "ice radius", "m", riceco2(:,:))
[2362]	1195
[2932]	1196	call write_output( "Tau3D1mic", " co2 ice opacities", " ", Qext1bins2(:,:))
[2362]	1197
[2932]	1198	call write_output( "tau1mic", "co2 ice opacity 1 micron", " ", tau1mic(:))
[2362]	1199
[2524]	1200	if (CLFvaryingCO2) then
[2932]	1201	call write_output( "co2cloudfrac", "co2 cloud fraction", " ", co2cloudfrac(:,:))
[2524]	1202	end if
[2362]	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