Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

improvedco2clouds_mod.F90

Last change on this file was 4019, checked in by aslmd, 5 days ago
Mars GCM: will the real double consistency please stand up
File size: 30.8 KB

Rev	Line
[2362]	1	!======================================================================================================================!
	2	! Module: Scheme of co2 cloud formation ===============================================================================!
	3	!----------------------------------------------------------------------------------------------------------------------!
	4	! Authors: Joaquim Audouard, Constantino Listowski, Anni Määttänen
	5	! Date: 09/2016
	6	!----------------------------------------------------------------------------------------------------------------------!
	7	! Contains subroutines:
	8	! - improvedco2clouds_mod: nucleation
	9	!======================================================================================================================!
	10	module improvedco2clouds_mod
	11
	12	implicit none
	13
	14	contains
	15	!======================================================================================================================!
	16	! SUBROUTINE: improvedco2clouds =======================================================================================!
	17	!======================================================================================================================!
	18	! Subject:
	19	!---------
	20	! This routine is used to form CO2 clouds when a parcel of the GCM is saturated.
	21	!----------------------------------------------------------------------------------------------------------------------!
	22	! Comments:
	23	!----------
	24	! Adaptation for CO2 clouds based on improvedclouds_mod.F
	25	!
	26	! It includes the ability to have supersaturation, a computation of the nucleation rates, growthrates and the
	27	! scavenging of dust particles by clouds. It is worth noting that the amount of dust is computed using the dust
	28	! optical depth computed in aeropacity.F.
	29	! That's why the variable called "tauscaling" is used to convert pq(dust_mass) and pq(dust_number), which are relative
	30	! quantities, to absolute and realistic quantities stored in zq. This has to be done to convert the inputs into
	31	! absolute values, but also to convert the outputs back into relative values which are then used by the sedimentation
	32	! and advection schemes.
	33	! CO2 ice particles can nucleate on both dust and water ice particles. When CO2 ice is deposited onto a water ice
	34	! particles, the particle is removed from the water tracers. Memory of the origin of the co2 particles is kept and
	35	! thus the water cycle shouldn't be modified by this.
	36	! There is an energy limit to how much co2 can sublimate/condensate. It is defined by the difference of the GCM
	37	! temperature with the co2 condensation temperature.
	38	!
	39	! /!\ WARNING /!\
	40	! No sedimentation of the water ice origin is performed in the microphysical timestep in co2cloud.F.
	41	!
	42	! If meteoritic particles are activated and turn into co2 ice particles, then they will be reversed in the dust
	43	! tracers if the cloud sublimates.
	44	!----------------------------------------------------------------------------------------------------------------------!
	45	! Paper:
	46	!-------
	47	! see co2cloud.F90
	48	!----------------------------------------------------------------------------------------------------------------------!
	49	! Algorithm:
	50	!-----------
	51	! 0. Firstcall
	52	! 0.1. Bonus: meteoritic component, extract data
	53	! 1. Initialization
	54	! 2. Compute saturation
	55	! 3. Bonus: additional meteoritic particles for nucleation
	56	! 4. Actual microphysics: Main loop over the GCM's grid
	57	! 4.1 Nucleation
	58	! 4.2. Ice growth: scheme for radius evolution
	59	! 4.3 Dust cores releasing if no more co2 ice
	60	! 5. Get cloud tendencies
	61	!======================================================================================================================!
	62	subroutine improvedCO2clouds(ngrid, nlay, microtimestep, pplay, pplev, pteff, sum_subpdt, pqeff, sum_subpdq, &
[2562]	63	subpdqcloudco2, subpdtcloudco2, nq, tauscaling, rb_cldco2, sigma_iceco2, dev2)
[2362]	64
	65	use comcstfi_h, only: pi, g, cpp
	66
	67	use updaterad, only: updaterice_micro, updaterice_microco2, updaterccnCO2
	68
	69	use tracer_mod, only: igcm_dust_mass, igcm_dust_number, rho_dust, igcm_h2o_ice, igcm_ccn_mass, igcm_ccn_number, &
[2562]	70	nuice_sed, igcm_co2, igcm_co2_ice, igcm_ccnco2_mass, igcm_ccnco2_number, &
	71	igcm_ccnco2_h2o_mass_ice, igcm_ccnco2_h2o_mass_ccn, igcm_ccnco2_h2o_number, nuiceco2_sed, &
[2589]	72	nuiceco2_ref, igcm_ccnco2_meteor_mass, igcm_ccnco2_meteor_number
[2362]	73
	74	use conc_mod, only: mmean
