!define erf
!define fourbyt
!define ifortc
!define in32
!define linux
!define ma28
!define selap
!define svnew
!define lic0
       module stgd2omod
!  Contains storage and subroutines for stgh2o which computes a table of
!  thermodynamic properties for heavy water as a function of temperature
!  and pressure.
!
       implicit none
       integer, parameter, private :: kr=selected_real_kind(12,200)
!
c      common /bend9/ a1,a2,a3,a4,a5
       real(kr) aa1,aa2,aa3,aa4,aa5
       common /check1/ ni
       real(kr) ni
c      common /checks/ dch1,dch2,pch1,pch2,pch3,tch1,tch2,tch3,dst,tst,
c    & hsch1,hsch2
       real(kr) dch1,dch2,pch1,pch2,pch3,tch1,tch2,tch3,dst,tst,
     & hsch1,hsch2
       data dch1,dch2,pch1,pch2,pch3,tch1,tch2,tch3,dst,tst,
     & hsch1,hsch2/0.d0,1.16d0,.000611d0,21.66d0,100.d0,.33333333d0,1.55
     #3d0,
     & 3.6108d0,.8d0,400.d0,0.d0,6052.2d0/
c      common /code/ messag
       character messag*128
       data messag /' Thermodynamic and transport properties for heavy w
     &ater pc=213.83atm,tc= 643.89 k,roc=.358 g/cc'/
c      common /cof/ a(10,7)
       real(kr) a(10,7)
       real(kr) a1(10),a2(10),a3(10),a4(10),a5(10),a6(10),a7(10)
       equivalence (a1(1),a(1,1)),(a2(1),a(1,2)),(a3(1),a(1,3)),
     & (a4(1),a(1,4)),(a5(1),a(1,5)),(a6(1),a(1,6)),(a7(1),a(1,7))
       data a1/73.13848592d0,-285.20415917d0,535.71659288d0,-649.8100061
     #4d0,
     & 574.63280680d0,-387.92157774d0,206.34569512d0,-79.89428513d0,
     & -996.36169097d0,-766.272900006d0/
       data a2 /24.74108348d0,-105.57317181d0,200.87302906d0,-235.187764
     #40d0,
     & 224.56976938d0,-40.09924297d0,128.77154771d0,-28.40907978d0,
     & -1389.08003142d0,-1672.09705556d0/
       data a3/11.64775625d0,-42.51820251d0,72.45541064d0,-82.55391089d0
     #,0.d0,
     & 0.d0,0.d0,0.d0,-267.85482520d0,-998.64982710d0/
       data a4/2.66566642d0,-9.19657655d0,15.13096920d0,-7.24860975d0,0.
     #d0,
     & 0.d0,0.d0,0.d0,-46.83904320d0,-227.34793319d0/
       data a5/-6.73408249d0,24.03602093d0,-41.08079830d0,45.39111005d0,
     #0.d0,
     & 0.d0,0.d0,0.d0,139.21659329d0,566.02305152d0/
       data a6/-5.24802962d0,18.52690633d0,-31.42397369d0,26.43208802d0,
     #0.d0,
     & 0.d0,0.d0,0.d0,96.31411481d0,453.20280933d0/
       data a7/-1.17583447d0,4.13816432d0,-6.55842224d0,4.75774631d0,0.d
     #0,0.d0,0.d0,
     & 0.d0,19.39184297d0,103.56819758d0/
c      common /consts/ tauc,rhoa,rhob,taua,e,r
       real(kr) tauc,rhoa,rhob,taua,e,r
       data tauc,rhoa,rhob,taua,e,r/1.553d0,0.7d0,1.1d0,2.53d0,4.3d0,0.4
     #1515d0/
c      common /conv1/ dconv(5)
       real(kr) dconv(5)
       data dconv/2*1.d0,62.4283d0,1000.d0,1.d0/
c      common /conv2/ tconv(5)
       real(kr) tconv(5)
       data tconv/1.d0,1.d0,.55555555d0,1.d0,1.d0/
c      common /conv3/ pconv(5)
       real(kr) pconv(5)
       data pconv/1.d0,9.8692327d0,145.038243d0,1.d6,1.d0/
c      common/conv4/ sconv(5)
       real(kr) sconv(5)
       data sconv/2*1.d0,.238849d0,1000.d0,1.d0/
c      common/conv5/ cconv(5)
       real(kr) cconv(5)
       data cconv/2*1.d0,.0328084d0,.01d0,1.d0/
c      common/conv6/ hconv(5)
       real(kr) hconv(5)
       data hconv/2*1.d0,.429929d0,1000.d0,1.d0/
c      common/conv7/ mconv(5)
       real(kr) mconv(5)
       data mconv/2*1.d0,.67196899d-1,.1d0,1.d0/
c      common/conv8/ kconv(5)
       real(kr) kconv(5)
       data kconv/2*1.d-2,1.60644d-4,100.d0,1.d0/
c      common/conv9/ stconv(5)
       real(kr) stconv(5)
       data stconv/2*1.d0,6.8521766d-5,.001d0,1.d0/
c      common/conv10/ alconv(5)
       real(kr) alconv(5)
       data alconv/2*.1d0,.328084d-2,.001d0,1.d0/
c      common /cosat/ cps1,cps2,cps3,cps4,cps5
       real(kr) cps1,cps2,cps3,cps4,cps5
       data cps1,cps2,cps3,cps4,cps5/
     & -7.81583d0,17.6012d0,-18.1747d0,-3.92488d0,4.19174d0/
!
c      common /crit/ rhocrt,pcrt,tcrt
       real(kr) rhocrt,pcrt,tcrt
       data rhocrt/0.358d0/,pcrt/21.66d0/,tcrt/643.89d0/
!
c      common /deriv/ pdt,ptv,pdtl,pdtv,ptvl,ptvv
       real(kr) pdt,ptv,pdtl,pdtv,ptvl,ptvv
!  pdt,pdtl,pdtv = derivative of pressure with respect to density
!                  at constant temperature.
!  ptv,ptvl,ptvv = derivative of pressure with respect to
!                  temperature at constant volume.
c      common /ierror/ irout
       integer irout
c      common /iterat/ tr,t0,t1,t2,t3,t4,t5,t6,t7,t8
       real(kr) tr,t0,t1,t2,t3,t4,t5,t6,t7,t8
c      common /minmax/ tmin,pmin,tmax,pmax
       real(kr) tmin,pmin,tmax,pmax
       data tmin/2.7695d+02/, pmin/6.60114d+02/
       data tmax/5.0d+03/, pmax/1.0d+08/
c      common /partls/ ptv,pdt
       real(kr) ptv1,pdt1
       real(kr) tdegf,tdegc
       data tdegf/0.d0/, tdegc/0.d0/
c      common /phcall/ ps,hs,ss
       real(kr) pszc,hszc,sszc
c      common /prhot/ ps,ds,tt
       real(kr) psz,dsz,ttz
c      common /prhott/ ps,ds,tt
       real(kr) pszz,dszz,tszz
c      common /propty/ c,cl,cp,cpl,cpv,cv,cvl,cvp,cvv,dl,dv,excesk,
c    & excl,excv,gamma,gammal,gammav,hl,hv,k,kl,kv,mu,mul,muv,s,sigma,
c    & sl,sv,ku
       real(kr) c,cl,cp,cpl,cpv,cv,cvl,cvp,cvv,dl,dv,excesk,excl,excv,
     & gamma,gammal,gammav,hl,hv,k,kl,kv,mu,mul,muv,s,sigma,sl,sv
       integer ku
