Context Navigation

co2cloud.F @ 1816

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

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

Download in other formats:

Original Format