[3902]	75	use tcondco2_mod, only: tcondco2
[2589]	76	use time_phylmdz_mod, only: daysec
[2562]	77	use nucleaco2_mod, only: nucleaco2
[2362]	78	use datafile_mod, only: datadir
[3008]	79	use massflowrateco2_mod, only: massflowrateco2
[2494]	80	use density_co2_ice_mod, only: density_co2_ice
[3008]	81	use microphys_h, only: nbinco2_cld, rad_cldco2, m0co2, mco2
	82	use microphys_h, only: mteta, mtetaco2
[3726]	83	use callkeys_mod, only: co2useh2o, meteo_flux
[2494]	84
[2362]	85	implicit none
	86
	87	!----------------------------------------------------------------------------------------------------------------------!
	88	! VARIABLES DECLARATION
	89	!----------------------------------------------------------------------------------------------------------------------!
	90	! Input arguments:
	91	!-----------------
	92	integer, intent(in) :: &
	93	nq, &! number of tracers
	94	ngrid, &! number of point grid
	95	nlay ! number of layer
	96
	97	real, intent(in) :: &
	98	microtimestep, &! physics time step (s)
	99	pplay(ngrid,nlay), &! mid-layer pressure (Pa)
	100	pplev(ngrid,nlay+1), &! inter-layer pressure (Pa)
	101	pteff(ngrid,nlay), &! temperature at the middle of the layers (K)
	102	sum_subpdt(ngrid,nlay), &! tendency on temperature from previous physical parametrizations
	103	pqeff(ngrid,nlay,nq), &! tracers (kg/kg)
	104	tauscaling(ngrid), &! convertion factor for qdust and Ndust
	105	sum_subpdq(ngrid,nlay,nq) ! tendencies on tracers before condensation (kg/kg.s-1)
	106
	107	real, intent(in) :: &
	108	sigma_iceco2 ! Variance of the co2 ice and CCN distributions
	109
	110	double precision, intent(in) :: &
	111	rb_cldco2(nbinco2_cld+1), & ! boundary values of each rad_cldco2
	112	dev2
	113	!----------------------------------------------------------------------------------------------------------------------!
	114	! Output arguments:
	115	!------------------
	116	real, intent(out) :: &
	117	subpdtcloudco2(ngrid,nlay), &! tendency on tracers due to CO2 condensation (K/s)
	118	subpdqcloudco2(ngrid,nlay,nq) ! tendency on tracers due to CO2 condensation (kg/kg.s-1)
	119	!----------------------------------------------------------------------------------------------------------------------!
	120	! Local:
	121	!-------
	122	!----1) Parameters:
	123	!------------------
	124	integer, parameter :: &
	125	! ---Meteoritic flux input file
	126	nbin_meteor = 100, &! Dimension 1
	127	nlev_meteor = 130, &! Dimension 2
	128	uMeteor = 666, &! File unit
	129	! ---Latent heat computation
	130	l0 = 595594d0, &! coeff from: Azreg-Aïnou (2005)
	131	l1 = 903.111d0, &! Title: "Low-temperature data for carbon dioxide"
	132	l2 = -11.5959d0, &! Pulication: eprint arXiv:1403.4403
	133	l3 = 0.0528288d0, &!
	134	l4 = -0.000103183d0 !
	135
	136	real, parameter :: &
[2660]	137	threshold = 1e-30, & ! limit value considering as zero
	138	threshold_2 = 1e-13 ! limit value considering the value is physical (below this value => computer noise for dble)
[2362]	139
[2660]	140	character(len=50), parameter:: &
[2592]	141	file_meteoritic_flux = 'Meteo_flux_Plane.dat'
[2362]	142	!----------------------------------------------------------------------------------------------------------------------!
	143	!----2) Saved:
	144	!-------------
	145	real, save :: &
	146	sigma_ice ! Variance of the h2o ice and CCN distributions
[2616]	147
	148	!$OMP THREADPRIVATE(sigma_ice)
[2362]	149
	150	double precision, save :: &
[2589]	151	pression_meteor(nlev_meteor), &! pressure from meteoritic flux input file
[2362]	152	meteor(nlev_meteor,nbin_meteor), &! Meteoritic flux read from file uMeteor
	153	dev3 ! 1. / ( sqrt(2.) * sigma_ice )
[2616]	154
	155	!$OMP THREADPRIVATE(meteor,dev3)
[2362]	156
	157	logical, save :: &
	158	firstcall = .true. ! Used to compute saved variables
[2616]	159
	160	!$OMP THREADPRIVATE(firstcall)
	161
[2362]	162	!----------------------------------------------------------------------------------------------------------------------!
	163	!----3) Variables:
	164	!-----------------
	165	integer :: &
	166	ig, &! loop on ngrid
	167	l, &! loop on nlay
	168	i, &! loop on nbinco2
	169	! ---Variables for meteoritic flux
	170	ibin, &! loop on nbin_meteor
