      Program Subsonic
C
C   This program solves the hydrodynamic escape equations for the subsonic
C   regime, following Roger Yelle's new paper (2003) and also Kasting and
C   Pollack (1983). The initial version attempts to look at a pure H2
C   atmosphere, similar to Watson et al. (1981).
C
      parameter(M=301,JJ=22)
      dimension r(M),rkm(M),z(M),zl(M),u(M),rho(M),T(M),g(M),r2(M),p(M),
     2  Ha(M),um(M),q(M),rhs(M),u02(M),ujeans(M),drp(M),drm(M),cp(M),
     3  flux(M),den(M),psi(M),alam(M),alamp(M),r0(M),delT(M)
      dimension a(M),b(M),c(M)
      dimension lwav(JJ),wav(JJ),phflux(JJ),sigh2(JJ),sigh(JJ),
     2  effh2(JJ),effh(JJ),eflux(JJ),solarf(M,JJ)
C
C   Define I/O files
      open(unit=1,file='subsonic.out')
      open(unit=2,file='solarf.out')
      open(unit=3,file='utden.out')
      open(unit=4,file='kasting2.dat',status='old')
C
C   Read the EUV data file
      read(4,400)
 400  format(/)
      do 401 j=1,JJ
      read(4,402) lwav(j),phflux(j),sigh2(j),sigh(j),effh2(j),
     2   effh(j),eflux(j)
 402  format(1x,i3,6(1x,e9.2))
      wav(j) = lwav(j)
C      write(1,402) lwav(j),phflux(j),sigh2(j),sigh(j),effh2(j),
C     2   effh(j),eflux(j)
 401  continue
C
C   Define numerical constants (in CGS units)
C
C   bigG = gravitational constant
C   sm = atomic hydrogen mass
C   bk = Boltzmann's constant
C   rE = radius of Earth
C   bigM = mass of Earth
C   aE = Earth-Sun distance
C   solarM = Sun's mass
C   appak0 = thermal conductivity constant (not used anymore for H2)
C   Rgas = universal gas constant
C   sigma = molecular collision cross section
C   umh2 = molecular wt. of H2 (the only gas in there right now)
C   dentot = assumed total number density at 100 km, including N2
C
      N = M
      bigG  = 6.67e-8
      sm = 1.67e-24
      bk = 1.38e-16
      rE = 6.371e8
      rEkm = rE/1.e5
      bigM = 5.98e27
      aE = 1.496e13
      solarM = 1.99e33
      Rgas = 8.314e7
      rootpi = sqrt(3.14159)
      sigma = 3.e-15
      umh2 = 2.
      dentot = 1.e13
C
C   Coefficients for H2 heat capacity (valid 298-1000 K)
      AH2 = 33.066178
      BH2 = -11.363417
      CH2 = 11.432816
      DH2 = -2.772874
      EH2 = -0.158558
C
C   Define an altitude grid (z = altitude in km, zl = log z)
      do 1 i=1,N
      zl(i) = 2. + (i-1)/100.
      z(i) = 10.**zl(i)
      rkm(i) = z(i) + rEkm
      r(i) = 1.e5*rkm(i)
      r2(i) = r(i)*r(i)
      g(i) = bigG*bigM/r2(i)
      um(i) = umh2
   1  continue
C
C      write(1,100)
 100  format(1x,'Altitudes in km')
C      write(1,101) z
 101  format(1x,1p10e10.3)
C      write(1,102)
 102  format(/1x,'Radial distance in cm')
C      write(1,101) r
C
C   Set lower boundary conditions
      T0 = 250.