c      common /qaux/ rbdif(8),radif(8),er,ed,tadif(7)
       real(kr) rbdif(8),radif(8),er,ed,tadif(7)
c      common /qs1/ sumi(7)
       real(kr) sumi(7)
c      common /sicof/ sic1,sic2,sic3,sic4,sic5
       real(kr) sic1,sic2,sic3,sic4,sic5
       data sic1,sic2,sic3,sic4,sic5/1834.5669d0,3.2661763d0,-0.00043593
     #99d0,
     & 15.4389d0,-0.9758409d0/
c      common /tparam/ ts
       real(kr) tdum
c      common /trip/ ttrip,ptrip,vtrip
       real(kr) ttrip,ptrip,vtrip
       data ttrip/2.7695d+2/, ptrip/6.60114d+2/, vtrip/9.04592d-4/
c      common /xminus/ xm1(7)
       real(kr) xm1(7)
       data xm1/1.d0,2.d0,3.d0,4.d0,5.d0,6.d0,7.d0/
!
c      character*11 p,state
c      character*10 day
c      character*8 timeod
c
c      data  next /1/, maxlpp /55/
!
       contains
!
       function check(ku,kz,t)
       implicit none
       integer, parameter :: kr=selected_real_kind(12,200)
!
       real(kr) check,dcheck,pcheck,tcheck
       integer ku,kz
       real(kr) d,p,t
!
       real*8 ch1,ch2,ch3
       integer iflag,kode
       character rout(11)*8
       data rout/'dens  ','press','temp','enth ','ent ','tempph',
     & 'tempps','cpprl','visc','therm','surf'/
!
!  Convert temperature t to degrees Kelvin and check for out of range.
!  Units are specified by ku. If kz is specified as 1, t is checked for
!  out of saturation range.
       entry tcheck(ku,kz,t)
       check = 1000.d0*tconv(ku)/t
       ch1 = 1000.d0/tch3
       ch2 = 1000.d0/tch2
       ch3 = 1000.d0/tch1
       kode = 1
    1  iflag = 1
       go to 10
!
!  Convert pressure to mn/m**2 and check for out of range. Units are
!  specified by ku. If kz is specified as 1, p is checked for out of
!  saturation range,
       entry  pcheck(ku,kz,p)
       check = p/pconv(ku)
       ch1 = pch1
       ch2 = pch2
       ch3 = pch3
       kode = 0
    2  iflag = 2
       go to 10
!
!  Convert density to g/cc and check for out of range. Units are
!  specified by ku.
       entry dcheck(ku,d)
       check = d/dconv(ku)
       ch1 = dch1
       ch3 = dch2
       kode = 0
    3  iflag = 3
       go to 20
   10  if (kz .eq. 1) go to 30
   20  if (check .lt. ch1) go to 40
       if (check .gt. ch3) go to 40
   25  if (kode.eq.1) check = t/tconv(ku)
       return
   30  if (check .lt. ch1) go to 40
       if (check .le. ch2) go to 25
   40  continue
       if (iflag .eq. 1) then
         write(*,800) check,'t',rout(irout)
       elseif (iflag .eq. 2) then
         write(*,800) check,'p',rout(irout)
       elseif (iflag .eq. 3) then
         write(*,800) check,'d',rout(irout)
       endif
  800  format (' ',g12.4,' is out of range for ',a,' in sub.-',a6)
       go to 25
       end function check
!
       subroutine cpprl (ku,t,d,cp,cv,gamma,c)
!  This subroutine returns the following to wasp in users units.
!    specific heat at constant pressure  =cp;
!    specific heat at constant volume    =cv;
!    specific heat ratio                 =gamma;
!    sonic velocity                      = c.
!  The partials ptv and pdt explained below are returned in common.
!
       implicit none
       integer ku
       real(kr) t,d,cp,cv,gamma,c
!
       real(kr) cq2dt,cq2d2t,cq2t2d,cq2t2t,cqcalc,cqdt,cqtd,cs,dhdd,
     & dhdt,ds,gammap,ts,tt
       integer kode,kz
!
       irout = 8
       ts = tcheck(ku,kz,t)
       tt = 1000.d0/ts
       ds = dcheck(ku,d)
       call qmust (ds)
       call qmust2 (ts)
       cq2t2d = q2t2d(ts)
       cqcalc = qcalc(ts)
       cqdt = qdt(ds,ts)
       cqtd = qtd(ts)
       cq2dt = q2dt(ds,ts)
       cq2d2t = q2d2t(ds,ts)
       cv = -2.d0*sic3*tt + sic4/tt - sic5 - r*ds*ts*ts*cq2t2d
!  ptv is partial of p by t (not tau).
!  pdt is partial of p by rho.
       ptv1 = r*ds*(1.d0+ ds*(cqcalc + ds*cqdt - ts*(cqtd + ds*cq2dt)))
       pdt1 = r*tt*(1.d0+ ds*(2.d0*cqcalc + ds*(4.d0*cqdt + ds*cq2d2t)))
       dhdt = -2.d0*sic3*tt + sic4/tt - sic5+r*(1.d0+
     & ds*(cqcalc + ds*cqdt - ts*(ds*cq2dt + cqtd + ts*cq2t2d)))
       dhdd = r*(tt*cqcalc + 1000.d0*cqtd +
     & ds*(tt*(3.d0*cqdt + ds*cq2d2t) +1000.d0*cq2dt))
       cp = dhdt - dhdd*(ptv1/pdt1)
       gamma = cp/cv
       cp = cp*sconv(ku)
       cv = cv*sconv(ku)
       gammap = gamma*10.d0*pdt1
       cs = 0.0d0
       if (gammap .gt. 0.0d0) cs = 1000.d0*sqrt(gammap)
       c = cs*cconv(ku)
       return
       end subroutine cpprl
!
       subroutine dens (ku,t,p,d,dl,dv,kz)
!  Compute density d given temperature t and pressure p.
!  Units are specified by ku. If kz is returned or specified as 1, the
!  saturated liquid and vapor densities, dl and dv respectively, are
!  computed as a function of t or p. The other value must be input as 0.
!
       implicit none
       integer ku,kz
       real(kr) t,p,d,dl,dv
!
       real(kr) dsincr,dsl,dsv,est,est1,est2,pss,rrr,test
!
       irout = 1
       if (kz .ne. 1) then
         tdum = tcheck(ku,kz,t)
         ttz = tdum
         call qmust2 (tdum)
       elseif (t .le. 0.0d0) then
         psz = pcheck(ku,kz,p)
         tdum = tss(psz)
         ttz = tdum
         if (t .le. 0.d0) t = tdum*tconv(ku)
         call qmust2 (tdum)
       else
         tdum = tcheck(ku,kz,t)
         ttz = tdum
         call qmust2 (tdum)
         call psss(psz)
         if (p. le.0 .d0) p = psz*pconv(ku)
       endif