[2589]	171	idx_min,&! index of min(diff_pressure)
[2362]	172	read_ok ! file uMeteor iostat
	173
	174	real :: &
[2562]	175	Nccnco2, &
	176	Qccnco2, &
	177	Nccnco2_h2o, &
	178	Qccnco2_h2o, &
[2362]	179	masse(ngrid,nlay), &! mass layer (kg.m-2)
	180	rice(ngrid,nlay), &! water ice mass mean radius (m): used for nucleation of CO2 on ice-coated ccns
	181	zq(ngrid,nlay,nq), &! local value of tracers (kg/kg)
	182	zq0(ngrid,nlay,nq), &! local init value of tracers (kg/kg)
	183	zt(ngrid,nlay), &! local value of temperature (K)
	184	zqsat(ngrid,nlay), &! saturation vapor pressure for CO2 (K)
	185	lw, &! Latent heat of sublimation (J.kg-1)
	186	rdust(ngrid,nlay), &! Dust geometric mean radius (m)
	187	rhocloud(ngrid,nlay), &! Cloud density (kg.m-3)
	188	rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
	189
[4019]	190	double precision :: &
	191	tcond(ngrid,nlay) ! condensation temperature of CO2 (K)
	192
[2362]	193	real(kind=8) :: &
	194	derf ! Error function
	195
	196	double precision :: &
	197	dMice, &! mass of condensed ice
	198	facteurmax, &! for energy limit on mass growth
	199	pco2, &! Co2 vapor partial pressure (Pa)
	200	satu, &! Co2 vapor saturation ratio over ice
[2660]	201	Mo, &! mass of dust particles
	202	No, &! number of dust particles
[2362]	203	Rn, &! logarithm of rdust/rice
	204	Rm, &! Rn * variance of ice and CCN distribution
	205	n_derf, &! derf( (rb_cldco2(1)+Rn) *dev3)
	206	m_derf, &! derf( (rb_cldco2(1)+Rm) *dev2)
[2660]	207	n_aer(nbinco2_cld), &! Radius used by the microphysical scheme (m) for dust
	208	m_aer(nbinco2_cld), &! number concentration V-1 of particle/each size bin for dust
	209	n_aer_meteor(nbinco2_cld), &! Radius used by the microphysical scheme (m) for meteor
	210	m_aer_meteor(nbinco2_cld), &! number concentration V-1 of particle/each size bin for meteor
	211	n_aer_h2oice(nbinco2_cld), &! mass mixing ratio of particle/each size bin for h2o ice
	212	m_aer_h2oice(nbinco2_cld), &! Same - for CO2 nucleation for h2o ice
