Context Navigation

co2cloud.F @ 1775

Last change on this file since 1775 was 1754, checked in by aslmd, 8 years ago
LMDZ.MARS : tcondco2 is not present, prevents the model to compile. commented the line.
Property svn:executable set to ``*
File size: 44.7 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,
[1720]	6	& rsedcloud,rhocloud,pzlev,pdqs_sedco2)
[1617]	7	! to use 'getin'
	8	use dimradmars_mod, only: naerkind
	9	USE comcstfi_h
	10	USE ioipsl_getincom
	11	USE updaterad
	12	use conc_mod, only: mmean
	13	use tracer_mod, only: nqmx, igcm_co2, igcm_co2_ice,
[1685]	14	& igcm_dust_mass, igcm_dust_number,igcm_h2o_ice,
	15	& igcm_ccn_mass,igcm_ccn_number,
[1617]	16	& igcm_ccnco2_mass, igcm_ccnco2_number,
	17	& rho_dust, nuiceco2_sed, nuiceco2_ref,
[1720]	18	& rho_ice_co2,r3n_q,rho_ice,nuice_sed
	19
[1651]	20
[1617]	21	IMPLICIT NONE
	22
	23
[1720]	24	#include "datafile.h"
	25	#include "callkeys.h"
	26	#include "microphys.h"
	27
	28
[1617]	29	c=======================================================================
	30	c CO2 clouds formation
	31	c
	32	c There is a time loop specific to cloud formation
	33	c due to timescales smaller than the GCM integration timestep.
	34	c microphysics subroutine is improvedCO2clouds.F
	35	c
	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
	45	c Several logicals can be set to .true. in callphys.def
	46	c co2clouds=.true. call this routine
	47	c co2useh2o=.true. allow the use of water ice particles as CCN for CO2
	48	c meteo_flux=.true. supply meteoritic particles
	49	c CLFvaryingCO2=.true. allows a subgrid temperature distribution
	50	c of amplitude spantCO2(=integer in callphys.def)
	51	c
[1617]	52	c=======================================================================
	53
	54	c-----------------------------------------------------------------------
	55	c declarations:
	56	c -------------
	57
	58	c Inputs:
	59	c ------
	60
	61	INTEGER ngrid,nlay
	62	INTEGER nq ! nombre de traceurs
	63	REAL ptimestep ! pas de temps physique (s)
	64	REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
	65	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
	66	REAL pdpsrf(ngrid) ! tendence surf pressure
	67	REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers
	68	REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K)
	69	REAL pdt(ngrid,nlay) ! tendence temperature des autres param.
[1720]	70	real,intent(in) :: pzlev(ngrid,nlay+1) ! altitude at the boundaries of the layers
[1617]	71
	72	real pq(ngrid,nlay,nq) ! traceur (kg/kg)
	73	real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1)
	74
	75	real rice(ngrid,nlay) ! Water Ice mass mean radius (m)
	76	! used for nucleation of CO2 on ice-coated ccns
[1685]	77	DOUBLE PRECISION rho_ice_co2T(ngrid,nlay)
[1617]	78	REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
	79	REAL tauscaling(ngrid) ! Convertion factor for dust amount
	80	real rdust(ngrid,nlay) ! Dust geometric mean radius (m)
	81
	82	c Outputs:
	83	c -------
	84
	85	real pdqcloudco2(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
	86	REAL pdtcloudco2(ngrid,nlay) ! tendence temperature due
	87	! a la chaleur latente
	88
[1629]	89	DOUBLE PRECISION riceco2(ngrid,nlay) ! Ice mass mean radius (m)
[1617]	90	! (r_c in montmessin_2004)
	91	REAL nuice(ngrid,nlay) ! Estimated effective variance
	92	! of the size distribution
	93	real rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius
	94	real rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
	95	real rhocloudco2t(ngrid,nlay) ! Cloud density (kg.m-3)
	96	real pdqs_sedco2(ngrid) ! CO2 flux at the surface
	97	c local:
	98	c ------
[1685]	99	!water
	100	real rsedcloud(ngrid,nlay) ! Cloud sedimentation radius
	101	real rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
[1617]	102	! for ice radius computation
[1720]	103
[1617]	104	REAl ccntyp
	105
	106	! for time loop
	107	INTEGER microstep ! time subsampling step variable
	108	INTEGER imicro ! time subsampling for coupled water microphysics & sedimentation
	109	SAVE imicro
	110	REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
	111	SAVE microtimestep
	112
	113	! tendency given by clouds (inside the micro loop)
	114	REAL subpdqcloudco2(ngrid,nlay,nq) ! cf. pdqcloud
	115	REAL subpdtcloudco2(ngrid,nlay) ! cf. pdtcloud
	116
	117	! global tendency (clouds+physics)
	118	REAL subpdq(ngrid,nlay,nq) ! cf. pdqcloud
	119	REAL subpdt(ngrid,nlay) ! cf. pdtcloud
	120	real wq(ngrid,nlay+1) ! ! displaced tracer mass (kg.m-2) during microtimestep because sedim (?/m-2)
	121
	122	REAL satuco2(ngrid,nlay) ! co2 satu ratio for output
	123	REAL zqsatco2(ngrid,nlay) ! saturation co2
	124
[1720]	125	INTEGER iq,ig,l,i
[1617]	126	LOGICAL,SAVE :: firstcall=.true.
[1720]	127	DOUBLE PRECISION Nccnco2, Niceco2,Nco2,mdustJA,ndustJA
[1617]	128	DOUBLE PRECISION Qccnco2
	129	real :: beta
	130
	131	real epaisseur (ngrid,nlay) ! Layer thickness (m)
	132	real masse (ngrid,nlay) ! Layer mass (kg.m-2)
[1651]	133	double precision diff,diff0
[1617]	134	real tempo_traceur_t(ngrid,nlay)
	135	real tempo_traceurs(ngrid,nlay,nq)
	136	real sav_trac(ngrid,nlay,nq)
	137	real pdqsed(ngrid,nlay,nq)
[1720]	138	REAL lw !Latent heat of sublimation (J.kg-1)
	139	REAL,save :: l0,l1,l2,l3,l4 !Coeffs for lw
[1685]	140
[1720]	141	DOUBLE PRECISION,allocatable,save :: memdMMccn(:,:) ! Nb particules intégré
[1685]	142	DOUBLE PRECISION,allocatable,save :: memdMMh2o(:,:)
	143	DOUBLE PRECISION,allocatable,save :: memdNNccn(:,:)
[1720]	144	DOUBLE PRECISION :: myT
	145	! What we need for Qext reading and tau computation
	146	DOUBLE PRECISION vrat_cld ! Volume ratio
	147	DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
	148	SAVE rb_cldco2
	149	DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-11 ! Minimum radius (m)
	150	DOUBLE PRECISION, PARAMETER :: rmax_cld = 3.e-6 ! Maximum radius (m)
	151	DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-12
	152	! Minimum boundary radius (m)
	153	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 5.e-6 ! Maximum boundary radius (m)
	154
	155	DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m)
	156	DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3)
	157	REAL sigma_iceco2 ! Variance of the ice and CCN distributions
	158
	159	logical :: file_ok
	160	double precision :: radv(10000),Qextv1mic(10000)
	161	double precision :: Qext1bins(100),Qtemp
	162	double precision :: ltemp1(10000),ltemp2(10000)
	163	integer :: nelem,lebon1,lebon2,uQext
	164	save Qext1bins
	165	character(len=100) scanline
	166	DOUBLE PRECISION n_aer(nbinco2_cld),Rn,No,n_derf,dev2
	167	DOUBLE PRECISION Qext1bins2(ngrid,nlay)
	168	DOUBLE PRECISION tau1mic(ngrid) !co2 ice column opacity at 1µm
	169
	170	! For sub grid T distribution
	171
	172	REAL zt(ngrid,nlay) ! local value of temperature
	173	REAL :: zq(ngrid, nlay,nq)
	174
	175	REAL :: tcond(ngrid,nlay)
	176	REAL :: zqvap(ngrid,nlay)
	177	REAL :: zqice(ngrid,nlay)
	178	REAL :: spant,zdelt ! delta T for the temperature distribution
	179	REAL :: zteff(ngrid, nlay)! effective temperature in the cloud,neb
	180	REAL :: pqeff(ngrid, nlay, nq)! effective tracers quantities in the cloud
	181	REAL :: cloudfrac(ngrid,nlay) ! cloud fraction
	182	REAL :: mincloud ! min cloud frac
	183	c logical :: CLFvaryingCO2
