Context Navigation

co2cloud.F @ 1713

Last change on this file since 1713 was 1685, checked in by jaudouard, 8 years ago
Further modifications on CO2 clouds scheme. Water ice clouds can now serve as CCN for CO2 clouds
Property svn:executable set to ``*
File size: 36.3 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,
	6	& rsedcloud,rhocloud,zlev,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,
[1685]	18	& rho_ice_co2,r3n_q,rho_ice,nuice_sed,
[1655]	19	& meteo_flux_mass,meteo_flux_number,
	20	& meteo_alt
[1651]	21
[1617]	22	IMPLICIT NONE
	23
	24
	25	c=======================================================================
	26	c CO2 clouds formation
	27	c
	28	c There is a time loop specific to cloud formation
	29	c due to timescales smaller than the GCM integration timestep.
	30	c microphysics subroutine is improvedCO2clouds.F
	31	c
	32	c The co2 clouds tracers (co2_ice, ccn mass and concentration) are
	33	c sedimented at each microtimestep. pdqs_sedco2 keeps track of the
	34	c CO2 flux at the surface
	35	c
	36	c Authors: 09/2016 Joachim Audouard & Constantino Listowski
	37	c Adaptation of the water ice clouds scheme (with specific microphysics)
	38	c of Montmessin, Navarro & al.
	39	c
	40	c
	41	c=======================================================================
	42
	43	c-----------------------------------------------------------------------
	44	c declarations:
	45	c -------------
	46
	47	!#include "dimensions.h"
	48	!#include "dimphys.h"
	49	#include "callkeys.h"
	50	!#include "tracer.h"
	51	!#include "comgeomfi.h"
	52	!#include "dimradmars.h"
	53	! naerkind is set in scatterers.h (built when compiling with makegcm -s #)
	54	!#include"scatterers.h"
	55	#include "microphys.h"
	56
	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.
	70	real,intent(in) :: zlev(ngrid,nlay+1) ! altitude at the boundaries of the layers
	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
	103	REAL Mo,No
	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
	125	INTEGER iq,ig,l
	126	LOGICAL,SAVE :: firstcall=.true.
	127	DOUBLE PRECISION Nccnco2, Niceco2,mdustJA,ndustJA
	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
	134	integer meteo_lvl
[1617]	135	real tempo_traceur_t(ngrid,nlay)
	136	real tempo_traceurs(ngrid,nlay,nq)
	137	real sav_trac(ngrid,nlay,nq)
	138	real pdqsed(ngrid,nlay,nq)
[1685]	139
	140	DOUBLE PRECISION,allocatable,save :: memdMMccn(:,:)
	141	DOUBLE PRECISION,allocatable,save :: memdMMh2o(:,:)
	142	DOUBLE PRECISION,allocatable,save :: memdNNccn(:,:)
[1617]	143	c ** un petit test de coherence
	144	c --------------------------
	145
	146	IF (firstcall) THEN
	147
	148	if (nq.gt.nqmx) then
	149	write(,) 'stop in co2cloud (nq.gt.nqmx)!'
	150	write(,) 'nq=',nq,' nqmx=',nqmx
	151	stop
	152	endif
[1685]	153	write(,) "co2cloud.F: rho_ice_co2 = ",rho_ice_co2
	154
[1617]	155	write(,) "co2cloud: igcm_co2=",igcm_co2
	156	write(,) " igcm_co2_ice=",igcm_co2_ice
	157
	158	write(,) "time subsampling for microphysic ?"
	159	#ifdef MESOSCALE
	160	imicro = 2
	161	#else
	162	imicro = 30
	163	#endif
	164	call getin("imicro",imicro)
	165	write(,)"imicro = ",imicro
	166
	167	microtimestep = ptimestep/real(imicro)
	168	write(,)"Physical timestep is",ptimestep
	169	write(,)"CO2 Microphysics timestep is",microtimestep
[1651]	170
[1685]	171
	172	if (meteo_flux_number .ne. 0) then
	173
[1651]	174	write(,) "Warning ! Meteoritic flux of dust is turned on"
	175	write(,) "Dust mass = ",meteo_flux_mass
	176	write(,) "Dust number = ",meteo_flux_number
	177	write(,) "Are added at the z-level (km)",meteo_alt
	178	write(,) "Every timestep in co2cloud.F"
[1685]	179	endif
	180	if (.not. allocated(memdMMccn)) allocate(memdMMccn(ngrid,nlay))
	181	if (.not. allocated(memdNNccn)) allocate(memdNNccn(ngrid,nlay))
	182	if (.not. allocated(memdMMh2o)) allocate(memdMMh2o(ngrid,nlay))
	183
	184	memdMMccn(:,:)=0.
	185	memdMMh2o(:,:)=0.
	186	memdNNccn(:,:)=0.
[1617]	187	firstcall=.false.
	188	ENDIF ! of IF (firstcall)
	189
	190	c-----Initialization
[1651]	191
[1685]	192	beta=0.85
	193	subpdq(1:ngrid,1:nlay,1:nq) = 0
	194	subpdt(1:ngrid,1:nlay) = 0
	195	subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0
	196	subpdtcloudco2(1:ngrid,1:nlay) = 0
	197
[1651]	198	c add meteoritic flux of dust
	199	!convert meteo_alt (in km) to z-level
	200	!pzlay altitudes of the layers
[1685]	201	if (meteo_flux_number .ne. 0) then
[1651]	202	do ig=1, ngrid
	203	diff0=1000
	204	meteo_lvl=1
	205	do l=1,nlay
	206	diff=pzlay(ig,l)/1000-meteo_alt
	207	if (abs(diff) .lt. diff0) then
	208	diff0=abs(diff)
	209	meteo_lvl=l
	210	endif
	211	enddo
	212	c write(,) "meteo_lvl=",meteo_lvl
[1685]	213	pdq(ig,meteo_lvl,igcm_dust_mass)=
	214	& pdq(ig,meteo_lvl,igcm_dust_mass)
	215	& +dble(meteo_flux_mass)/tauscaling(ig)
	216	pdq(ig,meteo_lvl,igcm_dust_number)=
	217	& pdq(ig,meteo_lvl,igcm_dust_number)
	218	& +dble(meteo_flux_number)/tauscaling(ig)
[1651]	219	enddo
[1685]	220	endif
	221	c call WRITEDIAGFI(ngrid,"pzlay","altitude","km",3,
	222	c & pzlay)
[1617]	223
	224
	225	wq(:,:)=0
	226	! default value if no ice
	227	rhocloudco2(1:ngrid,1:nlay) = rho_dust
	228	rhocloudco2t(1:ngrid,1:nlay) = rho_dust
	229	epaisseur(1:ngrid,1:nlay)=0
	230	masse(1:ngrid,1:nlay)=0
	231
	232	tempo_traceur_t(1:ngrid,1:nlay)=0
	233	tempo_traceurs(1:ngrid,1:nlay,1:nq)=0
	234	sav_trac(1:ngrid,1:nlay,1:nq)=0
	235	pdqsed(1:ngrid,1:nlay,1:nq)=0
	236
	237	do l=1,nlay
	238	do ig=1, ngrid
	239	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
	240	epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
	241	enddo
	242	enddo
	243
	244
	245
	246	c-------------------------------------------------------------------
	247	c 1. Tendencies:
	248	c------------------
	249
	250
	251
	252	c------------------------------------------------------------------
	253	c Time subsampling for microphysics
	254	c------------------------------------------------------------------
	255	DO microstep=1,imicro
	256	c------ Temperature tendency subpdt
	257	! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
	258	! If imicro=1 subpdt is the same as pdt
	259	DO l=1,nlay
	260	DO ig=1,ngrid
	261	c tempo_traceur_t(ig,l)=tempo_traceur_t(ig,l)
	262	c & + subpdtcloudco2(ig,l)
	263	!write(,) 'T micro= ', tempo_traceur_t(ig,l)
	264	c tempo_traceurs(ig,l,:)=tempo_traceurs(ig,l,:)
	265	c & +subpdqcloudco2(ig,l,:)
	266
	267	subpdt(ig,l) = subpdt(ig,l)
	268	& + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
	269
	270	subpdq(ig,l,igcm_dust_mass) =
	271	& subpdq(ig,l,igcm_dust_mass)
	272	& + pdq(ig,l,igcm_dust_mass)
	273
	274	subpdq(ig,l,igcm_dust_number) =
	275	& subpdq(ig,l,igcm_dust_number)
	276	& + pdq(ig,l,igcm_dust_number)
	277
	278	subpdq(ig,l,igcm_ccnco2_mass) =
	279	& subpdq(ig,l,igcm_ccnco2_mass)
	280	& + pdq(ig,l,igcm_ccnco2_mass)
	281
	282	subpdq(ig,l,igcm_ccnco2_number) =
	283	& subpdq(ig,l,igcm_ccnco2_number)
	284	& + pdq(ig,l,igcm_ccnco2_number)
	285
	286	subpdq(ig,l,igcm_co2_ice) =
	287	& subpdq(ig,l,igcm_co2_ice)
	288	& + pdq(ig,l,igcm_co2_ice)
[1685]	289
[1617]	290	subpdq(ig,l,igcm_co2) =
	291	& subpdq(ig,l,igcm_co2)
	292	& + pdq(ig,l,igcm_co2)
[1685]	293
	294	subpdq(ig,l,igcm_h2o_ice) =
	295	& subpdq(ig,l,igcm_h2o_ice)
	296	& + pdq(ig,l,igcm_h2o_ice)
[1617]	297
[1685]	298	subpdq(ig,l,igcm_ccn_mass) =
	299	& subpdq(ig,l,igcm_ccn_mass)
	300	& + pdq(ig,l,igcm_ccn_mass)
	301
	302	subpdq(ig,l,igcm_ccn_number) =
	303	& subpdq(ig,l,igcm_ccn_number)
	304	& + pdq(ig,l,igcm_ccn_number)
	305
[1617]	306	tempo_traceur_t(ig,l)= pt(ig,l)+subpdt(ig,l)*microtimestep
	307	tempo_traceurs(ig,l,:)= pq(ig,l,:)+subpdq(ig,l,:)
	308	& *microtimestep
	309	!Stepped entry for sedimentation
	310	ENDDO
	311	ENDDO
	312
	313	!RSEDCLOUD AND RICECO2 HERE
	314
	315	DO l=1, nlay
	316	DO ig=1,ngrid
[1685]	317
	318	rho_ice_co2T(ig,l)=1000.(1.72391-2.53e-4
	319	& tempo_traceur_t(ig,l)-2.87e-6*
	320	& tempo_traceur_t(ig,l)*tempo_traceur_t(ig,l))
	321
	322	rho_ice_co2=rho_ice_co2T(ig,l)
[1649]	323	Niceco2=max(tempo_traceurs(ig,l,igcm_co2_ice),1.e-30)
[1617]	324	Nccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_number),
	325	& 1.e-30)
	326	Qccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_mass),
	327	& 1.e-30)
	328	mdustJA= tempo_traceurs(ig,l,igcm_dust_mass)
	329	ndustJA=tempo_traceurs(ig,l,igcm_dust_number)