[2362]	213	Ic_rice, &! Mass transfer rate CO2 ice crystal
	214	ratioh2o_ccn, &! 1./(zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
	215	vo2co2, &! volume of co2 ice particle
	216	dN, &! number of particle of dust used as ccn
	217	dM, &! mass of dN
[2660]	218	dN_meteor, &! number of particle of meteoritic particles used as ccn
	219	dM_meteor, &! mass of dN_meteor
[2362]	220	dNh2o, &! number of particle of h2o ice used as ccn
	221	dMh2o, &! mass of dNh2o
	222	dNN, &! min(dN,zq(ig,l,igcm_dust_number))
	223	dMM, &! min(dM,zq(ig,l,igcm_dust_mass))
	224	dNNh2o, &! min(dNNh2o,zq(ig,l,igcm_ccn_number))
	225	dMh2o_ice, &! min(dMh2ozq(ig,l,igcm_h2o_ice)ratioh2o_ccn, zq(ig,l,igcm_h2o_ice))
	226	dMh2o_ccn, &! min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
	227	rate(nbinco2_cld), &! nucleation rate
[2660]	228	rate_meteor(nbinco2_cld), &! nucleation rate for meteor
[2362]	229	rateh2o(nbinco2_cld), &! nucleation rate for h2o
[2494]	230	rho_ice_co2T, &! density of co2 ice Temperature-dependent
[2362]	231	riceco2(ngrid,nlay), &! CO2 ice mean radius (m)
	232	vrat_cld, &! Volume ratio
[2456]	233	Proba, &! 1.0 - exp(-1.microtimesteprate(i))
	234	Probah2o, &! 1.0 - exp(-1.microtimesteprateh2o(i))
[2660]	235	Proba_meteor, &! 1.0 - exp(-1.microtimesteprate_meteor(i))
[2589]	236	diff_pressure(nlev_meteor), &! abs(pression_meteor(:)-pplev(ig,l))
[2660]	237	meteor_ccn(ngrid,nlay,nbinco2_cld) ! input flux of meteoritc particles
[2362]	238
	239	logical :: &
	240	file_ok ! test if meteoritic input file exists
	241	!======================================================================================================================!
	242	! BEGIN ===============================================================================================================!
	243	!======================================================================================================================!
	244	! 0. Firstcall
	245	!----------------------------------------------------------------------------------------------------------------------!
	246	if (firstcall) then
	247	firstcall = .false.
	248
	249	! Variance of the ice and CCN distributions
	250	sigma_ice = sqrt(log(1.+nuice_sed))
	251	dev3 = 1. / ( sqrt(2.) * sigma_ice )
	252	!----------------------------------------------------------------------------------------------------------------------!
	253	! 0.1. Bonus: meteoritic component, extract data
	254	!----------------------------------------------------------------------------------------------------------------------!
[2616]	255
[2362]	256	if (meteo_flux) then
	257	! Check if file exists
[2660]	258	inquire(file=trim(datadir)//'/'//trim(file_meteoritic_flux), exist=file_ok)
[2362]	259	if (.not. file_ok) then
[2660]	260	call abort_physic("CO2clouds", 'file '//trim(file_meteoritic_flux)//' should be in'//trim(datadir), 1)
[2362]	261	end if
[2660]	262	write(,)'Meteoritic flux file used: ', trim(file_meteoritic_flux)
[2362]	263
	264	! open file
[2660]	265	open(unit=uMeteor,file=trim(datadir)//'/'//trim(file_meteoritic_flux), FORM='formatted')
[2362]	266
	267	!skip 1 line
	268	read(uMeteor,*)
	269
	270	! extract pressure_meteor
	271	do i = 1, nlev_meteor
	272	read(uMeteor,*)pression_meteor(i)
	273	end do
	274
	275	!skip 1 line
	276	read(uMeteor,*)
	277
	278	! extract meteor flux
	279	do i = 1, nlev_meteor
	280	! les mêmes 100 bins size que la distri nuclea : on touche pas
	281	do ibin = 1, nbin_meteor
[2589]	282	read(uMeteor,'(F12.6)')meteor(i,ibin)
[2362]	283	end do
	284	end do
	285
	286	! close file
	287	close(uMeteor)
	288	end if ! of if meteo_flux
	289	end if ! firstcall
	290	!----------------------------------------------------------------------------------------------------------------------!
	291	! 1. Initialization
	292	!----------------------------------------------------------------------------------------------------------------------!
	293	rdust(:,:) = 0.
[2589]	294	meteor_ccn(:,:,:) = 0d0
[2362]	295	rice(:,:) = 1.e-8
	296	riceco2(:,:) = 1.e-11
	297
	298	! Initialize the tendencies
	299	subpdqcloudco2(:,:,:) = 0.
	300	subpdtcloudco2(:,:) = 0.
	301
	302	! pteff temperature layer; sum_subpdt dT.s-1
[2494]	303	zt(1:ngrid,1:nlay) = 0.
[2362]	304	zt(:,:) = pteff(:,:) + sum_subpdt(:,:) * microtimestep
	305
	306	! pqeff traceur kg/kg; sum_subpdq tendance idem .s-1
	307	zq(:,:,:) = pqeff(:,:,:) + sum_subpdq(:,:,:) * microtimestep
	308	where( zq(:,:,:) < threshold ) zq(:,:,:) = threshold
	309
	310	zq0(:,:,:) = zq(:,:,:)
[2494]	311	zqsat(:,:) = 0.
[2362]	312	!----------------------------------------------------------------------------------------------------------------------!
	313	! 2. Compute saturation
	314	!----------------------------------------------------------------------------------------------------------------------!
	315	call co2sat(ngrid*nlay,zt,zqsat)
	316	call tcondco2(ngrid,nlay,pplay, zq(:,:,igcm_co2) + zq(:,:,igcm_co2_ice),tcond)
	317	!----------------------------------------------------------------------------------------------------------------------!
	318	! 3. Bonus: additional meteoritic particles for nucleation
	319	!----------------------------------------------------------------------------------------------------------------------!
	320	if (meteo_flux) then
	321	do l = 1, nlay
	322	do ig = 1, ngrid
	323	masse(ig,l) = (pplev(ig,l) - pplev(ig,l+1)) / g
	324
[2589]	325	diff_pressure(:) = abs(pression_meteor(:)-pplev(ig,l))
[2362]	326
[2589]	327	idx_min = minloc(diff_pressure, DIM=1)
[2362]	328	! Par kg air csi par m carre, x epaisseur/masse pour par kg/air. Check original unit with J. Plane
[2589]	329	meteor_ccn(ig,l,:) = meteor(idx_min,:)/masse(ig,l)
[2362]	330	end do
	331	end do
	332	end if
	333	!----------------------------------------------------------------------------------------------------------------------!
	334	! 4. Actual microphysics: Main loop over the GCM's grid
	335	!----------------------------------------------------------------------------------------------------------------------!
	336	do l = 1, nlay
	337	do ig = 1, ngrid
	338
	339	! Get the partial pressure of co2 vapor and its saturation ratio
	340	pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/(mco21e3)) pplay(ig,l) ! mco2 (kg/mol) => mmean (g/mol)
	341	satu = pco2 / zqsat(ig,l)
	342
	343	!T-dependant CO2 ice density
[2494]	344	call density_co2_ice(dble(zt(ig,l)), rho_ice_co2T)
[2362]	345
[2494]	346	vo2co2 = m0co2 / rho_ice_co2T
[2362]	347	!----------------------------------------------------------------------------------------------------------------------!
	348	! 4.1 Nucleation
	349	!----------------------------------------------------------------------------------------------------------------------!
	350	! if there is condensation
	351	if ( satu >= 1 ) then
	352	call updaterccnCO2(zq(ig,l,igcm_dust_mass), zq(ig,l,igcm_dust_number), rdust(ig,l), tauscaling(ig))
	353
	354	! Expand the dust moments into a binned distribution
	355	n_aer(:) = 0d0 ! number of aerosol particles
	356	m_aer(:) = 0d0 ! mass of aerosol particles
	357
	358	No = zq(ig,l,igcm_dust_number) * tauscaling(ig)
	359
[2456]	360	Mo = (4./3.) * pi * rho_dust * No * rdust(ig,l)*3 exp(9.*nuiceco2_ref/2.)
[2362]	361
	362	if (No > threshold) then
[2456]	363	Rn = -log(rdust(ig,l))
[2362]	364
	365	Rm = Rn - 3. * sigma_iceco2 * sigma_iceco2
	366
	367	n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
	368	m_derf = derf( (rb_cldco2(1)+Rm) *dev2)
	369
	370	do i = 1, nbinco2_cld
	371	n_aer(i) = -0.5 * No * n_derf
	372	m_aer(i) = -0.5 * Mo * m_derf
	373
	374	n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
	375	m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
	376
	377	n_aer(i) = n_aer(i) + 0.5 * No * n_derf
	378	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
	379	end do
	380
[2562]	381	! Call to nucleation routine
	382	call nucleaco2(dble(pco2), zt(ig,l), dble(satu), n_aer, rate, vo2co2, mtetaco2)
	383	dN = 0.
	384	dM = 0.
	385
	386	do i = 1, nbinco2_cld
	387	Proba = 1.0 - exp(-1.microtimesteprate(i))
	388	dN = dN + n_aer(i) * Proba
	389	dM = dM + m_aer(i) * Proba
	390	end do
	391
	392	! Now increment CCN tracers and update dust tracers
[2589]	393	dNN = min(dN, zq(ig,l,igcm_dust_number)) ! dNN est devenu DN
	394	dMM = min(dM, zq(ig,l,igcm_dust_mass)) ! idem dans le min
[2562]	395
	396	zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dMM /tauscaling(ig)
	397	zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dNN /tauscaling(ig)
	398
	399	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - dMM /tauscaling(ig)
	400	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - dNN /tauscaling(ig)
[2589]	401	end if ! of if No > 1e-30
[2562]	402
[2589]	403	! Ajout meteor_ccn particles aux particules de poussière background
	404	if (meteo_flux) then
	405	n_aer_meteor(1:nbinco2_cld) = 0d0
	406	m_aer_meteor(1:nbinco2_cld) = 0d0
[2362]	407	do i = 1, nbinco2_cld
[2589]	408	n_aer_meteor(i) = meteor_ccn(ig,l,i)
	409	m_aer_meteor(i) = (4./3.) * pi * rho_dust * meteor_ccn(ig,l,i) * rad_cldco2(i)**3
	410	end do
	411	! Call to nucleation routine
	412	rate_meteor(1:nbinco2_cld) = 0d0
	413	call nucleaco2(dble(pco2), zt(ig,l), dble(satu), n_aer_meteor, rate_meteor, vo2co2, mtetaco2)
[2362]	414
[2589]	415	dN_meteor = 0.
	416	dM_meteor = 0.
	417	do i = 1, nbinco2_cld
	418	Proba_meteor = 1.0 - exp(-1.microtimesteprate_meteor(i))
	419	dN_meteor = dN_meteor + n_aer_meteor(i) * Proba_meteor
	420	dM_meteor = dM_meteor + m_aer_meteor(i) * Proba_meteor
	421	end do
	422	! Now increment CCN tracers and update dust tracers
	423	zq(ig,l,igcm_ccnco2_meteor_mass) = zq(ig,l,igcm_ccnco2_meteor_mass) + dM_meteor
	424	zq(ig,l,igcm_ccnco2_meteor_number) = zq(ig,l,igcm_ccnco2_meteor_number) + dN_meteor
	425
	426	zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dM_meteor
	427	zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dN_meteor
[2362]	428	end if
	429
[2562]	430	! Same but with h2o particles as CCN only if co2useh2o = .true.
	431	if (co2useh2o) then
	432	call updaterice_micro(zq(ig,l,igcm_h2o_ice), zq(ig,l,igcm_ccn_mass), zq(ig,l,igcm_ccn_number), &
	433	tauscaling(ig), rice(ig,l), rhocloud(ig,l))
[2362]	434
[2562]	435	! Total mass of H20 crystals,CCN included
	436	Mo = zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccn_mass) * tauscaling(ig) + threshold
[2362]	437
[2562]	438	No = zq(ig,l,igcm_ccn_number) * tauscaling(ig) + threshold
	439	if (No > threshold) then
[2456]	440	Rn = -log(rice(ig,l))
[2362]	441
	442	Rm = Rn - 3. * sigma_ice * sigma_ice
	443
	444	n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
	445	m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
	446
[2562]	447	n_aer_h2oice(:)=0.
	448	m_aer_h2oice(:)=0.
[2362]	449	do i = 1, nbinco2_cld
	450	n_aer_h2oice(i) = -0.5 * No * n_derf
	451	m_aer_h2oice(i) = -0.5 * Mo * m_derf
	452
	453	n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
	454	m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
	455
	456	n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf
	457	m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf
	458	end do
	459
[2562]	460	call nucleaco2(dble(pco2), zt(ig,l), dble(satu), n_aer_h2oice, rateh2o, vo2co2, mteta)
	461	dNh2o = 0.
	462	dMh2o = 0.
[2494]	463
[2362]	464	do i = 1, nbinco2_cld
[2456]	465	Probah2o = 1.0 - exp(-1.microtimesteprateh2o(i))
[2362]	466	dNh2o = dNh2o + n_aer_h2oice(i) * Probah2o
	467	dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o
	468	end do
	469
[2562]	470	! Update CCN for CO2 nucleating on H2O CCN : Warning: must keep memory of it
[2362]	471	dNNh2o = dNh2o/tauscaling(ig)
	472	dNNh2o = min(dNNh2o,zq(ig,l,igcm_ccn_number))
	473
	474	ratioh2o_ccn = 1./(zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
	475
	476	dMh2o_ccn = dMh2o * zq(ig,l,igcm_ccn_mass) * tauscaling(ig) * ratioh2o_ccn
	477	dMh2o_ccn = dMh2o_ccn/tauscaling(ig)
	478	dMh2o_ccn = min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
	479
	480	dMh2o_ice = dMh2o * zq(ig,l,igcm_h2o_ice) * ratioh2o_ccn
	481	dMh2o_ice = min(dMh2o_ice,zq(ig,l,igcm_h2o_ice))
	482
[2562]	483	zq(ig,l,igcm_ccnco2_h2o_mass_ice) = zq(ig,l,igcm_ccnco2_h2o_mass_ice) + dMh2o_ice
	484	zq(ig,l,igcm_ccnco2_h2o_mass_ccn) = zq(ig,l,igcm_ccnco2_h2o_mass_ccn) + dMh2o_ccn
	485	zq(ig,l,igcm_ccnco2_h2o_number) = zq(ig,l,igcm_ccnco2_h2o_number) + dNNh2o
[2362]	486
[2562]	487	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) - dNNh2o
[2362]	488	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) - dMh2o_ice
	489	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) - dMh2o_ccn
	490