[1617]	184	c ** un petit test de coherence
	185	c --------------------------
	186
	187	IF (firstcall) THEN
	188	if (nq.gt.nqmx) then
	189	write(,) 'stop in co2cloud (nq.gt.nqmx)!'
	190	write(,) 'nq=',nq,' nqmx=',nqmx
	191	stop
	192	endif
[1685]	193	write(,) "co2cloud.F: rho_ice_co2 = ",rho_ice_co2
[1617]	194	write(,) "co2cloud: igcm_co2=",igcm_co2
	195	write(,) " igcm_co2_ice=",igcm_co2_ice
	196
	197	write(,) "time subsampling for microphysic ?"
	198	#ifdef MESOSCALE
	199	imicro = 2
	200	#else
	201	imicro = 30
	202	#endif
	203	call getin("imicro",imicro)
	204	write(,)"imicro = ",imicro
	205
	206	microtimestep = ptimestep/real(imicro)
	207	write(,)"Physical timestep is",ptimestep
	208	write(,)"CO2 Microphysics timestep is",microtimestep
[1651]	209
[1720]	210	! Values for latent heat computation
	211	l0=595594d0
	212	l1=903.111d0
	213	l2=-11.5959d0
	214	l3=0.0528288d0
	215	l4=-0.000103183d0
	216	c$$$
	217	c$$$ if (meteo_flux_number .ne. 0) then
	218	c$$$
	219	c$$$ write(,) "Warning ! Meteoritic flux of dust is turned on"
	220	c$$$ write(,) "Dust mass = ",meteo_flux_mass
	221	c$$$ write(,) "Dust number = ",meteo_flux_number
	222	c$$$ write(,) "Are added at the z-level (km)",meteo_alt
	223	c$$$ write(,) "Every timestep in co2cloud.F"
	224	c$$$ endif
	225
[1685]	226	if (.not. allocated(memdMMccn)) allocate(memdMMccn(ngrid,nlay))
	227	if (.not. allocated(memdNNccn)) allocate(memdNNccn(ngrid,nlay))
	228	if (.not. allocated(memdMMh2o)) allocate(memdMMh2o(ngrid,nlay))
	229
	230	memdMMccn(:,:)=0.
	231	memdMMh2o(:,:)=0.
	232	memdNNccn(:,:)=0.
[1720]	233
	234	c Compute the size bins of the distribution of CO2 ice particles
	235	c --> used for opacity calculations
	236
	237	c rad_cldco2 is the primary radius grid used for microphysics computation.
	238	c The grid spacing is computed assuming a constant volume ratio
	239	c between two consecutive bins; i.e. vrat_cld.
	240	c vrat_cld is determined from the boundary values of the size grid:
	241	c rmin_cld and rmax_cld.
	242	c The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
	243	c dr_cld is the width of each rad_cldco2 bin.
	244	sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
	245	c Volume ratio between two adjacent bins
	246	! vrat_cld
	247	vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
	248	vrat_cld = exp(vrat_cld)
	249	rb_cldco2(1) = rbmin_cld
	250	rad_cldco2(1) = rmin_cld
	251	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
	252	do i=1,nbinco2_cld-1
	253	rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.)
	254	vol_cld(i+1) = vol_cld(i) * vrat_cld
	255	enddo
	256	do i=1,nbinco2_cld
	257	rb_cldco2(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
	258	& rad_cldco2(i)
	259	dr_cld(i) = rb_cldco2(i+1) - rb_cldco2(i)
	260	enddo
	261	rb_cldco2(nbinco2_cld+1) = rbmax_cld
	262	dr_cld(nbinco2_cld) = rb_cldco2(nbinco2_cld+1) -
	263	& rb_cldco2(nbinco2_cld)
	264
	265	c read the Qext values
	266	INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))//
	267	& '/optprop_co2ice_1mic.dat', EXIST=file_ok)
	268	IF (.not. file_ok) THEN
	269	write(,) 'file optprop_co2ice_1mic.dat should be in '
	270	& ,datafile
	271	STOP
	272	endif
	273	open(newunit=uQext,file=trim(datafile)//
	274	& '/optprop_co2ice_1mic.dat'
	275	& ,FORM='formatted')
	276	read(uQext,*) !skip 1 line
	277	do i=1,10000
	278	read(uQext,'(E11.5)') radv(i)
	279	enddo
	280	read(uQext,*) !skip 1 line
	281	do i=1,10000
	282	read(uQext,'(E11.5)') Qextv1mic(i)
	283	enddo
	284	close(uQext)
	285	c innterpol the Qext values
	286	!rice_out=rad_cldco2
	287	do i=1,nbinco2_cld
	288	ltemp1=abs(radv(:)-rb_cldco2(i))
	289	ltemp2=abs(radv(:)-rb_cldco2(i+1))
	290	lebon1=minloc(ltemp1,DIM=1)
	291	lebon2=minloc(ltemp2,DIM=1)
	292	nelem=lebon2-lebon1+1.
	293	Qtemp=0d0
	294	do l=0,nelem
	295	Qtemp=Qtemp+Qextv1mic(lebon1+l) !mean value in the interval
	296	enddo
	297	Qtemp=Qtemp/nelem
	298	Qext1bins(i)=Qtemp
	299	enddo
	300	Qext1bins(:)=Qext1bins(:)rad_cldco2(:)rad_cldco2(:)*pi
	301	! The actuall tau computation and output is performed in co2cloud.F
	302
	303	print*,'--------------------------------------------'
	304	print*,'Microphysics co2: size bin-Qext information:'
	305	print*,' i, rad_cldco2(i), Qext1bins(i)'
	306	do i=1,nbinco2_cld
	307	write(*,'(i3,3x,3(e12.6,4x))') i, rad_cldco2(i),
	308	& Qext1bins(i)
	309	enddo
	310	print*,'--------------------------------------------'
	311
	312
	313	do i=1,nbinco2_cld+1
	314	rb_cldco2(i) = log(rb_cldco2(i)) !! we save that so that it is not computed at each timestep and gridpoint
	315	enddo
	316	if (CLFvaryingCO2) then
	317	write(,)
	318	write(,) "CLFvaryingCO2 is set to true is callphys.def"
	319	write(,) "The temperature field is enlarged to +/-",spantCO2
	320	write(,) "for the CO2 microphysics "
	321	endif
