Context Navigation

co2cloud.F @ 1918

Last change on this file since 1918 was 1918, checked in by emillour, 7 years ago

Mars GCM:
Code cleanup:

remove "comorbit.h" since it is no longer used.
turn "datafile.h" into module datafile_mod.F90 (and rename variable "datafile" as "datadir" since it stores the path to the datafile directory).

Property svn:executable set to *

File size: 39.0 KB

Rev	Line
[1685]	1	SUBROUTINE co2cloud(ngrid,nlay,ptimestep,
[1617]	2	& pplev,pplay,pdpsrf,pzlay,pt,pdt,
	3	& pq,pdq,pdqcloudco2,pdtcloudco2,
	4	& nq,tau,tauscaling,rdust,rice,riceco2,nuice,
[1685]	5	& rsedcloudco2,rhocloudco2,
[1816]	6	& rsedcloud,rhocloud,pzlev,pdqs_sedco2,
	7	& pdu,pu)
[1820]	8	USE ioipsl_getincom, only: getin
[1617]	9	use dimradmars_mod, only: naerkind
[1820]	10	USE comcstfi_h, only: pi, g, cpp
	11	USE updaterad, only: updaterice_microco2, updaterice_micro,
	12	& updaterdust
[1816]	13	use conc_mod, only: mmean,rnew
[1617]	14	use tracer_mod, only: nqmx, igcm_co2, igcm_co2_ice,
[1685]	15	& igcm_dust_mass, igcm_dust_number,igcm_h2o_ice,
	16	& igcm_ccn_mass,igcm_ccn_number,
[1617]	17	& igcm_ccnco2_mass, igcm_ccnco2_number,
	18	& rho_dust, nuiceco2_sed, nuiceco2_ref,
[1720]	19	& rho_ice_co2,r3n_q,rho_ice,nuice_sed
[1913]	20	USE newsedim_mod, ONLY: newsedim
[1918]	21	USE datafile_mod, ONLY: datadir
[1617]	22	IMPLICIT NONE
	23
[1820]	24	include "callkeys.h"
	25	include "microphys.h"
[1720]	26
[1617]	27	c=======================================================================
	28	c CO2 clouds formation
	29	c
	30	c There is a time loop specific to cloud formation
	31	c due to timescales smaller than the GCM integration timestep.
	32	c microphysics subroutine is improvedCO2clouds.F
[1816]	33	c the microphysics time step is a fraction of the physical one
	34	c the integer imicroco2 must be set in callphys.def
	35	c
[1617]	36	c The co2 clouds tracers (co2_ice, ccn mass and concentration) are
	37	c sedimented at each microtimestep. pdqs_sedco2 keeps track of the
	38	c CO2 flux at the surface
	39	c
	40	c Authors: 09/2016 Joachim Audouard & Constantino Listowski
	41	c Adaptation of the water ice clouds scheme (with specific microphysics)
	42	c of Montmessin, Navarro & al.
	43	c
[1720]	44	c 07/2017 J.Audouard
[1816]	45	c Several logicals and integer must be set to .true. in callphys.def
[1818]	46	c if not, default values are .false.
[1720]	47	c co2clouds=.true. call this routine
	48	c co2useh2o=.true. allow the use of water ice particles as CCN for CO2
	49	c meteo_flux=.true. supply meteoritic particles
	50	c CLFvaryingCO2=.true. allows a subgrid temperature distribution
[1818]	51	c of amplitude spantCO2(=integer in callphys.def, typically 3)
	52	c satindexco2=.true. allows the filtering out of the sub-grid T distribution
	53	c if the GW saturates in the column. Based on Spiga et al
	54	c 2012
	55	c An index is computed for the column, and the sub-grid T
	56	c distribution is applied if the index remains < 0.1
	57	c setting to .false. applies the sub-grid T everywhere.
	58	c default value is .true., only applies if
	59	c CLFvaryingCO2=.true. anyway.
[1816]	60	c imicroco2=50
[1720]	61	c
[1816]	62	c The subgrid Temperature distribution is modulated (0 or 1) by Spiga et
	63	c al. (GRL 2012) Saturation Index to account for GW propagation or
	64	c dissipation upwards.
	65	c
	66	c 4D and column opacities are computed using Qext values at 1µm.
[1617]	67	c=======================================================================
	68
	69	c-----------------------------------------------------------------------
[1820]	70	c arguments:
[1617]	71	c -------------
	72
[1820]	73	INTEGER, INTENT(IN) :: ngrid,nlay
	74	REAL, INTENT(IN) :: ptimestep ! pas de temps physique (s)
	75	REAL, INTENT(IN) :: pplev(ngrid,nlay+1) ! Inter-layer pressures (Pa)
	76	REAL, INTENT(IN) :: pplay(ngrid,nlay) ! mid-layer pressures (Pa)
	77	REAL, INTENT(IN) :: pdpsrf(ngrid) ! tendency on surface pressure
	78	REAL, INTENT(IN) :: pzlay(ngrid,nlay) ! altitude at the middle of the layers
	79	REAL, INTENT(IN) :: pt(ngrid,nlay) ! temperature at the middle of the layers (K)
	80	REAL, INTENT(IN) :: pdt(ngrid,nlay) ! tendency on temperature from other parametrizations
	81	real, INTENT(IN) :: pq(ngrid,nlay,nq) ! tracers (kg/kg)
	82	real, INTENT(IN) :: pdq(ngrid,nlay,nq) ! tendencies before condensation (kg/kg.s-1)
[1911]	83	real, intent(OUT) :: pdqcloudco2(ngrid,nlay,nq) ! tendency due to CO2 condensation (kg/kg.s-1)
	84	real, intent(OUT) :: pdtcloudco2(ngrid,nlay) ! tendency on temperature due to latent heat
[1820]	85	INTEGER, INTENT(IN) :: nq ! number of tracers
	86	REAL, INTENT(IN) :: tau(ngrid,naerkind) ! Column dust optical depth at each point
	87	REAL, INTENT(IN) :: tauscaling(ngrid) ! Convertion factor for dust amount
	88	REAL, INTENT(OUT) :: rdust(ngrid,nlay) ! Dust geometric mean radius (m)
	89	real, intent(OUT) :: rice(ngrid,nlay) ! Water Ice mass mean radius (m)