[2562]	491	zq(ig,l,igcm_ccnco2_mass) = zq(ig,l,igcm_ccnco2_mass) + dMh2o_ice + dMh2o_ccn
	492	zq(ig,l,igcm_ccnco2_number) = zq(ig,l,igcm_ccnco2_number) + dNNh2o
	493	end if
	494	end if ! of if co2useh2o
[2362]	495	end if ! of is satu > 1
	496	!----------------------------------------------------------------------------------------------------------------------!
	497	! 4.2. Ice growth: scheme for radius evolution
	498	!----------------------------------------------------------------------------------------------------------------------!
	499	! We trigger crystal growth if and only if there is at least one nuclei (N>1). Indeed, if we are supersaturated
	500	! and still don't have at least one nuclei, we should better wait to avoid unrealistic value for nuclei radius
	501	! and so on for cases that remain negligible.
	502	!----------------------------------------------------------------------------------------------------------------------!
	503	! we trigger crystal growth
[2562]	504	if ((zq(ig,l,igcm_ccnco2_number)) * tauscaling(ig) + threshold >= 1) then
	505	Nccnco2 = dble(zq(ig,l,igcm_ccnco2_number))
	506	Qccnco2 = dble(zq(ig,l,igcm_ccnco2_mass))
	507	Nccnco2_h2o = 0.
	508	Qccnco2_h2o = 0.