[1617]	322	firstcall=.false.
	323	ENDIF ! of IF (firstcall)
	324
	325	c-----Initialization
[1720]	326	dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
[1685]	327	beta=0.85
	328	subpdq(1:ngrid,1:nlay,1:nq) = 0
	329	subpdt(1:ngrid,1:nlay) = 0
	330	subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0
	331	subpdtcloudco2(1:ngrid,1:nlay) = 0
[1720]	332	c$$$
	333	c$$$ call WRITEDIAGFI(ngrid,"pzlay","pzlay","km",3,
	334	c$$$ & pzlay)
	335	c$$$ call WRITEDIAGFI(ngrid,"pplay","pplay","Pa",3,
	336	c$$$ & pplay)
[1685]	337
[1617]	338	wq(:,:)=0
	339	! default value if no ice
	340	rhocloudco2(1:ngrid,1:nlay) = rho_dust
	341	rhocloudco2t(1:ngrid,1:nlay) = rho_dust
	342	epaisseur(1:ngrid,1:nlay)=0
	343	masse(1:ngrid,1:nlay)=0
	344
	345	sav_trac(1:ngrid,1:nlay,1:nq)=0
	346	pdqsed(1:ngrid,1:nlay,1:nq)=0
	347
	348	do l=1,nlay
	349	do ig=1, ngrid
	350	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
[1720]	351	epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
[1617]	352	enddo
	353	enddo
	354
[1720]	355	!CLFvaryingCO2=.true.
[1617]	356
	357	c-------------------------------------------------------------------
[1720]	358	c 0. Representation of sub-grid water ice clouds
	359	c------------------
	360	IF (CLFvaryingCO2) THEN
	361	spant=spantCO2 ! delta T for the temprature distribution
	362	mincloud=0.1 ! min cloudfrac when there is ice
	363	! IF (flagcloudco2) THEN
	364	! write(,) "CLFCO2 Delta T", spant
	365	! write(,) "CLFCO2 mincloud", mincloud
	366	! flagcloudco2=.false.
	367	! END IF
	368
	369	c-----Tendencies
	370	DO l=1,nlay
	371	DO ig=1,ngrid
	372	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	373	ENDDO
	374	ENDDO
	375	DO l=1,nlay
	376	DO ig=1,ngrid
	377	DO iq=1,nq
	378	zq(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
	379	ENDDO
	380	ENDDO
	381	ENDDO
	382	zqvap=zq(:,:,igcm_co2)
	383	zqice=zq(:,:,igcm_co2_ice)
[1754]	384	!! AS : this routine is not present?
	385	! CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
[1720]	386	! A tester: CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
	387	zdelt=spant
	388	DO l=1,nlay
	389	DO ig=1,ngrid
	390	IF (tcond(ig,l) .ge. (zt(ig,l)+zdelt)) THEN !Toute la fraction est saturée
	391	zteff(ig,l)=zt(ig,l)
	392	cloudfrac(ig,l)=1.
	393	ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN !Rien n'est saturé
	394	zteff(ig,l)=zt(ig,l)-zdelt
	395	cloudfrac(ig,l)=mincloud
	396	ELSE
	397	cloudfrac(ig,l)=(tcond(ig,l)-zt(ig,l)+zdelt)/
	398	& (2.0*zdelt)
	399	zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Ja Temperature moyenne de la fraction nuageuse
	400
	401	END IF
	402	zteff(ig,l)=zteff(ig,l)-pdt(ig,l)*ptimestep
	403	IF (cloudfrac(ig,l).le. mincloud) THEN
	404	cloudfrac(ig,l)=mincloud
	405	ELSE IF (cloudfrac(ig,l).ge. 1) THEN
	406	cloudfrac(ig,l)=1.
	407	END IF
	408	ENDDO
	409	ENDDO
	410	! Totalcloud frac of the column missing here
	411	c-----------------------
	412	c-----No sub-grid cloud representation (CLFvarying=false)
	413	ELSE
	414	DO l=1,nlay
	415	DO ig=1,ngrid
	416	zteff(ig,l)=pt(ig,l)
	417	END DO
	418	END DO
	419	END IF ! end if (CLFvaryingco2)--------------------------------------------
[1617]	420	c 1. Tendencies:
	421	c------------------
	422
[1720]	423	c------ Effective tracers quantities in the cloud fraction
	424	IF (CLFvaryingCO2) THEN
	425	pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier)
	426
	427	c pqeff(:,:,igcm_ccn_mass) =pq(:,:,igcm_ccn_mass)/
	428	c & cloudfrac(:,:)
	429	c pqeff(:,:,igcm_ccn_number)=
	430	c & pq(:,:,igcm_ccn_number)/cloudfrac(:,:)
	431	c pqeff(:,:,igcm_h2o_ice)= pq(:,:,igcm_h2o_ice)/
	432	c & cloudfrac(:,:)
	433	pqeff(:,:,igcm_ccnco2_mass) =pq(:,:,igcm_ccnco2_mass)/
	434	& cloudfrac(:,:)
	435	pqeff(:,:,igcm_ccnco2_number)=
	436	& pq(:,:,igcm_ccnco2_number)/cloudfrac(:,:)
	437	pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)/
	438	& cloudfrac(:,:)
[1617]	439
[1720]	440	ELSE
	441	pqeff(:,:,:)=pq(:,:,:)
	442	c pqeff(:,:,igcm_ccn_mass)= pq(:,:,igcm_ccn_mass)
	443	c pqeff(:,:,igcm_ccn_number)= pq(:,:,igcm_ccn_number)
	444	c pqeff(:,:,igcm_h2o_ice)= pq(:,:,igcm_h2o_ice)
	445	c pqeff(:,:,igcm_ccnco2_mass)= pq(:,:,igcm_ccnco2_mass)
	446	c pqeff(:,:,igcm_ccnco2_number)= pq(:,:,igcm_ccnco2_number)
	447	c pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)
	448	END IF
	449	tempo_traceur_t(1:ngrid,1:nlay)=zteff(1:ngrid,1:nlay)
	450	tempo_traceurs(1:ngrid,1:nlay,1:nq)=pqeff(1:ngrid,1:nlay,1:nq)