[1617]	90	! used for nucleation of CO2 on ice-coated ccns
[1820]	91	DOUBLE PRECISION, INTENT(out) :: riceco2(ngrid,nlay) ! Ice mass mean radius (m)
[1617]	92	! (r_c in montmessin_2004)
[1911]	93	REAL, INTENT(IN) :: nuice(ngrid,nlay) ! Estimated effective variance
[1617]	94	! of the size distribution
[1911]	95	real, intent(OUT) :: rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius
	96	real, intent(OUT) :: rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
	97	real, intent(OUT) :: rsedcloud(ngrid,nlay) ! Water Cloud sedimentation radius
	98	real, intent(OUT) :: rhocloud(ngrid,nlay) ! Water Cloud density (kg.m-3)
	99	real, intent(IN) :: pzlev(ngrid,nlay+1) ! altitude at the boundaries of the layers
	100	real, intent(OUT) :: pdqs_sedco2(ngrid) ! CO2 flux at the surface
[1820]	101	REAL, INTENT(IN) :: pdu(ngrid,nlay),pu(ngrid,nlay) !Zonal Wind: zu=pu+pdu*ptimestep
[1816]	102
[1617]	103	c local:
	104	c ------
[1816]	105
[1617]	106	! for time loop
	107	INTEGER microstep ! time subsampling step variable
[1820]	108	INTEGER, SAVE :: imicroco2 ! time subsampling for coupled water microphysics & sedimentation
	109	REAL, SAVE :: microtimestep ! integration timestep for coupled water microphysics & sedimentation
[1617]	110
	111	! tendency given by clouds (inside the micro loop)
	112	REAL subpdqcloudco2(ngrid,nlay,nq) ! cf. pdqcloud
	113	REAL subpdtcloudco2(ngrid,nlay) ! cf. pdtcloud
	114
	115	! global tendency (clouds+physics)
[1911]	116	REAL sum_subpdq(ngrid,nlay,nq) ! cf. pdqcloud
	117	REAL sum_subpdt(ngrid,nlay) ! cf. pdtcloud
[1617]	118	real wq(ngrid,nlay+1) ! ! displaced tracer mass (kg.m-2) during microtimestep because sedim (?/m-2)
	119
	120	REAL satuco2(ngrid,nlay) ! co2 satu ratio for output
	121	REAL zqsatco2(ngrid,nlay) ! saturation co2
	122
[1820]	123	DOUBLE PRECISION rho_ice_co2T(ngrid,nlay) !T-dependant CO2 ice density
	124	DOUBLE PRECISION :: myT ! temperature scalar for co2 density computation
	125
[1720]	126	INTEGER iq,ig,l,i
[1617]	127	LOGICAL,SAVE :: firstcall=.true.
[1816]	128	DOUBLE PRECISION Nccnco2, Niceco2,Nco2,Qccnco2
	129	real :: beta ! for sedimentation
[1617]	130
	131	real epaisseur (ngrid,nlay) ! Layer thickness (m)
	132	real masse (ngrid,nlay) ! Layer mass (kg.m-2)
[1911]	133	real ztsed(ngrid,nlay) ! tracers with real-time value in microtimeloop
	134	real zqsed(ngrid,nlay,nq)
	135	real zqsed0(ngrid,nlay,nq) !For sedimentation tendancy
	136	real subpdqsed(ngrid,nlay,nq)
	137	real sum_subpdqs_sedco2(ngrid) ! CO2 flux at the surface
[1685]	138
[1816]	139	DOUBLE PRECISION,allocatable,save :: memdMMccn(:,:) !memory of h2o particles
	140	DOUBLE PRECISION,allocatable,save :: memdMMh2o(:,:) !only if co2useh2o=.true.
	141	DOUBLE PRECISION,allocatable,save :: memdNNccn(:,:) !Nb particules H2O intégré
	142
	143	! What we need for Qext reading and tau computation : size distribution
[1720]	144	DOUBLE PRECISION vrat_cld ! Volume ratio
[1820]	145	DOUBLE PRECISION, SAVE :: rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
[1816]	146	DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m)
	147	DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6 ! Maximum radius (m)
	148	DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-10! Minimum boundary radius (m)
	149	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4 ! Maximum boundary radius (m)
[1720]	150	DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m)
	151	DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3)
[1820]	152	REAL, SAVE :: sigma_iceco2 ! Variance of the ice and CCN distributions
[1816]	153	logical :: file_ok !Qext file reading
[1720]	154	double precision :: radv(10000),Qextv1mic(10000)
[1820]	155	double precision, save :: Qext1bins(100)
	156	double precision :: Qtemp
[1720]	157	double precision :: ltemp1(10000),ltemp2(10000)
[1913]	158	integer :: nelem,lebon1,lebon2
	159	integer,parameter :: uQext=555
[1720]	160	DOUBLE PRECISION n_aer(nbinco2_cld),Rn,No,n_derf,dev2
	161	DOUBLE PRECISION Qext1bins2(ngrid,nlay)
	162	DOUBLE PRECISION tau1mic(ngrid) !co2 ice column opacity at 1µm
	163
	164	! For sub grid T distribution
	165
	166	REAL zt(ngrid,nlay) ! local value of temperature
	167	REAL :: zq(ngrid, nlay,nq)
	168
[1820]	169	real :: rhocloudco2t(ngrid,nlay) ! Cloud density (kg.m-3)
	170
[1816]	171	DOUBLE PRECISION :: tcond(ngrid,nlay) !CO2 condensation temperature
[1720]	172	REAL :: zqvap(ngrid,nlay)
[1816]	173	REAL :: zqice(ngrid,nlay)
[1720]	174	REAL :: spant,zdelt ! delta T for the temperature distribution
[1911]	175	REAL :: pteff(ngrid, nlay)! effective temperature in the cloud,neb
[1816]	176	REAL :: pqeff(ngrid, nlay, nq)! effective tracers quantities in the cloud
[1885]	177	REAL :: co2cloudfrac(ngrid,nlay) ! cloud fraction
[1816]	178	REAL :: mincloud ! min cloud frac
	179	DOUBLE PRECISION:: rho,zu,NN,gradT !For Saturation Index computation
	180	DOUBLE PRECISION :: SatIndex(ngrid,nlay),SatIndexmap(ngrid)
	181
[1720]	182	c logical :: CLFvaryingCO2
[1617]	183	c ** un petit test de coherence
	184	c --------------------------
	185
	186	IF (firstcall) THEN
	187	if (nq.gt.nqmx) then
	188	write(,) 'stop in co2cloud (nq.gt.nqmx)!'
	189	write(,) 'nq=',nq,' nqmx=',nqmx
	190	stop
	191	endif
[1685]	192	write(,) "co2cloud.F: rho_ice_co2 = ",rho_ice_co2
[1617]	193	write(,) "co2cloud: igcm_co2=",igcm_co2
	194	write(,) " igcm_co2_ice=",igcm_co2_ice
	195
	196	write(,) "time subsampling for microphysic ?"
	197	#ifdef MESOSCALE
