Context Navigation

improvedCO2clouds.F @ 1775

Last change on this file since 1775 was 1720, checked in by jaudouard, 8 years ago
Update on CO2 ice clouds scheme
Property svn:executable set to ``*
File size: 35.1 KB

Rev	Line
[1617]	1	subroutine improvedCO2clouds(ngrid,nlay,ptimestep,
[1720]	2	& pplay,pzlev,pt,pdt,
[1617]	3	& pq,pdq,pdqcloudco2,pdtcloudco2,
[1685]	4	& nq,tauscaling,
	5	& memdMMccn,memdMMh2o,memdNNccn)
[1617]	6	! to use 'getin'
	7	USE comcstfi_h
	8	USE ioipsl_getincom
	9	USE updaterad
	10	use tracer_mod
	11	!, only: rho_ice_co2, nuiceco2_sed, igcm_co2,
	12	! & rho_ice,igcm_h2o_ice, igcm_ccn_number,
	13	! & igcm_co2_ice, igcm_dust_mass,
	14	! & igcm_dust_number, igcm_ccnco2_mass,
	15	! & igcm_ccnco2_number
	16	use conc_mod, only: mmean
[1720]	17
[1617]	18	implicit none
	19
	20	c------------------------------------------------------------------
	21	c This routine is used to form CO2 clouds when a parcel of the GCM is
	22	c saturated. It includes the ability to have supersaturation, a
	23	c computation of the nucleation rates, growthrates and the
	24	c scavenging of dust particles by clouds.
	25	c It is worth noting that the amount of dust is computed using the
	26	c dust optical depth computed in aeropacity.F. That's why
	27	c the variable called "tauscaling" is used to convert
	28	c pq(dust_mass) and pq(dust_number), which are relative
	29	c quantities, to absolute and realistic quantities stored in zq.
	30	c This has to be done to convert the inputs into absolute
	31	c values, but also to convert the outputs back into relative
	32	c values which are then used by the sedimentation and advection
	33	c schemes.
[1685]	34	c CO2 ice particles can nucleate on both dust and on water ice particles
	35	c When CO2 ice is deposited onto a water ice particles, the particle is
	36	c removed from the water tracers.
	37	cWARNING: no sedimentation of the water ice origin is performed
	38	c in the microphysical timestep in co2cloud.F.
[1617]	39
	40	c Authors of the water ice clouds microphysics
	41	c J.-B. Madeleine, based on the work by Franck Montmessin
	42	c (October 2011)
	43	c T. Navarro, debug,correction, new scheme (October-April 2011)
	44	c A. Spiga, optimization (February 2012)
	45	c Adaptation for CO2 clouds by Joachim Audouard (09/16), based on the work
	46	c of Constantino Listowski
	47	c------------------------------------------------------------------
	48	!#include "dimensions.h"
	49	!#include "dimphys.h"
	50	#include "callkeys.h"
	51	!#include "tracer.h"
	52	!#include "comgeomfi.h"
	53	!#include "dimradmars.h"
	54	#include "microphys.h"
[1720]	55	#include "datafile.h"
[1617]	56	!#include "microphysCO2.h"
	57	!#include "conc.h"
	58	c------------------------------------------------------------------
	59	c Inputs:
	60
	61	INTEGER ngrid,nlay
	62	integer nq ! nombre de traceurs
	63	REAL ptimestep ! pas de temps physique (s)
	64	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
[1720]	65	REAL pzlev(ngrid,nlay) ! altitude au milieu des couches ()
	66
[1617]	67	REAL pt(ngrid,nlay) ! temperature at the middle of the
	68	! layers (K)
	69	REAL pdt(ngrid,nlay) ! tendance temperature des autres
	70	! param.
	71	REAL pq(ngrid,nlay,nq) ! traceur (kg/kg)
	72	REAL pdq(ngrid,nlay,nq) ! tendance avant condensation
	73	! (kg/kg.s-1)
	74	REAL tauscaling(ngrid) ! Convertion factor for qdust and Ndust
	75
[1685]	76	REAL rice(ngrid,nlay) ! Water Ice mass mean radius (m)
[1617]	77	! used for nucleation of CO2 on ice-coated ccns
[1720]	78	REAL rccnh2o(ngrid,nlay) ! Water Ice mass mean radius (m)
[1617]	79
	80	c Outputs:
	81	REAL pdqcloudco2(ngrid,nlay,nq) ! tendance de la condensation
	82	! CO2 (kg/kg.s-1)
	83	! condensation si igcm_co2_ice
	84	REAL pdtcloudco2(ngrid,nlay) ! tendance temperature due
	85	! a la chaleur latente
	86
	87	c------------------------------------------------------------------
	88	c Local variables:
	89
	90	LOGICAL firstcall
	91	DATA firstcall/.true./
	92	SAVE firstcall
	93
	94	REAL*8 derf ! Error function
	95	!external derf
	96
	97	!REAL*8 massflowrateCO2
	98	!external massflowrateCO2
	99
[1720]	100	INTEGER ig,l,i,flag_pourri
[1617]	101
	102	REAL zq(ngrid,nlay,nq) ! local value of tracers
	103	REAL zq0(ngrid,nlay,nq) ! local initial value of tracers
	104	REAL zt(ngrid,nlay) ! local value of temperature
	105	REAL zqsat(ngrid,nlay) ! saturation vapor pressure for CO2
	106	REAL lw !Latent heat of sublimation (J.kg-1)