[1617]	451
	452	c------------------------------------------------------------------
	453	c Time subsampling for microphysics
	454	c------------------------------------------------------------------
	455	DO microstep=1,imicro
	456	c------ Temperature tendency subpdt
	457	! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
	458	! If imicro=1 subpdt is the same as pdt
	459	DO l=1,nlay
	460	DO ig=1,ngrid
	461
	462	subpdt(ig,l) = subpdt(ig,l)
	463	& + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
	464	subpdq(ig,l,igcm_dust_mass) =
	465	& subpdq(ig,l,igcm_dust_mass)
	466	& + pdq(ig,l,igcm_dust_mass)
	467	subpdq(ig,l,igcm_dust_number) =
	468	& subpdq(ig,l,igcm_dust_number)
	469	& + pdq(ig,l,igcm_dust_number)
	470	subpdq(ig,l,igcm_ccnco2_mass) =
	471	& subpdq(ig,l,igcm_ccnco2_mass)
	472	& + pdq(ig,l,igcm_ccnco2_mass)
	473	subpdq(ig,l,igcm_ccnco2_number) =
	474	& subpdq(ig,l,igcm_ccnco2_number)
	475	& + pdq(ig,l,igcm_ccnco2_number)
	476	subpdq(ig,l,igcm_co2_ice) =
	477	& subpdq(ig,l,igcm_co2_ice)
	478	& + pdq(ig,l,igcm_co2_ice)
	479	subpdq(ig,l,igcm_co2) =
	480	& subpdq(ig,l,igcm_co2)
	481	& + pdq(ig,l,igcm_co2)
[1685]	482	subpdq(ig,l,igcm_h2o_ice) =
	483	& subpdq(ig,l,igcm_h2o_ice)
	484	& + pdq(ig,l,igcm_h2o_ice)
	485	subpdq(ig,l,igcm_ccn_mass) =
	486	& subpdq(ig,l,igcm_ccn_mass)
	487	& + pdq(ig,l,igcm_ccn_mass)
	488	subpdq(ig,l,igcm_ccn_number) =
	489	& subpdq(ig,l,igcm_ccn_number)
	490	& + pdq(ig,l,igcm_ccn_number)
[1617]	491	ENDDO
	492	ENDDO
[1720]	493
	494	c add meteoritic flux of dust (old version)
	495	cNew version (John Plane values) are added in improvedCO2clouds
	496	!convert meteo_alt (in km) to z-level
	497	!pzlay altitudes of the layers
	498	c$$$!zlev altitudes (nlayl+1) of the boundaries
	499	c$$$ if (meteo_flux_number .ge. 0) then
	500	c$$$ do ig=1, ngrid
	501	c$$$ l=1
	502	c$$$ do while ( (((meteo_alt .ge. pplev(ig,l)) .and.
	503	c$$$ & (meteo_alt .le. pplev(ig,l+1))) .eq. .false.)
	504	c$$$ & .and. (l .lt. nlay) )
	505	c$$$ l=l+1
	506	c$$$ enddo
	507	c$$$ meteo_lvl=l
	508	c$$$ subpdq(ig,meteo_lvl,igcm_dust_mass)=
	509	c$$$ & subpdq(ig,meteo_lvl,igcm_dust_mass)
	510	c$$$ & +meteo_flux_mass
	511	c$$$ subpdq(ig,meteo_lvl,igcm_dust_number)=
	512	c$$$ & subpdq(ig,meteo_lvl,igcm_dust_number)
	513	c$$$ & +meteo_flux_number
	514	c$$$ enddo
	515	c$$$ endif
	516	c-------------------------------------------------------------------
	517	c 2. Main call to the cloud schemes:
	518	c------------------------------------------------
	519	CALL improvedCO2clouds(ngrid,nlay,microtimestep,
	520	& pplay,pzlev,zteff,subpdt,
	521	& pqeff,subpdq,subpdqcloudco2,subpdtcloudco2,
	522	& nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn)
	523	c-------------------------------------------------------------------
	524	c 3. Updating tendencies after cloud scheme:
	525	c-----------------------------------------------
	526	DO l=1,nlay
	527	DO ig=1,ngrid
	528	subpdt(ig,l) =
	529	& subpdt(ig,l) + subpdtcloudco2(ig,l)
	530	subpdq(ig,l,igcm_dust_mass) =
	531	& subpdq(ig,l,igcm_dust_mass)
	532	& + subpdqcloudco2(ig,l,igcm_dust_mass)
	533	subpdq(ig,l,igcm_dust_number) =
	534	& subpdq(ig,l,igcm_dust_number)
	535	& + subpdqcloudco2(ig,l,igcm_dust_number)
	536	subpdq(ig,l,igcm_ccnco2_mass) =
	537	& subpdq(ig,l,igcm_ccnco2_mass)
	538	& + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
	539	subpdq(ig,l,igcm_ccnco2_number) =
	540	& subpdq(ig,l,igcm_ccnco2_number)
	541	& + subpdqcloudco2(ig,l,igcm_ccnco2_number)
	542	subpdq(ig,l,igcm_co2_ice) =
	543	& subpdq(ig,l,igcm_co2_ice)
	544	& + subpdqcloudco2(ig,l,igcm_co2_ice)
	545	subpdq(ig,l,igcm_co2) =
	546	& subpdq(ig,l,igcm_co2)
	547	& + subpdqcloudco2(ig,l,igcm_co2)
	548	subpdq(ig,l,igcm_h2o_ice) =
	549	& subpdq(ig,l,igcm_h2o_ice)
	550	& + subpdqcloudco2(ig,l,igcm_h2o_ice)
	551	subpdq(ig,l,igcm_ccn_mass) =
	552	& subpdq(ig,l,igcm_ccn_mass)
	553	& + subpdqcloudco2(ig,l,igcm_ccn_mass)
	554	subpdq(ig,l,igcm_ccn_number) =
	555	& subpdq(ig,l,igcm_ccn_number)
	556	& + subpdqcloudco2(ig,l,igcm_ccn_number)
	557	ENDDO
	558	ENDDO
	559
	560	!Sedimentation
	561	!Update traceurs, compute rsed
[1617]	562
	563	DO l=1, nlay
[1720]	564	DO ig=1,ngrid
	565	tempo_traceur_t(ig,l)=zteff(ig,l)+subpdt(ig,l)
	566	& *microtimestep
	567	tempo_traceurs(ig,l,:)=pqeff(ig,l,:)
	568	& +subpdq(ig,l,:)*microtimestep
	569
[1685]	570	rho_ice_co2T(ig,l)=1000.(1.72391-2.53e-4
	571	& tempo_traceur_t(ig,l)-2.87e-6*
	572	& tempo_traceur_t(ig,l)*tempo_traceur_t(ig,l))
[1720]	573
	574	rho_ice_co2=rho_ice_co2T(ig,l)
[1649]	575	Niceco2=max(tempo_traceurs(ig,l,igcm_co2_ice),1.e-30)
[1617]	576	Nccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_number),
	577	& 1.e-30)
	578	Qccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_mass),
	579	& 1.e-30)
	580	mdustJA= tempo_traceurs(ig,l,igcm_dust_mass)
	581	ndustJA=tempo_traceurs(ig,l,igcm_dust_number)
[1685]	582	if ((Nccnco2 .lt. tauscaling(ig)) .or. (Qccnco2 .lt.
[1617]	583	& 1.e-30 *tauscaling(ig))) then
	584	rdust(ig,l)=1.e-10
	585	else
[1685]	586	rdust(ig,l)=(3./4./pi/2500.Qccnco2/Nccnco2)*(1./3.)
	587	rdust(ig,l)=max(rdust(ig,l),1.e-10)
[1720]	588	c rdust(ig,l)=min(rdust(ig,l),5.e-6)
[1617]	589	end if
[1685]	590	rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2
	591	& + Qccnco2tauscaling(ig)rho_dust)
	592	& / (Niceco2 + Qccnco2*tauscaling(ig))
