Context Navigation

condense_co2.F90 @ 3562

Last change on this file since 3562 was 2972, checked in by emillour, 19 months ago
Generic PCM: Make number of scatterers fully dynamic (i.e. set in callphys.def and no longer by compilation option "-s #"). One should now specify naerkind = # in callphys.def (default is 0). EM
Property svn:executable set to ``*
File size: 18.4 KB

Rev	Line
[1477]	1	subroutine condense_co2(ngrid,nlayer,nq,ptimestep, &
[1482]	2	pcapcal,pplay,pplev,ptsrf,pt, &
[1485]	3	pdt,pdtsrf,pq,pdq, &
	4	pqsurf,pdqsurfc,albedo,pemisurf, &
[1482]	5	albedo_bareground,albedo_co2_ice_SPECTV, &
[1485]	6	pdtc,pdtsrfc,pdpsrfc,pdqc)
[305]	7
[2972]	8	use radinc_h, only : L_NSPECTV
[1216]	9	use gases_h, only: gfrac, igas_co2
[728]	10	use radii_mod, only : co2_reffrad
	11	use aerosol_mod, only : iaero_co2
[1482]	12	USE surfdat_h, only: emisice, emissiv
[1543]	13	USE geometry_mod, only: latitude ! in radians
[1216]	14	USE tracer_h, only: noms, rho_co2
[1384]	15	use comcstfi_mod, only: g, r, cpp
[305]	16
	17	implicit none
	18
	19	!==================================================================
	20	! Purpose
	21	! -------
[1485]	22	! Condense and/or sublime CO2 ice on the ground and in the atmosphere, and sediment the ice.
[305]	23	!
	24	! Inputs
	25	! ------
[1485]	26	! ngrid Number of vertical columns.
	27	! nlayer Number of vertical layers.
	28	! nq Number of tracers.
	29	! ptimestep Duration of the physical timestep (s).
	30	! pplay(ngrid,nlayer) Pressure layers (Pa).
	31	! pplev(ngrid,nlayer+1) Pressure levels (Pa).
	32	! pt(ngrid,nlayer) Atmospheric Temperatures (K).
	33	! ptsrf(ngrid) Surface temperatures (K).
	34	! pq(ngrid,nlayer,nq) Atmospheric tracers mixing ratios (kg/kg of air).
	35	! pqsurf(ngrid,nq) Surface tracers (kg/m2).
[305]	36	!
[1485]	37	! pdt(ngrid,nlayer) Time derivative before condensation/sublimation of pt.
	38	! pdtsrf(ngrid) Time derivative before condensation/sublimation of ptsrf.
	39	! pdq(ngrid,nlayer,nq) Time derivative before condensation/sublimation of
	40	!
	41	! albedo_bareground(ngrid) Albedo of the bare ground.
	42	! albedo_co2_ice_SPECTV(L_NSPECTV) Spectral albedo of CO2 ice.
[305]	43	!
	44	! Outputs
	45	! -------
[1485]	46	! pdpsrfc(ngrid) \ Contribution of condensation/sublimation
	47	! pdtc(ngrid,nlayer) \ to the time derivatives of
	48	! pdtsrfc(ngrid) / Surface Pressure, Atmospheric Temperatures,
	49	! pdqsurfc(ngrid) / Surface Temperatures, Surface Tracers,
	50	! pdqc(ngrid,nlayer,nq) / and Atmospheric Tracers.*
	51	!
	52	! pemisurf(ngrid) Emissivity of the surface.
[305]	53	!
	54	! Both
	55	! ----
[1485]	56	! albedo(ngrid,L_NSPECTV) Spectral albedo of the surface.
[305]	57	!
	58	! Authors
	59	! -------
	60	! Francois Forget (1996)
	61	! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
[526]	62	! Includes simplifed nucleation by J. Leconte (2011)
[305]	63	!
	64	!==================================================================
	65
[1485]	66	!--------------------------
	67	! Arguments
	68	!--------------------------