[2362]	509
[2562]	510	if (co2useh2o) then
	511	Nccnco2_h2o = zq(ig,l,igcm_ccnco2_h2o_number)
	512	Qccnco2_h2o = zq(ig,l,igcm_ccnco2_h2o_mass_ice) + zq(ig,l,igcm_ccnco2_h2o_mass_ccn)
	513	Nccnco2 = Nccnco2 - Nccnco2_h2o
	514	Qccnco2 = Qccnco2 - Qccnco2_h2o
[2660]	515	if (Nccnco2 <= 0.) then
	516	Nccnco2 = threshold
	517	Qccnco2 = threshold
	518	end if
[2562]	519	end if
	520
	521	call updaterice_microco2(dble(zq(ig,l,igcm_co2_ice)), dble(Qccnco2), dble(Nccnco2), dble(Qccnco2_h2o), &
	522	dble(Nccnco2_h2o), zt(ig,l), tauscaling(ig), riceco2(ig,l), rhocloudco2(ig,l))
[2362]	523	Ic_rice = 0.
	524
	525	! J.kg-1
	526	lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)*2 + l3 zt(ig,l)*3 + l4 zt(ig,l)**4
	527
	528	facteurmax = abs(tcond(ig,l)-zt(ig,l)) * (cpp/lw)
	529