!
!  Determine region.
       if (kz .eq. 1) go to 41
       psz = pcheck(ku,kz,p)
       if (psz .le. pch2) go to 110
       if (tdum .le. tch2) go to 130
       kz = 2
       est1 = dch2
       est2 = rhocrt
       test = 1000.d0/tdum - 273.15d0
       if (test .gt. 101.0d0) est1 = 1.122d0
       if (test .le. 180.0d0) then
         if (test .lt. 40.d0) test = 40.0d0
         est2 = (test*(test*(test*(test*.12476711d-09 - .52277795d-07) +
     &   0.54790571d-05) - 0.69617325d-03) + 1.0220277d0)
       endif
       call root (est1,est2,0.0d0,dsf,dsz)
       go to 150
  130  kz = 3
       est = rhocrt*3.d0
       call root (est,dch1,0.00d0,dsf,dsz)
       go to 150
  110  if (tdum .le. tch2) go to 50
       call psss (pss)
       if (abs((pss - psz)/pss) .lt. 1.d-4) go to 60
       if (psz .eq. pss) go to 60
       if (psz .lt. pss) go to 50
       kz = 2
   41  dsincr = 0.1d0
       dsz = dch2 + dsincr
       go to 81
   50  kz = 3
       dsz = psz*tdum/(1000.d0*r)
       go to 90
!
!  Region 1.
   60  kz = 1
       go to 41
   81  do
         dsz = dsz - dsincr
         dsl = solve(dsz,dsf,ddsf)
         if (dsl .gt. dch2) cycle
         if (dsl .ge. dch1) exit
         dsz = dsz + dsincr
         dsincr = 0.5d0*dsincr
       enddo
       if (kz .ne. 2) then
!        if (dsl.lt.rhocrt .and. 1000./tdum.ge.647.08)
!      & dsl =.118*sqrt(abs((647.36 - 1000./tdum)**2 - .0041)) + .317 +
!        2.64e-4
         dsz = psz*tdum/(1000.d0*r)
         if (tch2/tdum .gt. .985d0) dsz = .65d0*rhocrt
         if (tch2/tdum .gt. .995d0) dsz = .75d0*rhocrt
         if (tch2/tdum .gt. .999d0) dsz = .85d0*rhocrt
         if (tch2/tdum .gt. .9995d0) dsz = .90d0*rhocrt
         dsv = solve(dsz,dsf,ddsf)
!  The correction added to dsl can be added to dsv; however the vapor
!  locus displays a small cusp which matches current theory and so the
!  change will not be added at this time.
!  Previous coefficients for dl  and dv  were -.0036 instead of .0041
!  and -.098 and .0225 instead of -.098 and .0269.
!  The former were based on 647.3=tc and now on 647.286=tc.
!  The additions to dsl and dsv  (2.64e-4and1.2e-3) to match wasp
!  additions of dsl and dsv (2.64e-4&1.2e-3) to match wasp; .098 to .11.
         if ((1000.d0/tdum) .gt. 647.15d0) dsv = -.11d0*
     &   sqrt(abs((647.45d0- 1000.d0/tdum)**2 - .0269d0)) + .317d0+ 1.2d
     #-3
!
         dl = dsl*dconv(ku)
         dv = dsv*dconv(ku)
         rrr = .9d0*rhocrt
         if (dsv .lt. 0.d0) call root (rrr,.0100d0,0.00d0,dsf,dsv)
         return
       endif
!
!  Regions 2 and 3.
   90  if (kz .eq. 2) then
         dsz = dsl
       else
         dsz = solve(dsz,dsf,ddsf)
       endif
  150  d = dsz*dconv(ku)
       return
       end subroutine dens
!
       function dsf(d)
!  Function used to solve for density d  given temperature and pressure.
!
       implicit none
       real(kr) dsf,ddsf
       real(kr) d
!
       real(kr) pstate
!
       call qmust (d)
       pstate = 1000.d0*r*d/ttz*(1.d0+ d*(qcalc(ttz) + d*qdt(d,ttz)))
       dsf = pstate - psz
       return
!
       entry ddsf(d)
       call qmust (d)
       dsf = 1000.d0*r/ttz*(1.d0+ d*(2.0d0*qcalc(ttz) +
     & 4.0d0*d*qdt(d,ttz) + d*d*q2d2t(d,ttz)))
       return
       end function dsf
!
       function dtherm(xlm)
!  Function used to solve eq.(b52) for thermal conductivity.
!
       implicit none
       real(kr) dtherm,ddth
       real(kr) xlm
!
       real(kr) dth,tcalc
!
       tcalc = (xlm*(xlm*(xlm*(xlm*(xlm*(xlm*(xlm*t8 + t7) + t6) + t5) +
     & t4) + t3) + t2) + t1)*xlm + t0
       dtherm = tcalc - tr
       return
!
       entry ddth(xlm)
!  Derivative used to solve eq.(b52) in Newton Raphson iteration.
!
       dth = (xlm*(xlm*(xlm*(xlm*(xlm*(xlm*8.d0*t8 + 7.d0*t7) + 6.d0*t6)
     # +
     & 5.d0*t5) + 4.d0*t4) + 3.d0*t3) + 2.d0*t2)*xlm + t1
       dtherm = dth
       return
       end function dtherm
!
       subroutine ent (ku,tt,d,s)
!  This routine computes entropy given the temperature parameter tt
!  and the density d.  i/o units are specified by ku.
!  If saturation values are needed, this routine must be called twice
!  with dl and dv input as d.
!  Entropy is returned in s.
!
       implicit none
       integer ku
       real(kr) tt,d,s
!
       real(kr) ds,psit,sss,t,ts
       integer kz
!
       irout = 5
       ts = tcheck(ku,kz,tt)
       ds = dcheck(ku,d)
       call qmust (ds)
       call qmust2 (ts)
       t = 1000.d0/ts
       psit = 2.d0*sic3*t + sic2 + sic4/t + sic5*(1.d0+ log(t))
       sss= - r*(log(ds) + ds*(qcalc(ts) - ts*qtd(ts))) - psit
       s = sss*sconv(ku)
       return
       end subroutine ent
!
       subroutine enth (ku,tt,d,h)
!  This routine computes enthalpy given the temperature parameter tt
!  and the density d. i/o units are specified by ku.
!  If saturation values are needed,this routine must be called twice
!  with dl and dv input as d.
!  Enthalpy is returned in h.
!
       implicit none
       integer ku
       real(kr) ts,tt,d,h
!
       real(kr) ds,h1,h2,psio,psit,t,tcons
       integer ktr,kz
!
       irout = 4
       ts = tcheck(ku,kz,tt)
       ds = dcheck(ku,d)
       call qmust (ds)
       call qmust2 (ts)
       t = 1000.d0/ts
       psio = (sic3*t + sic2)*t + sic1 + (sic4 + sic5*t)*log(t)
       psit = 2.d0*sic3*t + sic2 + sic4/t + sic5*(1.d0+ log(t))
       h1 = psio - t*psit
       h2 = 1000.d0*r/ts*(1.d0+ ds*(qcalc(ts) + ts*qtd(ts) +
     & ds*qdt(ds,ts)))
       h = (h1 + h2)*hconv(ku)
       return
       end subroutine enth
!
       subroutine hwasp (ks,kp,tt,p,d,h,kz)