[1685]	330	if ((Nccnco2 .lt. tauscaling(ig)) .or. (Qccnco2 .lt.
[1617]	331	& 1.e-30 *tauscaling(ig))) then
	332	rdust(ig,l)=1.e-10
	333	else
[1685]	334	rdust(ig,l)=(3./4./pi/2500.Qccnco2/Nccnco2)*(1./3.)
	335	rdust(ig,l)=max(rdust(ig,l),1.e-10)
	336	! rdust(ig,l)=min(rdust(ig,l),5.e-6)
[1617]	337	end if
[1685]	338	rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2
	339	& + Qccnco2tauscaling(ig)rho_dust)
	340	& / (Niceco2 + Qccnco2*tauscaling(ig))
	341	c riceco2(ig,l)= Niceco2*3.0/
	342	c & (4.0rho_ice_co2pi*Nccnco2)
	343	c & +rdust(ig,l)rdust(ig,l)rdust(ig,l)
	344	c riceco2(ig,l)=riceco2(ig,l)**(1.0/3.0)
[1649]	345	c write(,) "in co2clouds, rice = ",riceco2(ig,l)
	346	c write(,) "in co2clouds, rho = ",rhocloudco2t(ig,l)
[1617]	347
[1685]	348	riceco2(ig,l)=(Niceco2*3.0/
	349	& (4.0rho_ice_co2pi*Nccnco2
	350	& tauscaling(ig)) +rdust(ig,l)rdust(ig,l)
	351	& rdust(ig,l))*(1.0/3.0)
	352
	353	riceco2(ig,l)=max(1.e-10,riceco2(ig,l))
	354	!riceco2(ig,l)=min(5.e-5,riceco2(ig,l))
	355
	356	c call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
	357	c & tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