[1720]	107	REAL,save :: l0,l1,l2,l3,l4
[1617]	108	REAL cste
[1685]	109	DOUBLE PRECISION dMice ! mass of condensed ice
[1617]	110	DOUBLE PRECISION sumcheck
	111	DOUBLE PRECISION pco2 ! Co2 vapor partial pressure (Pa)
	112	DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice
	113	DOUBLE PRECISION Mo,No
[1720]	114	DOUBLE PRECISION Rn, Rm, dev2,dev3, n_derf, m_derf
	115	DOUBLE PRECISION memdMMccn(ngrid,nlay)
	116	DOUBLE PRECISION memdMMh2o(ngrid,nlay)
[1685]	117	DOUBLE PRECISION memdNNccn(ngrid,nlay)
	118
[1617]	119	! Radius used by the microphysical scheme (m)
	120	DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin
	121	DOUBLE PRECISION m_aer(nbinco2_cld) ! mass mixing ratio of particle/each size bin
[1720]	122	DOUBLE PRECISION m_aer2(nbinco2_cld) ! mass mixing ratio of particle/each size bin
	123	DOUBLE PRECISION m_aer_h2oice2(nbinco2_cld) ! mass mixing ratio of particle/each size bin
[1617]	124
	125	DOUBLE PRECISION n_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
	126	DOUBLE PRECISION m_aer_h2oice(nbinco2_cld) ! Same - for CO2 nucleation
	127	DOUBLE PRECISION rad_h2oice(nbinco2_cld)
	128
	129	c REAL*8 sigco2 ! Co2-ice/air surface tension (N.m)
	130	c EXTERNAL sigco2
	131
[1685]	132	DOUBLE PRECISION dN,dM, dNh2o, dMh2o, dNN,dMM,dNNh2o,dMMh2o
[1720]	133	DOUBLE PRECISION dMh2o_ice,dMh2o_ccn
	134
[1617]	135	DOUBLE PRECISION rate(nbinco2_cld) ! nucleation rate
	136	DOUBLE PRECISION rateh2o(nbinco2_cld) ! nucleation rate
	137	REAL seq
[1685]	138	DOUBLE PRECISION rho_ice_co2T(ngrid,nlay)
[1617]	139	DOUBLE PRECISION riceco2(ngrid,nlay) ! CO2Ice mean radius (m)
[1685]	140	REAL rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
	141
[1617]	142	REAL rhocloudco2(ngrid,nlay) ! Cloud density (kg.m-3)
	143	REAL rdust(ngrid,nlay) ! Dust geometric mean radius (m)
	144
	145	c REAL res ! Resistance growth
	146	DOUBLE PRECISION Ic_rice ! Mass transfer rate CO2 ice crystal
[1685]	147	DOUBLE PRECISION ratioh2o_ccn
	148	DOUBLE PRECISION vo2co2
[1617]	149	c Parameters of the size discretization
	150	c used by the microphysical scheme
[1720]	151	DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-11 ! Minimum radius (m)
	152	DOUBLE PRECISION, PARAMETER :: rmax_cld = 3.e-6 ! Maximum radius (m)
	153	DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-12
	154	! Minimum bounary radius (m)
	155	DOUBLE PRECISION, PARAMETER :: rbmax_cld = 5.e-6 ! Maximum boundary radius (m)
[1617]	156	DOUBLE PRECISION vrat_cld ! Volume ratio
	157	DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
	158	SAVE rb_cldco2
	159
	160	DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m)
	161	DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3)
	162
[1685]	163	DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o
[1617]	164	REAL sigma_iceco2 ! Variance of the ice and CCN distributions
	165	SAVE sigma_iceco2
[1720]	166
[1617]	167
[1720]	168	REAL sigma_ice ! Variance of the ice and CCN distributions
	169	SAVE sigma_ice
	170
[1632]	171	DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2
[1617]	172
[1720]	173	integer coeffh2o
	174	!Variables for the meteoritic flux
	175
	176	double precision meteo_ccn(ngrid,nlay,100) !100=nbinco2_cld !!!
	177	double precision,save :: meteo(130,100)
	178	double precision mtemp(100)
	179	logical file_ok
	180	integer altitudes_meteo(130),nelem,lebon1,lebon2
	181	double precision :: ltemp1(130),ltemp2(130)
	182	integer ibin,uMeteo,j