!  Changed to hwasp, April 1988 J. E. Tolli, EG&g Idaho in accordance
!  with Hill-Macmillan-Lee fundamental equation for heavy water (same as
!  Keyes, et.al. equation of state cited above, except for constant
!  coefficients);  changes to original wasp routines are as follows:
!
       implicit none
       integer ks,kp,kz
       real(kr) tt,p,d,h
!
       real(kr) ds,t,tcons
       integer ktr,ktra
!
       real(kr) tcon(5)
       data tcon/1.d0,1.d0,1.8d0,1.d0,1.d0/
!
!        psss,tss,tssf:
!  Replace saturation equation with equation recommended by hill, et.al.
!  for heavy water.
!
!        blckdata:
!  Change fundamental equation coefficients, replace saturation equation
!  coefficients with values required for new saturation equation, change
!  other values (maximum and minimium values, critical point values,
!  etc.) as required.
!
!        dens:
!  Modify logic somewhat for solving liquid state density. Compute the
!  state relations and thermodynamic and transport properties of water
!  given temperature tt, pressure p, density d, or enthalpy h, or
!  entropy s. State relations are specified by ks. Thermodynamic and
!  transport properties are specified by kp. If kz is returned or
!  specified as 1, properties are computed at saturation.
!  tau is the temperature parameter used in the equation of state.
!  tau is equivalent to t in this subroutine.
!
       ktr = 0
       tcons = 0.0d0
       if (tdegf .gt. .9d0) tcons = 459.67d0
       if (tdegc .gt. .9d0) tcons = 273.15d0
       t = tt
       if (tt .gt. 0.d0)  t = 1000.d0/tt
       if (tcons.gt.0.0d0.and. tt.ne.0.0d0) t = 1000.d0/(tt+tcons)
       go to (10,20,30,40,45), ks
!  Compute density.
   10  call dens (ku,t,p,d,dl,dv,kz)
       if (tt .eq. 0.0d0) tt = 1000.d0/t - tcons
       go to 50
!  Compute pressure.
   20  call press (ku,t,d,p,kz)
       go to 50
!  Compute temperature.
   30  call temp (ku,p,d,t,kz)
       tt = 1000.d0/t -tcons
       go to 50
!  Compute temperature and density given pressure and enthalpy.
   40  call tempph (ku,p,h,t,d,dl,dv,kz)
       tt = 1000.d0/t - tcons
       go to 50
   45  call tempps (ku,p,s,t,d,dl,dv,kz)
       tt = 1000.d0/t - tcons
!
   50  if (kz.ne.1 .or. (ks.eq.1 .or. ks.gt.3)) go to 55
!  Obtain saturation densities dl and dv for ks=2 and ks=3 calls when
!      kz = 1
       call dens (ku,t,p,d,dl,dv,1)
!
 55    continue
       ktra = kp/2**ktr
       if (ktra .eq. 0) return
       ktr = ktr + 1
       if (mod(ktra,2) .eq. 0) goto 55
       go to (60,100,130,160,180,240), ktr
!
!  Compute enthalpy.
   60  if (kz .eq. 1) go to 65
       call enth (ku,t,d,h)
       go to 70
   65  call enth (ku,t,dl,hl)
       call enth (ku,t,dv,hv)
   70  go to 55
!
!  Compute entropy.
  100  if (kz .eq. 1) go to 105
       call ent (ku,t,d,s)
       go to 110
  105  call ent (ku,t,dl,sl)
       call ent (ku,t,dv,sv)
  110  go to 55
!
!  Compute specific heats and gamma and sonic velocity.
  130  if (kz .ne. 1) go to 135
       call cpprl (ku,t,dl,cpl,cvl,gammal,cl)
       ptvl = ptv1*pconv(ku)/tcon(ku)
       pdtl = pdt1*pconv(ku)/dconv(ku)
       call cpprl (ku,t,dv,cpv,cvv,gammav,cvp)
       pdtv = pdt1*pconv(ku)/dconv(ku)
       ptvv = ptv1*pconv(ku)/tcon (ku)
       go to 140
  135  call cpprl(ku,t,d,cp,cv,gamma,c)
       ptv = ptv1*pconv(ku)/tcon (ku)
       pdt = pdt1*pconv(ku)/dconv(ku)
 140   go to 55
!
!  Compute viscosity.
  160  if (kz .ne. 1) go to 165
       call visc (ku,kz,t,p,dl,mul)
       call visc (ku,kz,t,p,dv,muv)
       go to 170
  165  call visc (ku,kz,t,p,d,mu)
  170  go to 55
!
!  Compute thermal conductivity.
  180  if (kz .ne. 1) go to 220
       call therm (ku,kz,p,t,dl,excl,kl)
       call therm (ku,kz,p,t,dv,excv,kv)
       go to 190
  220  call therm (ku,kz,p,t,d,excesk,k)
  190  go to 55
!  Compute surface tension.
  240  call surf (ku,kz,t,sigma)
       return
       end subroutine hwasp
!
       subroutine press (ku,t,d,p,kz)
!  Compute pressure p given temperature t and density d.
!  Units are specified by ku.  if kz is returned or specified as 1, p is
!  computed at saturation as a function of t only.
!
       implicit none
       integer ku,kz
       real(kr) t,d,p
!
       real(kr) ds,dsl,dsv,ps,ts,ze
!
       irout = 2
       tdum = tcheck(ku,kz,t)
!
!  Determine region.
       if (kz - 1) 10,70,10
   10  ds = dcheck(ku,d)
       if ( kz .gt. 1) go to 80
       if (tdum - tch2) 50,50,20
   20  call dens (1,tdum,ze,ze,dsl,dsv,1)
       if (ds - dsl) 30,60,40
   30  if (ds - dsv) 50,60,60
   40  kz = 2
       go to 80
   50  kz = 3
       go to 80
!
!  Region 1.
   60  kz = 1
   70  call psss(ps)
       go to 90
!
!  Regions 2 and 3.
   80  call qmust (ds)
       call qmust2 (tdum)
       ps = 1000.d0*r*ds/tdum*(1.d0+ ds*(qcalc(ts) + ds*qdt(ds,tdum)))
   90  p = ps*pconv(ku)
       return
       end subroutine press
!
       subroutine psss (pss)
!
!  Compute saturation pressure pss in mn/m**2 as a function of
!  temperature in Kelvin degrees.
!
!  J. E. Tolli, EG&G Idaho Inc., 4/88.
!
!  Get temperature and adjusted tau value.
!
       implicit none
       real(kr) pss
!
       real(kr) atval,psrlog,tmp
!
       tmp = 1000.d0/tdum
       atval = 1.d0- tmp/tcrt
!
!  Calculate logarithm of psat/pcrt.
       if (abs(atval) .gt. 0.00001d0) then
         psrlog = tcrt/tmp*(cps1*atval + cps2*atval**1.9d0+
     &   cps3*atval**2 + cps4*atval**5.5d0+ cps5*atval**10)
       else
         psrlog = 0.0d0
       endif
!
!  Place saturation pressure in pss.
       pss = pcrt*exp(psrlog)
       return
       end subroutine psss
!
       function q2d2ta(tau)
!  Partial der of q  p2q/prho2
!
       implicit none
       real(kr) q2d2ta,q2d2t
       real(kr) d,tau
!
       real(kr) tsum,tsum1,tsum2
       integer i,j
