Context Navigation

improvedco2clouds_mod.F90 @ 2827

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

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

Download in other formats:

Original Format