[305]	69
[1485]	70
[858]	71	INTEGER,INTENT(IN) :: ngrid
	72	INTEGER,INTENT(IN) :: nlayer
	73	INTEGER,INTENT(IN) :: nq
	74	REAL,INTENT(IN) :: ptimestep
	75	REAL,INTENT(IN) :: pcapcal(ngrid)
	76	REAL,INTENT(IN) :: pplay(ngrid,nlayer)
	77	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)
	78	REAL,INTENT(IN) :: ptsrf(ngrid)
	79	REAL,INTENT(IN) :: pt(ngrid,nlayer)
	80	REAL,INTENT(IN) :: pdt(ngrid,nlayer)
	81	REAL,INTENT(IN) :: pdtsrf(ngrid)
	82	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
[1485]	83	REAL,INTENT(IN) :: pqsurf(ngrid,nq)
[858]	84	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
[1482]	85	REAL,INTENT(IN) :: albedo_bareground(ngrid)
	86	REAL,INTENT(IN) :: albedo_co2_ice_SPECTV(L_NSPECTV)
[1485]	87	REAL,INTENT(INOUT) :: albedo(ngrid,L_NSPECTV)
[858]	88	REAL,INTENT(OUT) :: pemisurf(ngrid)
	89	REAL,INTENT(OUT) :: pdtc(ngrid,nlayer)
	90	REAL,INTENT(OUT) :: pdtsrfc(ngrid)
[1485]	91	REAL,INTENT(OUT) :: pdpsrfc(ngrid)
[858]	92	REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq)
[1484]	93	REAL,INTENT(OUT) :: pdqsurfc(ngrid)
[305]	94
[1485]	95	!------------------------------
	96	! Local variables
	97	!------------------------------
[305]	98
[1485]	99	INTEGER l,ig,icap,ilay,iq,nw,igas,it
[305]	100
[1485]	101	REAL reffrad(ngrid,nlayer) ! Radius (m) of the CO2 ice particles.
	102	REAL*8 zt(ngrid,nlayer) ! Updated Atmospheric Temperatures (K).
	103	REAL ztsrf(ngrid) ! Updated Surface Temperatures (K).
	104	REAL zq(ngrid,nlayer,nq) ! Updated Atmospheric tracers mixing ratios (kg/kg of air).
	105	REAL piceco2(ngrid) ! Updated Surface Tracer (kg/m2).
	106	REAL ztcond (ngrid,nlayer) ! Atmospheric Temperatures of condensation of CO2.
	107	REAL ztnuc (ngrid,nlayer) ! Atmospheric Nucleation Temperatures.
	108	REAL ztcondsol(ngrid) ! Temperatures of condensation of CO2 at the surface.
	109	REAL zcondices(ngrid) ! Condensation rate on the ground (kg/m2/s).
	110	REAL zfallice(ngrid) ! Flux of ice falling on the surface (kg/m2/s).
	111	REAL Mfallice(ngrid) ! Total amount of ice fallen to the ground during the timestep (kg/m2).
	112	REAL wq(ngrid,nlayer+1) ! Total amount of ice fallen to the ground during the timestep (kg/m2).
	113	REAL subptimestep ! Duration of the subtimestep (s) for the sedimentation.
	114	Integer Ntime ! Number of subtimesteps.
	115	REAL masse (ngrid,nlayer) ! Mass of atmospheric layers (kg/m2)
	116	REAL w(ngrid,nlayer,nq) !
	117	REAL vstokes,reff !
	118	REAL ppco2 !
[305]	119
	120
[1485]	121	!------------------------------------------
	122	! Saved local variables
	123	!------------------------------------------
[305]	124
[1485]	125
[858]	126	REAL,SAVE :: latcond=5.9e5
	127	REAL,SAVE :: ccond
[787]	128	REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: emisref
[1485]	129	!$OMP THREADPRIVATE(latcond,ccond,emisref)
[858]	130	LOGICAL,SAVE :: firstcall=.true.
[1315]	131	!$OMP THREADPRIVATE(firstcall)
[305]	132	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
[1315]	133	!$OMP THREADPRIVATE(i_co2ice)
[305]	134	CHARACTER(LEN=20) :: tracername ! to temporarily store text
	135
	136