!
    1  tsum = 0.0d0
       do j = 2,7
         tsum = tsum + tadif(j)*sumi(j)
       enddo
       q2d2ta = sumi(1) + (tau -tauc)*tsum
       return
       entry q2d2t(d,tau)
       sumi(1) = 0.0d0
       do i = 3,8
         sumi(1) = sumi(1) + xm1(i-1)*xm1(i-2)*a(i,1)*radif(i-2)
       enddo
       sumi(1) = sumi(1) + er*(-e*a(10,1)*(2.0d0- ed) + e*e*a(9,1))
       do j = 2,7
         sumi(j) = 0.0d0
         do i = 3,8
           sumi(j) = sumi(j) + xm1(i-1)*xm1(i-2)*a(i,j)*rbdif(i-2)
         enddo
         sumi(j) = sumi(j) + er*(-e*a(10,j)*(2.0d0- ed)+e*e*a(9,j))
       enddo
       go to 1
       end function q2d2ta
!
       function q2dta(tau)
!  Partial der of q, p2q/prho-ptau.
!
       implicit none
       real(kr) q2dta,q2dt
       real(kr) d,tau
!
       real(kr) sum(6),tsum1,tsum2
       integer i,j
!
    1  tsum1 = 0.0d0
       tsum2 = 0.0d0
       do j = 3,7
         tsum1 = tsum1 + xm1(j-2)*tadif(j-1)*sum(j-1)
         tsum2 = tsum2 + tadif(j)*sum(j-1)
       enddo
       tsum2 = tsum2 + sum(1)
       tsum1 = tsum1*(tau - tauc)
       q2dta = tsum1 + tsum2
       return
       entry q2dt(d,tau)
       do j = 2,7
         sum(j-1) = 0.d0
         do i = 2,8
           sum(j-1) = sum(j-1) + xm1(i-1)*a(i,j)*rbdif(i-1)
         enddo
         sum(j-1) = sum(j-1) + er*(a(10,j) - e*(a(9,j) + a(10,j)*d))
       enddo
       go to 1
       end function q2dta
!
       function q2t2d(tau)
!  Partial der of q, p2q/ptau2.
!
       implicit none
       real(kr) q2t2d
       real(kr) tau
!
       real(kr) tsum1,tsum2
       integer j
!
       tsum1 = 0.0d0
       tsum2 = 0.0d0
       do j = 3,7
         tsum1 = tsum1 + xm1(j-2)*tadif(j-1)*sumi(j)
         if (j.eq.3) cycle
         tsum2 = tsum2 + xm1(j-2)*xm1(j-3)*tadif(j-2)*sumi(j)
       enddo
       q2t2d = 2.0d0*tsum1 + (tau - tauc)*tsum2
       return
       end function q2t2d
!
       function qcalc(tau)
!
       implicit none
       real(kr) qcalc
       real(kr) tau
!
       real(kr) tsum
       integer j
!
       tsum = 0.0d0
       do j = 2,7
         tsum = tsum + tadif(j)*sumi(j)
       enddo
       qcalc=sumi(1)+(tau-tch2)*tsum
       return
       end function qcalc
!
       function qdta(tau)
!  Partial der of q, pq/prho.
!
       implicit none
       real(kr) qdta,qdt
       real(kr) d,tau
!
       real(kr) sumi(7),sum,tsum
       equivalence (sumi(1),sum)
       integer i,j
!
    1  tsum = 0.0d0
       do j = 2,7
         tsum = tsum + tadif(j)*sumi(j)
       enddo
       qdta = sum +(tau - tauc)*tsum
       return
       entry qdt(d,tau)
       sum = 0.0d0
       do i = 2,8
         sum = sum + xm1(i-1)*a(i,1)*radif(i-1)
       enddo
       sum = sum +er*(a(10,1) - e*(a(9,1) + a(10,1)*d))
       do j = 2,7
         sumi(j) = 0.0d0
         do i = 2,8
           sumi(j) =sumi(j) +xm1(i-1)*a(i,j)*rbdif(i-1)
         enddo
       sumi(j) = sumi(j) + er*(a(10,j) - e*(a(9,j) + a(10,j)*d))
       enddo
       go to 1
       end function qdta
!
       subroutine qmust (d)
       implicit none
       real(kr) d,tau
!
       real(kr) tsum1,tsum2
       integer i,j
!
       radif(1) = 1.0d0
       radif(2) = d - rhoa
       rbdif(1) = 1.0d0
       rbdif(2) = d - rhob
       do i = 3,8
         rbdif(i) = rbdif(i-1)*rbdif(2)
         radif(i) = radif(i-1)*radif(2)
       enddo
       ed = e*d
       er = 1.0d0/exp(ed)
       sumi(1) = 0.0d0
       do i = 1,8
         sumi(1) = sumi(1) + a(i,1)*radif(i)
       enddo
       sumi(1) = sumi(1) + er*(a(9,1) + a(10,1)*d)
       do j = 2,7
         sumi(j) = 0.0d0
         do i = 1,8
           sumi(j) = sumi(j) + a(i,j)*rbdif(i)
         enddo
         sumi(j) = sumi(j) + er*(a(9,j) + a(10,j)*d)
       enddo
       return
       entry qmust2 (tau)
       tadif(1) = 0.0d0
       tadif(2) = 1.0d0
       tadif(3) = tau - taua
       do i = 4,7
         tadif(i) = tadif(i-1)*tadif(3)
       enddo
       return
       end subroutine qmust
!
       function qtd(tau)
!  Partial der of q, pq/ptau.
!
       implicit none
       real(kr) qtd
       real(kr) tau
!
       real(kr) tsum1,tsum2
       integer j
!
       tsum1 = 0.0d0
       tsum2 = 0.0d0
       do j = 3,7
         tsum1 = tsum1 + xm1(j-2)*tadif(j-1)*sumi(j)
         tsum2 = tsum2 + tadif(j)*sumi(j)
       enddo
       tsum2 = tsum2 + sumi(2)
       qtd = tsum2 + (tau - tauc)*tsum1
       return
       end function qtd
!
       subroutine root (x0,x2,fofx,func,x1)
!  Same as rootx--needed to prevent recursion.
!  solve for x1 such that func(x1) = fofx, where x1 lies between x0 and
!  x2.
!
       implicit none
       real(kr) x0,x2,fofx,x1
       real(kr) func
!
       real(kr) ans,b,f0,f2,fx,fxl,tol,x,xx0,xx2
       integer jump,kount