[1617]	183	c----------------------------------
	184	c TESTS
	185
	186
	187
	188	REAL satubf(ngrid,nlay),satuaf(ngrid,nlay)
	189	REAL res_out(ngrid,nlay)
	190
	191
	192	c------------------------------------------------------------------
	193
	194	IF (firstcall) THEN
	195	!=============================================================
	196	! 0. Definition of the size grid
	197	!=============================================================
	198	c rad_cldco2 is the primary radius grid used for microphysics computation.
	199	c The grid spacing is computed assuming a constant volume ratio
	200	c between two consecutive bins; i.e. vrat_cld.
	201	c vrat_cld is determined from the boundary values of the size grid:
	202	c rmin_cld and rmax_cld.
	203	c The rb_cldco2 array contains the boundary values of each rad_cldco2 bin.
	204	c dr_cld is the width of each rad_cldco2 bin.
	205
	206	c Volume ratio between two adjacent bins
	207	! vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
	208	! vrat_cld = exp(vrat_cld)
	209	vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
	210	vrat_cld = exp(vrat_cld)
	211	c write(,) "vrat_cld", vrat_cld
	212
	213	rb_cldco2(1) = rbmin_cld
	214	rad_cldco2(1) = rmin_cld
	215	vol_cld(1) = 4./3. * dble(pi) * rmin_cldrmin_cldrmin_cld
	216	! vol_cld(1) = 4./3. * pi * rmin_cldrmin_cldrmin_cld
	217
	218	do i=1,nbinco2_cld-1
	219	rad_cldco2(i+1) = rad_cldco2(i) * vrat_cld**(1./3.)
	220	vol_cld(i+1) = vol_cld(i) * vrat_cld
	221	enddo
	222
	223	do i=1,nbinco2_cld
	224	rb_cldco2(i+1)= ( (2.vrat_cld) / (vrat_cld+1.) )(1./3.)
	225	& rad_cldco2(i)
	226	dr_cld(i) = rb_cldco2(i+1) - rb_cldco2(i)
	227	enddo
	228	rb_cldco2(nbinco2_cld+1) = rbmax_cld
	229	dr_cld(nbinco2_cld) = rb_cldco2(nbinco2_cld+1) -
	230	& rb_cldco2(nbinco2_cld)
	231
	232	print*, ' '
	233	print*,'Microphysics co2: size bin information:'
	234	print*,'i,rb_cldco2(i), rad_cldco2(i),dr_cld(i)'
	235	print*,'-----------------------------------'
	236	do i=1,nbinco2_cld
	237	write(*,'(i3,3x,3(e12.6,4x))') i,rb_cldco2(i), rad_cldco2(i),
	238	& dr_cld(i)
	239	enddo
	240	write(*,'(i3,3x,e12.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1)
	241	print*,'-----------------------------------'
	242
	243	do i=1,nbinco2_cld+1
	244	rb_cldco2(i) = log(rb_cldco2(i)) !! we save that so that it is not computed
	245	!! at each timestep and gridpoint
	246	enddo
	247	c Contact parameter of co2 ice on dst ( m=cos(theta) )
	248	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	249	c mteta = 0.952
	250	c mtetaco2 = 0.952
[1685]	251	write(,) 'co2_param contact parameter:', mtetaco2
[1617]	252
	253	c Volume of a co2 molecule (m3)
	254	vo1 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
	255	vo1co2=vo1 ! AJOUT JA
	256	c Variance of the ice and CCN distributions
	257	sigma_iceco2 = sqrt(log(1.+nuiceco2_sed))
[1720]	258	sigma_ice = sqrt(log(1.+nuice_sed))
	259
[1617]	260
	261	c write(,) 'Variance of ice & CCN distribs :', sigma_iceco2
	262	c write(,) 'nuice for sedimentation:', nuiceco2_sed
	263	c write(,) 'Volume of a co2 molecule:', vo1co2
[1685]	264
	265	write(,) 'Variance of ice & CCN CO2 distribs :', sigma_iceco2
	266	write(,) 'nuice for co2 ice sedimentation:', nuiceco2_sed
	267	write(,) 'Volume of a co2 molecule:', vo1co2
[1720]	268
	269	coeffh2o=0
	270	if (co2useh2o) then
	271	write(,)
	272	write(,) "co2useh2o=.true. in callphys.def"
	273	write(,) "This means water ice particles can "
	274	write(,) "serve as CCN for CO2 microphysics"
	275	coeffh2o=1
	276	endif
	277	meteo_ccn(:,:,:)=0.
[1617]	278
[1720]	279	if (meteo_flux) then
	280	write(,)
	281	write(,) "meteo_flux=.true. in callphys.def"
	282	write(,) "meteoritic dust particles are available"
	283	write(,) "for co2 ice nucleation! "
	284	write(,) "Flux given by J. Plane (altitude,size bins)"
	285	! Initialisation of the flux: it is constant and is it saved
	286	!We must interpolate the table to the GCM altitudes
	287	INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))//
	288	& '/Meteo_flux_Plane.dat', EXIST=file_ok)
	289	IF (.not. file_ok) THEN
	290	write(,) 'file Meteo_flux_Plane.dat should be in '
	291	& ,datafile
	292	STOP
	293	endif
	294	!used Variables
	295	open(newunit=uMeteo,file=trim(datafile)//
	296	& '/Meteo_flux_Plane.dat'
	297	& ,FORM='formatted')
	298	!13000 records (130 altitudes x 100 bin sizes)
	299	do i=1,130
	300	do j=1,100 ! les mêmes 100 bins size que la distri nuclea: on touche pas
	301	read(uMeteo,'(F12.6)') meteo(i,j)
	302	enddo
	303	altitudes_meteo(i)=i ! On doit maintenant réinterpoler le tableau(130,100) sur les altitudes du GCM (nlay,100)
	304	enddo
	305	close(uMeteo)
	306	endif !of if meteo_flux
[1617]	307
[1720]	308	l0=595594d0
	309	l1=903.111d0
	310	l2=-11.5959d0
	311	l3=0.0528288d0
	312	l4=-0.000103183d0
[1617]	313	firstcall=.false.
[1720]	314
	315
[1617]	316	END IF
[1720]	317	c write(,) "max memdNN =",maxval(memdNNccn)
[1617]	318	!=============================================================
	319	! 1. Initialisation
	320	!=============================================================
	321	!cste = 4pirho_ice*ptimestep !not used for co2
[1720]	322	flag_pourri=0
[1617]	323
	324	res_out(:,:) = 0
