Context Navigation

co2cloud.F @ 1921

Last change on this file since 1921 was 1921, checked in by emillour, 8 years ago

Mars GCM:
CO2 code update:

nucleaCO2.F: code cleanup and use co2useh2o flag to handle cases where

h2o ccns have to be tracked and accounted for.

co2cloud.F & improvedco2clouds.F: code cleanup and use co2useh2o flag to

control updates on water variables if necessary.

physiq.F : cleanup on outputs & compute mtotco2 and icetotco2

Property svn:executable set to *

File size: 39.9 KB

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

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

Download in other formats:

Original Format