Context Navigation

improvedCO2clouds.F @ 1908

Last change on this file since 1908 was 1884, checked in by emillour, 8 years ago

Mars GCM:

updated massflowrateCO2 routine which now uses an analytical formula rather than an iterative solver
some code tidying in improvedCO2clouds.F

CL+EM

Property svn:executable set to *

File size: 28.5 KB

Line
1	subroutine improvedCO2clouds(ngrid,nlay,ptimestep,
2	& pplay,pplev,pt,pdt,
3	& pq,pdq,pdqcloudco2,pdtcloudco2,
4	& nq,tauscaling,
5	& memdMMccn,memdMMh2o,memdNNccn)
6	USE comcstfi_h, only: pi, g, cpp
7	USE updaterad, only: updaterice_micro, updaterice_microco2
8	use tracer_mod, only: igcm_dust_mass, igcm_dust_number, rho_dust,
9	& igcm_h2o_ice, igcm_ccn_mass,
10	& igcm_ccn_number, nuice_sed,
11	& igcm_co2, igcm_co2_ice, igcm_ccnco2_mass,
12	& igcm_ccnco2_number, nuiceco2_sed,
13	& rho_ice_co2
14	use conc_mod, only: mmean
15
16	implicit none
17
18	c------------------------------------------------------------------
19	c This routine is used to form CO2 clouds when a parcel of the GCM is
20	c saturated. It includes the ability to have supersaturation, a
21	c computation of the nucleation rates, growthrates and the
22	c scavenging of dust particles by clouds.
23	c It is worth noting that the amount of dust is computed using the
24	c dust optical depth computed in aeropacity.F. That's why
25	c the variable called "tauscaling" is used to convert
26	c pq(dust_mass) and pq(dust_number), which are relative
27	c quantities, to absolute and realistic quantities stored in zq.
28	c This has to be done to convert the inputs into absolute
29	c values, but also to convert the outputs back into relative
30	c values which are then used by the sedimentation and advection
31	c schemes.
32	c CO2 ice particles can nucleate on both dust and on water ice particles
33	c When CO2 ice is deposited onto a water ice particles, the particle is
34	c removed from the water tracers.
35	c Memory of the origin of the co2 particles is kept and thus the
36	c water cycle shouldn't be modified by this.
37	cWARNING: no sedimentation of the water ice origin is performed
38	c in the microphysical timestep in co2cloud.F.
39
40	c Authors of the water ice clouds microphysics
41	c J.-B. Madeleine, based on the work by Franck Montmessin
42	c (October 2011)
43	c T. Navarro, debug,correction, new scheme (October-April 2011)
44	c A. Spiga, optimization (February 2012)
45	c Adaptation for CO2 clouds by Joachim Audouard (09/16), based on the work
46	c of Constantino Listowski
47	c There is an energy limit to how much co2 can sublimate/condensate. It is
48	c defined by the difference of the GCM temperature with the co2 condensation
49	c temperature.
50	c Warning:
51	c If meteoritic particles are activated and turn into co2 ice particles,
52	c then they will be reversed in the dust tracers if the cloud sublimates
53
54	c------------------------------------------------------------------
55	include "callkeys.h"
56	include "microphys.h"
57	include "datafile.h"
58	c------------------------------------------------------------------
59	c Arguments:
60
61	INTEGER,INTENT(in) :: ngrid,nlay
62	integer,intent(in) :: nq ! number of tracers
63	REAL,INTENT(in) :: ptimestep ! physics time step (s)
64	REAL,INTENT(in) :: pplay(ngrid,nlay) ! mid-layer pressure (Pa)
65	REAL,INTENT(in) :: pplev(ngrid,nlay+1) ! inter-layer pressure (Pa)
66	REAL,INTENT(in) :: pt(ngrid,nlay) ! temperature at the middle of the
67	! layers (K)
68	REAL,INTENT(in) :: pdt(ngrid,nlay) ! tendency on temperature from
69	! previous physical parametrizations
70	REAL,INTENT(in) :: pq(ngrid,nlay,nq) ! tracers (kg/kg)
71	REAL,INTENT(in) :: pdq(ngrid,nlay,nq) ! tendencies on tracers
72	! before condensation (kg/kg.s-1)
73	REAL,INTENT(in) :: tauscaling(ngrid) ! Convertion factor for qdust and Ndust
74	c Outputs:
75	REAL,INTENT(out) :: pdqcloudco2(ngrid,nlay,nq) ! tendency on tracers
76	! due to CO2 condensation (kg/kg.s-1)
77	! condensation si igcm_co2_ice
78	REAL,INTENT(out) :: pdtcloudco2(ngrid,nlay) ! tendency on temperature due
79	! to latent heat
80
81	c------------------------------------------------------------------
82	c Local variables:
83	LOGICAL,SAVE :: firstcall=.true.
84	REAL*8 derf ! Error function
85	INTEGER ig,l,i,flag_pourri
86
87	real masse (ngrid,nlay) ! Layer mass (kg.m-2)
88	REAL rice(ngrid,nlay) ! Water Ice mass mean radius (m)
89	! used for nucleation of CO2 on ice-coated ccns
90	REAL rccnh2o(ngrid,nlay) ! Water Ice mass mean radius (m)
91	REAL zq(ngrid,nlay,nq) ! local value of tracers
92	REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
93	REAL zt(ngrid,nlay) ! local value of temperature
94	REAL zqsat(ngrid,nlay) ! saturation vapor pressure for CO2
95	real tcond(ngrid,nlay)
96	real zqco2(ngrid,nlay)
97	REAL lw !Latent heat of sublimation (J.kg-1)
98	REAL,save :: l0,l1,l2,l3,l4
99	DOUBLE PRECISION dMice ! mass of condensed ice
100	DOUBLE PRECISION sumcheck
101	DOUBLE PRECISION facteurmax!for energy limit on mass growth
102	DOUBLE PRECISION pco2,psat ! Co2 vapor partial pressure (Pa)
103	DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice
104	DOUBLE PRECISION Mo,No
105	DOUBLE PRECISION Rn, Rm, dev2,dev3, n_derf, m_derf
106	DOUBLE PRECISION memdMMccn(ngrid,nlay)
107	DOUBLE PRECISION memdMMh2o(ngrid,nlay)
108	DOUBLE PRECISION memdNNccn(ngrid,nlay)
109
110	! Radius used by the microphysical scheme (m)
111	DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin
112	DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin
113	DOUBLE PRECISION m_aer2(nbinco2_cld) ! mass mixing ratio of particle/each size bin
114	DOUBLE PRECISION m_aer_h2oice2(nbinco2_cld) ! mass mixing ratio of particle/each size bin
115
116	DOUBLE PRECISION n_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
117	DOUBLE PRECISION m_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
118	DOUBLE PRECISION rad_h2oice(nbinco2_cld)
119
120	c REAL*8 sigco2 ! Co2-ice/air surface tension (N.m)
121	c EXTERNAL sigco2
122
123	DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM,dNNh2o,dMMh2o
124	DOUBLE PRECISION dMh2o_ice,dMh2o_ccn
125
126	DOUBLE PRECISION rate(nbinco2_cld) ! nucleation rate
127	DOUBLE PRECISION rateh2o(nbinco2_cld) ! nucleation rate
128	REAL seq
129	DOUBLE PRECISION rho_ice_co2T(ngrid,nlay)
130	DOUBLE PRECISION riceco2(ngrid,nlay) ! CO2Ice mean radius (m)
131	REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
132
133	REAL rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
134	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
135
136	c REAL res ! Resistance growth
137	DOUBLE PRECISION Ic_rice ! Mass transfer rate CO2 ice crystal
138	DOUBLE PRECISION ratioh2o_ccn
139	DOUBLE PRECISION vo2co2
140	c Parameters of the size discretization
141	c used by the microphysical scheme
142	DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m)
143	DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6 ! Maximum radius (m)
144	DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-10! Minimum bounary radius (m)
145	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4 ! Maximum boundary radius (m)
146	DOUBLE PRECISION vrat_cld ! Volume ratio
147	DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
148	SAVE rb_cldco2
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
152	DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o
153	REAL sigma_iceco2 ! Variance of the co2 ice and CCN distributions
154	SAVE sigma_iceco2
155	REAL sigma_ice ! Variance of the h2o ice and CCN distributions
156	SAVE sigma_ice
157	DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2
158	integer,save :: coeffh2o !will be =0 is co2useh2o=.false.
159	! Variables for the meteoritic flux:
160	integer,parameter :: nbin_meteor=100
161	integer,parameter :: nlev_meteor=130
162	double precision meteor_ccn(ngrid,nlay,100) !100=nbinco2_cld !!!
163	double precision,save :: meteor(130,100)
164	double precision mtemp(100),pression_meteor(130)
165	logical file_ok
166	integer read_ok
167	integer nelem,lebon1,lebon2
168	double precision :: ltemp1(130),ltemp2(130)
169	integer ibin,uMeteor,j
170
171	IF (firstcall) THEN
172	!=============================================================
173	! 0. Definition of the size grid
174	!=============================================================
175	c rad_cldco2 is the primary radius grid used for microphysics computation.
176	c The grid spacing is computed assuming a constant volume ratio
177	c between two consecutive bins; i.e. vrat_cld.
178	c vrat_cld is determined from the boundary values of the size grid:
179	c rmin_cld and rmax_cld.
180	c The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
181	c dr_cld is the width of each rad_cldco2 bin.
182
183	vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
184	vrat_cld = dexp(vrat_cld)
185	rb_cldco2(1) = rbmin_cld
186	rad_cldco2(1) = rmin_cld
187	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
188	do i=1,nbinco2_cld-1
189	rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.)
190	vol_cld(i+1) = vol_cld(i) * vrat_cld
191	enddo
192	do i=1,nbinco2_cld
193	rb_cldco2(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
194	& rad_cldco2(i)
195	dr_cld(i) = rb_cldco2(i+1) - rb_cldco2(i)
196	enddo
197	rb_cldco2(nbinco2_cld+1) = rbmax_cld
198	dr_cld(nbinco2_cld) = rb_cldco2(nbinco2_cld+1) -
199	& rb_cldco2(nbinco2_cld)
200	print*, ' '
201	print*,'Microphysics co2: size bin information:'
202	print*,'i,rb_cldco2(i), rad_cldco2(i),dr_cld(i)'
203	print*,'-----------------------------------'
204	do i=1,nbinco2_cld
205	write(*,'(i3,3x,3(e13.6,4x))') i,rb_cldco2(i), rad_cldco2(i),
206	& dr_cld(i)
207	enddo
208	write(*,'(i3,3x,e13.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1)
209	print*,'-----------------------------------'
210	do i=1,nbinco2_cld+1
211	rb_cldco2(i) = dlog(rb_cldco2(i)) !! we save that so that it is not computed
212	!! at each timestep and gridpoint
213	enddo
214	c Contact parameter of co2 ice on dst ( m=cos(theta) )
215	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
216	c mteta = 0.952
217	c mtetaco2 = 0.952
218	write(,) 'co2_param contact parameter:', mtetaco2
219	c Volume of a co2 molecule (m3)
220	vo1 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
221	c below, vo1 will be recomputed using real rho ice co2 (T-dependant)
222	vo1co2=vo1 ! AJOUT JA
223	c Variance of the ice and CCN distributions
224	sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
225	sigma_ice = sqrt(log(1.+nuice_sed))
226
227
228	write(,) 'Variance of ice & CCN CO2 distribs :', sigma_iceco2
229	write(,) 'nuice for co2 ice sedimentation:', nuiceco2_sed
230	write(,) 'Volume of a co2 molecule:', vo1co2
231
232	coeffh2o=0
233	if (co2useh2o) then
234	write(,)
235	write(,) "co2useh2o=.true. in callphys.def"
236	write(,) "This means water ice particles can "
237	write(,) "serve as CCN for CO2 microphysics"
238	coeffh2o=1
239	endif
240	if (meteo_flux) then
241	write(,)
242	write(,) "meteo_flux=.true. in callphys.def"
243	write(,) "meteoritic dust particles are available"
244	write(,) "for co2 ice nucleation! "
245	write(,) "Flux given by J. Plane (pressions,size bins)"
246	! Initialisation of the flux: it is constant and is it saved
247	!We must interpolate the table to the GCM pressures
248	INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))//
249	& '/Meteo_flux_Plane.dat', EXIST=file_ok)
250	IF (.not. file_ok) THEN
251	write(,) 'file Meteo_flux_Plane.dat should be in '
252	& ,datafile
253	STOP
254	endif
255	!used Variables
256	open(newunit=uMeteor,file=trim(datafile)//
257	& '/Meteo_flux_Plane.dat'
258	& ,FORM='formatted')
259	!13000 records (130 pressions x 100 bin sizes)
260	read(uMeteor,*) !skip 1 line
261	do i=1,130
262	read(uMeteor,*,iostat=read_ok) pression_meteor(i)
263	if (read_ok==0) then
264	write(,) pression_meteor(i)
265	else
266	write(,) 'Error reading Meteo_flux_Plane.dat'
267	call abort_physic("CO2clouds",
268	& "Error reading Meteo_flux_Plane.dat",1)
269	endif
270	enddo
271	read(uMeteor,*) !skip 1 line
272	do i=1,130
273	do j=1,100 ! les mêmes 100 bins size que la distri nuclea: on touche pas
274	read(uMeteor,'(F12.6)',iostat=read_ok) meteor(i,j)
275	if (read_ok/=0) then
276	write(,) 'Error reading Meteo_flux_Plane.dat'
277	call abort_physic("CO2clouds",
278	& "Error reading Meteo_flux_Plane.dat",1)
279	endif
280	enddo
281	!On doit maintenant réinterpoler le tableau(130,100) sur les pressions du GCM (nlay,100)
282	enddo
283	close(uMeteor)
284	write(,) "File meteo_flux read, end of firstcall in co2 micro"
285	endif !of if meteo_flux
286
287	! Parameter values for Latent heat computation
288	l0=595594d0
289	l1=903.111d0
290	l2=-11.5959d0
291	l3=0.0528288d0
292	l4=-0.000103183d0
293	firstcall=.false.
294	END IF
295	!=============================================================
296	! 1. Initialisation
297	!=============================================================
298	flag_pourri=0
299	meteor_ccn(:,:,:)=0.
300	rice(:,:) = 1.e-8
301	riceco2(:,:) = 1.e-11
302
303	c Initialize the tendencies
304	pdqcloudco2(1:ngrid,1:nlay,1:nq)=0.
305	pdtcloudco2(1:ngrid,1:nlay)=0.
306
307	c pt temperature layer; pdt dT.s-1
308	c pq traceur kg/kg; pdq tendance idem .s-1
309	zt(1:ngrid,1:nlay) =
310	& pt(1:ngrid,1:nlay) +
311	& pdt(1:ngrid,1:nlay) * ptimestep
312	zq(1:ngrid,1:nlay,1:nq) =
313	& pq(1:ngrid,1:nlay,1:nq) +
314	& pdq(1:ngrid,1:nlay,1:nq) * ptimestep
315	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
316	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
317
318	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
319	!=============================================================
320	! 2. Compute saturation
321	!=============================================================
322	dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
323	dev3 = 1. / ( sqrt(2.) * sigma_ice )
324	call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2)
325	zqco2=zq(:,:,igcm_co2)+zq(:,:,igcm_co2_ice)
326	CALL tcondco2(ngrid,nlay,pplay,zqco2,tcond)
327	!=============================================================
328	! 3. Bonus: additional meteoritic particles for nucleation
329	!=============================================================
330	if (meteo_flux) then
331	!pression_meteo(130)
332	!pplev(ngrid,nlay+1)
333	!meteo(130,100)
334	!resultat: meteo_ccn(ngrid,nlay,100)
335	do l=1,nlay
336	do ig=1,ngrid
337	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
338	ltemp1=abs(pression_meteor(:)-pplev(ig,l))
339	ltemp2=abs(pression_meteor(:)-pplev(ig,l+1))
340	lebon1=minloc(ltemp1,DIM=1)
341	lebon2=minloc(ltemp2,DIM=1)
342	nelem=lebon2-lebon1+1.
343	mtemp(:)=0d0 !mtemp(100) : valeurs pour les 100bins
344	do ibin=1,100
345	mtemp(ibin)=sum(meteor(lebon1:lebon2,ibin))
346	enddo
347	meteor_ccn(ig,l,:)=mtemp(:)/nelem/masse(ig,l) !Par kg air
348	csi par m carre, x epaisseur/masse pour par kg/air
349	!write(,) "masse air ig l=",masse(ig,l)
350	!check original unit with J. Plane
351	enddo
352	enddo
353	endif
354	c ------------------------------------------------------------------------
355	c --------- Actual microphysics : Main loop over the GCM's grid ---------
356	c ------------------------------------------------------------------------
357	DO l=1,nlay
358	DO ig=1,ngrid
359	c Get the partial pressure of co2 vapor and its saturation ratio
360	pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/44.01) * pplay(ig,l)
361	satu = pco2 / zqsat(ig,l)
362
363	rho_ice_co2T(ig,l)=1000.(1.72391-2.53e-4zt(ig,l)
364	& -2.87e-6zt(ig,l)zt(ig,l)) !T-dependant CO2 ice density
365	vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l))
366	rho_ice_co2=rho_ice_co2T(ig,l)
367
368	!=============================================================
369	!4. Nucleation
370	!=============================================================
371	IF ( satu .ge. 1 ) THEN ! if there is condensation
372
373	c call updaterccn(zq(ig,l,igcm_dust_mass),
374	c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
375
376	!We do rdust computation "by hand" because we don't want to
377	! change the mininumum rdust value in updaterad... It would
378	! mess up various part of the GCM !
379
380	rdust(ig,l)= zq(ig,l,igcm_dust_mass)
381	& *0.75/pi/rho_dust
382	& / zq(ig,l,igcm_dust_number)
383	rdust(ig,l)= rdust(ig,l)**(1./3.)
384	if (zq(ig,l,igcm_dust_mass)*tauscaling(ig) .le. 1.e-20
385	& .or. zq(ig,l,igcm_dust_number)*tauscaling(ig) .le.1
386	& .or. rdust(ig,l) .le. 1.e-9) then
387	rdust(ig,l)=1.e-9
388	endif
389	rdust(ig,l)=min(5.e-4,rdust(ig,l))
390
391	c Expand the dust moments into a binned distribution
392	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)+1.e-30
393	No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30
394	Rn = rdust(ig,l)
395	Rn = -dlog(Rn)
396	Rm = Rn - 3. * sigma_iceco2*sigma_iceco2
397	n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
398	m_derf = derf( (rb_cldco2(1)+Rm) *dev2)
399
400	do i = 1, nbinco2_cld
401	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
402	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
403	n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
404	m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
405	n_aer(i) = n_aer(i) + 0.5 * No * n_derf +
406	& meteor_ccn(ig,l,i) !Ajout meteor_ccn particles
407	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
408	m_aer2(i)=4./3.pirho_dust
409	& meteor_ccn(ig,l,i)rad_cldco2(i)*rad_cldco2(i)
410	& *rad_cldco2(i)
411	enddo
412	if (meteo_flux) m_aer(:)=m_aer(:)+m_aer2(:)
413
414	!Same but with h2o particles as CCN
415	call updaterice_micro(zq(ig,l,igcm_h2o_ice),
416	& zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
417	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
418	Mo = zq(ig,l,igcm_h2o_ice) +
419	& zq(ig,l,igcm_ccn_mass)tauscaling(ig)+1.e-30!tauscaling(ig)
420	! & + 1.e-30 !Total mass of H20 crystals,CCN included
421	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
422	Rn = rice(ig,l)
423	Rn = -dlog(Rn)
424	Rm = Rn - 3. * sigma_ice*sigma_ice
425	n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
426	m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
427	do i = 1, nbinco2_cld
428	n_aer_h2oice(i) = -0.5 * No * n_derf
429	m_aer_h2oice(i) = -0.5 * Mo * m_derf
430	n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
431	m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
432	n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf
433	m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf
434	rad_h2oice(i) = rad_cldco2(i)
435	enddo
436	!Call to nucleation routine
437	call nucleaCO2(dble(pco2),zt(ig,l),dble(satu)
438	& ,n_aer,rate,n_aer_h2oice
439	& ,rad_h2oice,rateh2o,vo2co2)
440	dN = 0.
441	dM = 0.
442	dNh2o = 0.
443	dMh2o = 0.
444	do i = 1, nbinco2_cld
445	Proba=1.0-dexp(-1.ptimesteprate(i))
446	Probah2o=coeffh2o(1.0-dexp(-1.ptimestep*rateh2o(i))) !if co2useh2o=.false., this is =0
447	dNh2o = dNh2o + n_aer_h2oice(i) * Probah2o
448	dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o
449	dN = dN + n_aer(i) * Proba
450	dM = dM + m_aer(i) * Proba
451	enddo
452	! dM masse activée (kg) et dN nb particules par kg d'air
453	! Now increment CCN tracers and update dust tracers
454	dNN= dN/tauscaling(ig)
455	dMM= dM/tauscaling(ig)
456	dNN=min(dNN,zq(ig,l,igcm_dust_number))
457	dMM=min(dMM,zq(ig,l,igcm_dust_mass))
458	zq(ig,l,igcm_ccnco2_mass) =
459	& zq(ig,l,igcm_ccnco2_mass) + dMM
460	zq(ig,l,igcm_ccnco2_number) =
461	& zq(ig,l,igcm_ccnco2_number) + dNN
462	zq(ig,l,igcm_dust_mass)= zq(ig,l,igcm_dust_mass)-dMM
463	zq(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)-dNN
464
465	c Update CCN for CO2 nucleating on H2O CCN :
466	! Warning: must keep memory of it
467	if (co2useh2o) then
468	dNNh2o=dNh2o/tauscaling(ig)
469	dNNh2o=min(dNNh2o,zq(ig,l,igcm_ccn_number))
470	ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice)
471	& +zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
472	dMh2o_ice=dMh2ozq(ig,l,igcm_h2o_ice)ratioh2o_ccn
473	dMh2o_ccn=dMh2ozq(ig,l,igcm_ccn_mass)
474	& tauscaling(ig)*ratioh2o_ccn
475	dMh2o_ccn=dMh2o_ccn/tauscaling(ig)
476	dMh2o_ccn=min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
477	dMh2o_ice=min(dMh2o_ice,zq(ig,l,igcm_h2o_ice))
478	zq(ig,l,igcm_ccnco2_mass) =
479	& zq(ig,l,igcm_ccnco2_mass) + dMh2o_ice+dMh2o_ccn
480	zq(ig,l,igcm_ccnco2_number) =
481	& zq(ig,l,igcm_ccnco2_number) + dNNh2o
482	zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)-dNNh2o
483	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)-dMh2o_ice
484	zq(ig,l,igcm_ccn_mass)= zq(ig,l,igcm_ccn_mass)-dMh2o_ccn
485	memdMMh2o(ig,l)=memdMMh2o(ig,l)+dMh2o_ice
486	memdMMccn(ig,l)=memdMMccn(ig,l)+dMh2o_ccn
487	memdNNccn(ig,l)=memdNNccn(ig,l)+dNNh2o
488	endif !of if co2useh2o
489	ENDIF ! of is satu >1
490	!=============================================================
491	! 4. Ice growth: scheme for radius evolution
492	!=============================================================
493
494	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
495	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
496	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
497	No = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig)+1.e-30
498	IF (No .ge. 1)THEN ! we trigger crystal growth
499	call updaterice_microco2(zq(ig,l,igcm_co2_ice),
500	& zq(ig,l,igcm_ccnco2_mass),zq(ig,l,igcm_ccnco2_number),
501	& tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
502
503	Ic_rice=0.
504	flag_pourri=0
505	lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 +
506	& l3 * zt(ig,l)*3 + l4 zt(ig,l)**4 !J.kg-1
507	facteurmax=abs(Tcond(ig,l)-zt(ig,l))*cpp/lw
508	!specific heat of co2 ice = 1000 J.kg-1.K-1
509	!specific heat of atm cpp = 744.5 J.kg-1.K-1
510
511	ccccccc call scheme of microphys. mass growth for CO2
512	call massflowrateCO2(pplay(ig,l),zt(ig,l),
513	& satu,riceco2(ig,l),mmean(ig,l),Ic_rice)
514	c Ic_rice Mass transfer rate (kg/s) for a rice particle >0 si croissance !
515
516	if (isnan(Ic_rice) .or. Ic_rice .eq. 0.) then
517	Ic_rice=0.
518	flag_pourri=1
519	pdtcloudco2(ig,l)=-pdt(ig,l)
520	dMice=0
521
522	else
523	dMice=zq(ig,l,igcm_ccnco2_number)Ic_riceptimestep
524	& *tauscaling(ig) ! Kg par kg d'air, >0 si croissance !
525	!kg.s-1 par particule * nb particule par kg air*s
526	! = kg par kg air
527
528	dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
529	dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
530	!facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
531	! latent heat release >0 if growth i.e. if dMice >0
532	pdtcloudco2(ig,l)=dMice*lw/cpp/ptimestep
533	! kgco2/kgair* J/kgco2 * 1/(J.kgair-1.K-1)/s= K par seconde
534	!Now update tracers
535	zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice)+dMice
536	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)-dMice
537	endif
538
539	!=============================================================
540	! 5. Dust cores releasing if no more co2 ice :
541	!=============================================================
542
543	if (zq(ig,l,igcm_co2_ice).le. 1.e-25)THEN
544	!On sublime tout
545	if (co2useh2o) then
546	if (memdMMccn(ig,l) .gt. 0) then
547	zq(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
548	& +memdMMccn(ig,l)
549	endif
550	if (memdMMh2o(ig,l) .gt. 0) then
551	zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
552	& +memdMMh2o(ig,l)
553	endif
554
555	if (memdNNccn(ig,l) .gt. 0) then
556	zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
557	& +memdNNccn(ig,l)
558	endif
559	endif
560	zq(ig,l,igcm_dust_mass) =
561	& zq(ig,l,igcm_dust_mass)
562	& + zq(ig,l,igcm_ccnco2_mass)-
563	& (memdMMh2o(ig,l)+memdMMccn(ig,l))
564	zq(ig,l,igcm_dust_number) =
565	& zq(ig,l,igcm_dust_number)
566	& + zq(ig,l,igcm_ccnco2_number)-memdNNccn(ig,l)
567
568	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)
569	& + zq(ig,l,igcm_co2_ice)
570
571	zq(ig,l,igcm_ccnco2_mass)=0.
572	zq(ig,l,igcm_co2_ice)=0.
573	zq(ig,l,igcm_ccnco2_number)=0.
574	memdNNccn(ig,l)=0.
575	memdMMh2o(ig,l)=0.
576	memdMMccn(ig,l)=0.
577	riceco2(ig,l)=0.
578
579	endif !of if co2_ice <1e-25
580
581	ENDIF ! of if NCCN > 1
582	ENDDO ! of ig loop
583	ENDDO ! of nlayer loop
584
585	!=============================================================
586	! 6. END: get cloud tendencies
587	!=============================================================
588
589	! Get cloud tendencies
590	pdqcloudco2(1:ngrid,1:nlay,igcm_co2) =
591	& (zq(1:ngrid,1:nlay,igcm_co2) -
592	& zq0(1:ngrid,1:nlay,igcm_co2))/ptimestep
593	pdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) =
594	& (zq(1:ngrid,1:nlay,igcm_co2_ice) -
595	& zq0(1:ngrid,1:nlay,igcm_co2_ice))/ptimestep
596	pdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
597	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
598	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep
599	pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
600	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
601	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep
602	pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
603	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
604	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep
605	pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) =
606	& (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) -
607	& zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/ptimestep
608	pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) =
609	& (zq(1:ngrid,1:nlay,igcm_ccnco2_number) -
610	& zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/ptimestep
611	pdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) =
612	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
613	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep
614	pdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) =
615	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
616	& zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
617
618	end
619
620
621

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

Download in other formats:

Original Format