subroutine adsorption(p,T,fp,fT,nz)
C***********************************************************************
C   adsorption: calculates the amount adsorbed and returns the 
C               derivatives with respect temperature and pressure
C
C   p = partial pressure of H2O [Pa]
C   T = temperature [K]
C   f = mass density of adsorbed H20 [kg/m^3]
C   fp = partial derivative of f with respect to pressure
C   fT = partial derivative of T with respect to temperature
C   nz = number of grid points
C***********************************************************************
      implicit none
      integer NMAX
      parameter (NMAX=2000)

      real*8 p(NMAX), T(NMAX), f, fp(NMAX), fT(NMAX)
      integer nz, j

c     Jakosky et al. '97
      real*8 K, K0, e, nu, vm
      parameter (K0=1.57e-8, e=2573.9, nu=0.48) 
      parameter (vm=2000.*1.e5*2.84e-7) ! J/mol

c     Zent et al., basalt
!      real*8 rhor, gamma, B0, delta
!      parameter (rhor=2000., B0=0.51)  ! rhor = density of regolith
!      parameter (gamma=8.364e-16, delta=-2679.8)

      do j=1,nz
c        Jakosky et al. '97
         K = K0*exp(e/T(j))
         f = vm*(K*p(j)/(1.+K*p(j)))**nu
         fp(j) = f/p(j)*nu/(1.+K*p(j)) 
         fT(j) = -f*nu*e/T(j)**2/(1.+K*p(j))
c         Zent et al.
!         f = 1.e-5*rhor*(gamma*p(j))**B0*exp(-delta/T(j))
!         fp(j) = B0*f/p(j)      ! fails for p=0
!         fT(j) = delta*f/T(j)**2
      enddo
      end




      subroutine adsorption2(p,T,f,nz,psurf,Tsurf,fsurf)
C***********************************************************************
C   adsorption2: returns the amount of adsorbed H2O
C
C   (use exactly the same isotherms as in subroutine adsorption)
C
C   p = partial pressure of H2O [Pa]
C   T = temperature [K]
C   f = mass density of adsorbed H20 [kg/m^3]
C   nz = number of grid points
C   ?surf = quantities at surface
C***********************************************************************
      implicit none
      integer NMAX
      parameter (NMAX=2000)

      real*8 p(NMAX), T(NMAX), f(NMAX), psurf, Tsurf, fsurf
      integer nz, j

c     Jakosky et al. '97
      real*8 K, K0, e, nu, vm
      parameter (K0=1.57e-8, e=2573.9, nu=0.48) 
      parameter (vm=2000.*1.e5*2.84e-7) ! J/mol

c     Zent et al., basalt
!      real*8 rhor, gamma, B0, delta
!      parameter (rhor=2000., B0=0.51)
!      parameter (gamma=8.364e-16, delta=-2679.8)

      K = K0*exp(e/Tsurf)
      fsurf = vm*(K*psurf/(1.+K*psurf))**nu
      do j=1,nz
c        Jakosky et al. '97
         K = K0*exp(e/T(j))
         f(j) = vm*(K*p(j)/(1.+K*p(j)))**nu
c        Zent et al.
!        f(j) = 1.e-5*rhor*(gamma*p(j))**B0*exp(-delta/T(j))
      enddo
      end