Context Navigation

improvedCO2clouds.F @ 1758

Last change on this file since 1758 was 1720, checked in by jaudouard, 9 years ago
Update on CO2 ice clouds scheme
Property svn:executable set to ``*
File size: 35.1 KB

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

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

Download in other formats:

Original Format