[1485]	137	!------------------------------------------------
	138	! Initialization at the first call
	139	!------------------------------------------------
[305]	140
	141
	142	IF (firstcall) THEN
	143
[1485]	144	ALLOCATE(emisref(ngrid))
	145	! Find CO2 ice tracer.
[787]	146	do iq=1,nq
[305]	147	tracername=noms(iq)
	148	if (tracername.eq."co2_ice") then
	149	i_co2ice=iq
	150	endif
	151	enddo
	152
[1485]	153	write(,) "condense_co2: i_co2ice=",i_co2ice
[305]	154
[1485]	155	if((i_co2ice.lt.1))then
	156	print*,'In condens_cloud but no CO2 ice tracer, exiting.'
	157	print*,'Still need generalisation to arbitrary species!'
	158	stop
	159	endif
[305]	160
	161	ccond=cpp/(g*latcond)
[486]	162	print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond
[305]	163
[1485]	164	! Prepare special treatment if gas is not pure CO2
	165	! if (addn2) then
	166	! m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol)
	167	! m_noco2 = 28.02E-3 ! N2 molecular mass (kg/mol)
	168	! Compute A and B coefficient use to compute
	169	! mean molecular mass Mair defined by
	170	! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
	171	! 1/Mair = A*q(ico2) + B
	172	! A = (1/m_co2 - 1/m_noco2)
	173	! B = 1/m_noco2
	174	! endif
[305]	175
[1485]	176	! Minimum CO2 mixing ratio below which mixing occurs with layer above : qco2min =0.75
[305]	177
	178	firstcall=.false.
	179	ENDIF
	180
	181
[1485]	182	!------------------------------------------------
	183	! Tendencies initially set to 0
	184	!------------------------------------------------
[305]	185
[1485]	186
	187	pdqc(1:ngrid,1:nlayer,1:nq) = 0.
	188	pdtc(1:ngrid,1:nlayer) = 0.
	189	zq(1:ngrid,1:nlayer,1:nq) = 0.
	190	zt(1:ngrid,1:nlayer) = 0.
	191	Mfallice(1:ngrid) = 0.
	192	zfallice(1:ngrid) = 0.
	193	zcondices(1:ngrid) = 0.
	194	pdtsrfc(1:ngrid) = 0.
	195	pdpsrfc(1:ngrid) = 0.
	196	pdqsurfc(1:ngrid) = 0.
	197
	198
	199	!----------------------------------
[305]	200	! Atmospheric condensation
[1485]	201	!----------------------------------
[305]	202
	203
	204	! Compute CO2 Volume mixing ratio
	205	! -------------------------------
	206	! if (addn2) then
	207	! DO l=1,nlayer
	208	! DO ig=1,ngrid
	209	! qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep
[1485]	210	! Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2)
[305]	211	! mmean=1/(A*qco2 +B)
	212	! vmr_co2(ig,l) = qco2*mmean/m_co2
	213	! ENDDO
	214	! ENDDO
	215	! else
	216	! DO l=1,nlayer
	217	! DO ig=1,ngrid
	218	! vmr_co2(ig,l)=0.5
	219	! ENDDO
	220	! ENDDO
	221	! end if
	222
[1485]	223
	224	! Forecast the atmospheric frost temperature 'ztcond' and nucleation temperature 'ztnuc'.
[305]	225	DO l=1,nlayer
	226	DO ig=1,ngrid
[869]	227	ppco2=gfrac(igas_CO2)*pplay(ig,l)
[305]	228	call get_tcond_co2(ppco2,ztcond(ig,l))
[526]	229	call get_tnuc_co2(ppco2,ztnuc(ig,l))
[305]	230	ENDDO
	231	ENDDO
	232
[1485]	233	! Initialize zq and zt at the beginning of the sub-timestep loop and qsurf.
	234	DO ig=1,ngrid
	235	piceco2(ig)=pqsurf(ig,i_co2ice)
	236	DO l=1,nlayer