[1649]	358	c write(,) "in co2clouds, rice update = ",riceco2(ig,l)
	359	c write(,) "in co2clouds, rho update = "
	360	c & ,rhocloudco2t(ig,l)
[1617]	361
	362	rsedcloudco2(ig,l)=max(riceco2(ig,l)*
	363	& (1.+nuiceco2_sed)(1.+nuiceco2_sed)(1.+nuiceco2_sed),
	364	& rdust(ig,l))
[1685]	365	rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5)
[1649]	366	c write(,) 'Rsedcloud = ',rsedcloudco2(ig,l)
[1617]	367	!write(,) 'Rhocloudco2 = ',rhocloudco2t(ig,l)
	368
	369	ENDDO
	370	ENDDO
	371
	372	! Gravitational sedimentation
	373
	374	! sedimentation computed from radius computed from q in module radii_mod
	375	sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice)
	376	sav_trac(:,:,igcm_ccnco2_mass)=
	377	& tempo_traceurs(:,:,igcm_ccnco2_mass)
	378	sav_trac(:,:,igcm_ccnco2_number)=
	379	& tempo_traceurs(:,:,igcm_ccnco2_number)
	380
	381	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
	382	& microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
	383	& rsedcloudco2,rhocloudco2t,
	384	& tempo_traceurs(:,:,igcm_co2_ice),wq,beta) ! 3 traceurs
	385
	386	! sedim at the surface of co2 ice
	387	do ig=1,ngrid
	388	pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)
	389	end do
	390
	391	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
	392	& microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
	393	& rsedcloudco2,rhocloudco2t,
	394	& tempo_traceurs(:,:,igcm_ccnco2_mass),wq,beta)
	395
	396	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
	397	& microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
	398	& rsedcloudco2,rhocloudco2t,
	399	& tempo_traceurs(:,:,igcm_ccnco2_number),wq,beta)
	400
	401
	402	DO l = 1, nlay
	403	DO ig=1,ngrid
	404	pdqsed(ig,l,igcm_ccnco2_mass)=
	405	& (tempo_traceurs(ig,l,igcm_ccnco2_mass)-
	406	& sav_trac(ig,l,igcm_ccnco2_mass))/microtimestep
	407	pdqsed(ig,l,igcm_ccnco2_number)=
	408	& (tempo_traceurs(ig,l,igcm_ccnco2_number)-
	409	& sav_trac(ig,l,igcm_ccnco2_number))/microtimestep
	410	pdqsed(ig,l,igcm_co2_ice)=
	411	& (tempo_traceurs(ig,l,igcm_co2_ice)-
	412	& sav_trac(ig,l,igcm_co2_ice))/microtimestep
	413	ENDDO
	414	ENDDO
	415	!pdqsed est la tendance due a la sedimentation
	416
	417	DO l = 1, nlay
	418	DO ig=1,ngrid
	419	pdqsed(ig,l,igcm_ccnco2_mass)=
	420	& (tempo_traceurs(ig,l,igcm_ccnco2_mass)-
	421	& sav_trac(ig,l,igcm_ccnco2_mass))/microtimestep
	422	pdqsed(ig,l,igcm_ccnco2_number)=
	423	& (tempo_traceurs(ig,l,igcm_ccnco2_number)-
	424	& sav_trac(ig,l,igcm_ccnco2_number))/microtimestep
	425	pdqsed(ig,l,igcm_co2_ice)=
	426	& (tempo_traceurs(ig,l,igcm_co2_ice)-
	427	& sav_trac(ig,l,igcm_co2_ice))/microtimestep
	428	ENDDO
	429	ENDDO
	430	!pdqsed est la tendance due a la sedimentation
	431	DO l=1,nlay
	432	DO ig=1,ngrid
	433	subpdq(ig,l,igcm_ccnco2_mass) =
	434	& subpdq(ig,l,igcm_ccnco2_mass)
	435	& +pdqsed(ig,l,igcm_ccnco2_mass)
	436
	437	subpdq(ig,l,igcm_ccnco2_number) =
	438	& subpdq(ig,l,igcm_ccnco2_number)
	439	& +pdqsed(ig,l,igcm_ccnco2_number)
	440
	441	subpdq(ig,l,igcm_co2_ice) =
	442	& subpdq(ig,l,igcm_co2_ice)
	443	& +pdqsed(ig,l,igcm_co2_ice)
	444	ENDDO
	445	ENDDO
	446	c-------------------------------------------------------------------
	447	c 2. Main call to the different cloud schemes:
	448	c------------------------------------------------
	449	IF (microphysco2) THEN
	450	CALL improvedCO2clouds(ngrid,nlay,microtimestep,
	451	& pplay,pt,subpdt,
	452	& pq,subpdq,subpdqcloudco2,subpdtcloudco2,
[1685]	453	& nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn)
[1617]	454
	455	ELSE
	456
	457	write(,) ' no simpleCO2clouds procedure: STOP' ! listo