!
       tol = 1.d-5
       xx0 = x0
       xx2 = x2
       f0 = func(xx0)
       f2 = func(xx2)
       ans = (fofx - f0)/(f2 - f0)
       if (ans) 1007,120,120
  120  if (ans - 1.d0) 130,130,1008
  130  if (fofx - 0.d0) 80,70,80
   70  assign 100 to jump
       go to 90
   80  assign 110 to jump
   90  x = (xx0 + xx2)/2.d0
       kount = 0
  150  x1 = x
       kount = kount + 1
       ans = fofx - f2
       fx = func(x)
       fxl = f0 + (x - xx0)*(f2 - f0)/(xx2 - xx0)
       b = abs((fx - fxl)/(f2 - f0))
       if (ans*(fx - fofx) .lt. 0.d0) go to 1001
       xx0 = x
       f0 = fx
       if (b - 0.3d0) 10,20,20
   20  x = (x + xx2)/2.d0
       go to 40
 1001  xx2 = x
       f2 = fx
       if (b - .3d0) 10,30,30
   30  x = (xx0 + x)/2.d0
       go to 40
   10  x = xx0 + (fofx - f0)*(xx2 - xx0)/(f2 - f0)
   40  if (abs((x - x1)/x) - tol) 50,1000,1000
   50  go to jump,(100,110)
  100  if (abs(func(x)) - tol*10.0d0) 60,1000,1000
  110  if (abs((fofx - func(x))/fofx) - tol) 60,1000,1000
 1000  if (kount .gt. 40) tol = tol*10.d0
       if (kount .gt. 60) tol = tol*10.d0
       if (kount .gt. 80) tol = tol*10.d0
       if (kount .lt. 100) go to 150
  160  write (*,170) x1,x
  170  format ('0An iteration has been terminated at 100 iterations,',
     & ' the last two values were',3g15.5)
   60  x1 = x
       return
 1007  x1 = x0
       go to 140
 1008  x1 = x2
  140  write (*,141)
  141  format ('0Solution out of range')
       return
       end subroutine root
!
       subroutine rootx (x0,x2,fofx,func,x1)
!  Solve for x1 such that func(x1) = fofx, where x1 lies between x0 and
!  x2.
!
       implicit none
       real(kr) x0,x2,fofx,x1
       real(kr) func
       external func
!
       real(kr) ans,b,f0,f2,fx,fxl,tol,x,xx0,xx2
       integer jump,kount
!
       tol = 1.d-5
       xx0 = x0
       xx2 = x2
       f0 = func(xx0)
       f2 = func(xx2)
       ans = (fofx - f0)/(f2 - f0)
       if (ans) 1007,120,120
  120  if (ans - 1.d0) 130,130,1008
  130  if (fofx - 0.d0) 80,70,80
   70  assign 100 to jump
       go to 90
   80  assign 110 to jump
   90  x = (xx0 + xx2)/2.d0
       kount = 0
  150  x1 = x
       kount = kount + 1
       ans = fofx - f2
       fx = func(x)
       fxl = f0 + (x - xx0)*(f2 - f0)/(xx2 - xx0)
       b = abs((fx - fxl)/(f2 - f0))
       if (ans*(fx - fofx) .lt. 0.d0) go to 1001
       xx0 = x
       f0 = fx
       if (b - .3d0) 10,20,20
   20  x = (x + xx2)/2.d0
       go to 40
 1001  xx2 = x
       f2 = fx
       if (b - .3d0) 10,30,30
   30  x = (xx0 + x)/2.d0
       go to 40
   10  x = xx0 + (fofx - f0)*(xx2 - xx0)/(f2 - f0)
   40  if (abs((x - x1)/x) - tol) 50,1000,1000
   50  go to jump, (100,110)
  100  if (abs(func(x)) - tol*10.d0) 60,1000,1000
  110  if (abs((fofx - func(x))/fofx) -tol) 60,1000,1000
 1000  if (kount .gt. 40) tol = tol*10.d0
       if (kount .gt. 60) tol = tol*10.d0
       if (kount .gt. 80) tol = tol*10.d0
       if (kount .lt. 100) go to 150
  160  write (*,170) x1,x
  170  format ('0An iteration has been terminated at 100 iterations,',
     & ' the last two values were',3g15.5)
   60  x1 = x
       return
 1007  x1 = x0
       go to 140
 1008  x1 = x2
  140  write (*,141)
  141  format ('0 solution out of range')
       return
       end subroutine rootx
!
       function solve(xi,f,df)
!  Newton-Raphson iteration given an initial estimate xi and the
!  functions f and df.
!
       implicit none
       real(kr) solve
       real(kr) xi
       real(kr) df,f
       external df,f
!
       real(kr) ans,b,tol,xo,xoo,xn
!
       tol =1.d-5
       ni = 0
       xo = xi
       xn = xi
   10  xoo = xo
       xo = xn
c      xn = xo - f(xo)/df(xo)
       ans = f(xo)
       b = df(xo)
       xn = xo - ans/b
       ni = ni + 1
       if (abs((xn - xo)/xn) - tol) 70,20,20
   20  if (ni .gt. 40) tol = tol*10.d0
       if (ni .gt. 60) tol = tol*10.d0
       if (ni .gt. 80) tol = tol*10.d0
       if (ni - 100) 30,50,50
   30  if (abs((xn - xoo)/xn) -tol) 40,10,10
   40  xn = (xo + xn)/2.d0
       go to 10
   50  write (*,60) xoo,xo,xn
   60  format ('0An iteration has been terminated at 100 iterations.',
     & ' the last three values were',3g15.5)
   70  solve = xn
       return
       end function solve
!
       subroutine surf (ku,kz,tin,surft)
!  This routine calculates the surface tension of liquid water and the
!  Laplace constant.
!
       implicit none
       integer ku,kz
       real(kr) tin,surft
!
       real(kr) alc,t,tau,ze
! Units of  g - m/s**2.
       real(kr) g,tk
       data g,tk/9.80665d0,647.15d0/
!
       irout = 11
!  t is deg K.
       tau = tcheck(ku,kz,tin)
       t = 1000.d0/tau
       surft = 0.0d0
       alc = 0.0d0
       if (t .gt. tk) return
!  Units of surft must be dyne/cm.
!  Units of alc is mm.
       tr = t/tk
       r = abs(1.d0- tr)
       surft = 235.8d0*r**1.256d0*(1.d0- .625d0*r)
       if (tr .gt. .998d0) go to 2
       call dens (ku,tin,ze,ze,dl,dv,1)
       alc = sqrt(surft/(g*abs(dl - dv)))
!  Factor for results to be in mm as in the tables.
    2  alc = alc*alconv(ku)*sqrt(dconv(ku))
       surft = surft*stconv(ku)
       return
       end subroutine surf
!
       subroutine temp (ku,p,d,t,kz)
!  Compute tau=1000./temperature in users units given pressure and
!  density. If kz is specified as 1, tau will be a function of pressure
!  only.
!
       implicit none
       integer ku,kz
       real(kr) p,d,t
!
       real(kr) dsl,dsv,ze
!
       real(kr) dumvar
       irout = 3
       pszz = pcheck(ku,kz,p)
!
!  Determine region.
       if (kz - 1) 10,70,10
   10  dszz = dcheck(ku,d)
       if (pszz - pch2) 20,20,50
   20  tszz = tss(pszz)
       call dens (1,tszz,ze,ze,dsl,dsv,1)
       if (dszz - dsl) 30,60,40
   30  if (dszz - dsv) 50,60,60
   40  kz = 2
       tszz = tszz + .01d0
       go to 80
   50  kz = 3
       tszz = 1.2d0
       go to 80
!  Region 1.
   60  kz = 1
       go to 110
   70  tszz = tss(pszz)
       go to 110
!
!  Regions 2 and 3.
   80  call qmust (dszz)
       dumvar = qdt(dszz,tszz)
       dumvar = q2dt(dszz,tszz)
       tszz = solve(tszz,tsf,dtsf)
!  Verify region.
       if (pszz - pch2) 110,110,90
   90  if (tszz - tch2) 110,100,100
  100  kz = 2
  110  t = tszz*tconv(ku)
       return
       end subroutine temp
