Context Navigation

improvedCO2clouds.F @ 1918

Last change on this file since 1918 was 1918, checked in by emillour, 7 years ago

Mars GCM:
Code cleanup:

remove "comorbit.h" since it is no longer used.
turn "datafile.h" into module datafile_mod.F90 (and rename variable "datafile" as "datadir" since it stores the path to the datafile directory).

Property svn:executable set to *

File size: 28.8 KB

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

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

Download in other formats:

Original Format