[1685]	325	rice(:,:) = 1.e-8
[1617]	326	riceco2(:,:) = 1.e-11
	327
	328	c Initialize the tendencies
	329	pdqcloudco2(1:ngrid,1:nlay,1:nq)=0.
	330	pdtcloudco2(1:ngrid,1:nlay)=0.
	331
	332	c pt temperature layer; pdt dT.s-1
	333	c pq traceur kg/kg; pdq tendance idem .s-1
	334	zt(1:ngrid,1:nlay) =
	335	& pt(1:ngrid,1:nlay) +
	336	& pdt(1:ngrid,1:nlay) * ptimestep
[1685]	337	c call WRITEDIAGFI(ngrid,"Ztclouds",
	338	c & "Ztclouds",'K',3,zt)
	339	c call WRITEDIAGFI(ngrid,"pdtclouds",
	340	c & "pdtclouds",'K',3,pdt*ptimestep)
	341
[1617]	342	zq(1:ngrid,1:nlay,1:nq) =
	343	& pq(1:ngrid,1:nlay,1:nq) +
	344	& pdq(1:ngrid,1:nlay,1:nq) * ptimestep
	345
	346
	347	WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
	348	& zq(1:ngrid,1:nlay,1:nq) = 1.e-30
	349
	350	zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
	351
	352	!=============================================================
	353	! 2. Compute saturation
	354	!=============================================================
	355
[1720]	356
[1617]	357
	358	dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
[1720]	359	dev3 = 1. / ( sqrt(2.) * sigma_ice )
[1617]	360
	361	call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2)
	362
	363
[1720]	364	!=============================================================
	365	! Bonus: additional meteoritic particles for nucleation
	366	!=============================================================
	367	if (meteo_flux) then
	368	!altitude_meteo(130)
	369	!pzlev(ngrid,nlay)
	370	!meteo(130,100)
	371	!resultat: meteo_ccn(ngrid,nlay,100)
	372	do l=1,nlay
	373	do ig=1,ngrid
	374	ltemp1=abs(altitudes_meteo(:)-pzlev(ig,l))
	375	ltemp2=abs(altitudes_meteo(:)-pzlev(ig,l+1))
	376	lebon1=minloc(ltemp1,DIM=1)
	377	lebon2=minloc(ltemp2,DIM=1)
	378	nelem=lebon2-lebon1+1.
	379	mtemp(:)=0d0 !mtemp(100) : valeurs pour les 100bins
	380	do ibin=1,100
	381	mtemp(ibin)=sum(meteo(lebon1:lebon2,ibin))
	382	enddo
	383	meteo_ccn(ig,l,:)=mtemp(:)/nelem
	384	enddo
	385	enddo
	386	endif
	387
[1617]	388	c Main loop over the GCM's grid
[1685]	389	DO l=1,nlay
	390	DO ig=1,ngrid
[1617]	391	c Get the partial pressure of co2 vapor and its saturation ratio
	392	pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/44.01) * pplay(ig,l)
	393	c satu = zq(ig,l,igcm_co2) / zqsat(ig,l)
	394	satu = pco2 / zqsat(ig,l)
	395	!=============================================================
	396	! 3. Nucleation
	397	!=============================================================
[1685]	398	rho_ice_co2T(ig,l)=1000.(1.72391-2.53e-4zt(ig,l)
	399	& -2.87e-6zt(ig,l)zt(ig,l))
[1720]	400	vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l))
[1685]	401	rho_ice_co2=rho_ice_co2T(ig,l)
[1720]	402
[1617]	403	IF ( satu .ge. 1d0 ) THEN ! if there is condensation
	404
	405	rdust(ig,l)= zq(ig,l,igcm_dust_mass)
	406	& *0.75/pi/rho_dust
	407	& / zq(ig,l,igcm_dust_number)
	408	rdust(ig,l)= rdust(ig,l)**(1./3.)
[1720]	409	c write(,) "Improved2, l,Rdust = ",l,rdust(ig,l)
	410	rdust(ig,l)=max(1.e-9,rdust(ig,l))
	411	c rdust(ig,l)=min(1.e-5,rdust(ig,l))
[1685]	412	! write(,) "Improved3,Rdust = ",rdust(ig,l)
	413
[1617]	414	c Expand the dust moments into a binned distribution
[1685]	415	Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)+1.e-30
	416	No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30
[1617]	417	Rn = rdust(ig,l)
	418	Rn = -log(Rn)
	419	Rm = Rn - 3. * sigma_iceco2*sigma_iceco2
	420	n_derf = erf( (rb_cldco2(1)+Rn) *dev2)
	421	m_derf = erf( (rb_cldco2(1)+Rm) *dev2)
	422
	423	do i = 1, nbinco2_cld
	424	n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
	425	m_aer(i) = -0.5 * Mo * m_derf !! this ith previously computed
	426	n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
	427	m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
[1720]	428	n_aer(i) = n_aer(i) + 0.5 * No * n_derf +
	429	& meteo_ccn(ig,l,i) !Ajout meteo_ccn
[1617]	430	m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
[1720]	431	m_aer2(i)=4./3.pirho_dust
	432	& n_aer(i)rad_cldco2(i)*rad_cldco2(i)
	433	& *rad_cldco2(i)
	434	c write(,) "diff =",rad_cldco2(i),m_aer(i),m_aer2(i)
