Context Navigation

co2cloud.F @ 1919

Last change on this file since 1919 was 1918, checked in by emillour, 7 years ago

Mars GCM:
Code cleanup:

remove "comorbit.h" since it is no longer used.
turn "datafile.h" into module datafile_mod.F90 (and rename variable "datafile" as "datadir" since it stores the path to the datafile directory).

Property svn:executable set to *

File size: 39.0 KB

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

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

Download in other formats:

Original Format