[1720]	593
	594	rhocloudco2t(ig,l)=min(max(rhocloudco2t(ig,l)
	595	& ,rho_ice_co2),rho_dust)
[1685]	596	riceco2(ig,l)=(Niceco2*3.0/
	597	& (4.0rho_ice_co2pi*Nccnco2
	598	& tauscaling(ig)) +rdust(ig,l)rdust(ig,l)
	599	& rdust(ig,l))*(1.0/3.0)
	600	riceco2(ig,l)=max(1.e-10,riceco2(ig,l))
[1720]	601
	602	rsedcloudco2(ig,l)=max(riceco2(ig,l)*
[1617]	603	& (1.+nuiceco2_sed)(1.+nuiceco2_sed)(1.+nuiceco2_sed),
	604	& rdust(ig,l))
[1720]	605
[1617]	606	ENDDO
	607	ENDDO
	608
	609	! Gravitational sedimentation
[1720]	610
[1617]	611	sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice)
	612	sav_trac(:,:,igcm_ccnco2_mass)=
	613	& tempo_traceurs(:,:,igcm_ccnco2_mass)
	614	sav_trac(:,:,igcm_ccnco2_number)=
	615	& tempo_traceurs(:,:,igcm_ccnco2_number)
	616
	617	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
	618	& microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
	619	& rsedcloudco2,rhocloudco2t,
	620	& tempo_traceurs(:,:,igcm_co2_ice),wq,beta) ! 3 traceurs
	621
[1720]	622	! sedim at the surface of co2 ice : keep track of it for physiq_mod
[1617]	623	do ig=1,ngrid
[1720]	624	pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)/microtimestep
[1617]	625	end do
	626
	627	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
	628	& microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
	629	& rsedcloudco2,rhocloudco2t,
	630	& tempo_traceurs(:,:,igcm_ccnco2_mass),wq,beta)
	631
	632	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
	633	& microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
	634	& rsedcloudco2,rhocloudco2t,
	635	& tempo_traceurs(:,:,igcm_ccnco2_number),wq,beta)
[1720]	636
	637	DO l = 1, nlay !Compute tendencies
[1617]	638	DO ig=1,ngrid
	639	pdqsed(ig,l,igcm_ccnco2_mass)=
	640	& (tempo_traceurs(ig,l,igcm_ccnco2_mass)-
	641	& sav_trac(ig,l,igcm_ccnco2_mass))/microtimestep
	642	pdqsed(ig,l,igcm_ccnco2_number)=
	643	& (tempo_traceurs(ig,l,igcm_ccnco2_number)-
	644	& sav_trac(ig,l,igcm_ccnco2_number))/microtimestep
	645	pdqsed(ig,l,igcm_co2_ice)=
	646	& (tempo_traceurs(ig,l,igcm_co2_ice)-
	647	& sav_trac(ig,l,igcm_co2_ice))/microtimestep
	648	ENDDO
	649	ENDDO
	650	!pdqsed est la tendance due a la sedimentation
	651
	652	DO l=1,nlay
	653	DO ig=1,ngrid
	654	subpdq(ig,l,igcm_ccnco2_mass) =
	655	& subpdq(ig,l,igcm_ccnco2_mass)
	656	& +pdqsed(ig,l,igcm_ccnco2_mass)
	657	subpdq(ig,l,igcm_ccnco2_number) =
	658	& subpdq(ig,l,igcm_ccnco2_number)
	659	& +pdqsed(ig,l,igcm_ccnco2_number)
	660	subpdq(ig,l,igcm_co2_ice) =
	661	& subpdq(ig,l,igcm_co2_ice)
	662	& +pdqsed(ig,l,igcm_co2_ice)
	663	ENDDO
	664	ENDDO
	665
[1720]	666	ENDDO ! of DO microstep=1,imicro
[1617]	667
	668	c-------------------------------------------------------------------
	669	c 6. Compute final tendencies after time loop:
	670	c------------------------------------------------
	671	c CO2 flux at surface (kg.m-2.s-1)
	672	do ig=1,ngrid
[1720]	673	pdqs_sedco2(ig)=pdqs_sedco2(ig)/real(imicro)
[1617]	674	enddo
	675
	676	c------ Temperature tendency pdtcloud
	677	DO l=1,nlay
	678	DO ig=1,ngrid
	679	pdtcloudco2(ig,l) =
[1685]	680	& subpdt(ig,l)/real(imicro)-pdt(ig,l)
[1617]	681	ENDDO
	682	ENDDO
[1685]	683
[1617]	684	c------ Tracers tendencies pdqcloud
	685	DO l=1,nlay
	686	DO ig=1,ngrid
	687
	688	pdqcloudco2(ig,l,igcm_co2_ice) =
[1720]	689	& subpdq(ig,l,igcm_co2_ice)/real(imicro)
[1617]	690	& - pdq(ig,l,igcm_co2_ice)
	691	pdqcloudco2(ig,l,igcm_co2) =
[1720]	692	& subpdq(ig,l,igcm_co2)/real(imicro)
[1617]	693	& - pdq(ig,l,igcm_co2)
[1685]	694	pdqcloudco2(ig,l,igcm_h2o_ice) =
[1720]	695	& subpdq(ig,l,igcm_h2o_ice)/real(imicro)
[1685]	696	& - pdq(ig,l,igcm_h2o_ice)
[1617]	697	ENDDO
	698	ENDDO
	699
	700	DO l=1,nlay
	701	DO ig=1,ngrid
	702	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
[1685]	703	& subpdq(ig,l,igcm_ccnco2_mass)/real(imicro)
[1617]	704	& - pdq(ig,l,igcm_ccnco2_mass)
	705	pdqcloudco2(ig,l,igcm_ccnco2_number) =
[1685]	706	& subpdq(ig,l,igcm_ccnco2_number)/real(imicro)
[1617]	707	& - pdq(ig,l,igcm_ccnco2_number)
[1685]	708	pdqcloudco2(ig,l,igcm_ccn_mass) =
	709	& subpdq(ig,l,igcm_ccn_mass)/real(imicro)
	710	& - pdq(ig,l,igcm_ccn_mass)
	711	pdqcloudco2(ig,l,igcm_ccn_number) =
	712	& subpdq(ig,l,igcm_ccn_number)/real(imicro)
	713	& - pdq(ig,l,igcm_ccn_number)
[1617]	714	ENDDO
	715	ENDDO
	716
	717
	718	DO l=1,nlay
	719	DO ig=1,ngrid
	720	pdqcloudco2(ig,l,igcm_dust_mass) =
[1720]	721	& subpdq(ig,l,igcm_dust_mass)/real(imicro)
[1617]	722	& - pdq(ig,l,igcm_dust_mass)
	723	pdqcloudco2(ig,l,igcm_dust_number) =
	724	& subpdq(ig,l,igcm_dust_number)/real(imicro)
	725	& - pdq(ig,l,igcm_dust_number)
	726	ENDDO
	727	ENDDO
	728
	729	c------- Due to stepped entry, other processes tendencies can add up to negative values
	730	c------- Therefore, enforce positive values and conserve mass
	731
	732