[1617]	435	enddo
[1720]	436	c write(,) "Bilan =",sum(m_aer),sum(m_aer2)
[1617]	437
[1720]	438	c sumcheck = 0
	439	c do i = 1, nbinco2_cld
	440	c sumcheck = sumcheck + n_aer(i)
	441	c enddo
	442	c sumcheck = abs(sumcheck/No - 1)
	443	c if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
	444	c print*, "WARNING, No sumcheck PROBLEM"
	445	c print*, "sumcheck, No, rdust",sumcheck, No,rdust(ig,l)
	446	c print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
	447	c print*, "Dust binned distribution", n_aer
	448	c STOP
	449	c endif
	450
	451	c sumcheck = 0
	452	c do i = 1, nbinco2_cld
	453	c sumcheck = sumcheck + m_aer(i)
	454	c enddo
	455	c sumcheck = abs(sumcheck/Mo - 1)
	456	c if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld))
	457	c & then
	458	c print*, "WARNING, Mo sumcheck PROBLEM"
	459	c print*, "sumcheck, Mo",sumcheck, Mo
	460	c print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l
	461	c print*, "Dust binned distribution", m_aer
	462	c STOP
	463	c endif
	464	m_aer(:)=m_aer2(:)
	465
	466
	467	rccnh2o(ig,l)= zq(ig,l,igcm_ccn_mass)
	468	& *0.75/pi/rho_dust
	469	& / zq(ig,l,igcm_ccn_number)
	470
	471	rice(ig,l)=( zq(ig,l,igcm_h2o_ice)*3.0/
	472	& (4.0rho_ice_co2zq(ig,l,igcm_ccn_number)
	473	& pitauscaling(ig)) +rccnh2o(ig,l)*rccnh2o(ig,l)
	474	& rccnh2o(ig,l) )*(1.0/3.0)
	475	rhocloud(ig,l)=( zq(ig,l,igcm_h2o_ice)*rho_ice
	476	& +zq(ig,l,igcm_ccn_mass) tauscaling(ig)rho_dust)
	477	& / (zq(ig,l,igcm_h2o_ice)+zq(ig,l,igcm_ccn_mass)
	478	& *tauscaling(ig))
	479	c call updaterice_micro(
	480	c & zq(ig,l,igcm_h2o_ice), ! ice mass
	481	c & zq(ig,l,igcm_ccn_mass), ! ccn mass
	482	c & zq(ig,l,igcm_ccn_number), ! ccn number
	483	c & tauscaling(ig),rice(ig,l),rhocloud(ig,l))
[1685]	484	! rice radius of CCN + H20 crystal
[1720]	485	c write(,) "Improved1 Rice=",rice(ig,l)
[1685]	486	rice(ig,l)=max(1.e-10,rice(ig,l))
[1720]	487	c rice(ig,l)=min(1.e-5,rice(ig,l))
[1685]	488	!write(,) "Improved2 Rice=",rice(ig,l)
	489	Mo = zq(ig,l,igcm_h2o_ice) +
[1720]	490	& zq(ig,l,igcm_ccn_mass)tauscaling(ig)+1.e-30!tauscaling(ig)
	491	! & + 1.e-30 !Total mass of H20 crystals,CCN included
[1617]	492	No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
	493	Rn = rice(ig,l)
	494	Rn = -log(Rn)
[1720]	495	Rm = Rn - 3. * sigma_ice*sigma_ice
	496	n_derf = erf( (rb_cldco2(1)+Rn) *dev3)
	497	m_derf = erf( (rb_cldco2(1)+Rm) *dev3)
[1617]	498	do i = 1, nbinco2_cld
[1720]	499	n_aer_h2oice(i) = -0.5 * No * n_derf
	500	m_aer_h2oice(i) = -0.5 * Mo * m_derf
	501	n_derf = derf( (rb_cldco2(i+1)+Rn) *dev3)
	502	m_derf = derf( (rb_cldco2(i+1)+Rm) *dev3)
	503	n_aer_h2oice(i) = n_aer_h2oice(i) + 0.5 * No * n_derf
	504	m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf
[1685]	505	rad_h2oice(i) = rad_cldco2(i)
[1720]	506	m_aer_h2oice2(i)=4./3.pirhocloud(ig,l)
	507	& n_aer_h2oice(i)rad_h2oice(i)*rad_h2oice(i)
	508	& *rad_h2oice(i)
	509
	510
[1685]	511	c write(,) "before nuc, i,rad_h2o(i)= ",i,rad_cldco2(i)
	512	c & ,m_aer_h2oice(i),n_aer_h2oice(i)
[1617]	513	enddo
[1720]	514	c sumcheck = 0
	515	c do i = 1, nbinco2_cld
	516	c sumcheck = sumcheck + n_aer_h2oice(i)
	517	c enddo
	518	c sumcheck = abs(sumcheck/No - 1)
	519	c if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
	520	c print*, "WARNING, Noh2o sumcheck PROBLEM"
	521	c print*, "sumcheck, No, rice",sumcheck, No,rice(ig,l)
	522	c print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
	523	c print*, "Dust binned distribution", n_aer_h2oice
	524	c STOP
	525	c endif
[1617]	526
[1720]	527	c sumcheck = 0
	528	c do i = 1, nbinco2_cld
	529	c sumcheck = sumcheck + m_aer_h2oice(i)
	530	c enddo
	531	c sumcheck = abs(sumcheck/Mo - 1)
	532	c if ((sumcheck .gt. 1e-5) .and. (1./Rn .gt. rmin_cld))
	533	c & then
	534	c print*, "WARNING, Moh2o sumcheck PROBLEM"
	535	c print*, "sumcheck, Mo",sumcheck, Mo
	536	c print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l
	537	c print*, "Dust binned distribution", m_aer_h2oice
	538	c STOP
	539	c endif