[1685]	458	write(,) ' microphysco2 and co2clouds must be .true.' ! listo
	459	STOP
[1617]	460
	461	ENDIF
	462
	463
	464	c-------------------------------------------------------------------
	465	c 3. Updating tendencies after cloud scheme:
	466	c-----------------------------------------------
	467
	468	c IF (microphysco2) THEN
	469	DO l=1,nlay
	470	DO ig=1,ngrid
	471	subpdq(ig,l,igcm_dust_mass) =
	472	& subpdq(ig,l,igcm_dust_mass)
	473	& + subpdqcloudco2(ig,l,igcm_dust_mass)
	474
	475	subpdq(ig,l,igcm_dust_number) =
	476	& subpdq(ig,l,igcm_dust_number)
	477	& + subpdqcloudco2(ig,l,igcm_dust_number)
	478
	479	subpdq(ig,l,igcm_ccnco2_mass) =
	480	& subpdq(ig,l,igcm_ccnco2_mass)
	481	& + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
	482	c & +pdqsed(ig,l,igcm_ccnco2_mass)
	483
	484	subpdq(ig,l,igcm_ccnco2_number) =
	485	& subpdq(ig,l,igcm_ccnco2_number)
	486	& + subpdqcloudco2(ig,l,igcm_ccnco2_number)
	487	c & +pdqsed(ig,l,igcm_ccnco2_number)
	488
	489	subpdq(ig,l,igcm_co2_ice) =
	490	& subpdq(ig,l,igcm_co2_ice)
	491	& + subpdqcloudco2(ig,l,igcm_co2_ice)
	492	c & +pdqsed(ig,l,igcm_co2_ice)
	493
	494	subpdq(ig,l,igcm_co2) =
	495	& subpdq(ig,l,igcm_co2)
	496	& + subpdqcloudco2(ig,l,igcm_co2)
[1685]	497
	498	subpdq(ig,l,igcm_h2o_ice) =
	499	& subpdq(ig,l,igcm_h2o_ice)
	500	& + subpdqcloudco2(ig,l,igcm_h2o_ice)
	501
	502	subpdq(ig,l,igcm_ccn_mass) =
	503	& subpdq(ig,l,igcm_ccn_mass)
	504	& + subpdqcloudco2(ig,l,igcm_ccn_mass)
	505
	506	subpdq(ig,l,igcm_ccn_number) =
	507	& subpdq(ig,l,igcm_ccn_number)
	508	& + subpdqcloudco2(ig,l,igcm_ccn_number)
[1617]	509	ENDDO
	510	ENDDO
	511
[1685]	512
[1617]	513
	514
	515	IF (activice) THEN
	516	DO l=1,nlay
	517	DO ig=1,ngrid
	518	subpdt(ig,l) =
	519	& subpdt(ig,l) + subpdtcloudco2(ig,l)
	520	ENDDO
	521	ENDDO
	522	ENDIF
	523
	524
	525	ENDDO ! of DO microstep=1,imicro
	526
	527	c-------------------------------------------------------------------
	528	c 6. Compute final tendencies after time loop:
	529	c------------------------------------------------
	530	c CO2 flux at surface (kg.m-2.s-1)
	531	do ig=1,ngrid
	532	pdqs_sedco2(ig)=pdqs_sedco2(ig)/ptimestep
	533	enddo
	534
	535	c------ Temperature tendency pdtcloud
	536	DO l=1,nlay
	537	DO ig=1,ngrid
	538	pdtcloudco2(ig,l) =
[1685]	539	& subpdt(ig,l)/real(imicro)-pdt(ig,l)
[1617]	540	ENDDO
	541	ENDDO
[1685]	542
[1617]	543	c------ Tracers tendencies pdqcloud
	544	DO l=1,nlay
	545	DO ig=1,ngrid
	546
	547	pdqcloudco2(ig,l,igcm_co2_ice) =
[1685]	548	& subpdq(ig,l,igcm_co2_ice)/real(imicro)
[1617]	549	& - pdq(ig,l,igcm_co2_ice)
[1685]	550
[1617]	551	pdqcloudco2(ig,l,igcm_co2) =
[1685]	552	& subpdq(ig,l,igcm_co2)/real(imicro)
[1617]	553	& - pdq(ig,l,igcm_co2)
[1685]	554
	555	pdqcloudco2(ig,l,igcm_h2o_ice) =
	556	& subpdq(ig,l,igcm_h2o_ice)/real(imicro)
	557	& - pdq(ig,l,igcm_h2o_ice)
