Context Navigation

co2cloud.F @ 2156

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

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

Download in other formats:

Original Format