C   (This is Watson et al.'s lower B.C.)
      T(1) = T0
      u(1) = 0.
      den(1) = 1.e12
      rho(1) = sm*um(1)*den(1)
      p(1) = den(1)*bk*T(1)
      Ha(1) = bk*T(1)/(sm*um(1)*g(1))
C
C   Calculate diffusion-limited flux, assuming a fixed total number
C   density at the lower boundary
      fH2 = den(1)/dentot
      phidif = 2.5e13 * fH2
      write(1,120) den(1),fH2,phidif
 120  format(/1x,'nH2 =',1pe10.3,2x,'fH2 =',e10.3,2x,'phidif =',e10.3)
C
C   Set initial temperature at infinity. Then find the corresponding
C   isothermal critical point and sound speed
      Tinf = 600.
      u02iso = bk*Tinf/(um(1)*sm)
      u0iso = sqrt(u02iso)
      r0iso = bigG*bigM/(2.*u02iso)
      z0 = (r0iso-rE)/1.e5
      write(1,107) Tinf
 107  format(/1x,'Initial temperature at infinity: Tinf =',1pe10.3)
      write(1,106) r0iso,z0,u0iso
 106  format(/1x,'Isothermal critical point is at:',/1x,'r0 =',1pe10.3,
     2  2x,'z0 =',e10.3,' km',2x,'u0 =',e10.3/)
C
C   Set initial conditions over the rest of the grid
      do 2 i=2,N
      u(i) = 0.
      T(i) = T(1) + (Tinf-T0)*(z(i)-z(1))/2000.
      T(i) = amin1(T(i),Tinf)
      Ha(i) = bk*T(i)/(sm*um(i)*g(i))
      dr = r(i) - r(i-1)
      p(i) = p(i-1)*exp(-dr/Ha(i))
      rho(i) = sm*um(i)*p(i)/(bk*T(i))
      den(i) = rho(i)/(um(i)*sm)
   2  continue
C
C      write(1,103)
 103  format(/1x,'Scale height, Ha, in cm')
C      write(1,101) Ha
      write(1,104)
 104  format(/1x,'Initial temperature')
      write(1,101) T
C      write(1,105)
 105  format(/1x,'Initial (hydrostatic) density, rho, in gm/cm^3')
C      write(1,101) rho
C
      do 3 i=2,N-1
      drp(i) = r(i+1)-r(i)
      drm(i) = r(i) - r(i-1)
   3  continue
      drp(1) = r(2)-r(1)
      drm(N) = r(N)-r(N-1)
      imax = 0
      emax = 0.
      u(N) = 0.
      flux(N) = 0.
C
      time = 0.
      dt = 1.e-5
      tstop = 1.e10
      nsteps = 1
C
C *************************************************************
C   Start the time-stepping loop
      do 10 L=1,nsteps
      time = time + dt
      Nsave = N
      write(1,200) L,dt,N,z(N),u(N),T(N),flux(N),emax,imax
 200  format(1x,'L=',i4,2x,'dt=',1pe8.1,2x,'N=',i3,2x,'z=',1pe8.1,
     2  2x,'u=',e8.1,2x,'T=',e9.2,2x,'flux=',e8.1,2x,'emax=',e8.1,
     3  ' at i=',i3)
C
C   Calculate heating and cooling rates. Note the geometrical heating
C   effect: more energy is available at large distances.
      effEUV = 0.15
      angle = 60. * (3.14159/180.)
      angcos = cos(angle)
      column = den(N)*Ha(N)
      q(N) = 0.
C
      do 405 j=1,JJ
      solarf(N,j) = eflux(j)*exp(-sigh2(j)*column/angcos)
      q(N) = q(N) + solarf(N,j)*sigh2(j)*den(N)*effEUV*r2(N)/r2(1)
C
      do 30 i=N-1,1,-1
      denav = sqrt(den(i)*den(i+1))
      column = column + denav*drp(i)
      solarf(i,j) = eflux(j)*exp(-sigh2(j)*column/angcos)
      q(i) = q(i) + solarf(i,j)*sigh2(j)*den(i)*effEUV*r2(i)/r2(1)
  30  continue
  405 continue
C      write(1,406)
  406 format(/1x,'q')
C      write(1,101) q
C
C   Set H2 heat capacity
      do 4 i=1,N
      tl = T(i)/1000.
      tl2 = tl*tl
      tl3 = tl2*tl
   4  cp(i) = (AH2 + BH2*tl + CH2*tl2 + DH2*tl3 + EH2/tl2)/2. * 1.e7
C      write(1,212)
 212  format(/1x,'cp')
C      write(1,101) cp
C
C   Advance the temperature equation using FTCS differencing
      rhs(1) = 0.
      do 5 i=2,N-1
      rip = (r(i)+r(i+1))/2.
      rip2 = rip*rip
      rim = (r(i)+r(i-1))/2.
      rim2 = rim*rim
      Tp = (T(i)+T(i+1))/2.
      Tm = (T(i)+T(i-1))/2.
      appakp = thercon(Tp)
      appakm = thercon(Tm)
      drp2 = drp(i)*drp(i)
      drm2 = drm(i)*drm(i)
      dr2 = drp(i) + drm(i)
      denom1 = drp(i)*drm(i)*dr2
      denom2 = rho(i)*cp(i)*r2(i)*dr2
C
C   Set things up in tridiagonal form
      a(i) = u(i)*rip2/denom1 + 2.*rim2*appakm/drm(i)/denom2
      b(i) = - u(i)*(rip2-rim2)/denom1 - (rim2*appakm/drm(i)
     2   + rip2*appakp/drp(i))*2./denom2
      c(i) = - u(i)*rim2/denom1 + 2.*rip2*appakp/drp(i)/denom2
      d = q(i)/(rho(i)*cp(i)) - u(i)*g(i)/cp(i)
      rhs(i) = a(i)*T(i-1) + b(i)*T(i) + c(i)*T(i+1) + d
C      write(1,3000) i,q(i),a(i),b(i),c(i),d,rhs(i)
 3000 format(1x,i3,2x,1p6e10.3)
   5  continue
C
C   Set boundary conditions
      a(1) = 0.
      b(1) = b(2)
      c(1) = 0.
      rhs(1) = 0.
      a(N) = - b(N-1)
      b(N) = b(N-1)
      c(N) = 0.
      rhs(N) = b(N)*(T(N)-T(N-1))
C
C   Compute the matrix (I/dt-J) and store as vectors a, b, and c
      do 25 i=1,N
      a(i) = - a(i)
      b(i) = 1./dt - b(i)
      c(i) = -c(i)
  25  continue
C
C   Call the matrix routine
      call tridag(a,b,c,rhs,delT,N)
C      write(1,301)
 301  format(/1x,'delT')
C      write(1,101) delT
C
C   Compute new temperatures and maximum change in T at any height
      emax = 0.
      do 26 i=1,N
      chg = abs(delT(i)/T(i))
      if(chg.gt.emax) imax = i
      emax = amax1(emax,chg)
      T(i) = T(i) + delT(i)
  26  continue
C
C   Recompute T-dependent quantities
      do 6 i=1,N
      u02(i) = bk*T(i)/(um(i)*sm)
      Ha(i) = bk*T(i)/(um(i)*sm*g(i))
      r0(i) = bigG*bigM/(2.*u02(i))
      alam(i) = r(i)/Ha(i)
      ve = sqrt(2.*bigG*bigM/r(i))
      veu = ve - u(i)
      veu2 = veu*veu
      alamp(i) = veu2/(2.*u02(i))
      ujeans(i) = rootpi/2. * exp(-alamp(i))*sqrt(2.*u02(i))*
     2  (1.+alamp(i))
C      ujeans(i) = rootpi/2. * exp(-alam(i))*sqrt(2.*u02(i))*(1.+alam(i))
C   I have commented out the standard Jeans velocity. The one being used is
C   my modified Jeans velocity from a moving bulk atmosphere
   6  continue
C
C      write(1,201)
 201  format(/1x,'Temperature, T, in Kelvins')
C      write(1,101) T
C
C   Now solve the wind equation from the outside in using Jeans escape 
C   velocity conditions at the upper boundary. Jeans' formula is
C   modified to account for non-zero flow velocity at the exobase
      u(N) = ujeans(N)
C
      do 7 i=N,2,-1
      r0av = 0.5*(r0(i)+r0(i-1))
      rhsu = alog(T(i)/T(i-1)) - alog(um(i)/um(i-1)) + 2.*r0av*
     2  (1./r(i-1)-1./r(i)) -2.*alog(r(i)/r(i-1))
      ui2 = u(i)*u(i)
      rhsu1 = rhsu/(1.-ui2/u02(i))
      x = u(i)*exp(-rhsu1)
C   Use 1 step of Newton's method to refine estimate of u(i-1)
      f = alog(u(i)) - alog(x) + (x*x-ui2)/(2.*u02(i)) - rhsu
      fp = x/u02(i) - 1./x
      u(i-1) = x - f/fp
   7  continue
C
C      write(1,202)
 202  format(/1x,'Velocity, u, in cm/s')
C      write(1,101) u
C
C   Finally, solve the density equation upward from the lower boundary
      do 8 i=2,N
      u02av = 0.5*(u02(i)+u02(i-1))
      r0av = 0.5*(r0(i)+r0(i-1))
      rhsr = - (u(i)*u(i)-u(i-1)*u(i-1))/(2.*u02av) + 2.*r0av*(1./r(i)-
     2  1./r(i-1)) - alog(T(i)/T(i-1)) + alog(um(i)/um(i-1))
      rho(i) = rho(i-1)*exp(rhsr)
   8  continue
C
C      write(1,203)
 203  format(/1x,'Density, rho, in gm/cm^3')
C      write(1,101) rho
C
C   Calculate number densities and flux conservation and print out results
      do 9 i=1,N
      den(i) = rho(i)/(um(i)*sm)
      flux(i) = den(i)*u(i)*r2(i)/r2(1)
      psi(i) = u(i)*u(i)/u02(i)
      p(i) = den(i)*bk*T(i)
   9  continue
C
      if(L.lt.nsteps) go to 21
      write(1,208) time
 208  format(/1x,'time =',1pe10.3)
      write(3,204)
 204  format(3x,'i',4x,'z(km)',5x,'rho',7x,'den',7x,'u',8x,'ujeans',
     2  5x,'psi',7x,'alam',6x,'T',8x,'delT',6x,'flux')
      do 11 i=1,N
      write(3,205) i,z(i),rho(i),den(i),u(i),ujeans(i),psi(i),alam(i),
     2  T(i),delT(i),flux(i)
 205  format(1x,i3,1p10e10.3)
  11  continue
      write(3,204)
      write(3,207)
 207  format(1x,i2)
C
      write(2,304)
 304  format(3x,'i',4x,'z(km)',4x,'solarf',5x,'q',9x,'den',7x,'p',9x,
     2  'Ha',7x,'r(cm)')
      do 111 i=1,N
      write(2,205) i,z(i),solarf(i,1),q(i),den(i),p(i),Ha(i),r(i)
 111  continue
      write(2,304)
C   Estimate the pressure at infinity assuming hydrostatic equilibrium
      pinf = p(N)*exp(-r(N)/Ha(N))
      write(2,305) pinf
 305  format(/1x,'pinf =',1pe10.3)
  21  continue
C
C   Calculate the exobase height
      siginv = 1./sigma
      do 12 i=1,N
      is = i
      prod = den(i)*Ha(i)
      if (prod.lt.siginv) go to 13
  12  continue
  13  continue
C
C   Readjust the top of the grid
      if(is.lt.N) go to 24
C   If this statement executed, the exobase was below the old height,
C     so it should move down
C
C   The exobase is at or above the top of the grid, so need to expand
C   the grid to see if it should move up. To do this, first calculate
C   T, u, and rho going outwards from the old exobase
C
C      write(1,1000)
 1000 format(/1x,'Grid moves up')
C
C   Assume that T is constant
      do 16 i=N+1,M
      T(i) = T(N)
      um(i) = um(N)
      u02(i) = bk*T(i)/(um(i)*sm)
      Ha(i) = bk*T(i)/(um(i)*sm*g(i))
      r0(i) = bigG*bigM/(2.*u02(i))
      alam(i) = r(i)/Ha(i)
      ve = sqrt(2.*bigG*bigM/r(i))
      veu = ve - u(i)
      veu2 = veu*veu
      alamp(i) = veu2/(2.*u02(i))
      ujeans(i) = rootpi/2. * exp(-alamp(i))*sqrt(2.*u02(i))*
     2  (1.+alamp(i))
C      ujeans(i) = rootpi/2. * exp(-alam(i))*sqrt(2.*u02(i))*(1.+alam(i))
  16  continue
C
C      write(1,1001) T(M)
 1001 format(5x,'T(M) =',1pe10.3)
C
C   Recalculate u
      do 17 i=N+1,M
      r0av = 0.5*(r0(i)+r0(i-1))
      rhsu = alog(T(i)/T(i-1)) - alog(um(i)/um(i-1)) + 2.*r0av*
     2  (1./r(i-1)-1./r(i)) -2.*alog(r(i)/r(i-1))
      ui2 = u(i-1)*u(i-1)
      rhsu1 = rhsu/(1.-ui2/u02(i-1))
      u(i) = u(i-1)*exp(rhsu1)
  17  continue
C
C      write(1,1002) u(M)
 1002 format(5x,'u(M) =',1pe10.3)
C
C   Now solve the density equation upward from the previous upper
C     boundary
      do 18 i=N+1,M
      u02av = 0.5*(u02(i)+u02(i-1))
      r0av = 0.5*(r0(i)+r0(i-1))
      rhsr = - (u(i)*u(i)-u(i-1)*u(i-1))/(2.*u02av) + 2.*r0av*(1./r(i)-
     2  1./r(i-1)) - alog(T(i)/T(i-1)) + alog(um(i)/um(i-1))
      rho(i) = rho(i-1)*exp(rhsr)
      den(i) = rho(i)/(um(i)*sm)
      flux(i) = den(i)*u(i)*r2(i)/r2(1)
  18  continue
C
C      write(1,1003) rho(M)
 1003 format(5x,'rho(M) =',1pe10.3)
C
C   Finally, recalculate the exobase height starting from the old value
      do 22 i=N,M
      is = i
      prod = den(i)*Ha(i)
      if (prod.lt.siginv) go to 24
  22  continue
  24  N = is
      Nm5 = Nsave - 5.
      if(N.lt.Nm5) N = Nm5
C
C   Adjust the time step based on how fast the solution is changing
      dt0 = dt
      if(emax.lt.1.e-3) dt = dt0*1.1
      if(emax.lt.1.e-4) dt = dt0*2.
      if(emax.lt.1.e-5) dt = dt0*10.
      if(emax.gt.2.e-3) dt = dt0*0.9
      if(emax.gt.3.e-3) dt = dt0*0.5
C
      if (time.gt.tstop) stop
  10  continue
C   End the time-stepping loop
C *******************************************************************
C
      stop
      end
C
      Function thercon(T)
C   This subroutine calculates the thermal conductivity of H2 from the
C   formula given by Saxena and Saxena, J.Phys. A.: Gern. Phys. 3, 309-
C   320, 1970.
C
      appakH2 = 20.37 + 8.2e-2*T + 3.56e-6*T*T
      thercon = appakH2 * 4.184 * 100.
      return
      end
C
      SUBROUTINE TRIDAG(A,B,C,R,U,N)
      PARAMETER (M=301)
      DIMENSION GAM(M),A(M),B(M),C(M),R(M),U(M)
      IF(B(1).EQ.0.)PAUSE
      BET=B(1)
      U(1)=R(1)/BET
      DO 11 J=2,N
        GAM(J)=C(J-1)/BET
        BET=B(J)-A(J)*GAM(J)
        IF(BET.EQ.0.)PAUSE
        U(J)=(R(J)-A(J)*U(J-1))/BET
11    CONTINUE
      DO 12 J=N-1,1,-1
        U(J)=U(J)-GAM(J+1)*U(J+1)
12    CONTINUE
      RETURN
      END
