Context Navigation

co2cloud.F @ 1962

Last change on this file since 1962 was 1938, checked in by aslmd, 7 years ago
cancel commit 1935 will all due apologies. for Mars mesoscale applications do not use versions later than LMDZ.MARS rev 1779. to be fixed in the future.
Property svn:executable set to ``*
File size: 40.6 KB

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

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

Download in other formats:

Original Format