[1617]	540	c Get the rates of nucleation
[1720]	541
[1617]	542	call nucleaCO2(dble(pco2),zt(ig,l),dble(satu)
	543	& ,n_aer,rate,n_aer_h2oice
[1685]	544	& ,rad_h2oice,rateh2o,vo2co2)
[1720]	545	m_aer_h2oice(:)=m_aer_h2oice2(:)
[1617]	546	dN = 0.
	547	dM = 0.
	548	dNh2o = 0.
	549	dMh2o = 0.
	550	do i = 1, nbinco2_cld
	551	Proba=1.0-dexp(-1.ptimesteprate(i))
[1720]	552	Probah2o=coeffh2o(1.0-dexp(-1.ptimestep*rateh2o(i)))
[1685]	553	dNh2o = dNh2o + n_aer_h2oice(i) * Probah2o
	554	dMh2o = dMh2o + m_aer_h2oice(i) * Probah2o
[1617]	555	dN = dN + n_aer(i) * Proba
	556	dM = dM + m_aer(i) * Proba
[1720]	557
[1617]	558	enddo
	559	! dM masse activée (kg) et dN nb particules par kg d'air
	560
	561	dNN= dN/tauscaling(ig)
	562	dMM= dM/tauscaling(ig)
[1720]	563	dNN=min(dNN,zq(ig,l,igcm_dust_number))
	564	dMM=min(dMM,zq(ig,l,igcm_dust_mass))
	565	! if (dNN .gt. 0) then
	566	! write(,) "Nuclea dNN crees=",dNN
	567	! write(,) "Nuclea dMM crees=",dMM
	568	! endif
[1685]	569	dNNh2o=dNh2o/tauscaling(ig)
[1720]	570	dNNh2o=min(dNNh2o,zq(ig,l,igcm_ccn_number))
[1617]	571
[1720]	572	ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice)
	573	& +zq(ig,l,igcm_ccn_mass)*tauscaling(ig))
	574
	575	dMh2o_ice=dMh2ozq(ig,l,igcm_h2o_ice)ratioh2o_ccn
	576	dMh2o_ccn=dMh2ozq(ig,l,igcm_ccn_mass)
	577	& tauscaling(ig)*ratioh2o_ccn
	578
	579	dMh2o_ccn=dMh2o_ccn/tauscaling(ig)
	580	dMh2o_ccn=min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
	581	dMh2o_ice=min(dMh2o_ice,zq(ig,l,igcm_h2o_ice))
[1617]	582
	583	zq(ig,l,igcm_ccnco2_mass) =
	584	& zq(ig,l,igcm_ccnco2_mass) + dMM
[1720]	585	zq(ig,l,igcm_ccnco2_number) =
[1617]	586	& zq(ig,l,igcm_ccnco2_number) + dNN
	587
[1720]	588	zq(ig,l,igcm_dust_mass)= zq(ig,l,igcm_dust_mass)-dMM
	589	zq(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)-dNN
	590
[1685]	591	c Update CCN for CO2 nucleating on H2O CCN :
	592	! Warning: must keep memory of it
	593	zq(ig,l,igcm_ccnco2_mass) =
[1720]	594	& zq(ig,l,igcm_ccnco2_mass) + dMh2o_ice+dMh2o_ccn
[1685]	595	zq(ig,l,igcm_ccnco2_number) =
	596	& zq(ig,l,igcm_ccnco2_number) + dNNh2o
[1617]	597
[1720]	598	zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)-dNNh2o
	599	zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)-dMh2o_ice
	600	zq(ig,l,igcm_ccn_mass)= zq(ig,l,igcm_ccn_mass)-dMh2o_ccn
	601
[1617]	602
[1720]	603	memdMMh2o(ig,l)=memdMMh2o(ig,l)+dMh2o_ice
	604	memdMMccn(ig,l)=memdMMccn(ig,l)+dMh2o_ccn
[1685]	605	memdNNccn(ig,l)=memdNNccn(ig,l)+dNNh2o
[1617]	606	ENDIF ! of is satu >1
	607	!=============================================================
	608	! 4. Ice growth: scheme for radius evolution
	609	!=============================================================
	610
	611	c We trigger crystal growth if and only if there is at least one nuclei (N>1).
	612	c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
	613	c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
[1720]	614	IF (zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) .ge. 1)THEN
[1685]	615	! we trigger crystal growth
	616	c
[1720]	617	c write(,) "ccn number mass=",zq(ig,l,igcm_ccnco2_number),
	618	c & zq(ig,l,igcm_ccnco2_mass)
[1617]	619
[1685]	620	c Niceco2 = zq(ig,l,igcm_co2_ice)
	621	c Qccnco2 = zq(ig,l,igcm_ccnco2_mass)
	622	c Nccnco2 = zq(ig,l,igcm_ccnco2_number)
	623	c call updaterice_microco2(Niceco2,Qccnco2,Nccnco2,
	624	c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
	625	c write(,) "updater rice=",riceco2(ig,l)
	626
[1617]	627	rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass)
	628	& *0.75/pi/rho_dust
	629	& / zq(ig,l,igcm_ccnco2_number)
[1685]	630	rdust(ig,l)= rdust(ig,l)**(1./3.)
[1617]	631	rdust(ig,l)=max(1.e-10,rdust(ig,l))
[1720]	632	c rdust(ig,l)=min(5.e-6,rdust(ig,l))
	633
