Context Navigation

improvedco2clouds_mod.F @ 2316

Last change on this file since 2316 was 2261, checked in by emillour, 5 years ago
Mars GCM: Revert co2cloud.F and improvedco2clouds_mod.F to what they were in r2257 (they were accidentaly reverted to earlier version by r2260). EM
Property svn:executable set to ``*
File size: 29.6 KB

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

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

Download in other formats:

Original Format