[2388]	530	! call scheme of microphys. mass growth for CO2 (evaporation/condensation)
[2362]	531	call massflowrateCO2(pplay(ig,l), zt(ig,l), satu, riceco2(ig,l), mmean(ig,l), Ic_rice)
	532
	533	! Ic_rice Mass transfer rate (kg/s) for a rice particle > 0 si croissance !
	534	if (isnan(Ic_rice) .or. Ic_rice == 0.) then
	535	Ic_rice = 0.
	536	subpdtcloudco2(ig,l) = -sum_subpdt(ig,l)
	537	dMice = 0
	538	else
	539	! Kg par kg d'air, >0 si croissance !
	540	! kg.s-1 par particule * nb particule par kg air*s = kg par kg air
	541	dMice = zq(ig,l,igcm_ccnco2_number) * Ic_rice * microtimestep * tauscaling(ig)
	542
	543	! facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
	544	! latent heat release > 0 if growth i.e. if dMice > 0
	545	dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
	546	dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
	547
	548	! kgco2/kgair* J/kgco2 * 1/(J.kgair-1.K-1)/s = K /s
	549	subpdtcloudco2(ig,l) = dMice * lw / cpp / microtimestep
	550
	551	!Now update tracers
	552	zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice) + dMice
	553	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) - dMice
	554	end if
[2394]	555	end if ! if zq(ccnco2_number) >= 1
[2362]	556	!----------------------------------------------------------------------------------------------------------------------!
	557	! 4.3 Dust cores releasing if no more co2 ice
	558	!----------------------------------------------------------------------------------------------------------------------!
[2394]	559	! On sublime tout
[2660]	560	if ((zq(ig,l,igcm_co2_ice) < threshold_2).or.(zq(ig,l,igcm_ccnco2_number) * tauscaling(ig) < 1.)) then
[2562]	561	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) + zq(ig,l,igcm_ccnco2_mass)
	562	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) + zq(ig,l,igcm_ccnco2_number)
	563	zq(ig,l,igcm_ccnco2_mass) = 0.
	564	zq(ig,l,igcm_ccnco2_number) = 0.
[2589]	565	if (meteo_flux) then
	566	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - zq(ig,l,igcm_ccnco2_meteor_mass)
	567	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - zq(ig,l,igcm_ccnco2_meteor_number)
	568	zq(ig,l,igcm_ccnco2_meteor_mass) = 0.
	569	zq(ig,l,igcm_ccnco2_meteor_number) = 0.
	570	end if
[2562]	571	if (co2useh2o) then
	572	zq(ig,l,igcm_dust_mass) = zq(ig,l,igcm_dust_mass) - zq(ig,l,igcm_ccnco2_h2o_mass_ccn) - &
	573	zq(ig,l,igcm_ccnco2_h2o_mass_ice)
	574	zq(ig,l,igcm_dust_number) = zq(ig,l,igcm_dust_number) - zq(ig,l,igcm_ccnco2_h2o_number)