[1617]	634	riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/
[1685]	635	& (4.0rho_ice_co2zq(ig,l,igcm_ccnco2_number)
	636	& pitauscaling(ig)) +rdust(ig,l)*rdust(ig,l)
[1617]	637	& rdust(ig,l) )*(1.0/3.0)
	638
[1685]	639	c riceco2(ig,l)=max(1.e-10,riceco2(ig,l))
	640	c riceco2(ig,l)=min(1.e-5,riceco2(ig,l))
[1617]	641	! WATCH OUT: CO2 nuclei is supposed to be dust
	642	! only when deriving rhocloud (otherwise would need to keep info on water embedded in co2) - listo
[1649]	643	c write(,) "Rdust before growth = ",rdust(ig,l)
	644	c write(,) "Riceco2 before growth = ",riceco2(ig,l)
[1617]	645
[1632]	646	!! Niceco2,Qccnco2,Nccnco2
[1685]	647	c Niceco2 = zq(ig,l,igcm_co2_ice)
	648	c Qccnco2 = zq(ig,l,igcm_ccnco2_mass)
	649	c Nccnco2 = zq(ig,l,igcm_ccnco2_number)
	650	c call updaterice_microco2(Niceco2,Qccnco2,Nccnco2,
	651	c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
[1720]	652	c write(,) "Riceco2 before growth = ",riceco2(ig,l)
[1685]	653	c write(,) "rdust before growth = ",rdust(ig,l)
[1720]	654	c write(,) "co2ice before growth=",zq(ig,l,igcm_co2_ice)
	655	c write(,) "co2 before growth=",zq(ig,l,igcm_co2)
[1685]	656	c write(,) "pplay before growth=",pplay(ig,l)
	657	c write(,) "zt before growth =",zt(ig,l)
[1720]	658	c write(,) "Satu before growth=",satu
[1685]	659	c riceco2(ig,l)=max(riceco2(ig,l),rdust(ig,l))
	660	No = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig)+1.e-30
[1617]	661	! No nb de particules de poussieres mis à l'échelle pour donner une opacité optique
	662
	663	c saturation at equilibrium
	664	c rice should not be too small, otherwise seq value is not valid
	665	c seq = exp(2.sigco2mco2 / (rho_ice_co2rgpzt(ig,l)*
	666	c & max(riceco2(ig,l),1.e-7))) !Exponant sans unité OK
	667
	668	ccccccc Scheme of microphys. mass growth for CO2
	669
	670	call massflowrateCO2(pplay(ig,l),zt(ig,l),
	671	& satu,riceco2(ig,l),mmean(ig,l),Ic_rice) ! Mass transfer rate (kg/s) for a rice particle
	672	! Ic_rice mass flux kg.s-1 <0 si croissance !
[1720]	673	if (isnan(Ic_rice)) then
	674	Ic_rice=0.
	675	flag_pourri=1
	676	endif
	677	c drsurdt=-1.0/(4.0piriceco2(ig,l)*
	678	c & riceco2(ig,l)rho_ice_co2)Ic_rice
[1685]	679	dMice = No * Ic_rice*ptimestep ! Kg par kg d'air, <0 si croissance !
[1649]	680	c write(,) "dMicev0 in improved = " , dMice
[1617]	681
[1685]	682	if (dMice .ge. 0d0) then
	683	dMice = min(dMice,abs(zq(ig,l,igcm_co2_ice)))
	684	else
	685	dMice =-1.* min(abs(dMice),abs(zq(ig,l,igcm_co2)))
	686	endif
[1720]	687	c riceco2(ig,l)=riceco2(ig,l)+drsurdt*ptimestep
[1617]	688	c write(,) "riceco2+dr/dt = ", riceco2(ig,l)
[1720]	689	!write(,) "dMice in improved = " , dMice
[1617]	690
[1720]	691	zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice)-dMice
[1685]	692	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)+dMice
[1617]	693
[1720]	694	! latent heat release >0 if growth i.e. if dMice <0
[1617]	695
	696
	697	lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 +
	698	& l3 * zt(ig,l)*3 + l4 zt(ig,l)**4 !J.kg-1
	699	c write(,) "CPP= ",cpp ! = 744.5
[1720]	700	!Peut etre un probleme de signe ici!
	701	pdtcloudco2(ig,l)= -1.dMicelw/cpp/ptimestep ! kg par couche * J par kg /J par K / s = K par seconde
[1617]	702
[1720]	703	!write(,) "pdtcloudco2 after growth = ",pdtcloudco2(ig,l)
[1617]	704
[1720]	705	!write(,) "co2 after growth = ",zq(ig,l,igcm_co2)
	706	!write(,) "co2_ice after growth = ",zq(ig,l,igcm_co2_ice)
[1617]	707
	708	!deltaT par condens/subli. qui remplace le dT du CO2 du newcondens pré-constantino
	709	!PDT should be in K/s.
	710	!=============================================================
	711	! 5. Dust cores released, tendancies, latent heat, etc ...
	712	!=============================================================
	713
[1720]	714	c if (abs(dMice) .ge. 1.e-10) then
	715
	716	c write(,) "DIAG zq pdt",(zq(ig,l,igcm_co2_ice)-
	717	c & zq0(ig,l,igcm_co2_ice))/ptimestep,pdtcloudco2(ig,l)
	718	c endif
	719	ENDIF ! of if NCCN > 1