!
       subroutine tempph (ku,p,h,t,d,dl,dv,kz)
       implicit none
       integer ku,kz
       real(kr) p,h,t,d,dl,dv
!
       real(kr) ds,dsl,dsv,hsl,hsv,ts,ts1,ts2,ze
!
       pszc = pcheck(ku,kz,p)
       irout = 6
       hszc = h/hconv(ku)
       if (hszc - hsch1) 20,10,10
   10  if (hszc - hsch2) 40,40,20
!
!  Input h - out of range tag.
   20  write (*,301) hszc,hsch1,hsch2
  301  format ('0input h =  ',g14.6,'j/g  is out of range of hmin=',
     & f4.1, 'and hmax = ',f7.1,  'j/g')
!
   40  if (pszc - pch2) 140,140,130
  130  ts1 = tch1
       ts2 = tch3
       go to 110
  140  ts = 0.0d0
       call dens (1,ts,pszc,ze,dl,dv,1)
       if (kz - 1) 50,70,50
   50  call enth (1,ts,dl,hsl)
       call enth (1,ts,dv,hsv)
       if (hszc - hsl) 90,70,60
   60  if (hszc - hsv) 70,70,100
!  Region 1.
   70  kz = 1
   80  call dens (1,ts,ze,ze,dsl,dsv,1)
       dl = dsl*dconv(ku)
       dv = dsv*dconv(ku)
       go to 120
!  Region 2.
   90  kz = 2
       ts1 = tch3
       pszc = pszc*1.00011d0
       ts2 = 1000.d0/(1000.d0/ts - 1.d-5)
       go to 110
!  Region 3.
  100  kz = 3
       ts1 = 1000.d0/(1000.d0/ts + 1.d-5)
       pszc = pszc*.99988d0
       ts2 = tch1
!  Regions 2 and 3.
  110  call rootx (ts1,ts2,hszc,tshf,ts)
       call dens (1,ts,pszc,ds,ze,ze,kz)
       d = ds*dconv(ku)
!  Verify region.
       if (pszc - pch2) 120,120,150
  150  if (ts -tch2) 170,170,160
  160  kz = 2
       go to 120
  170  kz = 3
  120  t = ts*tconv(ku)
       return
       end subroutine tempph
!
       subroutine tempps (ku,p,s,t,d,dl,dv,kz)
       implicit none
       integer ku,kz
       real(kr) p,s,t,d,dl,dv
!
       real(kr) ds,dsl,dsv,smax,ss,ssl,ssv,ts,ts1,ts2,ze
!
       irout = 7
       smax = 13.26d0
       pszc = pcheck(ku,kz,p)
       sszc = s/sconv(ku)
       if (sszc .lt. 0.0d0) go to 20
       if (sszc .le. smax) go to 40
!
!  Input s out of range tag.
   20  write (*,301) sszc,smax
  301  format ('0input s = ',g14.6,'j/g-k is out of range of  smin=0.',
     & '0 and smax=   ',f7.1,'j/g-k')
   40  if (pszc - pch2) 140,140,130
  130  ts1 = tch1
       ts2 = tch3
       go to 110
  140  ts = 0.0d0
       call dens (1,ts,pszc,ze,dl,dv,1)
       if (kz - 1) 50,70,50
   50  call ent(1,ts,dl,ssl)
       call ent(1,ts,dv,ssv)
       if (sszc - ssl)  90,70,60
   60  if (sszc - ssv) 70,70,100
!  Region 1.
   70  kz = 1
   80  call dens (1,ts,ze,ze,dsl,dsv,1)
       dl = dsl*dconv(ku)
       dv = dsv*dconv(ku)
       go to 120
!  Region 2.
   90  kz = 2
       ts1 = tch3
       pszc = pszc*1.00011d0
       ts2 = ts*1.00001d0
       go to 110
! Region 3.
  100  kz = 3
       ts1 = ts*.99999d0
       pszc = pszc*.99988d0
       ts2 = tch1
!  Regions 2 and 3.
  110  call rootx (ts1,ts2,sszc,tpsf,ts)
       call dens (1,ts,pszc,ds,ze,ze,kz)
       d = ds*dconv(ku)
!  Verify region.
       if (pszc - pch2) 120,120,150
  150  if (ts - tch2) 170,170,160
  160  kz = 2
       go to 120
  170  kz = 3
  120  t = ts*tconv(ku)
       return
       end subroutine tempps
!
       subroutine therm (ku,kz,pin,tin,din,excesk,tcond)
!  The yata minamiyama form for therm cond is not necessairly the form
!  acceptable to iaps but added here because of its simplicity and
!  accuracy   rch    5  july   1977.
!  Subroutine calculates the thermal conductivity in internal units of
!  w/cm-k and converts to users units.
!  Equations are the internationally agreed upon ones in regions where
!  same are available and are proposed equations in other regions.
!  the near subcritical region is the authors fit.
!
       implicit none
       integer ku,kz
       real(kr) ds,pin,tin,din,excesk,tcond
!
       real(kr) az,ans,b,d,d2,delt2,dtsf1,pbar,pr,pmn,pstate,sum1,
     & sum2,sum4,t,tk
!
       real(kr) za0,za1,za2,za3
       data za0,za1,za2,za3/1.02811d-2,2.99621d-2,1.56146d-2,
     & -4.22464d-3/
       real(kr) tc,dc
       data tc,dc/647.3d0,.3177d0/
       real(kr) c1,c2,c3,c4,c5,c6
       data c1,c2,c3,c4,c5,c6/.642857d0,-4.11717d0,-6.17937d0,3.08976d-3
     #,
     & 8.22994d-2,10.0932d0/
       real(kr) q1,q2,q3,q4
       data q1,q2,q3,q4/7.01309d-2,1.1852d-2,1.69937d-3,-1.02d0/
!
       irout = 10
       pmn = pcheck(ku,kz,pin)
       tdum = tcheck(ku,kz,tin)
       ds = dcheck(ku,din)
!  Convert tdum and pmn to various units.
       tk = 1000.d0/tdum
       tr = tk/tcrt
       t = tk - 273.15d0
       pbar = pmn*10.d0
       pr = pmn/pcrt