[1720]	733
[1617]	734	DO l=1,nlay
[1720]	735	DO ig=1,ngrid
	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.e-20)
	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-30)) THEN
	744
	745	pdqcloudco2(ig,l,igcm_ccnco2_number) =
	746	& - pqeff(ig,l,igcm_ccnco2_number)/ptimestep
	747	& - pdq(ig,l,igcm_ccnco2_number)+1.e-20
	748	pdqcloudco2(ig,l,igcm_dust_number) =
	749	& -pdqcloudco2(ig,l,igcm_ccnco2_number)
	750	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
	751	& - pqeff(ig,l,igcm_ccnco2_mass)/ptimestep
	752	& - pdq(ig,l,igcm_ccnco2_mass)+1.e-30
	753	pdqcloudco2(ig,l,igcm_dust_mass) =
	754	& -pdqcloudco2(ig,l,igcm_ccnco2_mass)
[1617]	755
[1720]	756	ENDIF
	757	ENDDO
[1617]	758	ENDDO
[1720]	759
[1617]	760
	761
[1720]	762	DO l=1,nlay
	763	DO ig=1,ngrid
	764	IF ( (pqeff(ig,l,igcm_dust_number) +
	765	& ptimestep* (pdq(ig,l,igcm_dust_number) +
	766	& pdqcloudco2(ig,l,igcm_dust_number)) .le. 1.e-30)
	767	& .or. (pqeff(ig,l,igcm_dust_mass)+
	768	& ptimestep* (pdq(ig,l,igcm_dust_mass) +
	769	& pdqcloudco2(ig,l,igcm_dust_mass))
	770	& .le. 1.e-30)) then
	771
	772	pdqcloudco2(ig,l,igcm_dust_number) =
	773	& - pqeff(ig,l,igcm_dust_number)/ptimestep
	774	& - pdq(ig,l,igcm_dust_number)+1.e-30
	775	pdqcloudco2(ig,l,igcm_ccnco2_number) =
[1617]	776	& -pdqcloudco2(ig,l,igcm_dust_number)
[1720]	777	pdqcloudco2(ig,l,igcm_dust_mass) =
	778	& - pqeff(ig,l,igcm_dust_mass)/ptimestep
[1685]	779	& - pdq(ig,l,igcm_dust_mass) +1.e-30
[1720]	780	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
	781	& -pdqcloudco2(ig,l,igcm_dust_mass)
[1617]	782
[1720]	783	ENDIF
	784	ENDDO
	785	ENDDO
	786	!pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq
	787	c$$$
	788	c$$$
	789	c$$$ DO l=1,nlay
	790	c$$$ DO ig=1,ngrid
	791	c$$$ IF (pq(ig,l,igcm_co2_ice) + ptimestep*
	792	c$$$ & (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice))
	793	c$$$ & .lt. 1.e-30) THEN
	794	c$$$ pdqcloudco2(ig,l,igcm_co2_ice) =
	795	c$$$ & - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice)
	796	c$$$ pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice)
	797	c$$$ !write(,) "WARNING CO2 ICE in co2cloud.F"
	798	c$$$
	799	c$$$ ENDIF
	800	c$$$ IF (pq(ig,l,igcm_co2) + ptimestep*
	801	c$$$ & (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2))
	802	c$$$ & .lt. 0.5) THEN
	803	c$$$ pdqcloudco2(ig,l,igcm_co2) =
	804	c$$$ & - pdq(ig,l,igcm_co2_ice) !- pdq(ig,l,igcm_co2)
	805	c$$$c pdqcloudco2(ig,l,igcm_co2_ice) = -pdqcloudco2(ig,l,igcm_co2)
	806	c$$$ pdqcloudco2(ig,l,igcm_co2_ice) =
	807	c$$$ & - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice)
	808	c$$$ ! write(,) "WARNING CO2 in co2cloud.F"
	809	c$$$ ENDIF
	810	c$$$ ENDDO
	811	c$$$ ENDDO
	812	c$$$
[1617]	813
[1720]	814	DO l=1, nlay
	815	DO ig=1,ngrid
