Context Navigation

improvedCO2clouds.F @ 1711

Last change on this file since 1711 was 1685, checked in by jaudouard, 8 years ago
Further modifications on CO2 clouds scheme. Water ice clouds can now serve as CCN for CO2 clouds
Property svn:executable set to ``*
File size: 32.3 KB

Line
1	subroutine improvedCO2clouds(ngrid,nlay,ptimestep,
2	& pplay,pt,pdt,
3	& pq,pdq,pdqcloudco2,pdtcloudco2,
4	& nq,tauscaling,
5	& memdMMccn,memdMMh2o,memdNNccn)
6	! to use 'getin'
7	USE comcstfi_h
8	USE ioipsl_getincom
9	USE updaterad
10	use tracer_mod
11	!, only: rho_ice_co2, nuiceco2_sed, igcm_co2,
12	! & rho_ice,igcm_h2o_ice, igcm_ccn_number,
13	! & igcm_co2_ice, igcm_dust_mass,
14	! & igcm_dust_number, igcm_ccnco2_mass,
15	! & igcm_ccnco2_number
16	use conc_mod, only: mmean
17	implicit none
18
19
20	c------------------------------------------------------------------
21	c This routine is used to form CO2 clouds when a parcel of the GCM is
22	c saturated. It includes the ability to have supersaturation, a
23	c computation of the nucleation rates, growthrates and the
24	c scavenging of dust particles by clouds.
25	c It is worth noting that the amount of dust is computed using the
26	c dust optical depth computed in aeropacity.F. That's why
27	c the variable called "tauscaling" is used to convert
28	c pq(dust_mass) and pq(dust_number), which are relative
29	c quantities, to absolute and realistic quantities stored in zq.
30	c This has to be done to convert the inputs into absolute
31	c values, but also to convert the outputs back into relative
32	c values which are then used by the sedimentation and advection
33	c schemes.
34	c CO2 ice particles can nucleate on both dust and on water ice particles
35	c When CO2 ice is deposited onto a water ice particles, the particle is
36	c removed from the water tracers.
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------------------------------------------------------------------
48	!#include "dimensions.h"
49	!#include "dimphys.h"
50	#include "callkeys.h"
51	!#include "tracer.h"
52	!#include "comgeomfi.h"
53	!#include "dimradmars.h"
54	#include "microphys.h"
55	!#include "microphysCO2.h"
56	!#include "conc.h"
57	c------------------------------------------------------------------
58	c Inputs:
59
60	INTEGER ngrid,nlay
61	integer nq ! nombre de traceurs
62	REAL ptimestep ! pas de temps physique (s)
63	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
64
65	REAL pt(ngrid,nlay) ! temperature at the middle of the
66	! layers (K)
67	REAL pdt(ngrid,nlay) ! tendance temperature des autres
68	! param.
69	REAL pq(ngrid,nlay,nq) ! traceur (kg/kg)
70	REAL pdq(ngrid,nlay,nq) ! tendance avant condensation
71	! (kg/kg.s-1)
72	REAL tauscaling(ngrid) ! Convertion factor for qdust and Ndust
73
74	REAL rice(ngrid,nlay) ! Water Ice mass mean radius (m)
75	! used for nucleation of CO2 on ice-coated ccns
76
77	c Outputs:
78	REAL pdqcloudco2(ngrid,nlay,nq) ! tendance de la condensation
79	! CO2 (kg/kg.s-1)
80	! condensation si igcm_co2_ice
81	REAL pdtcloudco2(ngrid,nlay) ! tendance temperature due
82	! a la chaleur latente
83
84	c------------------------------------------------------------------
85	c Local variables:
86
87	LOGICAL firstcall
88	DATA firstcall/.true./
89	SAVE firstcall
90
91	REAL*8 derf ! Error function
92	!external derf
93
94	!REAL*8 massflowrateCO2
95	!external massflowrateCO2
96
97	INTEGER ig,l,i
98
99	REAL zq(ngrid,nlay,nq) ! local value of tracers
100	REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
101	REAL zt(ngrid,nlay) ! local value of temperature
102	REAL zqsat(ngrid,nlay) ! saturation vapor pressure for CO2
103	REAL lw !Latent heat of sublimation (J.kg-1)
104	REAL l0,l1,l2,l3,l4
105	REAL cste
106	DOUBLE PRECISION dMice ! mass of condensed ice
107	DOUBLE PRECISION sumcheck
108	DOUBLE PRECISION pco2 ! Co2 vapor partial pressure (Pa)
109	DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice
110	DOUBLE PRECISION Mo,No
111	DOUBLE PRECISION Rn, Rm, dev2, n_derf, m_derf
112	DOUBLE PRECISION memdMMccn(ngrid,nlay)
113	DOUBLE PRECISION memdMMh2o(ngrid,nlay)
114	DOUBLE PRECISION memdNNccn(ngrid,nlay)
115
116	! Radius used by the microphysical scheme (m)
117	DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin
118	DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin
119
120	DOUBLE PRECISION n_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
121	DOUBLE PRECISION m_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
122	DOUBLE PRECISION rad_h2oice(nbinco2_cld)
123
124	c REAL*8 sigco2 ! Co2-ice/air surface tension (N.m)
125	c EXTERNAL sigco2
126
127	DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM,dNNh2o,dMMh2o
128	DOUBLE PRECISION rate(nbinco2_cld) ! nucleation rate
129	DOUBLE PRECISION rateh2o(nbinco2_cld) ! nucleation rate
130	REAL seq
131	DOUBLE PRECISION rho_ice_co2T(ngrid,nlay)
132	DOUBLE PRECISION riceco2(ngrid,nlay) ! CO2Ice mean radius (m)
133	REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
134
135	REAL rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
136	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
137
138	c REAL res ! Resistance growth
139	DOUBLE PRECISION Ic_rice ! Mass transfer rate CO2 ice crystal
140	DOUBLE PRECISION ratioh2o_ccn
141	DOUBLE PRECISION vo2co2
142	c Parameters of the size discretization
143	c used by the microphysical scheme
144	DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-10 ! Minimum radius (m)
145	DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-4 ! Maximum radius (m)
146	DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-11
147	! Minimum boundary radius (m)
148	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-3 ! Maximum boundary radius (m)
149	DOUBLE PRECISION vrat_cld ! Volume ratio
150	DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
151	SAVE rb_cldco2
152
153	DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m)
154	DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3)
155
156	DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o
157	REAL sigma_iceco2 ! Variance of the ice and CCN distributions
158	SAVE sigma_iceco2
159
160	DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2
161
162	c----------------------------------
163	c TESTS
164
165	INTEGER countcells
166
167	LOGICAL test_flag ! flag for test/debuging outputs
168	SAVE test_flag
169
170
171	REAL satubf(ngrid,nlay),satuaf(ngrid,nlay)
172	REAL res_out(ngrid,nlay)
173
174
175	c------------------------------------------------------------------
176
177	IF (firstcall) THEN
178	!=============================================================
179	! 0. Definition of the size grid
180	!=============================================================
181	c rad_cldco2 is the primary radius grid used for microphysics computation.
182	c The grid spacing is computed assuming a constant volume ratio
183	c between two consecutive bins; i.e. vrat_cld.
184	c vrat_cld is determined from the boundary values of the size grid:
185	c rmin_cld and rmax_cld.
186	c The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
187	c dr_cld is the width of each rad_cldco2 bin.
188
189	c Volume ratio between two adjacent bins
190	! vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
191	! vrat_cld = exp(vrat_cld)
192	vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
193	vrat_cld = exp(vrat_cld)
194	c write(,) "vrat_cld", vrat_cld
195
196	rb_cldco2(1) = rbmin_cld
197	rad_cldco2(1) = rmin_cld
198	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
199	! vol_cld(1) = 4./3. * pi * rmin_cldrmin_cldrmin_cld
200
201	do i=1,nbinco2_cld-1
202	rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.)
203	vol_cld(i+1) = vol_cld(i) * vrat_cld
204	enddo
205
206	do i=1,nbinco2_cld
207	rb_cldco2(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
208	& rad_cldco2(i)
209	dr_cld(i) = rb_cldco2(i+1) - rb_cldco2(i)
210	enddo
211	rb_cldco2(nbinco2_cld+1) = rbmax_cld
212	dr_cld(nbinco2_cld) = rb_cldco2(nbinco2_cld+1) -
213	& rb_cldco2(nbinco2_cld)
214
215	print*, ' '
216	print*,'Microphysics co2: size bin information:'
217	print*,'i,rb_cldco2(i), rad_cldco2(i),dr_cld(i)'
218	print*,'-----------------------------------'
219	do i=1,nbinco2_cld
220	write(*,'(i3,3x,3(e12.6,4x))') i,rb_cldco2(i), rad_cldco2(i),
221	& dr_cld(i)
222	enddo
223	write(*,'(i3,3x,e12.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1)
224	print*,'-----------------------------------'
225
226	do i=1,nbinco2_cld+1
227	rb_cldco2(i) = log(rb_cldco2(i)) !! we save that so that it is not computed
228	!! at each timestep and gridpoint
229	enddo
230
231	c Contact parameter of co2 ice on dst ( m=cos(theta) )
232	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
233	c mteta = 0.952
234	c mtetaco2 = 0.952
235	write(,) 'co2_param contact parameter:', mtetaco2
236
237	c Volume of a co2 molecule (m3)
238	vo1 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
239	vo1co2=vo1 ! AJOUT JA
240	c Variance of the ice and CCN distributions
241	sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
242
243
244	c write(,) 'Variance of ice & CCN distribs :', sigma_iceco2
245	c write(,) 'nuice for sedimentation:', nuiceco2_sed
246	c write(,) 'Volume of a co2 molecule:', vo1co2
247
248	write(,) 'Variance of ice & CCN CO2 distribs :', sigma_iceco2
249	write(,) 'nuice for co2 ice sedimentation:', nuiceco2_sed
250	write(,) 'Volume of a co2 molecule:', vo1co2
251
252
253
254	test_flag = .false.
255	firstcall=.false.
256
257	END IF
258
259
260	!=============================================================
261	! 1. Initialisation
262	!=============================================================
263	!cste = 4pirho_ice*ptimestep !not used for co2
264
265	res_out(:,:) = 0
266	rice(:,:) = 1.e-8
267	riceco2(:,:) = 1.e-11
268
269	c Initialize the tendencies
270	pdqcloudco2(1:ngrid,1:nlay,1:nq)=0.
271	pdtcloudco2(1:ngrid,1:nlay)=0.
272
273	c pt temperature layer; pdt dT.s-1
274	c pq traceur kg/kg; pdq tendance idem .s-1
275	zt(1:ngrid,1:nlay) =
276	& pt(1:ngrid,1:nlay) +
277	& pdt(1:ngrid,1:nlay) * ptimestep
278	c call WRITEDIAGFI(ngrid,"Ztclouds",
279	c & "Ztclouds",'K',3,zt)
280	c call WRITEDIAGFI(ngrid,"pdtclouds",
281	c & "pdtclouds",'K',3,pdt*ptimestep)
282
283	zq(1:ngrid,1:nlay,1:nq) =
284	& pq(1:ngrid,1:nlay,1:nq) +
285	& pdq(1:ngrid,1:nlay,1:nq) * ptimestep
286
287
288	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
289	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
290
291	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
292
293	!=============================================================
294	! 2. Compute saturation
295	!=============================================================
296
297
298	dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
299
300	call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2)
301
302	countcells = 0
303
304	c Faire rice co2 update en n-1 puis a chaque microdt, mettre a jour riceco2
305
306	c Main loop over the GCM's grid
307	DO l=1,nlay
308	DO ig=1,ngrid
309
310
311	c Get the partial pressure of co2 vapor and its saturation ratio
312	pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/44.01) * pplay(ig,l)
313	c satu = zq(ig,l,igcm_co2) / zqsat(ig,l)
314	satu = pco2 / zqsat(ig,l)
315	!=============================================================
316	! 3. Nucleation
317	!=============================================================
318	rho_ice_co2T(ig,l)=1000.(1.72391-2.53e-4zt(ig,l)
319	& -2.87e-6zt(ig,l)zt(ig,l))
320	vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l)) ! m0co2
321	rho_ice_co2=rho_ice_co2T(ig,l)
322	c call updaterccn(zq(ig,l,igcm_dust_mass),
323	c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
324	c write(,) "l, pco2, satu= ",l,pco2,satu
325	IF ( satu .ge. 1d0 ) THEN ! if there is condensation
326	c write(,)
327
328	!write(,) "Zt, rho=",zt(ig,l),rho_ice_co2
329	c Masse_atm=mmean(ig,l)1.e-3pplay(ig,l)/rgp/zt(ig,l) !Kg par couche
330
331
332	c call updaterccn(zq(ig,l,igcm_dust_mass),
333	c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
334
335	c call updaterccn(zq(ig,l,igcm_dust_mass),
336	c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
337	c write(,) "Improved, l,Rdust = ",l,rdust(ig,l)
338
339	rdust(ig,l)= zq(ig,l,igcm_dust_mass)
340	& *0.75/pi/rho_dust
341	& / zq(ig,l,igcm_dust_number)
342	rdust(ig,l)= rdust(ig,l)**(1./3.)
343	!write(,) "Improved2, l,Rdust = ",l,rdust(ig,l)
344	rdust(ig,l)=max(1.e-10,rdust(ig,l))
345	rdust(ig,l)=min(5.e-5,rdust(ig,l))
346	! write(,) "Improved3,Rdust = ",rdust(ig,l)
347
348	c Expand the dust moments into a binned distribution
349	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)+1.e-30
350	No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30
351	c write(,) "Improved dust number, mass = ",
352	c & zq(ig,l,igcm_dust_number)* tauscaling(ig),
353	c & zq(ig,l,igcm_dust_mass)* tauscaling(ig)
354	c write(,) "No, Mo = ",No, Mo
355	Rn = rdust(ig,l)
356	Rn = -log(Rn)
357	Rm = Rn - 3. * sigma_iceco2*sigma_iceco2
358	n_derf = erf( (rb_cldco2(1)+Rn) *dev2)
359	m_derf = erf( (rb_cldco2(1)+Rm) *dev2)
360
361	do i = 1, nbinco2_cld
362	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
363	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
364	n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
365	m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
366	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
367	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
368	c write(,) "i, rad_cldco2(i) = ",i, rad_cldco2(i),
369	c & n_aer(i)
370	enddo
371
372
373	sumcheck = 0
374	do i = 1, nbinco2_cld
375	sumcheck = sumcheck + n_aer(i)
376	enddo
377	sumcheck = abs(sumcheck/No - 1)
378	if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
379	print*, "WARNING, No sumcheck PROBLEM"
380	print*, "sumcheck, No",sumcheck, No
381	print*, "rdust =",rdust(ig,l)
382	print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
383	print*, "Dust binned distribution", n_aer
384	STOP
385	endif
386
387	sumcheck = 0
388	do i = 1, nbinco2_cld
389	sumcheck = sumcheck + m_aer(i)
390	enddo
391	sumcheck = abs(sumcheck/Mo - 1)
392	if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld))
393	& then
394	print*, "WARNING, Mo sumcheck PROBLEM"
395	print*, "sumcheck, Mo",sumcheck, Mo
396	print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l
397	print*, "Dust binned distribution", m_aer
398	STOP
399	endif
400
401	call updaterice_micro(
402	& zq(ig,l,igcm_h2o_ice), ! ice mass
403	& zq(ig,l,igcm_ccn_mass), ! ccn mass
404	& zq(ig,l,igcm_ccn_number), ! ccn number
405	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
406	! rice radius of CCN + H20 crystal
407	!write(,) "Improved1 Rice=",rice(ig,l)
408	rice(ig,l)=max(1.e-10,rice(ig,l))
409	rice(ig,l)=min(5.e-5,rice(ig,l))
410	!write(,) "Improved2 Rice=",rice(ig,l)
411	Mo = zq(ig,l,igcm_h2o_ice) +
412	& zq(ig,l,igcm_ccn_mass)*tauscaling(ig)
413	& + 1.e-30 !Total mass of H20 crystals,CCN included
414	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
415	Rn = rice(ig,l)
416	Rn = -log(Rn)
417	Rm = Rn - 3. * sigma_iceco2*sigma_iceco2
418	n_derf = erf( (rb_cldco2(1)+Rn) *dev2)
419	m_derf = erf( (rb_cldco2(1)+Rm) *dev2)
420	do i = 1, nbinco2_cld
421	n_aer_h2oice(i) = -0.5 * No * n_derf !! this ith previously computed
422	m_aer_h2oice(i) = -0.5 * Mo * m_derf !! this ith previously computed
423	n_derf = derf( (rb_cldco2(i+1)+Rn) *dev2)
424	m_derf = derf( (rb_cldco2(i+1)+Rm) *dev2)
425	n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf ! vector not really needed - temp var - listo
426	m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf ! vector not really needed - temp var
427	rad_h2oice(i) = rad_cldco2(i)
428	c write(,) "before nuc, i,rad_h2o(i)= ",i,rad_cldco2(i)
429	c & ,m_aer_h2oice(i),n_aer_h2oice(i)
430	enddo
431	sumcheck = 0
432	do i = 1, nbinco2_cld
433	sumcheck = sumcheck + n_aer_h2oice(i)
434	enddo
435	sumcheck = abs(sumcheck/No - 1)
436	if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
437	print*, "WARNING, Noh2o sumcheck PROBLEM"
438	print*, "sumcheck, No",sumcheck, No
439	print*, "rice =",rice(ig,l)
440	print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
441	print*, "Dust binned distribution", n_aer_h2oice
442	STOP
443	endif
444
445
446	c Get the rates of nucleation
447	call nucleaCO2(dble(pco2),zt(ig,l),dble(satu)
448	& ,n_aer,rate,n_aer_h2oice
449	& ,rad_h2oice,rateh2o,vo2co2)
450	! regarder rateh20, et mettre = 0 si non nul pour le moment
451	dN = 0.
452	dM = 0.
453	dNh2o = 0.
454	dMh2o = 0.
455	do i = 1, nbinco2_cld
456	Proba=1.0-dexp(-1.ptimesteprate(i))
457	Probah2o=1.0-dexp(-1.ptimesteprateh2o(i))
458
459	dNh2o = dNh2o + n_aer_h2oice(i) * Probah2o
460	dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o
461
462	dN = dN + n_aer(i) * Proba
463	dM = dM + m_aer(i) * Proba
464	c write(,) "i, dNi, dN= ",i,n_aer(i)*Proba,dN
465	enddo
466
467	! dM masse activée (kg) et dN nb particules par kg d'air
468
469	c write(,) " nuclea dM = ",dM/tauscaling(ig),
470	c & " nuclea dN = ", dN/tauscaling(ig)
471
472	dNN= dN/tauscaling(ig)
473	dMM= dM/tauscaling(ig)
474	dNNh2o=dNh2o/tauscaling(ig)
475	dMMh2o=dMh2o/tauscaling(ig)
476
477	dNN=min(dNN,abs(zq(ig,l,igcm_dust_number)))
478	dMM=min(dMM,abs(zq(ig,l,igcm_dust_mass)))
479	c
480	c write(,) "Nuclea dNN crees=",dNN
481	dNNh2o=min(dNNh2o,abs(zq(ig,l,igcm_ccn_number)))
482	dMMh2o=min(dMMh2o,abs(zq(ig,l,igcm_h2o_ice)/tauscaling(ig)
483	& +zq(ig,l,igcm_ccn_mass))) !Total mass of H2O crystals available
484
485	c write(,) "Nuclea dNNh2o crees=",dNNh2o
486
487	c Update CCNs for CO2 crystals
488	! WARNING dM dMh2o, interaction nuages eau-co2 -- h20 set to 0 for now
489	zq(ig,l,igcm_ccnco2_mass) =
490	& zq(ig,l,igcm_ccnco2_mass) + dMM
491	zq(ig,l,igcm_ccnco2_number) =
492	& zq(ig,l,igcm_ccnco2_number) + dNN
493	zq(ig,l,igcm_dust_mass) =
494	& zq(ig,l,igcm_dust_mass) - dMM
495	zq(ig,l,igcm_dust_number) =
496	& zq(ig,l,igcm_dust_number) - dNN
497
498	c Update CCN for CO2 nucleating on H2O CCN :
499	! Warning: must keep memory of it
500	zq(ig,l,igcm_ccnco2_mass) =
501	& zq(ig,l,igcm_ccnco2_mass) + dMMh2o
502	zq(ig,l,igcm_ccnco2_number) =
503	& zq(ig,l,igcm_ccnco2_number) + dNNh2o
504
505
506	zq(ig,l,igcm_ccn_number) =
507	& zq(ig,l,igcm_ccn_number) - dNNh2o
508
509	ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice)
510	& +zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
511
512
513	memdMMh2o(ig,l)= memdMMh2o(ig,l)+zq(ig,l,igcm_h2o_ice)*
514	& dMMh2o*ratioh2o_ccn
515	memdMMccn(ig,l)= memdMMccn(ig,l)+zq(ig,l,igcm_ccn_mass)*
516	& tauscaling(ig)dMMh2oratioh2o_ccn
517	memdNNccn(ig,l)=memdNNccn(ig,l)+dNNh2o
518	c if (dMMh2o .gt. 0) then
519	c write(,) 'test h2o'
520	c write(,) "dMMh2o=",dMMh2o
521	c write(,) "2 =",zq(ig,l,igcm_ccn_mass)tauscaling(ig)
522	c & dMMh2oratioh2o_ccn+zq(ig,l,igcm_h2o_ice)
523	c & dMMh2o*ratioh2o_ccn
524	c write(,) "3=",zq(ig,l,igcm_ccn_mass)tauscaling(ig)
525	c & dMMh2o*ratioh2o_ccn
526	c write(,) "4=",zq(ig,l,igcm_h2o_ice)*
527	c & dMMh2o*ratioh2o_ccn
528	c write(,) "tauscaling=",tauscaling(ig)
529	c endif
530	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)*
531	& (1.-dMMh2o*ratioh2o_ccn)
532	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass)*
533	& tauscaling(ig)(1.-dMMh2oratioh2o_ccn)
534
535
536	ENDIF ! of is satu >1
537	!=============================================================
538	! 4. Ice growth: scheme for radius evolution
539	!=============================================================
540
541	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
542	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
543	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
544	c IF ( zq(ig,l,igcm_ccnco2_number)*tauscaling(ig).ge. 1.0) THEN
545
546	IF (zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) .ge. 1.)
547	& THEN
548	! we trigger crystal growth
549	c
550
551	c Niceco2 = zq(ig,l,igcm_co2_ice)
552	c Qccnco2 = zq(ig,l,igcm_ccnco2_mass)
553	c Nccnco2 = zq(ig,l,igcm_ccnco2_number)
554	c call updaterice_microco2(Niceco2,Qccnco2,Nccnco2,
555	c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
556	c write(,) "updater rice=",riceco2(ig,l)
557
558	rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass)
559	& *0.75/pi/rho_dust
560	& / zq(ig,l,igcm_ccnco2_number)
561	rdust(ig,l)= rdust(ig,l)**(1./3.)
562	rdust(ig,l)=max(1.e-10,rdust(ig,l))
563	! rdust(ig,l)=min(5.e-6,rdust(ig,l))
564
565	riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/
566	& (4.0rho_ice_co2zq(ig,l,igcm_ccnco2_number)
567	& pitauscaling(ig)) +rdust(ig,l)*rdust(ig,l)
568	& rdust(ig,l) )*(1.0/3.0)
569
570	c riceco2(ig,l)=max(1.e-10,riceco2(ig,l))
571	c riceco2(ig,l)=min(1.e-5,riceco2(ig,l))
572	! WATCH OUT: CO2 nuclei is supposed to be dust
573	! only when deriving rhocloud (otherwise would need to keep info on water embedded in co2) - listo
574	c write(,) "Rdust before growth = ",rdust(ig,l)
575	c write(,) "Riceco2 before growth = ",riceco2(ig,l)
576
577	!! Niceco2,Qccnco2,Nccnco2
578	c Niceco2 = zq(ig,l,igcm_co2_ice)
579	c Qccnco2 = zq(ig,l,igcm_ccnco2_mass)
580	c Nccnco2 = zq(ig,l,igcm_ccnco2_number)
581	c call updaterice_microco2(Niceco2,Qccnco2,Nccnco2,
582	c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
583	!write(,) "Riceco2 update before growth = ",riceco2(ig,l)
584	c write(,) "rdust before growth = ",rdust(ig,l)
585	c write(,) "co2 before growth=",zq(ig,l,igcm_co2)
586	c write(,) "pplay before growth=",pplay(ig,l)
587	c write(,) "zt before growth =",zt(ig,l)
588	! write(,) "Satu before growth=",satu
589	c riceco2(ig,l)=max(riceco2(ig,l),rdust(ig,l))
590	No = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig)+1.e-30
591	! No nb de particules de poussieres mis à l'échelle pour donner une opacité optique
592
593	c saturation at equilibrium
594	c rice should not be too small, otherwise seq value is not valid
595	c seq = exp(2.sigco2mco2 / (rho_ice_co2rgpzt(ig,l)*
596	c & max(riceco2(ig,l),1.e-7))) !Exponant sans unité OK
597
598	ccccccc Scheme of microphys. mass growth for CO2
599
600	call massflowrateCO2(pplay(ig,l),zt(ig,l),
601	& satu,riceco2(ig,l),mmean(ig,l),Ic_rice) ! Mass transfer rate (kg/s) for a rice particle
602	! Ic_rice mass flux kg.s-1 <0 si croissance !
603	drsurdt=-1.0/(4.0piriceco2(ig,l)*
604	& riceco2(ig,l)rho_ice_co2)Ic_rice
605	dMice = No * Ic_rice*ptimestep ! Kg par kg d'air, <0 si croissance !
606	c write(,) "dMicev0 in improved = " , dMice
607
608	if (dMice .ge. 0d0) then
609	dMice = min(dMice,abs(zq(ig,l,igcm_co2_ice)))
610	else
611	dMice =-1.* min(abs(dMice),abs(zq(ig,l,igcm_co2)))
612	endif
613	riceco2(ig,l)=riceco2(ig,l)+drsurdt*ptimestep
614	c write(,) "riceco2+dr/dt = ", riceco2(ig,l)
615	c write(,) "dMice in improved = " , dMice
616
617	zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice)
618	& -dMice
619	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)+dMice
620	c write(,) "Improved zq co2 ice = ", zq(ig,l,igcm_co2_ice)
621	! countcells = countcells + 1
622
623	c riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/
624	c & (4.0rho_ice_co2pi*zq(ig,l,igcm_ccnco2_number)
625	c & tauscaling(ig)) +rdust(ig,l)rdust(ig,l)
626	c & rdust(ig,l) )*(1.0/3.0)
627	c write(,) "Improved new riceco2 = ",riceco2(ig,l)
628
629	c write(,) "new riceco2 improvedupdaterad= ",riceco2(ig,l)
630
631	! latent heat release
632
633	l0=595594.
634	l1=903.111
635	l2=-11.5959
636	l3=0.0528288
637	l4=-0.000103183
638
639	lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 +
640	& l3 * zt(ig,l)*3 + l4 zt(ig,l)**4 !J.kg-1
641	c write(,) "CPP= ",cpp ! = 744.5
642
643	pdtcloudco2(ig,l)= dMicelw/cpp/ptimestep ! kg par couche J par kg /J par K / s = K par seconde
644
645	c write(,) "pdtcloudco2 after growth = ",pdtcloudco2(ig,l)
646
647
648	!deltaT par condens/subli. qui remplace le dT du CO2 du newcondens pré-constantino
649	!PDT should be in K/s.
650	!=============================================================
651	! 5. Dust cores released, tendancies, latent heat, etc ...
652	!=============================================================
653
654	c If all the ice particles sublimate, all the condensation
655	c nuclei are released:
656
657	c !!! with CO2 ice nuclei: dust/H2O nuclei are not released because
658	c they were not subtracted to dust_number
659	c Their counter is just set to "0".
660	c (see end of section 3.) : On peut les enlever à dust
661
662	c interaction ho-co2 ici, dans la mise a jour des traceurs WARNING reflechir
663	!! Niceco2,Qccnco2,Nccnco2
664
665
666
667	rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass)
668	& *0.75/pi/rho_dust
669	& / zq(ig,l,igcm_ccnco2_number)
670	rdust(ig,l)= rdust(ig,l)**(1./3.)
671	rdust(ig,l)=max(1.e-10,rdust(ig,l))
672	!rdust(ig,l)=min(5.e-6,rdust(ig,l))
673
674	riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/
675	& (4.0rho_ice_co2zq(ig,l,igcm_ccnco2_number)
676	& pitauscaling(ig)) +rdust(ig,l)*rdust(ig,l)
677	& rdust(ig,l) )*(1.0/3.0)
678	!Niceco2 = zq(ig,l,igcm_co2_ice)
679	!Qccnco2 = zq(ig,l,igcm_ccnco2_mass)
680	!Nccnco2 = zq(ig,l,igcm_ccnco2_number)
681	c
682	c call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
683	c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
684
685	if ((zq(ig,l,igcm_co2_ice).le. 1.e-23)
686	& .or.(riceco2(ig,l) .le. rdust(ig,l))) then
687
688	c write(,) "Riceco2 improved before reset=",riceco2(ig,l)
689	c write(,) "Niceco2 improved before reset=",
690	c & zq(ig,l,igcm_co2_ice)
691	c write(,) "Rdust improved before reset=",rdust(ig,l)
692
693	if (memdMMccn(ig,l) .gt. 0) then
694	zq(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
695	& +memdMMccn(ig,l)
696	endif
697	if (memdMMh2o(ig,l) .gt. 0) then
698	zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
699	& +memdMMh2o(ig,l)
700	endif
701
702	if (memdNNccn(ig,l) .gt. 0) then
703	zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
704	& +memdNNccn(ig,l)
705	endif
706
707	if (zq(ig,l,igcm_ccnco2_mass) .gt. 1.e-30) then
708	zq(ig,l,igcm_dust_mass) =
709	& zq(ig,l,igcm_dust_mass)
710	& + zq(ig,l,igcm_ccnco2_mass)-
711	& (memdMMh2o(ig,l)+memdMMccn(ig,l))
712	endif
713	if (zq(ig,l,igcm_ccnco2_number) .gt. 1.e-30) then
714	zq(ig,l,igcm_dust_number) =
715	& zq(ig,l,igcm_dust_number)
716	& + zq(ig,l,igcm_ccnco2_number)-memdNNccn(ig,l)
717	endif
718
719	if (zq(ig,l,igcm_co2_ice) .gt. 1.e-30) then
720	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)
721	& + zq(ig,l,igcm_co2_ice)
722	endif
723
724	zq(ig,l,igcm_ccnco2_mass)=0.
725	zq(ig,l,igcm_co2_ice)=0.
726	zq(ig,l,igcm_ccnco2_number)=0.
727	memdNNccn(ig,l)=0.
728	memdMMh2o(ig,l)=0.
729	memdMMccn(ig,l)=0.
730	riceco2(ig,l)=0.
731	pdtcloudco2(ig,l)=0.
732	endif
733
734	ENDIF ! of if NCCN > 1
735	ENDDO ! of ig loop
736	ENDDO ! of nlayer loop
737
738
739	! Get cloud tendencies
740	pdqcloudco2(1:ngrid,1:nlay,igcm_co2) =
741	& (zq(1:ngrid,1:nlay,igcm_co2) -
742	& zq0(1:ngrid,1:nlay,igcm_co2))/ptimestep
743
744	pdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) =
745	& (zq(1:ngrid,1:nlay,igcm_co2_ice) -
746	& zq0(1:ngrid,1:nlay,igcm_co2_ice))/ptimestep
747
748	pdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
749	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
750	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep
751
752	pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
753	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
754	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep
755
756	pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
757	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
758	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep
759
760	pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) =
761	& (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) -
762	& zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/ptimestep
763
764	pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) =
765	& (zq(1:ngrid,1:nlay,igcm_ccnco2_number) -
766	& zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/ptimestep
767
768
769	c if (scavenging) then
770
771	pdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) =
772	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
773	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep
774
775	pdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) =
776	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
777	& zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
778	c endif
779
780	return
781	end
782
783
784

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

Download in other formats:

Original Format