[2362]	575
[2562]	576	zq(ig,l,igcm_ccn_mass) = zq(ig,l,igcm_ccn_mass) + zq(ig,l,igcm_ccnco2_h2o_mass_ccn)
	577	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice) + zq(ig,l,igcm_ccnco2_h2o_mass_ice)
	578	zq(ig,l,igcm_ccn_number) = zq(ig,l,igcm_ccn_number) + zq(ig,l,igcm_ccnco2_h2o_number)
	579	zq(ig,l,igcm_ccnco2_h2o_number) = 0.
	580	zq(ig,l,igcm_ccnco2_h2o_mass_ccn) = 0.
	581	zq(ig,l,igcm_ccnco2_h2o_mass_ice) = 0.
[2394]	582	end if
[2562]	583	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) + zq(ig,l,igcm_co2_ice)
	584	zq(ig,l,igcm_co2_ice) = 0.
	585	riceco2(ig,l) = 0.
[2394]	586	end if !of if co2_ice < threshold or zq(ccnco2_number) < 1
[2362]	587	end do ! of ig loop
	588	end do ! of nlayer loop
	589	!----------------------------------------------------------------------------------------------------------------------!
	590	! 5. Get cloud tendencies
	591	!----------------------------------------------------------------------------------------------------------------------!
	592	subpdqcloudco2(:,:,igcm_co2) = ( zq(:,:,igcm_co2) - zq0(:,:,igcm_co2) ) / microtimestep
	593
	594	subpdqcloudco2(:,:,igcm_co2_ice) = ( zq(:,:,igcm_co2_ice) - zq0(:,:,igcm_co2_ice) ) / microtimestep
	595
	596	subpdqcloudco2(:,:,igcm_ccnco2_mass) = ( zq(:,:,igcm_ccnco2_mass) - zq0(:,:,igcm_ccnco2_mass) ) / microtimestep
	597
	598	subpdqcloudco2(:,:,igcm_ccnco2_number) = ( zq(:,:,igcm_ccnco2_number) - zq0(:,:,igcm_ccnco2_number))/microtimestep
	599
	600	subpdqcloudco2(:,:,igcm_dust_mass) = ( zq(:,:,igcm_dust_mass) - zq0(:,:,igcm_dust_mass) ) / microtimestep
	601
	602	subpdqcloudco2(:,:,igcm_dust_number) = ( zq(:,:,igcm_dust_number) - zq0(:,:,igcm_dust_number) ) / microtimestep
	603
[2589]	604	if (meteo_flux) then
	605	subpdqcloudco2(:,:,igcm_ccnco2_meteor_mass) = (zq(:,:,igcm_ccnco2_meteor_mass)-zq0(:,:,igcm_ccnco2_meteor_mass) &
	606	)/microtimestep
	607
	608	subpdqcloudco2(:,:,igcm_ccnco2_meteor_number) = ( zq(:,:,igcm_ccnco2_meteor_number) - &
	609	zq0(:,:,igcm_ccnco2_meteor_number) )/microtimestep
	610	end if
	611
[2362]	612	if (co2useh2o) then
	613	subpdqcloudco2(:,:,igcm_h2o_ice) = ( zq(:,:,igcm_h2o_ice) - zq0(:,:,igcm_h2o_ice) ) / microtimestep
	614
	615	subpdqcloudco2(:,:,igcm_ccn_mass) = ( zq(:,:,igcm_ccn_mass) - zq0(:,:,igcm_ccn_mass) ) / microtimestep
	616
	617	subpdqcloudco2(:,:,igcm_ccn_number) = ( zq(:,:,igcm_ccn_number) - zq0(:,:,igcm_ccn_number) ) / microtimestep
[2562]	618
	619	subpdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ice) = (zq(:,:,igcm_ccnco2_h2o_mass_ice)-zq0(:,:,igcm_ccnco2_h2o_mass_ice)&
	620	)/microtimestep
	621
	622	subpdqcloudco2(:,:,igcm_ccnco2_h2o_mass_ccn) = (zq(:,:,igcm_ccnco2_h2o_mass_ccn)-zq0(:,:,igcm_ccnco2_h2o_mass_ccn)&
	623	)/microtimestep
	624
	625	subpdqcloudco2(:,:,igcm_ccnco2_h2o_number) = ( zq(:,:,igcm_ccnco2_h2o_number) - zq0(:,:,igcm_ccnco2_h2o_number) &
	626	)/microtimestep
[2362]	627	end if
	628	!======================================================================================================================!
	629	! END =================================================================================================================!
	630	!======================================================================================================================!
	631	end subroutine improvedCO2clouds
	632
	633	end module improvedCO2clouds_mod
	634

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

Download in other formats:

Original Format