Context Navigation

improvedCO2clouds.F @ 1818

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

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

Download in other formats:

Original Format