[1617]	558	ENDDO
	559	ENDDO
	560
	561
	562	! call WRITEdiagfi(ngrid,"co2cloud00","co2 traceur","kg/kg",1,
	563	! & pq(1,:,igcm_co2_ice) + ptimestep*
	564	! & (pdq(1,:,igcm_co2_ice) + pdqcloudco2(1,:,igcm_co2_ice)))
	565
	566
[1685]	567	c IF(microphysco2) THEN
[1617]	568	DO l=1,nlay
	569	DO ig=1,ngrid
	570	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
[1685]	571	& subpdq(ig,l,igcm_ccnco2_mass)/real(imicro)
[1617]	572	& - pdq(ig,l,igcm_ccnco2_mass)
	573	pdqcloudco2(ig,l,igcm_ccnco2_number) =
[1685]	574	& subpdq(ig,l,igcm_ccnco2_number)/real(imicro)
[1617]	575	& - pdq(ig,l,igcm_ccnco2_number)
[1685]	576	pdqcloudco2(ig,l,igcm_ccn_mass) =
	577	& subpdq(ig,l,igcm_ccn_mass)/real(imicro)
	578	& - pdq(ig,l,igcm_ccn_mass)
	579	pdqcloudco2(ig,l,igcm_ccn_number) =
	580	& subpdq(ig,l,igcm_ccn_number)/real(imicro)
	581	& - pdq(ig,l,igcm_ccn_number)
[1617]	582	ENDDO
	583	ENDDO
[1685]	584	c ENDIF
[1617]	585
	586
	587	IF(scavenging) THEN
	588	DO l=1,nlay
	589	DO ig=1,ngrid
	590	pdqcloudco2(ig,l,igcm_dust_mass) =
	591	& subpdq(ig,l,igcm_dust_mass)/real(imicro)
	592	& - pdq(ig,l,igcm_dust_mass)
	593	pdqcloudco2(ig,l,igcm_dust_number) =
	594	& subpdq(ig,l,igcm_dust_number)/real(imicro)
	595	& - pdq(ig,l,igcm_dust_number)
	596	ENDDO
	597	ENDDO
	598	ENDIF
	599
	600	c ENDIF
	601	c------- Due to stepped entry, other processes tendencies can add up to negative values
	602	c------- Therefore, enforce positive values and conserve mass
	603
	604
	605	IF(microphysco2) THEN
	606	DO l=1,nlay
	607	DO ig=1,ngrid
	608	IF ((pq(ig,l,igcm_ccnco2_number) +
	609	& ptimestep* (pdq(ig,l,igcm_ccnco2_number) +
	610	& pdqcloudco2(ig,l,igcm_ccnco2_number))
[1685]	611	& .lt. 1.e-15)
	612	& .or. (pq(ig,l,igcm_ccnco2_mass) +
[1617]	613	& ptimestep* (pdq(ig,l,igcm_ccnco2_mass) +
	614	& pdqcloudco2(ig,l,igcm_ccnco2_mass))
[1685]	615	& .lt. 1.e-30)) THEN
[1617]	616
	617	pdqcloudco2(ig,l,igcm_ccnco2_number) =
	618	& - pq(ig,l,igcm_ccnco2_number)/ptimestep
[1685]	619	& - pdq(ig,l,igcm_ccnco2_number)+1.e-15
[1617]	620
	621	pdqcloudco2(ig,l,igcm_dust_number) =
	622	& -pdqcloudco2(ig,l,igcm_ccnco2_number)
	623
	624	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
	625	& - pq(ig,l,igcm_ccnco2_mass)/ptimestep
[1685]	626	& - pdq(ig,l,igcm_ccnco2_mass)+1.e-30
[1617]	627
	628	pdqcloudco2(ig,l,igcm_dust_mass) =
	629	& -pdqcloudco2(ig,l,igcm_ccnco2_mass)
	630
	631	ENDIF
	632	ENDDO
	633	ENDDO
	634	ENDIF
	635
	636
	637	IF(scavenging) THEN
	638	DO l=1,nlay
	639	DO ig=1,ngrid
	640	IF ( (pq(ig,l,igcm_dust_number) +
	641	& ptimestep* (pdq(ig,l,igcm_dust_number) +
[1685]	642	& pdqcloudco2(ig,l,igcm_dust_number)) .le. 1.e-15)
[1617]	643	& .or. (pq(ig,l,igcm_dust_mass)+
	644	& ptimestep* (pdq(ig,l,igcm_dust_mass) +
	645	& pdqcloudco2(ig,l,igcm_dust_mass))
[1685]	646	& .le. 1.e-30)) then
[1617]	647
	648	pdqcloudco2(ig,l,igcm_dust_number) =
	649	& - pq(ig,l,igcm_dust_number)/ptimestep
[1685]	650	& - pdq(ig,l,igcm_dust_number)+1.e-15
[1617]	651
	652	pdqcloudco2(ig,l,igcm_ccnco2_number) =
	653	& -pdqcloudco2(ig,l,igcm_dust_number)
	654
	655	pdqcloudco2(ig,l,igcm_dust_mass) =
	656	& - pq(ig,l,igcm_dust_mass)/ptimestep