[305]	237	zt(ig,l)=pt(ig,l)
	238	zq(ig,l,i_co2ice)=pq(ig,l,i_co2ice)
	239	IF( zq(ig,l,i_co2ice).lt.-1.e-6 ) THEN
	240	print*,'Uh-oh, zq = ',zq(ig,l,i_co2ice),'at ig,l=',ig,l
	241	if(l.eq.1)then
	242	print*,'Perhaps the atmosphere is collapsing on surface...?'
	243	endif
	244	END IF
	245	ENDDO
	246	ENDDO
	247
[1485]	248	! Calculate the mass of each atmospheric layer (kg.m-2)
[305]	249	do ilay=1,nlayer
[787]	250	DO ig=1,ngrid
[305]	251	masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g
	252	end do
	253	end do
	254
[1485]	255
	256	!-----------------------------------------------------------
[305]	257	! START CONDENSATION/SEDIMENTATION SUB-TIME LOOP
[1485]	258	!-----------------------------------------------------------
	259
	260
	261	Ntime = 20 ! number of sub-timestep
[305]	262	subptimestep = ptimestep/float(Ntime)
	263
[1485]	264	! Add the tendencies from other physical processes at each subtimstep.
[305]	265	DO it=1,Ntime
	266	DO l=1,nlayer
	267	DO ig=1,ngrid
	268	zt(ig,l) = zt(ig,l) + pdt(ig,l) * subptimestep
	269	zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdq(ig,l,i_co2ice) * subptimestep
	270	END DO
	271	END DO
	272
[1485]	273	! Gravitational sedimentation starts.
[305]	274
[1485]	275	! Sedimentation computed from radius computed from q in module radii_mod.
[1308]	276	call co2_reffrad(ngrid,nlayer,nq,zq,reffrad)
[728]	277
[1485]	278	DO ilay=1,nlayer
[787]	279	DO ig=1,ngrid
[305]	280
[858]	281	reff = reffrad(ig,ilay)
[305]	282
	283	call stokes &
	284	(pplev(ig,ilay),pt(ig,ilay), &
	285	reff,vstokes,rho_co2)
	286
	287	!w(ig,ilay,i_co2ice) = 0.0
	288	w(ig,ilay,i_co2ice) = vstokes * subptimestep * &
	289	pplev(ig,ilay)/(r*pt(ig,ilay))
	290
[1485]	291	END DO
	292	END DO
[305]	293
[1485]	294	! Computing q after sedimentation
[1308]	295	call vlz_fi(ngrid,nlayer,zq(1,1,i_co2ice),2.,masse,w(1,1,i_co2ice),wq)
[305]	296
	297
[1485]	298	! Progressively accumulating the flux to the ground.
	299	! Mfallice is the total amount of ice fallen to the ground.
[787]	300	DO ig=1,ngrid
[305]	301	Mfallice(ig) = Mfallice(ig) + wq(ig,i_co2ice)
[1485]	302	END DO
[305]	303
[1485]	304	!----------------------------------------------------------
	305	! Condensation / sublimation in the atmosphere
	306	!----------------------------------------------------------
	307	! (MODIFICATIONS FOR EARLY MARS: falling heat neglected, condensation of CO2 into tracer i_co2ice)
[305]	308
	309
	310	DO l=nlayer , 1, -1
	311	DO ig=1,ngrid
	312	pdtc(ig,l)=0.
	313
[1485]	314	! ztcond-> ztnuc in test beneath to nucleate only when super saturation occurs(JL 2011)
[526]	315	IF ((zt(ig,l).LT.ztnuc(ig,l)).or.(zq(ig,l,i_co2ice).gt.1.E-10)) THEN
[305]	316	pdtc(ig,l) = (ztcond(ig,l) - zt(ig,l))/subptimestep
	317	pdqc(ig,l,i_co2ice) = pdtc(ig,l)ccondg
	318
[1485]	319	! Case when the ice from above sublimes entirely
[305]	320	IF ((zq(ig,l,i_co2ice).lt.-pdqc(ig,l,i_co2ice)*subptimestep) &
	321	.AND. (zq(ig,l,i_co2ice).gt.0)) THEN
	322
	323	pdqc(ig,l,i_co2ice) = -zq(ig,l,i_co2ice)/subptimestep
	324	pdtc(ig,l) =-zq(ig,l,i_co2ice)/(ccondgsubptimestep)
	325
	326	END IF
	327
