Context Navigation

co2cloud.F @ 2228

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

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

Download in other formats:

Original Format