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

Rev	Line
[1617]	1	subroutine improvedCO2clouds(ngrid,nlay,ptimestep,
	2	& pplay,pt,pdt,
	3	& pq,pdq,pdqcloudco2,pdtcloudco2,
[1685]	4	& nq,tauscaling,
	5	& memdMMccn,memdMMh2o,memdNNccn)
[1617]	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.
[1685]	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.
[1617]	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
[1685]	74	REAL rice(ngrid,nlay) ! Water Ice mass mean radius (m)
[1617]	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
[1685]	106	DOUBLE PRECISION dMice ! mass of condensed ice
[1617]	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
[1685]	112	DOUBLE PRECISION memdMMccn(ngrid,nlay)
	113	DOUBLE PRECISION memdMMh2o(ngrid,nlay)
	114	DOUBLE PRECISION memdNNccn(ngrid,nlay)
	115
[1617]	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
[1685]	127	DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM,dNNh2o,dMMh2o
[1617]	128	DOUBLE PRECISION rate(nbinco2_cld) ! nucleation rate
	129	DOUBLE PRECISION rateh2o(nbinco2_cld) ! nucleation rate
	130	REAL seq
[1685]	131	DOUBLE PRECISION rho_ice_co2T(ngrid,nlay)
[1617]	132	DOUBLE PRECISION riceco2(ngrid,nlay) ! CO2Ice mean radius (m)
[1685]	133	REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
	134
[1617]	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
[1685]	140	DOUBLE PRECISION ratioh2o_ccn
	141	DOUBLE PRECISION vo2co2
[1617]	142	c Parameters of the size discretization
	143	c used by the microphysical scheme
[1685]	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
[1617]	147	! Minimum boundary radius (m)
[1685]	148	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 1.e-3 ! Maximum boundary radius (m)
[1617]	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
[1685]	156	DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o
[1617]	157	REAL sigma_iceco2 ! Variance of the ice and CCN distributions
	158	SAVE sigma_iceco2
	159
[1632]	160	DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2
[1617]	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
[1685]	235	write(,) 'co2_param contact parameter:', mtetaco2
[1617]	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
[1685]	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
[1617]	251
	252
[1685]	253
[1617]	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
[1685]	266	rice(:,:) = 1.e-8
[1617]	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
[1685]	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
[1617]	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
[1685]	307	DO l=1,nlay
	308	DO ig=1,ngrid
[1617]	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	!=============================================================
[1685]	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)
[1617]	322	c call updaterccn(zq(ig,l,igcm_dust_mass),
	323	c & zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
[1685]	324	c write(,) "l, pco2, satu= ",l,pco2,satu
[1617]	325	IF ( satu .ge. 1d0 ) THEN ! if there is condensation
[1649]	326	c write(,)
[1685]	327
	328	!write(,) "Zt, rho=",zt(ig,l),rho_ice_co2
[1617]	329	c Masse_atm=mmean(ig,l)1.e-3pplay(ig,l)/rgp/zt(ig,l) !Kg par couche
	330
	331
[1685]	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))
[1649]	337	c write(,) "Improved, l,Rdust = ",l,rdust(ig,l)
[1617]	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.)
[1685]	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
[1617]	348	c Expand the dust moments into a binned distribution
[1685]	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 = ",
[1649]	352	c & zq(ig,l,igcm_dust_number)* tauscaling(ig),
	353	c & zq(ig,l,igcm_dust_mass)* tauscaling(ig)
[1617]	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
[1685]	368	c write(,) "i, rad_cldco2(i) = ",i, rad_cldco2(i),
	369	c & n_aer(i)
[1617]	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
[1685]	381	print*, "rdust =",rdust(ig,l)
[1617]	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
[1685]	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
[1617]	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)
[1685]	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)
[1617]	430	enddo
[1685]	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
[1617]	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
[1685]	449	& ,rad_h2oice,rateh2o,vo2co2)
[1617]	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))
[1685]	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
[1617]	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)
[1685]	474	dNNh2o=dNh2o/tauscaling(ig)
	475	dMMh2o=dMh2o/tauscaling(ig)
[1617]	476
[1685]	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
[1617]	484
[1685]	485	c write(,) "Nuclea dNNh2o crees=",dNNh2o
	486
[1617]	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
[1685]	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
[1617]	504
[1685]	505
	506	zq(ig,l,igcm_ccn_number) =
	507	& zq(ig,l,igcm_ccn_number) - dNNh2o
[1617]	508
[1685]	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
[1617]	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
[1685]	546	IF (zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) .ge. 1.)
	547	& THEN
	548	! we trigger crystal growth
	549	c
[1617]	550
[1685]	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
[1617]	558	rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass)
	559	& *0.75/pi/rho_dust
	560	& / zq(ig,l,igcm_ccnco2_number)
[1685]	561	rdust(ig,l)= rdust(ig,l)**(1./3.)
[1617]	562	rdust(ig,l)=max(1.e-10,rdust(ig,l))
[1685]	563	! rdust(ig,l)=min(5.e-6,rdust(ig,l))
	564
[1617]	565	riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/
[1685]	566	& (4.0rho_ice_co2zq(ig,l,igcm_ccnco2_number)
	567	& pitauscaling(ig)) +rdust(ig,l)*rdust(ig,l)
[1617]	568	& rdust(ig,l) )*(1.0/3.0)
	569
[1685]	570	c riceco2(ig,l)=max(1.e-10,riceco2(ig,l))
	571	c riceco2(ig,l)=min(1.e-5,riceco2(ig,l))
[1617]	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
[1649]	574	c write(,) "Rdust before growth = ",rdust(ig,l)
	575	c write(,) "Riceco2 before growth = ",riceco2(ig,l)
[1617]	576
[1632]	577	!! Niceco2,Qccnco2,Nccnco2
[1685]	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
[1617]	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
[1685]	605	dMice = No * Ic_rice*ptimestep ! Kg par kg d'air, <0 si croissance !
[1649]	606	c write(,) "dMicev0 in improved = " , dMice
[1617]	607
[1685]	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
[1617]	614	c write(,) "riceco2+dr/dt = ", riceco2(ig,l)
[1685]	615	c write(,) "dMice in improved = " , dMice
[1617]	616
[1685]	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
[1617]	622
[1685]	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)
[1617]	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
[1685]	645	c write(,) "pdtcloudco2 after growth = ",pdtcloudco2(ig,l)
[1617]	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
[1685]	663	!! Niceco2,Qccnco2,Nccnco2
	664
	665
[1617]	666
[1685]	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))
[1617]	684
[1685]	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
[1617]	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
[1685]	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
[1617]	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
[1685]	767
[1617]	768
[1685]	769	c if (scavenging) then
[1617]	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
[1685]	774
[1617]	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
[1685]	778	c endif
[1617]	779
	780	return
	781	end
	782
	783
	784

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

Download in other formats:

Original Format