[1685]	657	& - pdq(ig,l,igcm_dust_mass) +1.e-30
[1617]	658
	659	pdqcloudco2(ig,l,igcm_ccnco2_mass) =
	660	& -pdqcloudco2(ig,l,igcm_dust_mass)
	661	ENDIF
	662	ENDDO
	663	ENDDO
	664	ENDIF !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq
	665
	666
[1685]	667	c DO l=1,nlay
	668	c DO ig=1,ngrid
	669	c IF (pq(ig,l,igcm_co2_ice) + ptimestep*
	670	c & (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice))
	671	c & .lt. 1.e-25) THEN
	672	c pdqcloudco2(ig,l,igcm_co2_ice) =
	673	c & - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice)
	674	c pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice)
	675	c write(,) "WARNING CO2 ICE in co2cloud.F"
	676
	677	c ENDIF
	678	c IF (pq(ig,l,igcm_co2) + ptimestep*
	679	c & (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2))
	680	c & .lt. 1.e-10) THEN
	681	c pdqcloudco2(ig,l,igcm_co2) =
	682	c & - pq(ig,l,igcm_co2)/ptimestep - pdq(ig,l,igcm_co2)
	683	c pdqcloudco2(ig,l,igcm_co2_ice) = -pdqcloudco2(ig,l,igcm_co2)
	684	c write(,) "WARNING CO2 in co2cloud.F"
	685	c ENDIF
	686	c ENDDO
	687	c ENDDO
[1617]	688
	689
	690
	691
	692	c------Update the ice and dust particle size "riceco2" for output or photochemistry
	693	c------Only rsedcloudco2 is used for the co2 (cloud) cycle
	694
[1685]	695	c IF(scavenging) THEN
	696	c DO l=1, nlay
	697	c DO ig=1,ngrid
