Context Navigation

co2cloud.F @ 1944

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

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

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

Download in other formats:

Original Format