 
|
Project On-line Deliverables: D05.0
Page under construction;please bear with us.
Model Calibration and Uncertainty AnalysisTechnical Report (DRAFT)
 A APPENDIXA.1 Original Routine
subroutine evapor (t,rpar,evap)
c
c This subroutine computes the evaporation rate per unit area
c at time t after the spill of a (volatile) liquid dependent on
c geographical, climatological and meteorological data of the
c environment as well as the physico-chemical properties of the
c liquid, the surrounding air and the underlying ground.
c Ref.: P.I. Kawamura and D. Mackay,
c J. Hazardous Materials, 15 (1987) 343-364.
c
c INPUT:
c
c t ...... time elapsed since the spill [s]
c
c sa = rpar(1) [rad]
c cloud = rpar(2) [%]
c ubar10 = rpar(3) [m/s]
c tempa = rpar(4) [K]
c diamp = rpar(5) [m]
c depthp = rpar(6) [m]
c wtmol = rpar(7) [g/mol]
c dab = rpar(8) [m2/s]
c conliq = rpar(9) [J/msK]
c boilp = rpar(10) [oC]
c vapht = rpar(11) [J/mol]
c a = rpar(12)
c b = rpar(13)
c c = rpar(14)
c
c sa ....... solar altitude [rad]
c cloud .... cloud cover [%]
c ubar10 ... wind speed measured at 10 m [m/s]
c tempa .... ambient air temperature [K]
c diamp .... pool diameter or downwind length of pool [m]
c depthp ... initial depth of the pool [m]
c wtmol .... molecular weight of the liquid [g/mol]
c dab ...... diffusivity of the liquid in air [m2/s]
c conliq ... thermal conductivity of the liquid [J/msK]
c boilp .... boiling point of the liquid [oC]
c vapht .... heat of vaporization (latent heat of evaporation)
[J/mol]
c a ...... } coefficients for the vapor pressure equation of the
c b ...... } form press = exp(a-b/(T+c)) [Pa]
c c ...... } where T ... temperature [K]
c
c OUTPUT:
c
c evap ..... evaporation rate [g/m2s]
c
c ENDDOC-------------------------------last revision:1988-01-13 mp
implicit undefined (a-z)
c
real t, rpar(*), evap
c
real sa, cloud, ubar10, tempa
real diamp, depthp, a, b, c
real wtmol, dab, conliq, boilp, vapht
real rgas, visair, congrd, difgrd, pi
real qsol, schmdt, trmass, hgrd
real hliqi, ugrd, er, press, ea, z,
s, beta
c
c rgas ..... universal gas constant [Pam3/molK]
c visair ... kinematic viscosity of air [m2/s]
c congrd ... thermal conductivity of saturated sand [J/msK]
c difgrd ... thermal diffusivity of sand [m2/s]
c
data rgas /8.314/
data visair /0.000015/
data congrd /2.083333/
data difgrd /0.7e-6/
data pi /3.1415927/
c
sa = rpar(1)
cloud = rpar(2)
ubar10 = rpar(3)
tempa = rpar(4)
diamp = rpar(5)
depthp = rpar(6)
wtmol = rpar(7)
dab = rpar(8)
conliq = rpar(9)
boilp = rpar(10)
vapht = rpar(11)
a = rpar(12)
b = rpar(13)
c = rpar(14)
c
c *** Preliminary calculations ***
c
c Compute net solar insulation [J/m2s]:
c Ref.: J.M. Raphael,
c Prediction of temperature in rivers and reservoirs,
c Proc. Amer. Soc. Civ. Eng., J.Power Div., 88 P02, 1962.
c
qsol = 1111.1111*(1.-0.000071*cloud*cloud)*(sin(sa)-0.1)
if (qsol .lt. 0.) qsol = 0.
c
c Compute mass transfer coefficient [m/s]:
c Ref.: D. Mackay and R.S. Matsugu, Evaporation rate of hydrocarbon
c spills on water and land, Can. J. Chem. Eng., 5 (1973) 434.
c
schmdt = visair/dab
trmass = 0.0047864*ubar10**0.78/(diamp**0.11*schmdt**0.67)
c
c Compute ground heat transfer coefficients:
c
c 1) time dependent heat transfer coefficient at the
c liquid-ground interface [J/m2sK]:
c Ref.: J.P. Holman, Heat Transfer, McGraw Hill, New York, 1976.
c
hgrd = 2.*congrd/sqrt(pi*difgrd*t)
c
c 2) (invers of) liquid heat transfer coefficient [J/m2sK]:
c
hliqi = depthp/(conliq*(1.+exp(-0.06*(boilp-70.))))
c
c 3) overall heat transfer coefficient [J/m2sK]:
c
ugrd = 1./(1./hgrd+hliqi)
c
c *** Direct evaporation model ***
c
c Evaporation due to solar insulation:
er = qsol*wtmol/vapht
c
c Liquid evaporation rate:
c
press = exp(a-b/(tempa-273+c))
ea = trmass*wtmol*press/(rgas*tempa)
c
c Compute factor z and slope of vapor pressure curve s [Pa/K]:
c
c
z = (3650120.*schmdt**0.67+ugrd*rgas*tempa/trmass)/vapht
s = press*b/(tempa-273+c)**2
c
c Finally compute evaporation rate [g/m2s]:
c
beta = s/(s+z)
evap = er*beta+ea*(1.-beta)
return
end
A.2 Changed Declarations
c
subroutine evapor (t,upar,rpar,evap)
c
c INPUT:
c
c t ...... time elapsed since the spill [s]
c
sa = upar(1)
cloud = upar(2)
ubar10 = upar(3)
tempa = upar(4)
diamp = upar(5)
depthp = upar(6)
wtmol = rpar(1)
dab = rpar(2)
conliq = rpar(3)
boilp = rpar(4)
vapht = rpar(5)
a = rpar(6)
b = rpar(7)
c = rpar(8)
A.3 Routine created by ADIFOR
C DISCLAIMER
C
C This file was generated on 11/17/97 by the version of
C ADIFOR compiled on Aug 21 1995.
C
C ADIFOR was prepared as an account of work sponsored by an
C agency of the United States Government, Rice University, and
C the University of Chicago. NEITHER THE AUTHOR(S), THE UNITED
C STATES GOVERNMENT NOR ANY AGENCY THEREOF, NOR RICE UNIVERSITY,
C NOR THE UNIVERSITY OF CHICAGO, INCLUDING ANY OF THEIR EMPLOYEES
C OR OFFICERS, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES
C ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY,
COMPLETE-
C NESS, OR USEFULNESS OF ANY INFORMATION OR PROCESS DISCLOSED, OR
C REPRESENTS THAT ITS USE WOULD NOT INFRINGE PRIVATELY OWNED
RIGHTS.
C
C
subroutine g_evapor(g_p_, t, upar, g_upar, rpar, evap,
g_evap)
C
C
C This subroutine computes the evaporation rate per unit area
C at time t after the spill of a (volatile) liquid dependent on
C geographical, climatological and meteorological data of the
C environment as well as the physico-chemical properties of the
C liquid, the surrounding air and the underlying ground.
C Ref.: P.I. Kawamura and D. Mackay,
C J. Hazardous Materials, 15 (1987) 343-364.
C
C INPUT:
C
C t ...... time elapsed since the spill [s]
C
C sa = upar(1) [rad]
C cloud = upar(2) [%]
C ubar10 = upar(3) [m/s]
C tempa = upar(4) [K]
C diamp = upar(5) [m]
C depthp = upar(6) [m]
C wtmol = rpar(1) [g/mol]
C dab = rpar(2) [m2/s]
C conliq = rpar(3) [J/msK]
C boilp = rpar(4) [oC]
C vapht = rpar(5) [J/mol]
C a = rpar(6)
C b = rpar(7)
C c = rpar(8)
C
C sa ...... solar altitude [rad]
C cloud ... cloud cover [%]
C ubar10 .. wind speed measured at 10 m [m/s]
C tempa ... ambient air temperature [K]
C diamp ... pool diameter or downwind length of pool [m]
C depthp .. initial depth of the pool [m]
C wtmol ... molecular weight of the liquid [g/mol]
C dab ..... diffusivity of the liquid in air [m2/s]
C conliq .. thermal conductivity of the liquid [J/msK]
C boilp ... boiling point of the liquid [\373C]
C vapht ... heat of vaporization (=latent heat of evaporation)
[J/mol]
C a ..... } coefficients for the vapor pressure equation of the
C b ......} form press = exp(a-b/(T+c)) [Pa]
C c ..... } where T ... temperature [K]
C
C OUTPUT:
C
C evap ..... evaporation rate [g/m2s]
C
C ENDDOC-------------------------------last revision:1988-01-13
mp
C implicit undefined (a-z)
C
real t, upar(*), rpar(*), evap
C
real sa, cloud, ubar10, tempa
real diamp, depthp, a, b, c
real wtmol, dab, conliq, boilp, vapht
real rgas, visair, congrd, difgrd, pi
real qsol, schmdt, trmass, hgrd
real hliqi, ugrd, er, press, ea, z, s, beta
C
C rgas ..... universal gas constant [Pam3/molK]
C visair ... kinematic viscosity of air [m2/s]
C congrd ... thermal conductivity of saturated sand [J/msK]
C difgrd ... thermal diffusivity of sand [m2/s]
C
integer g_pmax_
parameter (g_pmax_ = 6)
integer g_i_, g_p_
real r1_w, r2_p, r1_p, r9_b, r8_b, r7_b, r6_b, r5_b, r2_v,
r3_v
real r4_v, r5_v, r6_v, r7_v, r8_v, r4_b, r3_b, r2_b,
g_sa(g_pmax_), g_upar(g_pmax_, *)
real g_cloud(g_pmax_), g_ubar10(g_pmax_),
g_tempa(g_pmax_), g_diamp(g_pmax_), g_depthp(g_pmax_),
g_qsol(g_pmax_), g_trmass(g_pmax_), g_hliqi(g_pmax_),
g_ugrd(g_pmax_),g_er(g_pmax_)
real g_r1_w(g_pmax_), g_press(g_pmax_), g_ea(g_pmax_),
g_z(g_pmax_), g_s(g_pmax_), g_beta(g_pmax_),
g_evap(g_pmax_)
save g_er, g_r1_w, g_press, g_ea, g_z, g_s, g_beta
save g_sa, g_cloud, g_ubar10, g_tempa, g_diamp, g_depthp,
g_qsol, g_trmass, g_hliqi, g_ugrd
intrinsic real
data rgas /8.314/
data visair /0.000015/
data congrd /2.083333/
data difgrd /0.7e-6/
data pi /3.1415927/
C
integer g_ehfid
data g_ehfid /0/
C
call ehsfid(g_ehfid, \324evapor\325,\325g_evapor.f\325)
C
if (g_p_ .gt. g_pmax_) then
print *, \324Parameter g_p_ is greater than g_pmax_\325
stop
endif
do g_i_ = 1, g_p_
g_sa(g_i_) = g_upar(g_i_, 1)
enddo
sa = upar(1)
C--------
do g_i_ = 1, g_p_
g_cloud(g_i_) = g_upar(g_i_, 2)
enddo
cloud = upar(2)
C--------
do g_i_ = 1, g_p_
g_ubar10(g_i_) = g_upar(g_i_, 3)
enddo
ubar10 = upar(3)
C--------
do g_i_ = 1, g_p_
g_tempa(g_i_) = g_upar(g_i_, 4)
enddo
tempa = upar(4)
C--------
do g_i_ = 1, g_p_
g_diamp(g_i_) = g_upar(g_i_, 5)
enddo
diamp = upar(5)
C--------
do g_i_ = 1, g_p_
g_depthp(g_i_) = g_upar(g_i_, 6)
enddo
depthp = upar(6)
C--------
wtmol = rpar(1)
dab = rpar(2)
conliq = rpar(3)
boilp = rpar(4)
vapht = rpar(5)
a = rpar(6)
b = rpar(7)
c = rpar(8)
C
C *** Preliminary calculations ***
C
C Compute net solar insulation [J/m2s]:
C Ref.: J.M. Raphael,
C Prediction of temperature in rivers and reservoirs,
C Proc. Amer. Soc. Civ. Eng., J.Power Div., 88 P02, 1962.
C
r2_v = 0.000071 * cloud
r5_v = 1111.1111 * (1. - r2_v * cloud)
r7_v = sin(sa)
r1_p = cos(sa)
r8_v = r7_v - 0.1
r5_b = r5_v * r1_p
r7_b = -(r8_v * 1111.1111)
r9_b = r7_b * r2_v + r7_b * cloud * 0.000071
do g_i_ = 1, g_p_
g_qsol(g_i_) = r9_b * g_cloud(g_i_) + r5_b * g_sa(g_i_)
enddo
qsol = r5_v * r8_v
C--------
if (qsol .lt. 0.) then
do g_i_ = 1, g_p_
g_qsol(g_i_) = 0.0
enddo
qsol = 0.
C--------
endif
C
C Compute mass transfer coefficient [m/s]:
C Ref.: D. Mackay and R.S. Matsugu,
C Evaporation rate of hydrocarbon spills on water and land,
C Can. J. Chem. Eng., 5 (1973) 434.
C
schmdt = visair / dab
if ( ubar10 .ne. 0.0e0 ) then
r2_v = ubar10 ** ( 0.78 - 2.0e0)
r2_v = r2_v * ubar10
r2_p = 0.78 * r2_v
r2_v = r2_v * ubar10
else
r2_v = ubar10 ** 0.78
if ( (0.0e0 .lt. 0.78) .and. ( 0.78 .lt. 1.0e0) ) then
call ehbfSO (10,ubar10, 0.78, r2_v, r2_p, 0.0, +g_ehfid,+197)
else
r2_p = 0.0e0
endif
endif
if ( diamp .ne. 0.0e0 ) then
r5_v = diamp ** ( 0.11 - 2.0e0)
r5_v = r5_v * diamp
r1_p = 0.11 * r5_v
r5_v = r5_v * diamp
else
r5_v = diamp ** 0.11
if ( (0.0e0 .lt. 0.11) .and. ( 0.11 .lt. 1.0e0) )
then
call ehbfSO (10,diamp, 0.11, r5_v, r1_p, 0.0,
+g_ehfid, +212)
else
r1_p = 0.0e0
endif
endif
r6_v = r5_v * schmdt ** 0.67
r7_v = 0.0047864 * r2_v / r6_v
r5_b = -r7_v / r6_v * schmdt ** 0.67 * r1_p
r7_b = 1.0 / r6_v * 0.0047864 * r2_p
do g_i_ = 1, g_p_
g_trmass(g_i_) = r7_b * g_ubar10(g_i_) + r5_b
* g_diamp(g_i_)
enddo
trmass = r7_v
C--------
C
C Compute ground heat transfer coefficients:
C
C 1) time dependent heat transfer coefficient at the
C liquid-ground interface [J/m2sK]:
C Ref.: J.P. Holman, Heat Transfer, McGraw Hill, New York, 1976.
C
hgrd = 2. * congrd / sqrt(pi * difgrd * t)
C
C 2) (invers of) liquid heat transfer coefficient [J/m2sK]:
C
r2_b = 1.0 / (conliq * (1. + exp(-0.06 * (boilp - 70.))))
do g_i_ = 1, g_p_
g_hliqi(g_i_) = r2_b * g_depthp(g_i_)
enddo
hliqi = depthp / (conliq * (1. + exp(-0.06 * (boilp - 70.))))
C--------
C
C 3) overall heat transfer coefficient [J/m2sK]:
C
r2_v = 1. / hgrd + hliqi
r3_v = 1. / r2_v
r3_b = -r3_v / r2_v
do g_i_ = 1, g_p_
g_ugrd(g_i_) = r3_b * g_hliqi(g_i_)
enddo
ugrd = r3_v
C--------
C
C *** Direct evaporation model ***
C
C Evaporation due to solar insulation:
C
r3_b = 1.0 / vapht * wtmol
do g_i_ = 1, g_p_
g_er(g_i_) = r3_b * g_qsol(g_i_)
enddo
er = qsol * wtmol / vapht
C--------
C
C Liquid evaporation rate:
C
r3_v = tempa - real(273) + c
r4_v = b / r3_v
r5_b = -(-r4_v / r3_v)
do g_i_ = 1, g_p_
g_r1_w(g_i_) = r5_b * g_tempa(g_i_)
enddo
r1_w = a - r4_v
r2_v = exp(r1_w)
r1_p = r2_v
do g_i_ = 1, g_p_
g_press(g_i_) = r1_p * g_r1_w(g_i_)
enddo
press = r2_v
C--------
r2_v = trmass * wtmol
r6_v = rgas * tempa
r7_v = r2_v * press / r6_v
r2_b = 1.0 / r6_v
r4_b = -r7_v / r6_v * rgas
r6_b = r2_b * r2_v
r7_b = r2_b * press * wtmol
do g_i_ = 1, g_p_
g_ea(g_i_) = r7_b * g_trmass(g_i_) + r6_b * g_press(g_i_)
+ r4 *_b * g_tempa(g_i_)
enddo
ea = r7_v
C--------
C
C Compute factor z and slope of vapor pressure curve s [Pa/K]:
C
r2_v = ugrd * rgas
r6_v = r2_v * tempa / trmass
r3_b = 1.0 / vapht
r4_b = r3_b * (1.0 / trmass)
r5_b = r3_b * (-r6_v / trmass)
r7_b = r4_b * r2_v
r8_b = r4_b * tempa * rgas
do g_i_ = 1, g_p_
g_z(g_i_) = r8_b * g_ugrd(g_i_) + r7_b * g_tempa(g_i_)
+ r5_b ** g_trmass(g_i_)
enddo
z = (3650120. * schmdt ** 0.67 + r6_v) / vapht
C--------
r6_v = (tempa - real(273) + c) * (tempa - real(273) + c)
r1_p = 2.0e0 * (tempa - real(273) + c)
r7_v = press * b / r6_v
r6_b = -r7_v / r6_v * r1_p
r7_b = 1.0 / r6_v * b
do g_i_ = 1, g_p_
g_s(g_i_) = r7_b * g_press(g_i_) + r6_b * g_tempa(g_i_)
enddo
s = r7_v
C--------
C
C Finally compute evaporation rate [g/m2s]:
C
r3_v = s + z
r4_v = s / r3_v
r3_b = -r4_v / r3_v
r2_b = 1.0 / r3_v + r3_b
do g_i_ = 1, g_p_
g_beta(g_i_) = r2_b * g_s(g_i_) + r3_b * g_z(g_i_)
enddo
beta = r4_v
C--------
r5_v = 1. - beta
r6_b = -ea + er
do g_i_ = 1, g_p_
g_evap(g_i_) = beta * g_er(g_i_) + r6_b * g_beta(g_i_)
+ r5_v ** g_ea(g_i_)
enddo
evap = er * beta + ea * r5_v
C--------
return
end
|
© Copyright 1995-2002 by: ESS Environmental Software and Services GmbH AUSTRIA