Context Navigation

co2cloud.F @ 1880

Last change on this file since 1880 was 1820, checked in by emillour, 7 years ago

Mars GCM:

Bug fix in co2cloud.F (in the computation of the saturation index) along with some code tidying.

Property svn:executable set to *

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

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

Download in other formats:

Original Format