[1617]	698
	699	c call updaterdust(
	700	c & pq(ig,l,igcm_dust_mass) + ! dust mass
	701	c & (pdq(ig,l,igcm_dust_mass) + ! dust mass
	702	c & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, ! dust mass
	703	c & pq(ig,l,igcm_dust_number) + ! dust number
	704	c & (pdq(ig,l,igcm_dust_number) + ! dust number
	705	c & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number
	706	c & rdust(ig,l))
	707	c write(,) "in co2clouds, rdust(ig,l)= ",rdust(ig,l)
[1685]	708	c mdustJA= pq(ig,l,igcm_dust_mass) +
	709	c & (pdq(ig,l,igcm_dust_mass) +
	710	c & pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep
	711	c ndustJA=pq(ig,l,igcm_dust_number) +
	712	c & (pdq(ig,l,igcm_dust_number) +
	713	c & pdqcloudco2(ig,l,igcm_dust_number))*ptimestep
	714	c if ((ndustJA .lt. tauscaling(ig)) .or. (mdustJA .lt.
	715	c & 1.e-30 *tauscaling(ig))) then
	716	c rdust(ig,l)=1.e-10
	717	c else
	718	c rdust(ig,l)=(3./4./pi/2500.mdustJA/ndustJA)*(1./3.)
	719	c rdust(ig,l)=min(rdust(ig,l),5.e-6)
	720	c rdust(ig,l)=max(rdust(ig,l),1.e-10)
	721	c endif
	722	c ENDDO
	723	c ENDDO
	724	c ENDIF
[1617]	725
	726
	727	IF(microphysco2) THEN
	728
	729	DO l=1, nlay
	730	DO ig=1,ngrid
	731
	732	c call updaterice_microco2(
	733	c & pq(ig,l,igcm_co2_ice) + ! ice mass
	734	c & (pdq(ig,l,igcm_co2_ice) + ! ice mass
	735	c & pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep, ! ice mass
	736	c & pq(ig,l,igcm_ccnco2_mass) + ! ccn mass
	737	c & (pdq(ig,l,igcm_ccnco2_mass) + ! ccn mass
	738	c & pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep, ! ccn mass
	739	c & pq(ig,l,igcm_ccnco2_number) + ! ccn number
	740	c & (pdq(ig,l,igcm_ccnco2_number) + ! ccn number
	741	c & pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep, ! ccn number
	742	c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
	743	c write(,) "in co2clouds, riceco2(ig,l)= ",riceco2(ig,l)
	744
	745
	746	Niceco2=pq(ig,l,igcm_co2_ice) +
	747	& (pdq(ig,l,igcm_co2_ice) +
	748	& pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep
	749	Nccnco2=max((pq(ig,l,igcm_ccnco2_number) +
	750	& (pdq(ig,l,igcm_ccnco2_number) +
	751	& pdqcloudco2(ig,l,igcm_ccnco2_number))ptimestep)
	752	& tauscaling(ig),1.e-30)
	753	Qccnco2=max((pq(ig,l,igcm_ccnco2_mass) +
	754	& (pdq(ig,l,igcm_ccnco2_mass) +
	755	& pdqcloudco2(ig,l,igcm_ccnco2_mass))ptimestep)
	756	& tauscaling(ig),1.e-30)
	757
	758
[1685]	759	if ((Nccnco2 .lt. 1) .or. (Qccnco2 .lt.
	760	& 1.e-30)) then
	761	rdust(ig,l)=1.e-10
	762	else
	763
	764	rdust(ig,l)= Qccnco2
	765	& *0.75/pi/rho_dust
	766	& / Nccnco2
	767	rdust(ig,l)= rdust(ig,l)**(1./3.)
	768	rdust(ig,l)=max(1.e-10,rdust(ig,l))
	769	!rdust(ig,l)=min(5.e-5,rdust(ig,l))
	770	endif
	771
	772	rho_ice_co2T(ig,l)=1000.(1.72391-2.53e-4
	773	& (pt(ig,l)+(subpdt(ig,l)+pdtcloudco2(ig,l))*ptimestep)
	774	& -2.87e-6*
	775	& (pt(ig,l)+(subpdt(ig,l)+pdtcloudco2(ig,l))ptimestep)
	776	& (pt(ig,l)+(subpdt(ig,l)+pdtcloudco2(ig,l))*ptimestep))
	777	rho_ice_co2=rho_ice_co2T(ig,l)
[1617]	778
[1685]	779	riceco2(ig,l)=(Niceco2*3.0/
[1617]	780	& (4.0rho_ice_co2pi*Nccnco2)
[1685]	781	& +rdust(ig,l)*rdust(ig,l)
	782	& rdust(ig,l) )*(1.0/3.0)
[1617]	783
[1685]	784	if ( (Niceco2 .le. 1.e-23) .or.
	785	& (riceco2(ig,l) .le. rdust(ig,l))
	786	& .or. (Nccnco2 .le. 0.01))then
[1617]	787
[1685]	788	c write(,) "Riceco2 co2cloud before reset=",riceco2(ig,l)
	789	c write(,) "Niceco2 co2cloud before reset=",Niceco2
	790	c write(,) "Nccnco2 co2cloud before reset=",Nccnco2
	791	c write(,) "Rdust co2cloud before reset=",rdust(ig,l)
[1617]	792
[1685]	793	c$$$ !NO CLOUD : RESET TRACERS AND CONSERVE MASS
	794	if (pq(ig,l,igcm_co2_ice)+(pdq(ig,l,igcm_co2_ice)+
	795	& pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep .le. 0.) then
	796	pdqcloudco2(ig,l,igcm_co2)=0.
	797	pdqcloudco2(ig,l,igcm_co2_ice)=0.
	798	else
	799	pdqcloudco2(ig,l,igcm_co2)= pq(ig,l,igcm_co2_ice)
	800	& /ptimestep+pdq(ig,l,igcm_co2_ice)
	801	pdqcloudco2(ig,l,igcm_co2_ice)=-pq(ig,l,igcm_co2_ice)
	802	& /ptimestep-pdq(ig,l,igcm_co2_ice)
	803	endif
	804	!Reverse h2o ccn and ices into h2o tracers
	805	if (memdMMccn(ig,l) .gt. 0) then
	806	pdqcloudco2(ig,l,igcm_ccn_mass)=memdMMccn(ig,l)
	807	else
	808	memdMMccn(ig,l)=0.
	809	pdqcloudco2(ig,l,igcm_ccn_mass)=0.
	810	endif
	811	if (memdNNccn(ig,l) .gt. 0) then
	812	pdqcloudco2(ig,l,igcm_ccn_number)=memdNNccn(ig,l)
	813	else
	814	memdNNccn(ig,l)=0.
	815	pdqcloudco2(ig,l,igcm_ccn_number)=0.
	816	endif
	817	if (memdMMh2o(ig,l) .gt. 0) then
	818	pdqcloudco2(ig,l,igcm_h2o_ice)=memdMMh2o(ig,l)
	819	else
	820	memdMMh2o(ig,l)=0.
	821	pdqcloudco2(ig,l,igcm_h2o_ice)=0.
	822	endif
[1617]	823
[1685]	824	if (pq(ig,l,igcm_ccnco2_mass)+(pdq(ig,l,igcm_ccnco2_mass)+
	825	& pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep
	826	& .le. 1.e-30) then
	827	pdqcloudco2(ig,l,igcm_ccnco2_mass)=0.
	828	pdqcloudco2(ig,l,igcm_ccnco2_number)=0.
	829	pdqcloudco2(ig,l,igcm_co2)=0.
	830	pdqcloudco2(ig,l,igcm_co2_ice)=0.
	831	else
	832	pdqcloudco2(ig,l,igcm_ccnco2_mass)=
	833	& -pq(ig,l,igcm_ccnco2_mass)
	834	& /ptimestep-pdq(ig,l,igcm_ccnco2_mass)
	835	endif
[1617]	836
[1685]	837	if (pq(ig,l,igcm_ccnco2_number)+
	838	& (pdq(ig,l,igcm_ccnco2_number)+
	839	& pdqcloudco2(ig,l,igcm_ccnco2_number))
	840	& *ptimestep.le. 1.e-30)
	841	& then
	842	pdqcloudco2(ig,l,igcm_ccnco2_mass)=0.
	843	pdqcloudco2(ig,l,igcm_ccnco2_number)=0.
	844	pdqcloudco2(ig,l,igcm_co2)=0.
	845	pdqcloudco2(ig,l,igcm_co2_ice)=0.
	846	else
	847	pdqcloudco2(ig,l,igcm_ccnco2_number)=
	848	& -pq(ig,l,igcm_ccnco2_number)
	849	& /ptimestep-pdq(ig,l,igcm_ccnco2_number)
	850	endif
[1617]	851
[1685]	852	if (pq(ig,l,igcm_dust_number)+
	853	& (pdq(ig,l,igcm_dust_number)+
	854	& pdqcloudco2(ig,l,igcm_dust_number))
	855	& *ptimestep.le. 1.e-30)
	856	& then
	857	pdqcloudco2(ig,l,igcm_dust_number)=0.
	858	pdqcloudco2(ig,l,igcm_dust_mass)=0.
	859	else
	860	pdqcloudco2(ig,l,igcm_dust_number)=
	861	& pq(ig,l,igcm_ccnco2_number)
	862	& /ptimestep+pdq(ig,l,igcm_ccnco2_number)
	863	& -memdNNccn(ig,l)
	864	endif
	865
	866	if (pq(ig,l,igcm_dust_mass)+
	867	& (pdq(ig,l,igcm_dust_mass)+
	868	& pdqcloudco2(ig,l,igcm_dust_mass))
	869	& *ptimestep .le. 1.e-30)
	870	& then
	871	pdqcloudco2(ig,l,igcm_dust_number)=0.
	872	pdqcloudco2(ig,l,igcm_dust_mass)=0.
	873	else
	874	pdqcloudco2(ig,l,igcm_dust_mass)=
	875	& pq(ig,l,igcm_ccnco2_mass)
	876	& /ptimestep+pdq(ig,l,igcm_ccnco2_mass)
	877	& -(memdMMh2o(ig,l)+memdMMccn(ig,l))
	878	endif
	879	memdMMccn(ig,l)=0.
	880	memdMMh2o(ig,l)=0.
	881	memdNNccn(ig,l)=0.
	882	riceco2(ig,l)=0.
	883	pdtcloudco2(ig,l)=0.
	884	endif
[1617]	885
[1685]	886
[1617]	887
[1685]	888	!update rice water
	889	call updaterice_micro(
	890	& pq(ig,l,igcm_h2o_ice) + ! ice mass
	891	& (pdq(ig,l,igcm_h2o_ice) + ! ice mass
	892	& pdqcloudco2(ig,l,igcm_h2o_ice))*ptimestep, ! ice mass
	893	& pq(ig,l,igcm_ccn_mass) + ! ccn mass
	894	& (pdq(ig,l,igcm_ccn_mass) + ! ccn mass
	895	& pdqcloudco2(ig,l,igcm_ccn_mass))*ptimestep, ! ccn mass
	896	& pq(ig,l,igcm_ccn_number) + ! ccn number
	897	& (pdq(ig,l,igcm_ccn_number) + ! ccn number
	898	& pdqcloudco2(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
	899	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
	900
[1617]	901
[1685]	902	call updaterdust(
	903	& pq(ig,l,igcm_dust_mass) + ! dust mass
	904	& (pdq(ig,l,igcm_dust_mass) + ! dust mass
	905	& pdqcloudco2(ig,l,igcm_dust_mass))*ptimestep, ! dust mass
	906	& pq(ig,l,igcm_dust_number) + ! dust number
	907	& (pdq(ig,l,igcm_dust_number) + ! dust number
	908	& pdqcloudco2(ig,l,igcm_dust_number))*ptimestep, ! dust number
	909	& rdust(ig,l))
	910
	911	ENDDO
	912	ENDDO
	913
	914
	915
[1617]	916	ELSE ! no microphys ! not of concern for co2 clouds - listo
	917
	918	ENDIF ! of IF(microphysco2)
	919
	920
	921	c TO CHEK for relevancy - listo
	922
	923	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
	924	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	925	c Then that should not affect the ice particle radius
[1685]	926	!do ig=1,ngrid
	927	! if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
	928	! if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
	929	&! riceco2(ig,2)=riceco2(ig,3)
	930	! riceco2(ig,1)=riceco2(ig,2)
	931	! end if
	932	!end do
	933
	934	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
	935	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	936	c Then that should not affect the ice particle radius
[1617]	937	do ig=1,ngrid
	938	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
	939	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
	940	& riceco2(ig,2)=riceco2(ig,3)
	941	riceco2(ig,1)=riceco2(ig,2)
	942	end if
	943	end do
	944
	945
	946	DO l=1,nlay
	947	DO ig=1,ngrid
[1685]	948	rsedcloud(ig,l)=max(rice(ig,l)*
	949	& (1.+nuice_sed)(1.+nuice_sed)(1.+nuice_sed),
	950	& rdust(ig,l))
	951	! rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
	952	ENDDO
	953	ENDDO
	954
	955	DO l=1,nlay
	956	DO ig=1,ngrid
[1617]	957	rsedcloudco2(ig,l)=max(riceco2(ig,l)*
	958	& (1.+nuiceco2_sed)(1.+nuiceco2_sed)(1.+nuiceco2_sed),
	959	& rdust(ig,l))
[1685]	960	rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5)
[1617]	961	ENDDO
	962	ENDDO
	963
	964	call co2sat(ngrid*nlay,pt,pplay,zqsatco2)
	965	do ig=1,ngrid
	966	do l=1,nlay
[1685]	967	satuco2(ig,l) = (pq(ig,l,igcm_co2) +
	968	& (pdq(ig,l,igcm_co2) +
	969	& pdqcloudco2(ig,l,igcm_co2))ptimestep)
[1617]	970	& (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)
	971
[1649]	972	c write(,) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l)
[1617]	973	enddo
	974	enddo
[1649]	975	call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
[1629]	976	& satuco2)
	977	call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
[1649]	978	& ,3,riceco2)
[1617]	979	! or output in diagfi.nc (for testphys1d)
	980	c call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
	981	c call WRITEDIAGFI(ngrid,'temp','Temperature ',
	982	c & 'K JA',1,pt)
	983
[1649]	984	call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2","m",3,
[1617]	985	& rsedcloudco2)
	986
	987	! used for rad. transfer calculations
	988	! nuice is constant because a lognormal distribution is prescribed
	989	c nuice(1:ngrid,1:nlay)=nuice_ref
	990
	991
	992
	993	c=======================================================================
	994
	995	END
	996

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

Download in other formats:

Original Format