!  Out of range check on press and temp.
       if (pbar.lt..01d0.or. pbar.gt.1000.d0) write(*,151) tin,pin
       if (t.lt.0.0d0.or. t.gt.800.d0) write (*,151) tin,pin
  151  format ('0  t =',f12.4,' or  p =',f12.4,' is out of range,',
     & ' returned thermal conductivity is extrapolated')
       t = tk
       tr = t/tc
       d = ds
       r = d/dc
       sum1 = sqrt(tr)*(((za3*tr + za2)*tr + za1)*tr + za0)
       sum2 = -.39707d0+ .400302d0*r + 1.06d0*exp(-.171587d0*(r + 2.3921
     #9d0)**2)
       delt2 = abs(tr - 1.d0) + c4
       az = 2.d0+ c5/delt2**.6d0
       b = 1.d0+ az
       d2 = 1.d0/delt2
       if (tr .lt. 1.d0) d2 = c6/delt2**.6d0
       sum4 = (q1/tr**10 + q2)*r**1.8d0*exp(c1*(1.d0- r**2.8d0)) +
     & q3*d2*r**az*exp(az*(1.d0- r**b)/b) + q4*exp(c2*tr**1.5d0+ c3/r**5
     #)
       ans = sum1 + sum2 + sum4
       tcond = ans*kconv(ku)
       return
       end subroutine therm
!
       function tsf(ts)
!  Function used to solve for temperature ts given pressure and density.
!
       implicit none
       real(kr) tsf,dtsf
       real(kr) ts
!
       real(kr) dtsf1,pstate
!
       call qmust2 (ts)
       pstate = 1000.d0*r*dszz/ts*(1.d0+ dszz*(qcalc(ts) + dszz*qdta(ts)
     #))
       tsf = pstate - pszz
       return
c
       entry dtsf(ts)
       call qmust2 (ts)
       dtsf1 = r*dszz*((1.0d0+ dszz*dszz*qdta(ts) + dszz*qcalc(ts)) -
     & ts*dszz*(dszz*q2dta(ts) + qtd(ts)))
       tsf = dtsf1*(-1000.d0/(ts*ts))
       return
       end function tsf
!
       function tshf(ts)
       implicit none
       real(kr) tshf,tpsf
       real(kr) ts
!
       real(kr) ds,hsc,ssc,ze
       integer kz
!
       kz = 0
       call dens (1,ts,pszc,ds,ze,ze,kz)
       call enth (1,ts,ds,hsc)
       tshf = hsc
       return
c
       entry tpsf(ts)
       kz = 0
       call dens (1,ts,pszc,ds,ze,ze,kz)
       call ent (1,ts,ds,ssc)
       tshf = ssc
       return
       end function tshf
!
       function tss(ps)
!  Compute saturation temperature (k) for given saturation pressure ps
!  (mpa) and return tss as corresponding tau value (1/k).
!
!  J. E. Tolli, EG&G Idaho, Inc., 4/88.
!
       implicit none
       real(kr) tss
       real(kr) ps
!
       real(kr) atestm,tmpest
!
!  Check for critical pressure.
       if (abs(ps - pcrt) .le. 1.0d-5) then
         tss = 0.0d0
         go to 10
       endif
!  Calculate constant term for saturation equation.
       aa1 = -log (ps/pcrt)
!  Get initial temperature estimate.
       tmpest = 278.d0
       if (ps .gt. 0.001d0) tmpest = 323.d0
       if (ps .gt. 0.1d0) tmpest = 423.d0
       if (ps .gt. 1.0d0) tmpest = 473.d0
!  Calculate adjusted tau estimated value.
       atestm = 1.d0- tmpest/tcrt
!  Solve saturation equation.
       tss = solve (atestm,tssf,dtssf)
!  Calculate corresponding temperature.
   10  tss = (1.d0- tss)*tcrt
c--calculate tau value
       tss = 1000.d0/tss
       return
       end function tss
!
       function tssf(tss)
!  Function used to solve for saturation temperature tss given pressure.
!
       implicit none
       real(kr) tssf,dtssf
       real(kr) tss
!
       tssf = (cps1*tss + cps2*tss**1.9d0+ cps3*tss**2 + cps4*tss**5.5d0
     #+
     & cps5*tss**10)/(1.d0- tss) + aa1
       return
       entry dtssf(tss)
!  Derivative of function used to solve for saturation temperature tss
!  given pressure.
!
       tssf = ((cps1 + 1.9d0*cps2*tss**0.9d0+ 2.0d0*cps3*tss +
     & 5.5d0*cps4*tss**4.5d0+ 10.d0*cps5*tss**9) +
     & ((cps1*tss + cps2*tss**1.9d0+ cps3*tss**2 + cps4*tss**5.5d0+
     & cps5*tss**10)/(1.d0- tss)))/(1.d0- tss)
       return
       end function tssf
!
       subroutine visc (ku,kz,tin,pin,din,svisc)
!  Calculate the viscosity given tau,p,and d in user-s units ku.
!  Answer returned in user-units in svisc.
!
       implicit none
       integer ku,kz
       real(kr) ds,tin,pin,din,svisc
!
       real(kr) amuo,deld,delmu,dr,p,pr,pmn,spvr,sum,t,tk,tsr
       real(kr) avo(4),dt(6),dd(5),bv(6,5)
       real(kr) bv1(6),bv2(6),bv3(6),bv4(6),bv5(6)
       equivalence (bv1(1),bv(1,1)),(bv2(1),bv(1,2)),(bv3(1),bv(1,3)),
     & (bv4(1),bv(1,4)),(bv5(1),bv(1,5))
       real(kr) rhos,tstr
       integer i,j
!
       data avo/1.81583d-2,1.77624d-2,1.05287d-2,-3.6744d-3/
       data rhos,tstr/.317763d0,647.27d0/
       data bv1/0.501938d0,0.162888d0,-0.130356d0,0.907919d0,-0.551119d0
     #,0.146543d0/
       data bv2/0.235622d0,0.789393d0, 0.673665d0,1.207552d0, 0.0670665d
     #0,
     & -0.0843370d0/
       data bv3/-0.274637d0,-0.743539d0,-0.959456d0,-0.687343d0,-0.49708
     #9d0,
     & 0.195286d0/
       data bv4/0.145831d0,0.263129d0,0.347247d0,0.213486d0,0.100754d0,-
     #0.032932d0/
       data bv5/-0.0270448d0,-0.0253093d0,-0.0267758d0,-0.0822904d0,0.06
     #02253d0,
     & -0.0 202595d0/
!
       irout = 9
!
!  tk is deg k, t is deg c, tr is reduced temp.
       tk = 1000.d0/tcheck(ku,kz,tin)
       t = tk - 273.15d0
       tr = tk/tcrt
!  p is bars, pmn is mega newtons/m*m ,pr is reduced pressure.
       pmn = pcheck(ku,kz,pin)
       p = 10.0d0*pmn
       pr=pmn/pcrt
!  ds is g/cc, spvr is reduced specific volume..
       ds = dcheck(ku,din)
       spvr = rhocrt/ds
!  Check for out of range on p and t.
       if (p.lt..99d0.or. p.gt.1000.01d0) write(*,100) t,p
       if (t.gt.900.0d0.or. t.lt.0.0d0) write(*,100) t,p
  100  format (' Out of range; viscosity extrapolated  for t=',
     & f12.4,' p= ',f12.4)
!  Calculate viscosity for all regions to 900 cent and 1000 bar.
       tst = tstr/tk -1.d0
       dt(1) = 1.d0
       dt(2) = tst
       dt(3) = tst*dt(2)
       dt(4) = tst*dt(3)
       dt(5)  =tst*dt(4)
       dt(6) = tst*dt(5)
       dr = ds/rhos
       deld = dr - 1.d0
       dd(1) = 1.d0
       dd(2) = deld
       dd(3) = deld*dd(2)
       dd(4) = deld*dd(3)
       dd(5) = deld*dd(4)
       sum = 0.d0
  101  do i = 1,6
         do j = 1,5
           sum = sum + bv(i,j)*dt(i)*dd(j)
         enddo
       enddo
       delmu = exp(dr*sum)
  102  tsr = tstr/tk
       amuo = sqrt(tk/tstr)/
     & (((tsr*avo(4) + avo(3))*tsr + avo(2))*tsr +avo(1))
       svisc = delmu*amuo*1.d-5 *mconv(ku)
       return
       end subroutine visc
!
       end module stgd2omod