[1485]	328	! Temperature and q after condensation
[305]	329	zt(ig,l) = zt(ig,l) + pdtc(ig,l) * subptimestep
	330	zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdqc(ig,l,i_co2ice) * subptimestep
	331	END IF
	332
	333	ENDDO
	334	ENDDO
[1485]	335
	336	ENDDO! end of subtimestep loop.
[305]	337
[1485]	338	! Computing global tendencies after the subtimestep.
[305]	339	DO l=1,nlayer
	340	DO ig=1,ngrid
	341	pdtc(ig,l) = &
	342	(zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep
	343	pdqc(ig,l,i_co2ice) = &
	344	(zq(ig,l,i_co2ice)-(pq(ig,l,i_co2ice)+pdq(ig,l,i_co2ice)*ptimestep))/ptimestep
	345	END DO
	346	END DO
	347	DO ig=1,ngrid
	348	zfallice(ig) = Mfallice(ig)/ptimestep
	349	END DO
	350
	351
	352	!-----------------------------------------------------------------------
[1485]	353	! Condensation/sublimation on the ground
	354	!-----------------------------------------------------------------------
[305]	355
[1485]	356
	357	! Forecast of ground temperature ztsrf and frost temperature ztcondsol.
[305]	358	DO ig=1,ngrid
[869]	359	ppco2=gfrac(igas_CO2)*pplay(ig,1)
[305]	360	call get_tcond_co2(ppco2,ztcondsol(ig))
	361
	362	ztsrf(ig) = ptsrf(ig)
	363
	364	if((ztsrf(ig).le.ztcondsol(ig)+2.0).and.(ngrid.eq.1))then
	365	print*,'CO2 is condensing on the surface in 1D. This atmosphere is doomed.'
	366	print*,'T_surf = ',ztsrf,'K'
	367	print*,'T_cond = ',ztcondsol,'K'
	368	open(116,file='surf_vals.out')
	369	write(116,*) 0.0, pplev(1,1), 0.0, 0.0
	370	close(116)
	371	call abort
	372	endif
	373
	374	ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep
	375
	376	ENDDO
	377
	378	DO ig=1,ngrid
[1485]	379
[305]	380	IF(ig.GT.ngrid/2+1) THEN
	381	icap=2
	382	ELSE
	383	icap=1
	384	ENDIF
	385
[1485]	386	! Loop over where we have condensation / sublimation
[305]	387	IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. & ! ground condensation
	388	(zfallice(ig).NE.0.) .OR. & ! falling snow
	389	((ztsrf(ig) .GT. ztcondsol(ig)) .AND. & ! ground sublimation
	390	((piceco2(ig)+zfallice(ig)*ptimestep) .NE. 0.))) THEN
	391
[1485]	392
	393	! Condensation or partial sublimation of CO2 ice
[305]	394	zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
	395	/(latcond*ptimestep)
	396	pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
	397
[1485]	398	! If the entire CO_2 ice layer sublimes
	399	! (including what has just condensed in the atmosphere)
[305]	400	IF((piceco2(ig)/ptimestep+zfallice(ig)).LE. &
	401	-zcondices(ig))THEN
	402	zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig)
	403	pdtsrfc(ig)=(latcond/pcapcal(ig))* &
	404	(zcondices(ig))
	405	END IF
	406
[1485]	407	! Changing CO2 ice amount and pressure
	408	piceco2(ig) = piceco2(ig) + pdqsurfc(ig)*ptimestep
	409	pdqsurfc(ig) = zcondices(ig) + zfallice(ig)
	410	pdpsrfc(ig) = -pdqsurfc(ig)*g
[305]	411
[1485]	412	IF(ABS(pdpsrfc(ig)*ptimestep).GT.pplev(ig,1)) THEN
[1526]	413	PRINT*,'STOP in condens in condense_co2'
[305]	414	PRINT*,'condensing more than total mass'
	415	PRINT*,'Grid point ',ig
	416	PRINT*,'Ps = ',pplev(ig,1)
