Context Navigation

co2cloud.F @ 2349

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

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

Download in other formats:

Original Format