[1617]	816
	817	c call updaterice_microco2(
	818	c & pq(ig,l,igcm_co2_ice) + ! ice mass
	819	c & (pdq(ig,l,igcm_co2_ice) + ! ice mass
	820	c & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep, ! ice mass
	821	c & pq(ig,l,igcm_ccnco2_mass) + ! ccn mass
	822	c & (pdq(ig,l,igcm_ccnco2_mass) + ! ccn mass
	823	c & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep, ! ccn mass
	824	c & pq(ig,l,igcm_ccnco2_number) + ! ccn number
	825	c & (pdq(ig,l,igcm_ccnco2_number) + ! ccn number
	826	c & pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep, ! ccn number
	827	c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
	828	c write(,) "in co2clouds, riceco2(ig,l)= ",riceco2(ig,l)
[1720]	829	Niceco2=pqeff(ig,l,igcm_co2_ice) +
	830	& (pdq(ig,l,igcm_co2_ice) +
	831	& pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep
	832	Nco2=pqeff(ig,l,igcm_co2) +
	833	& (pdq(ig,l,igcm_co2) +
	834	& pdqcloudco2(ig,l,igcm_co2))*ptimestep
	835	Nccnco2=max((pqeff(ig,l,igcm_ccnco2_number) +
	836	& (pdq(ig,l,igcm_ccnco2_number) +
	837	& pdqcloudco2(ig,l,igcm_ccnco2_number))ptimestep)
	838	& tauscaling(ig),1.e-30)
	839	Qccnco2=max((pqeff(ig,l,igcm_ccnco2_mass) +
	840	& (pdq(ig,l,igcm_ccnco2_mass) +
	841	& pdqcloudco2(ig,l,igcm_ccnco2_mass))ptimestep)
	842	& tauscaling(ig),1.e-30)
	843
	844	if (Nccnco2 .lt. 0.1) then
	845	rdust(ig,l)=1.e-10
	846	else
[1617]	847
[1720]	848	rdust(ig,l)= Qccnco2
	849	& *0.75/pi/rho_dust
	850	& / Nccnco2
	851	rdust(ig,l)= rdust(ig,l)**(1./3.)
	852	rdust(ig,l)=max(1.e-10,rdust(ig,l))
	853	c rdust(ig,l)=min(5.e-6,rdust(ig,l))
	854	endif
	855	myT=zteff(ig,l)+(pdt(ig,l)+pdtcloudco2(ig,l))*ptimestep
	856	rho_ice_co2T(ig,l)=1000.(1.72391-2.53e-4
	857	& myT-2.87e-6* myT* myT)
	858	rho_ice_co2=rho_ice_co2T(ig,l)
[1617]	859
[1720]	860	lw = l0 + l1 * myT + l2 myT*2 +
	861	& l3 * myT*3 + l4 myT**4 !J.kg-1
[1617]	862
[1720]	863	riceco2(ig,l)=(Niceco2*3.0/
	864	& (4.0rho_ice_co2pi*Nccnco2)
	865	& +rdust(ig,l)*rdust(ig,l)
	866	& rdust(ig,l) )*(1.0/3.0)
	867	c & .or. (riceco2(ig,l) .le. rdust(ig,l))
	868	if ( (Niceco2 .le. 1.e-25).or. (Nccnco2 .le. 0.1) )THEN !anciennement 200
	869	c riceco2(ig,l)=0.
[1617]	870
[1720]	871	c & .or. (Nccnco2 .le. 1.)
	872	c endif
	873	!Flux chaleur latente <0 quand sublimation
[1617]	874
[1720]	875	pdtcloudco2(ig,l)= pdtcloudco2(ig,l)-Niceco2*lw/cpp/ptimestep
	876	c$$$ !NO CLOUD : RESET TRACERS AND CONSERVE MASS
	877	c if (pq(ig,l,igcm_co2_ice)+(pdq(ig,l,igcm_co2_ice)+
	878	c & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep .le. 0.) then
	879	c pdqcloudco2(ig,l,igcm_co2)=0.
	880	c pdqcloudco2(ig,l,igcm_co2_ice)=0.
	881	c else
	882	pdqcloudco2(ig,l,igcm_co2_ice)=-1.*pqeff(ig,l,igcm_co2_ice)
[1685]	883	& /ptimestep-pdq(ig,l,igcm_co2_ice)
[1720]	884	pdqcloudco2(ig,l,igcm_co2)=-1.*
	885	& pdqcloudco2(ig,l,igcm_co2_ice)
	886	c endif
	887	! Reverse h2o ccn and ices into h2o tracers
[1685]	888	if (memdMMccn(ig,l) .gt. 0) then
[1720]	889	pdqcloudco2(ig,l,igcm_ccn_mass)=memdMMccn(ig,l)/ptimestep
[1685]	890	else
	891	memdMMccn(ig,l)=0.
	892	pdqcloudco2(ig,l,igcm_ccn_mass)=0.
	893	endif
	894	if (memdNNccn(ig,l) .gt. 0) then
[1720]	895	pdqcloudco2(ig,l,igcm_ccn_number)=memdNNccn(ig,l)/ptimestep
[1685]	896	else
	897	memdNNccn(ig,l)=0.
	898	pdqcloudco2(ig,l,igcm_ccn_number)=0.
	899	endif
	900	if (memdMMh2o(ig,l) .gt. 0) then
[1720]	901	pdqcloudco2(ig,l,igcm_h2o_ice)=memdMMh2o(ig,l)/ptimestep
[1685]	902	else
	903	memdMMh2o(ig,l)=0.
	904	pdqcloudco2(ig,l,igcm_h2o_ice)=0.
	905	endif
[1720]	906	c if (pq(ig,l,igcm_ccnco2_mass)+(pdq(ig,l,igcm_ccnco2_mass)+
	907	c & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep
	908	c & .le. 1.e-30) then
	909	c pdqcloudco2(ig,l,igcm_ccnco2_mass)=0.
	910	c pdqcloudco2(ig,l,igcm_ccnco2_number)=0.
	911	c pdqcloudco2(ig,l,igcm_co2)=0.
	912	c pdqcloudco2(ig,l,igcm_co2_ice)=0.
	913	c else
[1685]	914	pdqcloudco2(ig,l,igcm_ccnco2_mass)=
[1720]	915	& -1.*pqeff(ig,l,igcm_ccnco2_mass)
[1685]	916	& /ptimestep-pdq(ig,l,igcm_ccnco2_mass)
[1720]	917	c endif
	918	c if (pq(ig,l,igcm_ccnco2_number)+
	919	c & (pdq(ig,l,igcm_ccnco2_number)+
	920	c & pdqcloudco2(ig,l,igcm_ccnco2_number))
	921	c & *ptimestep.le. 1.e-30)
	922	c & then
	923	c pdqcloudco2(ig,l,igcm_ccnco2_mass)=0.
	924	c pdqcloudco2(ig,l,igcm_ccnco2_number)=0.
	925	c pdqcloudco2(ig,l,igcm_co2)=0.
	926	c pdqcloudco2(ig,l,igcm_co2_ice)=0.
	927	c else
[1685]	928	pdqcloudco2(ig,l,igcm_ccnco2_number)=
[1720]	929	& -1.*pqeff(ig,l,igcm_ccnco2_number)
[1685]	930	& /ptimestep-pdq(ig,l,igcm_ccnco2_number)
[1720]	931	c endif
	932	c if (pq(ig,l,igcm_dust_number)+
	933	c & (pdq(ig,l,igcm_dust_number)+
	934	c & pdqcloudco2(ig,l,igcm_dust_number))
	935	c & *ptimestep.le. 1.e-30)
	936	c & then
	937	c pdqcloudco2(ig,l,igcm_dust_number)=0.
	938	c pdqcloudco2(ig,l,igcm_dust_mass)=0.
	939	c else
[1685]	940	pdqcloudco2(ig,l,igcm_dust_number)=
[1720]	941	& pqeff(ig,l,igcm_ccnco2_number)
[1685]	942	& /ptimestep+pdq(ig,l,igcm_ccnco2_number)
[1720]	943	& -memdNNccn(ig,l)/ptimestep
	944	c endif
	945	c if (pq(ig,l,igcm_dust_mass)+
	946	c & (pdq(ig,l,igcm_dust_mass)+
	947	c & pdqcloudco2(ig,l,igcm_dust_mass))
	948	c & *ptimestep .le. 1.e-30)
	949	c & then
	950	c pdqcloudco2(ig,l,igcm_dust_number)=0.
	951	c pdqcloudco2(ig,l,igcm_dust_mass)=0.
	952	c else
[1685]	953	pdqcloudco2(ig,l,igcm_dust_mass)=
[1720]	954	& pqeff(ig,l,igcm_ccnco2_mass)
[1685]	955	& /ptimestep+pdq(ig,l,igcm_ccnco2_mass)
[1720]	956	& -(memdMMh2o(ig,l)+memdMMccn(ig,l))/ptimestep
	957	c endif
[1685]	958	memdMMccn(ig,l)=0.
	959	memdMMh2o(ig,l)=0.
	960	memdNNccn(ig,l)=0.
	961	riceco2(ig,l)=0.
	962	endif
[1720]	963	c Compute opacities
	964	No=Nccnco2
	965	Rn=-log(riceco2(ig,l))
	966	n_derf = erf( (rb_cldco2(1)+Rn) *dev2)
	967	Qext1bins2(ig,l)=0.
	968	do i = 1, nbinco2_cld !Qext below 50 is negligible
	969	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
	970	n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
	971	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
	972	Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i)
	973	enddo
	974	!l'opacité de la case ig est la somme sur l de Qext1bins2
	975
	976
[1685]	977	!update rice water
	978	call updaterice_micro(
[1720]	979	& pqeff(ig,l,igcm_h2o_ice) + ! ice mass
[1685]	980	& (pdq(ig,l,igcm_h2o_ice) + ! ice mass
	981	& pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep, ! ice mass
[1720]	982	& pqeff(ig,l,igcm_ccn_mass) + ! ccn mass
[1685]	983	& (pdq(ig,l,igcm_ccn_mass) + ! ccn mass
	984	& pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep, ! ccn mass
[1720]	985	& pqeff(ig,l,igcm_ccn_number) + ! ccn number
[1685]	986	& (pdq(ig,l,igcm_ccn_number) + ! ccn number
	987	& pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
	988	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
	989
[1617]	990
[1685]	991	call updaterdust(
[1720]	992	& pqeff(ig,l,igcm_dust_mass) + ! dust mass
[1685]	993	& (pdq(ig,l,igcm_dust_mass) + ! dust mass
	994	& pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, ! dust mass
[1720]	995	& pqeff(ig,l,igcm_dust_number) + ! dust number
[1685]	996	& (pdq(ig,l,igcm_dust_number) + ! dust number
	997	& pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number
	998	& rdust(ig,l))
	999
	1000	ENDDO
	1001	ENDDO
[1720]	1002	tau1mic(:)=0.
	1003	do l=1,nlay
	1004	do ig=1,ngrid
	1005	tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l)
	1006	enddo
	1007	enddo
[1685]	1008
[1720]	1009	c$$$
	1010	c$$$ if (riceco2(725,22) .ge. 1.e-10) then
	1011	c$$$
	1012	c$$$ write(,) 'DIAGJA co2 ',pqeff(725,22,igcm_co2) +
	1013	c$$$ & (pdq(725,22,igcm_co2) +
	1014	c$$$ & pdqcloudco2(725,22,igcm_co2))*ptimestep
	1015	c$$$ write(,) 'DIAGJA co2_ice',pqeff(725,22,igcm_co2_ice) +
	1016	c$$$ & (pdq(725,22,igcm_co2_ice) +
	1017	c$$$ & pdqcloudco2(725,22,igcm_co2_ice))*ptimestep
	1018	c$$$
	1019	c$$$ write(,) 'DIAGJA riceco2',riceco2(725,22)
	1020	c$$$ write(,) 'DIAGJA T',zteff(725,22) +
	1021	c$$$ & (pdt(725,22) + pdtcloudco2(725,22))*ptimestep
	1022	c$$$ write(,) 'DIAG pdtcloud',pdtcloudco2(725,22)*ptimestep
	1023	c$$$ write(,) 'DIAGJA ccnNco2',pqeff(725,22,igcm_ccnco2_number)+
	1024	c$$$ & (pdq(725,22,igcm_ccnco2_number) +
	1025	c$$$ & pdqcloudco2(725,22,igcm_ccnco2_number))*ptimestep
	1026	c$$$
	1027	c$$$ write(,) 'DIAGJA dustN',pqeff(725,22,igcm_dust_number) +
	1028	c$$$ & (pdq(725,22,igcm_dust_number) +
	1029	c$$$ & pdqcloudco2(725,22,igcm_dust_number))*ptimestep
	1030	c$$$ ENDIF
	1031	c$$$
[1685]	1032
[1617]	1033
[1720]	1034
[1617]	1035	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
	1036	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1037	c Then that should not affect the ice particle radius
	1038	do ig=1,ngrid
	1039	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
	1040	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
	1041	& riceco2(ig,2)=riceco2(ig,3)
	1042	riceco2(ig,1)=riceco2(ig,2)
	1043	end if
	1044	end do
	1045
	1046
	1047	DO l=1,nlay
	1048	DO ig=1,ngrid
[1685]	1049	rsedcloud(ig,l)=max(rice(ig,l)*
	1050	& (1.+nuice_sed)(1.+nuice_sed)(1.+nuice_sed),
	1051	& rdust(ig,l))
	1052	! rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
	1053	ENDDO
	1054	ENDDO
	1055
	1056	DO l=1,nlay
	1057	DO ig=1,ngrid
[1617]	1058	rsedcloudco2(ig,l)=max(riceco2(ig,l)*
	1059	& (1.+nuiceco2_sed)(1.+nuiceco2_sed)(1.+nuiceco2_sed),
	1060	& rdust(ig,l))
[1720]	1061	c rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5)
[1617]	1062	ENDDO
	1063	ENDDO
	1064
[1720]	1065	call co2sat(ngridnlay,zteff+(pdt+pdtcloudco2)ptimestep
	1066	& ,pplay,zqsatco2)
	1067	do l=1,nlay
	1068	do ig=1,ngrid
	1069	satuco2(ig,l) = (pqeff(ig,l,igcm_co2) +
	1070	& (pdq(ig,l,igcm_co2) +
	1071	& pdqcloudco2(ig,l,igcm_co2))ptimestep)
	1072	& (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)
	1073	!write(,) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l)
[1617]	1074	enddo
	1075	enddo
[1720]	1076	!Tout ce qui est modifié par la microphysique de CO2 doit être rapporté a cloudfrac
	1077	IF (CLFvaryingCO2) THEN
	1078	DO l=1,nlay
	1079	DO ig=1,ngrid
	1080	pdqcloudco2(ig,l,igcm_ccn_mass)=
	1081	& pdqcloudco2(ig,l,igcm_ccn_mass)*cloudfrac(ig,l)
	1082	pdqcloudco2(ig,l,igcm_ccnco2_mass)=
	1083	& pdqcloudco2(ig,l,igcm_ccnco2_mass)*cloudfrac(ig,l)
	1084	pdqcloudco2(ig,l,igcm_ccn_number)=
	1085	& pdqcloudco2(ig,l,igcm_ccn_number)*cloudfrac(ig,l)
	1086	pdqcloudco2(ig,l,igcm_ccnco2_number)=
	1087	& pdqcloudco2(ig,l,igcm_ccnco2_number)*cloudfrac(ig,l)
	1088	pdqcloudco2(ig,l,igcm_dust_mass)=
	1089	& pdqcloudco2(ig,l,igcm_dust_mass)*cloudfrac(ig,l)
	1090	pdqcloudco2(ig,l,igcm_dust_number)=
	1091	& pdqcloudco2(ig,l,igcm_dust_number)*cloudfrac(ig,l)
	1092	pdqcloudco2(ig,l,igcm_h2o_ice)=
	1093	& pdqcloudco2(ig,l,igcm_h2o_ice)*cloudfrac(ig,l)
	1094	pdqcloudco2(ig,l,igcm_co2_ice)=
	1095	& pdqcloudco2(ig,l,igcm_co2_ice)*cloudfrac(ig,l)
	1096	pdqcloudco2(ig,l,igcm_co2)=
	1097	& pdqcloudco2(ig,l,igcm_co2)*cloudfrac(ig,l)
	1098	pdtcloudco2(ig,l)=pdtcloudco2(ig,l)*cloudfrac(ig,l)
	1099	ENDDO
	1100	ENDDO
	1101	ENDIF
	1102
	1103	call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
[1629]	1104	& satuco2)
[1720]	1105	call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
[1649]	1106	& ,3,riceco2)
[1720]	1107	! or output in diagfi.nc (for testphys1d)
[1617]	1108	c call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
	1109	c call WRITEDIAGFI(ngrid,'temp','Temperature ',
	1110	c & 'K JA',1,pt)
	1111
[1649]	1112	call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2","m",3,
[1617]	1113	& rsedcloudco2)
	1114
[1720]	1115	call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"," ",2,
	1116	& tau1mic)
	1117	call WRITEdiagfi(ngrid,"cloudfrac","co2 cloud fraction"," ",3,
	1118	& cloudfrac)
[1617]	1119	! used for rad. transfer calculations
	1120	! nuice is constant because a lognormal distribution is prescribed
	1121	c nuice(1:ngrid,1:nlay)=nuice_ref
	1122
	1123
	1124
	1125	c=======================================================================
	1126
	1127	END
	1128

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

Download in other formats:

Original Format