[1816]	198	imicroco2 = 2
[1617]	199	#else
[1816]	200	imicroco2 = 30
[1617]	201	#endif
[1816]	202	call getin("imicroco2",imicroco2)
	203	write(,)"imicroco2 = ",imicroco2
[1617]	204
[1816]	205	microtimestep = ptimestep/real(imicroco2)
[1617]	206	write(,)"Physical timestep is",ptimestep
	207	write(,)"CO2 Microphysics timestep is",microtimestep
[1651]	208
[1816]	209
[1685]	210	if (.not. allocated(memdMMccn)) allocate(memdMMccn(ngrid,nlay))
	211	if (.not. allocated(memdNNccn)) allocate(memdNNccn(ngrid,nlay))
	212	if (.not. allocated(memdMMh2o)) allocate(memdMMh2o(ngrid,nlay))
	213
	214	memdMMccn(:,:)=0.
	215	memdMMh2o(:,:)=0.
	216	memdNNccn(:,:)=0.
[1720]	217
	218	c Compute the size bins of the distribution of CO2 ice particles
	219	c --> used for opacity calculations
	220
	221	c rad_cldco2 is the primary radius grid used for microphysics computation.
	222	c The grid spacing is computed assuming a constant volume ratio
	223	c between two consecutive bins; i.e. vrat_cld.
	224	c vrat_cld is determined from the boundary values of the size grid:
	225	c rmin_cld and rmax_cld.
	226	c The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
	227	c dr_cld is the width of each rad_cldco2 bin.
	228	sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
	229	c Volume ratio between two adjacent bins
	230	! vrat_cld
	231	vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
	232	vrat_cld = exp(vrat_cld)
	233	rb_cldco2(1) = rbmin_cld
	234	rad_cldco2(1) = rmin_cld
	235	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
	236	do i=1,nbinco2_cld-1
	237	rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.)
	238	vol_cld(i+1) = vol_cld(i) * vrat_cld
	239	enddo
	240	do i=1,nbinco2_cld
	241	rb_cldco2(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
	242	& rad_cldco2(i)
	243	dr_cld(i) = rb_cldco2(i+1) - rb_cldco2(i)
	244	enddo
	245	rb_cldco2(nbinco2_cld+1) = rbmax_cld
	246	dr_cld(nbinco2_cld) = rb_cldco2(nbinco2_cld+1) -
	247	& rb_cldco2(nbinco2_cld)
	248
	249	c read the Qext values
[1918]	250	INQUIRE(FILE=TRIM(datadir)//
[1720]	251	& '/optprop_co2ice_1mic.dat', EXIST=file_ok)
	252	IF (.not. file_ok) THEN
	253	write(,) 'file optprop_co2ice_1mic.dat should be in '
[1918]	254	& ,trim(datadir)
[1720]	255	STOP
	256	endif
[1918]	257	! open(newunit=uQext,file=trim(datadir)//
	258	open(unit=uQext,file=trim(datadir)//
[1720]	259	& '/optprop_co2ice_1mic.dat'
	260	& ,FORM='formatted')
	261	read(uQext,*) !skip 1 line
	262	do i=1,10000
	263	read(uQext,'(E11.5)') radv(i)
	264	enddo
	265	read(uQext,*) !skip 1 line
	266	do i=1,10000
	267	read(uQext,'(E11.5)') Qextv1mic(i)
	268	enddo
	269	close(uQext)
	270	c innterpol the Qext values
	271	!rice_out=rad_cldco2
	272	do i=1,nbinco2_cld
	273	ltemp1=abs(radv(:)-rb_cldco2(i))
	274	ltemp2=abs(radv(:)-rb_cldco2(i+1))
	275	lebon1=minloc(ltemp1,DIM=1)
[1816]	276	lebon2=min(minloc(ltemp2,DIM=1),10000)
[1720]	277	nelem=lebon2-lebon1+1.
	278	Qtemp=0d0
	279	do l=0,nelem
[1816]	280	Qtemp=Qtemp+Qextv1mic(min(lebon1+l,10000)) !mean value in the interval
[1720]	281	enddo
	282	Qtemp=Qtemp/nelem
	283	Qext1bins(i)=Qtemp
	284	enddo
	285	Qext1bins(:)=Qext1bins(:)rad_cldco2(:)rad_cldco2(:)*pi
	286	! The actuall tau computation and output is performed in co2cloud.F
	287
	288	print*,'--------------------------------------------'
	289	print*,'Microphysics co2: size bin-Qext information:'
	290	print*,' i, rad_cldco2(i), Qext1bins(i)'
	291	do i=1,nbinco2_cld
[1885]	292	write(*,'(i3,3x,3(e13.6,4x))') i, rad_cldco2(i),
[1720]	293	& Qext1bins(i)
	294	enddo
	295	print*,'--------------------------------------------'
	296
	297
	298	do i=1,nbinco2_cld+1
	299	rb_cldco2(i) = log(rb_cldco2(i)) !! we save that so that it is not computed at each timestep and gridpoint
	300	enddo
	301	if (CLFvaryingCO2) then
	302	write(,)
	303	write(,) "CLFvaryingCO2 is set to true is callphys.def"
	304	write(,) "The temperature field is enlarged to +/-",spantCO2
	305	write(,) "for the CO2 microphysics "
	306	endif
[1816]	307
[1617]	308	firstcall=.false.
	309	ENDIF ! of IF (firstcall)
	310
	311	c-----Initialization
[1720]	312	dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
[1685]	313	beta=0.85
[1911]	314	sum_subpdq(1:ngrid,1:nlay,1:nq) = 0
	315	sum_subpdt(1:ngrid,1:nlay) = 0
[1685]	316	subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0
	317	subpdtcloudco2(1:ngrid,1:nlay) = 0
	318
[1617]	319	wq(:,:)=0
	320	! default value if no ice
	321	rhocloudco2(1:ngrid,1:nlay) = rho_dust
	322	rhocloudco2t(1:ngrid,1:nlay) = rho_dust
	323	epaisseur(1:ngrid,1:nlay)=0
	324	masse(1:ngrid,1:nlay)=0
	325
[1911]	326	zqsed0(1:ngrid,1:nlay,1:nq)=0
	327	sum_subpdqs_sedco2(1:ngrid)=0
	328	subpdqsed(1:ngrid,1:nlay,1:nq)=0
[1617]	329
	330	do l=1,nlay
	331	do ig=1, ngrid
	332	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
[1720]	333	epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
[1617]	334	enddo
	335	enddo
	336
	337	c-------------------------------------------------------------------
[1720]	338	c 0. Representation of sub-grid water ice clouds
[1816]	339	c-------------------------------------------------------------------
[1720]	340	IF (CLFvaryingCO2) THEN
[1816]	341
[1720]	342	spant=spantCO2 ! delta T for the temprature distribution
[1885]	343	mincloud=0.1 ! min co2cloudfrac when there is ice
[1911]	344	pteff(:,:)=pt(:,:)
[1885]	345	co2cloudfrac(:,:)=mincloud
[1720]	346
	347	c-----Tendencies
	348	DO l=1,nlay
	349	DO ig=1,ngrid
	350	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	351	ENDDO
	352	ENDDO
	353	DO l=1,nlay
	354	DO ig=1,ngrid
	355	DO iq=1,nq
	356	zq(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
	357	ENDDO
	358	ENDDO
	359	ENDDO
	360	zqvap=zq(:,:,igcm_co2)
	361	zqice=zq(:,:,igcm_co2_ice)
[1816]	362
	363
[1885]	364	call WRITEDIAGFI(ngrid,"co2cloud_pzlev","pzlev","km",3,
[1816]	365	& pzlev)
[1885]	366	call WRITEDIAGFI(ngrid,"co2cloud_pzlay","pzlay","km",3,
[1816]	367	& pzlay)
[1885]	368	call WRITEDIAGFI(ngrid,"co2cloud_pplay","pplay","Pa",3,
[1816]	369	& pplay)
	370
[1818]	371	if (satindexco2) then !logical in callphys.def
[1820]	372	DO l=12,26
	373	! layers 12 --> 26 ~ 12->85 km
	374	DO ig=1,ngrid
	375	! compute N^2 static stability
	376	gradT=(zt(ig,l+1)-zt(ig,l))/(pzlev(ig,l+1)-pzlev(ig,l))
	377	NN=sqrt(g/zt(iq,l)*(max(gradT,-g/cpp)+g/cpp))
	378	! compute absolute value of zonal wind field
	379	zu=abs(pu(ig,l) + pdu(ig,l)*ptimestep)
	380	! compute background density
	381	rho=pplay(ig,l)/(rnew(ig,l)*zt(ig,l))
	382	!saturation index
	383	SatIndex(ig,l)=sqrt(7.5e-7150.e3/(2.pi)*NN/
	384	& (rhozuzu*zu))
	385	ENDDO
	386	ENDDO
[1816]	387	!Then compute Satindex map
	388	! layers 12 --> 26 ~ 12->85 km
[1820]	389	DO ig=1,ngrid
	390	SatIndexmap(ig)=maxval(SatIndex(ig,12:26))
	391	ENDDO
[1816]	392
[1820]	393	call WRITEDIAGFI(ngrid,"SatIndexmap","SatIndexmap","km",2,
	394	& SatIndexmap)
[1818]	395	else
[1820]	396	do ig=1,ngrid
	397	SatIndexmap(ig)=0.05 !maxval(SatIndex(ig,12:26))
	398	enddo
	399	endif ! of if (satindexco2)
[1816]	400
	401	!Modulate the DeltaT by GW propagation index :
	402	! Saturation index S in Spiga 2012 paper
	403	!Assuming like in the paper,
	404	!GW phase speed (stationary waves) c=0 m.s-1
	405	!lambdaH =150 km
	406	!Fo=7.5e-7 J.m-3
	407
	408	CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
[1720]	409	zdelt=spant
[1820]	410	DO ig=1,ngrid
[1816]	411
[1820]	412	IF (SatIndexmap(ig) .le. 0.1) THEN
	413	DO l=1,nlay-1
[1816]	414
	415	IF (tcond(ig,l) .ge. (zt(ig,l)+zdelt)
[1820]	416	& .or. tcond(ig,l) .le. 0 ) THEN !The entire fraction is saturated
[1911]	417	pteff(ig,l)=zt(ig,l)
[1885]	418	co2cloudfrac(ig,l)=1.
[1820]	419	ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN ! No saturation at all
[1911]	420	pteff(ig,l)=zt(ig,l)-zdelt
[1885]	421	co2cloudfrac(ig,l)=mincloud
[1720]	422	ELSE
[1885]	423	co2cloudfrac(ig,l)=(tcond(ig,l)-zt(ig,l)+zdelt)/
[1720]	424	& (2.0*zdelt)
[1911]	425	pteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Mean temperature of the cloud fraction
[1816]	426	END IF !ig if (tcond(ig,l) ...
[1911]	427	pteff(ig,l)=pteff(ig,l)-pdt(ig,l)*ptimestep
[1885]	428	IF (co2cloudfrac(ig,l).le. mincloud) THEN
	429	co2cloudfrac(ig,l)=mincloud
	430	ELSE IF (co2cloudfrac(ig,l).gt. 1) THEN
	431	co2cloudfrac(ig,l)=1.
[1720]	432	END IF
[1820]	433	ENDDO
	434	ELSE
[1816]	435	!SatIndex not favorable for GW : leave pt untouched
[1911]	436	pteff(ig,l)=pt(ig,l)
[1885]	437	co2cloudfrac(ig,l)=mincloud
[1820]	438	ENDIF ! of if(SatIndexmap...
	439	ENDDO ! of DO ig=1,ngrid
[1720]	440	! Totalcloud frac of the column missing here
	441	c-----------------------
	442	c-----No sub-grid cloud representation (CLFvarying=false)
	443	ELSE
	444	DO l=1,nlay
	445	DO ig=1,ngrid
[1911]	446	pteff(ig,l)=pt(ig,l)
[1720]	447	END DO
	448	END DO
[1816]	449	END IF ! end if (CLFvaryingco2)
[1617]	450	c------------------------------------------------------------------
[1816]	451	c microtimestep timeloop for microphysics:
	452	c 0.Stepped entry for tendancies
	453	c 1.Compute sedimentation and update tendancies
	454	c 2.Call co2clouds microphysics
	455	c 3.Update tendancies
[1617]	456	c------------------------------------------------------------------
[1816]	457	DO microstep=1,imicroco2
[1617]	458	c------ Temperature tendency subpdt
	459	! If imicro=1 subpdt is the same as pdt
[1820]	460	DO l=1,nlay
	461	DO ig=1,ngrid
[1911]	462	sum_subpdt(ig,l) = sum_subpdt(ig,l)
[1816]	463	& + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
[1911]	464	sum_subpdq(ig,l,igcm_dust_mass) =
	465	& sum_subpdq(ig,l,igcm_dust_mass)
[1816]	466	& + pdq(ig,l,igcm_dust_mass)
[1911]	467	sum_subpdq(ig,l,igcm_dust_number) =
	468	& sum_subpdq(ig,l,igcm_dust_number)
[1816]	469	& + pdq(ig,l,igcm_dust_number)
	470
[1911]	471	sum_subpdq(ig,l,igcm_ccnco2_mass) =
	472	& sum_subpdq(ig,l,igcm_ccnco2_mass)
[1816]	473	& + pdq(ig,l,igcm_ccnco2_mass)
[1911]	474	sum_subpdq(ig,l,igcm_ccnco2_number) =
	475	& sum_subpdq(ig,l,igcm_ccnco2_number)
[1816]	476	& + pdq(ig,l,igcm_ccnco2_number)
	477
[1911]	478	sum_subpdq(ig,l,igcm_co2_ice) =
	479	& sum_subpdq(ig,l,igcm_co2_ice)
[1816]	480	& + pdq(ig,l,igcm_co2_ice)
[1911]	481	sum_subpdq(ig,l,igcm_co2) =
	482	& sum_subpdq(ig,l,igcm_co2)
[1816]	483	& + pdq(ig,l,igcm_co2)
	484
[1911]	485	sum_subpdq(ig,l,igcm_h2o_ice) =
	486	& sum_subpdq(ig,l,igcm_h2o_ice)
[1816]	487	& + pdq(ig,l,igcm_h2o_ice)
[1911]	488	sum_subpdq(ig,l,igcm_ccn_mass) =
	489	& sum_subpdq(ig,l,igcm_ccn_mass)
[1816]	490	& + pdq(ig,l,igcm_ccn_mass)
[1911]	491	sum_subpdq(ig,l,igcm_ccn_number) =
	492	& sum_subpdq(ig,l,igcm_ccn_number)
[1816]	493	& + pdq(ig,l,igcm_ccn_number)
[1820]	494	ENDDO
	495	ENDDO
[1816]	496	c- Effective tracers quantities in the cloud fraction
[1820]	497	IF (CLFvaryingCO2) THEN
[1816]	498	pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier)
	499	pqeff(:,:,igcm_ccnco2_mass) =pq(:,:,igcm_ccnco2_mass)/
[1885]	500	& co2cloudfrac(:,:)
[1816]	501	pqeff(:,:,igcm_ccnco2_number)=
[1885]	502	& pq(:,:,igcm_ccnco2_number)/co2cloudfrac(:,:)
[1816]	503	pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)/
[1885]	504	& co2cloudfrac(:,:)
[1820]	505	ELSE
[1816]	506	pqeff(:,:,:)=pq(:,:,:)
[1820]	507	END IF
[1617]	508
[1816]	509	c------------------------------------------------------
	510	c 1.SEDIMENTATION : update tracers, compute parameters,
	511	c call to sedimentation routine, update tendancies
	512	c------------------------------------------------------
[1911]	513	IF (sedimentation) THEN
	514
[1820]	515	DO l=1, nlay
[1720]	516	DO ig=1,ngrid
[1911]	517	ztsed(ig,l)=pteff(ig,l)
	518	& +sum_subpdt(ig,l)*microtimestep
	519	zqsed(ig,l,:)=pqeff(ig,l,:)
	520	& +sum_subpdq(ig,l,:)*microtimestep
[1685]	521	rho_ice_co2T(ig,l)=1000.(1.72391-2.53e-4
[1911]	522	& ztsed(ig,l)-2.87e-6*
	523	& ztsed(ig,l)*ztsed(ig,l))
[1720]	524
	525	rho_ice_co2=rho_ice_co2T(ig,l)
[1911]	526	Niceco2=max(zqsed(ig,l,igcm_co2_ice),1.e-30)
	527	Nccnco2=max(zqsed(ig,l,igcm_ccnco2_number),
[1617]	528	& 1.e-30)
[1911]	529	Qccnco2=max(zqsed(ig,l,igcm_ccnco2_mass),
[1617]	530	& 1.e-30)
[1816]	531	call updaterice_microco2(Niceco2,
	532	& Qccnco2,Nccnco2,
	533	& tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
	534	if (Niceco2 .le. 1.e-25
	535	& .or. Nccnco2*tauscaling(ig) .le. 1) THEN
	536	riceco2(ig,l)=1.e-9
	537	endif
[1720]	538	rhocloudco2t(ig,l)=min(max(rhocloudco2t(ig,l)
	539	& ,rho_ice_co2),rho_dust)
	540	rsedcloudco2(ig,l)=max(riceco2(ig,l)*
[1617]	541	& (1.+nuiceco2_sed)(1.+nuiceco2_sed)(1.+nuiceco2_sed),
[1816]	542	& riceco2(ig,l))
[1617]	543	ENDDO
[1820]	544	ENDDO
[1816]	545	! Gravitational sedimentation
[1911]	546	zqsed0(:,:,igcm_co2_ice)=zqsed(:,:,igcm_co2_ice)
	547	zqsed0(:,:,igcm_ccnco2_mass)=zqsed(:,:,igcm_ccnco2_mass)
	548	zqsed0(:,:,igcm_ccnco2_number)=zqsed(:,:,igcm_ccnco2_number)
[1816]	549	!We save actualized tracer values to compute sedimentation tendancies
[1820]	550	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
[1911]	551	& microtimestep,pplev,masse,epaisseur,ztsed,
[1617]	552	& rsedcloudco2,rhocloudco2t,
[1911]	553	& zqsed(:,:,igcm_co2_ice),wq,beta) ! 3 traceurs
[1720]	554	! sedim at the surface of co2 ice : keep track of it for physiq_mod
[1820]	555	do ig=1,ngrid
[1911]	556	sum_subpdqs_sedco2(ig)=
	557	& sum_subpdqs_sedco2(ig)+ wq(ig,1)/microtimestep
[1820]	558	end do
	559	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
[1911]	560	& microtimestep,pplev,masse,epaisseur,ztsed,
[1617]	561	& rsedcloudco2,rhocloudco2t,
[1911]	562	& zqsed(:,:,igcm_ccnco2_mass),wq,beta)
[1820]	563	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
[1911]	564	& microtimestep,pplev,masse,epaisseur,ztsed,
[1617]	565	& rsedcloudco2,rhocloudco2t,
[1911]	566	& zqsed(:,:,igcm_ccnco2_number),wq,beta)
[1820]	567	DO l = 1, nlay !Compute tendencies
	568	DO ig=1,ngrid
[1911]	569	subpdqsed(ig,l,igcm_ccnco2_mass)=
	570	& (zqsed(ig,l,igcm_ccnco2_mass)-
	571	& zqsed0(ig,l,igcm_ccnco2_mass))/microtimestep
	572	subpdqsed(ig,l,igcm_ccnco2_number)=
	573	& (zqsed(ig,l,igcm_ccnco2_number)-
	574	& zqsed0(ig,l,igcm_ccnco2_number))/microtimestep
	575	subpdqsed(ig,l,igcm_co2_ice)=
	576	& (zqsed(ig,l,igcm_co2_ice)-
	577	& zqsed0(ig,l,igcm_co2_ice))/microtimestep
[1820]	578	ENDDO
	579	ENDDO
[1816]	580	!update subtimestep tendencies with sedimentation input
	581	DO l=1,nlay
[1617]	582	DO ig=1,ngrid
[1911]	583	sum_subpdq(ig,l,igcm_ccnco2_mass) =
	584	& sum_subpdq(ig,l,igcm_ccnco2_mass)
	585	& +subpdqsed(ig,l,igcm_ccnco2_mass)
	586	sum_subpdq(ig,l,igcm_ccnco2_number) =
	587	& sum_subpdq(ig,l,igcm_ccnco2_number)
	588	& +subpdqsed(ig,l,igcm_ccnco2_number)
	589	sum_subpdq(ig,l,igcm_co2_ice) =
	590	& sum_subpdq(ig,l,igcm_co2_ice)
	591	& +subpdqsed(ig,l,igcm_co2_ice)
[1617]	592	ENDDO
[1911]	593	ENDDO
	594
	595	END IF !(end if sedimentation)
	596
[1816]	597	c------------------------------------------------------
	598	c 2. Main call to the cloud schemes:
	599	c------------------------------------------------------
[1820]	600	CALL improvedCO2clouds(ngrid,nlay,microtimestep,
[1911]	601	& pplay,pplev,pteff,sum_subpdt,
	602	& pqeff,sum_subpdq,subpdqcloudco2,subpdtcloudco2,
[1816]	603	& nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn)
	604	c-----------------------------------------------------
	605	c 3. Updating tendencies after cloud scheme:
	606	c-----------------------------------------------------
[1820]	607	DO l=1,nlay
	608	DO ig=1,ngrid
[1911]	609	sum_subpdt(ig,l) =
	610	& sum_subpdt(ig,l) + subpdtcloudco2(ig,l)
[1816]	611
[1911]	612	sum_subpdq(ig,l,igcm_dust_mass) =
	613	& sum_subpdq(ig,l,igcm_dust_mass)
[1816]	614	& + subpdqcloudco2(ig,l,igcm_dust_mass)
[1911]	615	sum_subpdq(ig,l,igcm_dust_number) =
	616	& sum_subpdq(ig,l,igcm_dust_number)
[1816]	617	& + subpdqcloudco2(ig,l,igcm_dust_number)
	618
[1911]	619	sum_subpdq(ig,l,igcm_ccnco2_mass) =
	620	& sum_subpdq(ig,l,igcm_ccnco2_mass)
[1816]	621	& + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
[1911]	622	sum_subpdq(ig,l,igcm_ccnco2_number) =
	623	& sum_subpdq(ig,l,igcm_ccnco2_number)
[1816]	624	& + subpdqcloudco2(ig,l,igcm_ccnco2_number)
	625
[1911]	626	sum_subpdq(ig,l,igcm_co2_ice) =
	627	& sum_subpdq(ig,l,igcm_co2_ice)
[1816]	628	& + subpdqcloudco2(ig,l,igcm_co2_ice)
[1911]	629	sum_subpdq(ig,l,igcm_co2) =
	630	& sum_subpdq(ig,l,igcm_co2)
[1816]	631	& + subpdqcloudco2(ig,l,igcm_co2)
	632
[1911]	633	sum_subpdq(ig,l,igcm_h2o_ice) =
	634	& sum_subpdq(ig,l,igcm_h2o_ice)
[1816]	635	& + subpdqcloudco2(ig,l,igcm_h2o_ice)
[1911]	636	sum_subpdq(ig,l,igcm_ccn_mass) =
	637	& sum_subpdq(ig,l,igcm_ccn_mass)
[1816]	638	& + subpdqcloudco2(ig,l,igcm_ccn_mass)
[1911]	639	sum_subpdq(ig,l,igcm_ccn_number) =
	640	& sum_subpdq(ig,l,igcm_ccn_number)
[1816]	641	& + subpdqcloudco2(ig,l,igcm_ccn_number)
[1820]	642	ENDDO
	643	ENDDO
[1720]	644	ENDDO ! of DO microstep=1,imicro
[1617]	645
	646	c------------------------------------------------
[1816]	647	c Compute final tendencies after time loop:
	648	c------------------------------------------------
[1617]	649	c CO2 flux at surface (kg.m-2.s-1)
	650	do ig=1,ngrid
[1911]	651	pdqs_sedco2(ig)=sum_subpdqs_sedco2(ig)/real(imicroco2)
[1617]	652	enddo
	653	c------ Temperature tendency pdtcloud
[1820]	654	DO l=1,nlay
	655	DO ig=1,ngrid
[1617]	656	pdtcloudco2(ig,l) =
[1911]	657	& sum_subpdt(ig,l)/real(imicroco2)-pdt(ig,l)
[1820]	658	ENDDO
	659	ENDDO
[1617]	660	c------ Tracers tendencies pdqcloud
[1820]	661	DO l=1,nlay
	662	DO ig=1,ngrid
[1816]	663	pdqcloudco2(ig,l,igcm_co2_ice) =
[1911]	664	& sum_subpdq(ig,l,igcm_co2_ice)/real(imicroco2)
[1816]	665	& - pdq(ig,l,igcm_co2_ice)
	666	pdqcloudco2(ig,l,igcm_co2) =
[1911]	667	& sum_subpdq(ig,l,igcm_co2)/real(imicroco2)
[1816]	668	& - pdq(ig,l,igcm_co2)
	669	pdqcloudco2(ig,l,igcm_h2o_ice) =
[1911]	670	& sum_subpdq(ig,l,igcm_h2o_ice)/real(imicroco2)
[1816]	671	& - pdq(ig,l,igcm_h2o_ice)
[1820]	672	ENDDO
	673	ENDDO
	674	DO l=1,nlay
	675	DO ig=1,ngrid
[1816]	676	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
[1911]	677	& sum_subpdq(ig,l,igcm_ccnco2_mass)/real(imicroco2)
[1816]	678	& - pdq(ig,l,igcm_ccnco2_mass)
	679	pdqcloudco2(ig,l,igcm_ccnco2_number) =
[1911]	680	& sum_subpdq(ig,l,igcm_ccnco2_number)/real(imicroco2)
[1816]	681	& - pdq(ig,l,igcm_ccnco2_number)
	682	pdqcloudco2(ig,l,igcm_ccn_mass) =
[1911]	683	& sum_subpdq(ig,l,igcm_ccn_mass)/real(imicroco2)
[1816]	684	& - pdq(ig,l,igcm_ccn_mass)
	685	pdqcloudco2(ig,l,igcm_ccn_number) =
[1911]	686	& sum_subpdq(ig,l,igcm_ccn_number)/real(imicroco2)
[1816]	687	& - pdq(ig,l,igcm_ccn_number)
[1820]	688	ENDDO
	689	ENDDO
	690	DO l=1,nlay
	691	DO ig=1,ngrid
[1816]	692	pdqcloudco2(ig,l,igcm_dust_mass) =
[1911]	693	& sum_subpdq(ig,l,igcm_dust_mass)/real(imicroco2)
[1816]	694	& - pdq(ig,l,igcm_dust_mass)
	695	pdqcloudco2(ig,l,igcm_dust_number) =
[1911]	696	& sum_subpdq(ig,l,igcm_dust_number)/real(imicroco2)
[1816]	697	& - pdq(ig,l,igcm_dust_number)
[1820]	698	ENDDO
	699	ENDDO
[1816]	700	c-------Due to stepped entry, other processes tendencies can add up to negative values
	701	c-------Therefore, enforce positive values and conserve mass
[1820]	702	DO l=1,nlay
	703	DO ig=1,ngrid
[1816]	704	IF ((pqeff(ig,l,igcm_ccnco2_number) +
	705	& ptimestep* (pdq(ig,l,igcm_ccnco2_number) +
	706	& pdqcloudco2(ig,l,igcm_ccnco2_number))
	707	& .lt. 1.)
	708	& .or. (pqeff(ig,l,igcm_ccnco2_mass) +
	709	& ptimestep* (pdq(ig,l,igcm_ccnco2_mass) +
	710	& pdqcloudco2(ig,l,igcm_ccnco2_mass))
	711	& .lt. 1.e-20)) THEN
	712	pdqcloudco2(ig,l,igcm_ccnco2_number) =
	713	& - pqeff(ig,l,igcm_ccnco2_number)/ptimestep
	714	& - pdq(ig,l,igcm_ccnco2_number)+1.
	715	pdqcloudco2(ig,l,igcm_dust_number) =
	716	& -pdqcloudco2(ig,l,igcm_ccnco2_number)
	717	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
	718	& - pqeff(ig,l,igcm_ccnco2_mass)/ptimestep
	719	& - pdq(ig,l,igcm_ccnco2_mass)+1.e-20
	720	pdqcloudco2(ig,l,igcm_dust_mass) =
	721	& -pdqcloudco2(ig,l,igcm_ccnco2_mass)
	722	ENDIF
[1820]	723	ENDDO
	724	ENDDO
	725	DO l=1,nlay
	726	DO ig=1,ngrid
[1816]	727	IF ( (pqeff(ig,l,igcm_dust_number) +
	728	& ptimestep* (pdq(ig,l,igcm_dust_number) +
	729	& pdqcloudco2(ig,l,igcm_dust_number)) .le. 1.)
	730	& .or. (pqeff(ig,l,igcm_dust_mass)+
	731	& ptimestep* (pdq(ig,l,igcm_dust_mass) +
	732	& pdqcloudco2(ig,l,igcm_dust_mass))
	733	& .le. 1.e-20)) then
	734	pdqcloudco2(ig,l,igcm_dust_number) =
	735	& - pqeff(ig,l,igcm_dust_number)/ptimestep
	736	& - pdq(ig,l,igcm_dust_number)+1.
	737	pdqcloudco2(ig,l,igcm_ccnco2_number) =
	738	& -pdqcloudco2(ig,l,igcm_dust_number)
	739	pdqcloudco2(ig,l,igcm_dust_mass) =
	740	& - pqeff(ig,l,igcm_dust_mass)/ptimestep
	741	& - pdq(ig,l,igcm_dust_mass) +1.e-20
	742	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
	743	& -pdqcloudco2(ig,l,igcm_dust_mass)
	744	ENDIF
[1820]	745	ENDDO
	746	ENDDO
[1816]	747	!pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq
[1820]	748	DO l=1,nlay
	749	DO ig=1,ngrid
[1816]	750	IF (pqeff(ig,l,igcm_co2_ice) + ptimestep*
	751	& (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice))
	752	& .lt. 1.e-15) THEN
	753	pdqcloudco2(ig,l,igcm_co2_ice) =
	754	& - pqeff(ig,l,igcm_co2_ice)/ptimestep-pdq(ig,l,igcm_co2_ice)
	755	pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice)
	756	ENDIF
	757	IF (pqeff(ig,l,igcm_co2) + ptimestep*
	758	& (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2))
	759	& .lt. 0.1) THEN
	760	pdqcloudco2(ig,l,igcm_co2) =
	761	& - pqeff(ig,l,igcm_co2)/ptimestep - pdq(ig,l,igcm_co2)
	762	pdqcloudco2(ig,l,igcm_co2_ice)= -pdqcloudco2(ig,l,igcm_co2)
	763	ENDIF