[1617]	720
[1685]	721	rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass)
	722	& *0.75/pi/rho_dust
	723	& / zq(ig,l,igcm_ccnco2_number)
	724	rdust(ig,l)= rdust(ig,l)**(1./3.)
[1720]	725	rdust(ig,l)=max(1.e-9,rdust(ig,l))
	726	c rdust(ig,l)=min(5.e-6,rdust(ig,l))
[1685]	727
	728	riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/
	729	& (4.0rho_ice_co2zq(ig,l,igcm_ccnco2_number)
	730	& pitauscaling(ig)) +rdust(ig,l)*rdust(ig,l)
	731	& rdust(ig,l) )*(1.0/3.0)
	732	!Niceco2 = zq(ig,l,igcm_co2_ice)
	733	!Qccnco2 = zq(ig,l,igcm_ccnco2_mass)
	734	!Nccnco2 = zq(ig,l,igcm_ccnco2_number)
	735	c
	736	c call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
	737	c & tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
[1720]	738
	739	c If there is no more co2 ice, all the ice particles sublimate, all the condensation nuclei are released:
	740
	741	if ((zq(ig,l,igcm_co2_ice).le. 1.e-25)) then
	742	c & .or. flag_pourri .eq. 1
	743	c & .or.(riceco2(ig,l) .le. rdust(ig,l))
	744	c & .or. (l .ge.1 .and. l .le. 5)
	745	c & .or. (zq(ig,l,igcm_co2_ice) .ge. 0.1)
	746	lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 +
	747	& l3 * zt(ig,l)*3 + l4 zt(ig,l)**4 !J.kg-1
	748	c write(,) "CPP= ",cpp ! = 744.5
	749
	750	pdtcloudco2(ig,l)=pdtcloudco2(ig,l)-1.
	751	& zq(ig,l,igcm_co2_ice)lw/cpp/ptimestep !
	752	!On sublime tout !
[1685]	753	c write(,) "Riceco2 improved before reset=",riceco2(ig,l)
	754	c write(,) "Niceco2 improved before reset=",
	755	c & zq(ig,l,igcm_co2_ice)
	756	c write(,) "Rdust improved before reset=",rdust(ig,l)
	757
	758	if (memdMMccn(ig,l) .gt. 0) then
	759	zq(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
	760	& +memdMMccn(ig,l)
	761	endif
	762	if (memdMMh2o(ig,l) .gt. 0) then
	763	zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)
	764	& +memdMMh2o(ig,l)
	765	endif
	766
	767	if (memdNNccn(ig,l) .gt. 0) then
	768	zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)
	769	& +memdNNccn(ig,l)
	770	endif
	771
[1720]	772	c if (zq(ig,l,igcm_ccnco2_mass) .gt. 1.e-30) then
[1685]	773	zq(ig,l,igcm_dust_mass) =
	774	& zq(ig,l,igcm_dust_mass)
	775	& + zq(ig,l,igcm_ccnco2_mass)-
	776	& (memdMMh2o(ig,l)+memdMMccn(ig,l))
[1720]	777	c endif
	778	c if (zq(ig,l,igcm_ccnco2_number) .gt. 1.e-30) then
[1685]	779	zq(ig,l,igcm_dust_number) =
	780	& zq(ig,l,igcm_dust_number)
	781	& + zq(ig,l,igcm_ccnco2_number)-memdNNccn(ig,l)
[1720]	782	c endif
[1685]	783
[1720]	784	c if (zq(ig,l,igcm_co2_ice) .gt. 1.e-30) then
[1685]	785	zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)
	786	& + zq(ig,l,igcm_co2_ice)
[1720]	787	c endif
[1685]	788
	789	zq(ig,l,igcm_ccnco2_mass)=0.
	790	zq(ig,l,igcm_co2_ice)=0.
	791	zq(ig,l,igcm_ccnco2_number)=0.
	792	memdNNccn(ig,l)=0.
	793	memdMMh2o(ig,l)=0.
	794	memdMMccn(ig,l)=0.
	795	riceco2(ig,l)=0.
[1720]	796
[1685]	797	endif
[1720]	798	ENDDO ! of ig loop
	799	ENDDO ! of nlayer loop
	800	! Get cloud tendencies
[1617]	801	pdqcloudco2(1:ngrid,1:nlay,igcm_co2) =
	802	& (zq(1:ngrid,1:nlay,igcm_co2) -
	803	& zq0(1:ngrid,1:nlay,igcm_co2))/ptimestep
	804	pdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) =
	805	& (zq(1:ngrid,1:nlay,igcm_co2_ice) -
	806	& zq0(1:ngrid,1:nlay,igcm_co2_ice))/ptimestep
[1685]	807	pdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
	808	& (zq(1:ngrid,1:nlay,igcm_h2o_ice) -
	809	& zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep
	810	pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
	811	& (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
	812	& zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep
	813	pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
	814	& (zq(1:ngrid,1:nlay,igcm_ccn_number) -
	815	& zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep
[1617]	816	pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) =
	817	& (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) -
	818	& zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/ptimestep
	819	pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) =
	820	& (zq(1:ngrid,1:nlay,igcm_ccnco2_number) -
	821	& zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/ptimestep
[1720]	822	pdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) =
	823	& (zq(1:ngrid,1:nlay,igcm_dust_mass) -
	824	& zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep
	825	pdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) =
	826	& (zq(1:ngrid,1:nlay,igcm_dust_number) -
	827	& zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
[1617]	828	return
	829	end
	830
	831
	832

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

Download in other formats:

Original Format