Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

improvedco2clouds_mod.F90

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

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

Download in other formats:

Original Format