[1485]	417	PRINT*,'d Ps = ',pdpsrfc(ig)
[305]	418	STOP
	419	ENDIF
	420	END IF
[1485]	421
	422	ENDDO ! end of ngrid loop.
[305]	423
[1485]	424
	425	!---------------------------------------------------------------------------------------------
[305]	426	! Surface albedo and emissivity of the ground below the snow (emisref)
[1485]	427	!---------------------------------------------------------------------------------------------
	428
	429
[787]	430	DO ig=1,ngrid
[1485]	431
[1542]	432	IF(latitude(ig).LT.0.) THEN
[1216]	433	icap=2 ! Southern Hemisphere
[305]	434	ELSE
[1216]	435	icap=1 ! Nortnern hemisphere
[305]	436	ENDIF
	437
	438	if(.not.piceco2(ig).ge.0.) THEN
	439	if(piceco2(ig).le.-1.e-8) print*, &
[1484]	440	'WARNING : in condense_co2cloud: piceco2(',ig,')=', piceco2(ig)
[305]	441	piceco2(ig)=0.
	442	endif
[1484]	443	if (piceco2(ig) .gt. 1.) then ! CO2 Albedo condition changed to ~1 mm coverage. Change by MT2015.
[1482]	444	DO nw=1,L_NSPECTV
	445	albedo(ig,nw) = albedo_co2_ice_SPECTV(nw)
	446	ENDDO
[305]	447	emisref(ig) = emisice(icap)
	448	else
[1482]	449	DO nw=1,L_NSPECTV
	450	albedo(ig,nw) = albedo_bareground(ig) ! Note : If you have some water, it will be taken into account in the "hydrol" routine.
	451	ENDDO
[305]	452	emisref(ig) = emissiv
	453	pemisurf(ig) = emissiv
	454	end if
[1484]	455
[1485]	456	END DO
[305]	457
	458	return
[1485]	459
	460	end subroutine condense_co2
[305]	461
[1485]	462
	463
[305]	464	!-------------------------------------------------------------------------
[1485]	465	!-------------------------------------------------------------------------
	466	!-------------------------------------------------------------------------
	467	!-------------------------------------------------------------------------
	468	!-------------------------------------------------------------------------
	469	!-------------------------------------------------------------------------
[305]	470
[1485]	471
	472
	473	subroutine get_tcond_co2(p,tcond) ! Calculates the condensation temperature for CO2
	474
	475
[305]	476	implicit none
	477
	478	real p, peff, tcond
	479	real, parameter :: ptriple=518000.0
	480
[526]	481	peff=p
[305]	482
[1485]	483	if(peff.lt.ptriple) then
	484	tcond = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula.
[305]	485	else
[1485]	486	tcond = 684.2-92.3log(peff)+4.32log(peff)**2 ! liquid-vapour transition (based on CRC handbook 2003 data)
[305]	487	endif
	488	return
	489
[1485]	490	end subroutine get_tcond_co2
	491
	492
	493
[526]	494	!-------------------------------------------------------------------------
[1485]	495	!-------------------------------------------------------------------------
	496	!-------------------------------------------------------------------------
	497	!-------------------------------------------------------------------------
	498	!-------------------------------------------------------------------------
	499	!-------------------------------------------------------------------------
[526]	500
[1485]	501
	502
	503	subroutine get_tnuc_co2(p,tnuc)
	504	! Calculates the nucleation temperature for CO2, based on a simple super saturation criterion. JL 2011.
	505
[1397]	506	use callkeys_mod, only: co2supsat
	507
[526]	508	implicit none
	509
	510	real p, peff, tnuc
	511	real, parameter :: ptriple=518000.0
	512
	513	peff=p/co2supsat
	514
[1485]	515	if(peff.lt.ptriple) then
[526]	516	tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula
	517	else
	518	tnuc = 684.2-92.3log(peff)+4.32log(peff)**2
	519	! liquid-vapour transition (based on CRC handbook 2003 data)
	520	endif
[1485]	521
[526]	522	return
	523
[1485]	524	end subroutine get_tnuc_co2

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

Download in other formats:

Original Format