[1617]	764	ENDDO
[1820]	765	ENDDO
[1617]	766
[1816]	767	c Update clouds parameters values in the cloud fraction (for output)
[1820]	768	DO l=1, nlay
	769	DO ig=1,ngrid
[1617]	770
[1720]	771	Niceco2=pqeff(ig,l,igcm_co2_ice) +
	772	& (pdq(ig,l,igcm_co2_ice) +
	773	& pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep
	774	Nco2=pqeff(ig,l,igcm_co2) +
	775	& (pdq(ig,l,igcm_co2) +
	776	& pdqcloudco2(ig,l,igcm_co2))*ptimestep
	777	Nccnco2=max((pqeff(ig,l,igcm_ccnco2_number) +
	778	& (pdq(ig,l,igcm_ccnco2_number) +
[1816]	779	& pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep)
	780	& ,1.e-30)
[1720]	781	Qccnco2=max((pqeff(ig,l,igcm_ccnco2_mass) +
	782	& (pdq(ig,l,igcm_ccnco2_mass) +
[1816]	783	& pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep)
	784	& ,1.e-30)
[1720]	785
[1911]	786	myT=pteff(ig,l)+(pdt(ig,l)+pdtcloudco2(ig,l))*ptimestep
[1720]	787	rho_ice_co2T(ig,l)=1000.(1.72391-2.53e-4
	788	& myT-2.87e-6* myT* myT)
	789	rho_ice_co2=rho_ice_co2T(ig,l)
[1816]	790	c rho_ice_co2 is shared by tracer_mod and used in updaterice
	791	c Compute particle size
	792	call updaterice_microco2(Niceco2,
	793	& Qccnco2,Nccnco2,
	794	& tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
	795
	796	if ( (Niceco2 .le. 1.e-25 .or.
	797	& Nccnco2*tauscaling(ig) .le. 1.) )THEN
	798	riceco2(ig,l)=0.
[1820]	799	Qext1bins2(ig,l)=0.
	800	else
	801	c Compute opacities
	802	No=Nccnco2*tauscaling(ig)
	803	Rn=-dlog(riceco2(ig,l))
	804	n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
	805	Qext1bins2(ig,l)=0.
	806	do i = 1, nbinco2_cld
	807	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
	808	n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
	809	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
	810	Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i)
	811	enddo
[1816]	812	endif
	813
[1685]	814	!update rice water
[1820]	815	call updaterice_micro(
[1720]	816	& pqeff(ig,l,igcm_h2o_ice) + ! ice mass
[1685]	817	& (pdq(ig,l,igcm_h2o_ice) + ! ice mass
	818	& pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep, ! ice mass
[1720]	819	& pqeff(ig,l,igcm_ccn_mass) + ! ccn mass
[1685]	820	& (pdq(ig,l,igcm_ccn_mass) + ! ccn mass
	821	& pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep, ! ccn mass
[1720]	822	& pqeff(ig,l,igcm_ccn_number) + ! ccn number
[1685]	823	& (pdq(ig,l,igcm_ccn_number) + ! ccn number
	824	& pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
	825	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
	826
[1820]	827	call updaterdust(
[1720]	828	& pqeff(ig,l,igcm_dust_mass) + ! dust mass
[1685]	829	& (pdq(ig,l,igcm_dust_mass) + ! dust mass
	830	& pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, ! dust mass
[1720]	831	& pqeff(ig,l,igcm_dust_number) + ! dust number
[1685]	832	& (pdq(ig,l,igcm_dust_number) + ! dust number
	833	& pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number
	834	& rdust(ig,l))
	835
[1820]	836	ENDDO
	837	ENDDO
[1720]	838
[1617]	839	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
	840	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	841	c Then that should not affect the ice particle radius
[1816]	842
[1820]	843	do ig=1,ngrid
	844	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
[1816]	845	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
	846	& riceco2(ig,2)=riceco2(ig,3)
	847	riceco2(ig,1)=riceco2(ig,2)
[1820]	848	endif
[1617]	849	end do
	850
[1816]	851	DO l=1,nlay
[1617]	852	DO ig=1,ngrid
[1685]	853	rsedcloud(ig,l)=max(rice(ig,l)*
	854	& (1.+nuice_sed)(1.+nuice_sed)(1.+nuice_sed),
	855	& rdust(ig,l))
	856	! rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
	857	ENDDO
	858	ENDDO
	859
[1820]	860	DO l=1,nlay
[1685]	861	DO ig=1,ngrid
[1617]	862	rsedcloudco2(ig,l)=max(riceco2(ig,l)*
	863	& (1.+nuiceco2_sed)(1.+nuiceco2_sed)(1.+nuiceco2_sed),
	864	& rdust(ig,l))
[1720]	865	c rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5)
[1617]	866	ENDDO
[1820]	867	ENDDO
[1617]	868
[1911]	869	call co2sat(ngridnlay,pteff+(pdt+pdtcloudco2)ptimestep
[1720]	870	& ,pplay,zqsatco2)
[1820]	871	do l=1,nlay
	872	do ig=1,ngrid
[1720]	873	satuco2(ig,l) = (pqeff(ig,l,igcm_co2) +
	874	& (pdq(ig,l,igcm_co2) +
	875	& pdqcloudco2(ig,l,igcm_co2))ptimestep)
	876	& (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)
[1820]	877	enddo
	878	enddo
[1885]	879	!Everything modified by CO2 microphysics must be wrt co2cloudfrac
[1820]	880	IF (CLFvaryingCO2) THEN
[1885]	881	DO l=1,nlay
	882	DO ig=1,ngrid
	883	pdqcloudco2(ig,l,igcm_ccn_mass)=
	884	& pdqcloudco2(ig,l,igcm_ccn_mass)*co2cloudfrac(ig,l)
	885	pdqcloudco2(ig,l,igcm_ccnco2_mass)=
	886	& pdqcloudco2(ig,l,igcm_ccnco2_mass)*co2cloudfrac(ig,l)
	887	pdqcloudco2(ig,l,igcm_ccn_number)=
	888	& pdqcloudco2(ig,l,igcm_ccn_number)*co2cloudfrac(ig,l)
	889	pdqcloudco2(ig,l,igcm_ccnco2_number)=
	890	& pdqcloudco2(ig,l,igcm_ccnco2_number)*co2cloudfrac(ig,l)
	891	pdqcloudco2(ig,l,igcm_dust_mass)=
	892	& pdqcloudco2(ig,l,igcm_dust_mass)*co2cloudfrac(ig,l)
	893	pdqcloudco2(ig,l,igcm_dust_number)=
	894	& pdqcloudco2(ig,l,igcm_dust_number)*co2cloudfrac(ig,l)
	895	pdqcloudco2(ig,l,igcm_h2o_ice)=
	896	& pdqcloudco2(ig,l,igcm_h2o_ice)*co2cloudfrac(ig,l)
	897	pdqcloudco2(ig,l,igcm_co2_ice)=
	898	& pdqcloudco2(ig,l,igcm_co2_ice)*co2cloudfrac(ig,l)
	899	pdqcloudco2(ig,l,igcm_co2)=
	900	& pdqcloudco2(ig,l,igcm_co2)*co2cloudfrac(ig,l)
	901	pdtcloudco2(ig,l)=pdtcloudco2(ig,l)*co2cloudfrac(ig,l)
	902	Qext1bins2(ig,l)=Qext1bins2(ig,l)*co2cloudfrac(ig,l)
	903	ENDDO
	904	ENDDO
[1820]	905	ENDIF
	906	! opacity in mesh ig is the sum over l of Qext1bins2: Is this true ?
	907	tau1mic(:)=0.
	908	do l=1,nlay
	909	do ig=1,ngrid
	910	tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l)
	911	enddo
	912	enddo
[1816]	913	!Outputs:
[1820]	914	call WRITEDIAGFI(ngrid,"SatIndex","SatIndex"," ",3,
[1816]	915	& SatIndex)
[1820]	916	call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
[1629]	917	& satuco2)
[1820]	918	call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
[1816]	919	& ,3,riceco2)
[1885]	920	call WRITEdiagfi(ngrid,"co2cloudfrac","co2 cloud fraction"
	921	& ," ",3,co2cloudfrac)
[1820]	922	call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2"
[1816]	923	& ,"m",3,rsedcloudco2)
[1820]	924	call WRITEdiagfi(ngrid,"Tau3D1mic"," co2 ice opacities"
[1816]	925	& ," ",3,Qext1bins2)
[1820]	926	call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"
[1816]	927	& ," ",2,tau1mic)
[1617]	928
[1820]	929	END

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

Download in other formats:

Original Format