      subroutine st_msa(ipar,istna,stn_adj_file,ivelm,vel_mod_file,
     1                  elatd,elond,stn,mcode,vpa,vpb,vpc,
     2                  vsa,vsb,vsc,vh,psa,ssa)
c
c ...... extract model and station adjustments
c        elatd & elond are in geographical degrees
c
c     parameters:
c         ipar = 1 -> extract model
c         ipar = 2 -> extract station adjustment
c
c	Modified 21 nov 94 to acquire file names from config.h
c
c	Region	Name
c
c	a	Northern Coast
c	b	Northern California
c	c	Ukiah-Willits
c	d	Central Coast
c	e	Sierra Nevada
c	g	Markleeville
c	h	Halls Valley
c	k	Coalinga
c	m	Mammoth
c	s	Bishop
c	z	Santa Cruz Mtns
c	o	Outside
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
	include "config.h"
      character*80 safile,stn_adj_file,vel_mod_file
      character*1 mcode
      character*4 stn
      integer*4 ipar,istna,ivelm
      real*8 elatd,elond,vpa,vpb,vpc,vsa,vsb,vsc,vh,psa,ssa
c
      if(ipar.eq.1) then
        go to 100
      elseif(ipar.eq.2) then
        go to 200
      endif
      stop
c
  100 continue
c
c ...... extract appropriate model parameters
c
      if(ivelm.eq.1) then
c
c ...... read in external velocity model info
c
        open(27,file=vel_mod_file(1:lnblnk(vel_mod_file)),
     1       status='old')
        read(27,'(7f10.0)') vpa,vpb,vpc,vsa,vsb,vsc,vd
        close(27)
        mcode='x'
        return
      endif
      if(elatd.ge.40.d0) then
        if(elond.le.-123.d0) then
          mcode = 'a'
          vpa = 5.24d0
          vpb = 0.068d0
          vpc = 7.98d0
          vsa = 3.03d0
          vsb = 0.039d0
          vsc = 4.61d0
          vh  = 25.d0
        else
          mcode = 'b'
          vpa = 6.00d0
          vpb = 0.024d0
          vpc = 7.85d0
          vsa = 3.46d0
          vsb = 0.014d0
          vsc = 4.53d0
          vh  = 45.d0
        endif
      elseif(elatd.gt.39.0d0) then
        if(elond.le.-121.75d0) then
          mcode = 'c'
          vpa = 5.28d0
          vpb = 0.075d0
          vpc = 7.70d0
          vsa = 2.98d0
          vsb = 0.043d0
          vsc = 4.36d0
          vh  = 25.d0
        else
          mcode = 'e'
          vpa = 6.00d0
          vpb = 0.024d0
          vpc = 7.85d0
          vsa = 3.46d0
          vsb = 0.014d0
          vsc = 4.53d0
          vh  = 45.d0
        endif
      elseif(elatd.gt.38.5d0) then
        if(elond.le.-121.75d0) then
          mcode = 'c'
          vpa = 5.28d0
          vpb = 0.075d0
          vpc = 7.70d0
          vsa = 2.98d0
          vsb = 0.043d0
          vsc = 4.36d0
          vh  = 25.d0
        else
          if(elond.lt.-120.d0) then
            mcode = 'e'
            vpa = 6.00d0
            vpb = 0.024d0
            vpc = 7.85d0
            vsa = 3.46d0
            vsb = 0.014d0
            vsc = 4.53d0
            vh  = 45.d0
          elseif(elond.le.-119.5d0) then
            mcode = 'g'
            vpa = 6.00d0
            vpb = 0.024d0
            vpc = 7.85d0
            vsa = 3.46d0
            vsb = 0.014d0
            vsc = 4.53d0
            vh  = 45.d0
          else
            mcode = 'e'
            vpa = 6.00d0
            vpb = 0.024d0
            vpc = 7.85d0
            vsa = 3.46d0
            vsb = 0.014d0
            vsc = 4.53d0
            vh  = 45.d0
          endif
        endif
      else
        x=-119.5d0-0.9d0*(elatd-36.d0)
        if(elond.le.x) then
          vpa = 5.28d0
          vpb = 0.075d0
          vpc = 7.70d0
          vsa = 2.98d0
          vsb = 0.043d0
          vsc = 4.36d0
          vh  = 25.d0
          if(elatd.ge.37.53d0) then
            mcode='d'
          elseif(elatd.ge.36.70d0) then
            xhr=-121.49d0-0.76471d0*(elatd-37.19d0)
            xhl=-121.66d0-0.76471d0*(elatd-37.10d0)
            xht=-121.92d0+1.88889d0*(elatd-37.44d0)
            xhb=-121.66d0+1.88889d0*(elatd-37.10d0)
            xzr=-121.38d0-1.39623d0*(elatd-36.87d0)
            xzl=-121.55d0-1.36364d0*(elatd-36.70d0)
            xzt=-122.30d0+1.20000d0*(elatd-37.25d0)
            xzb=-121.55d0+1.00000d0*(elatd-36.70d0)
            if(elatd.gt.37.44d0) then
              if(elond.lt.xht) then
                mcode='d'
              elseif(elond.le.xhr) then
                mcode='h'
              else
                mcode='d'
              endif
            elseif(elatd.gt.37.40d0) then
              if(elond.lt.xhl) then
                mcode='d'
              elseif(elond.le.xhr) then
                mcode='h'
              else
                mcode='d'
              endif
            elseif(elatd.ge.37.25d0) then
              if(elond.lt.xzt) then
                mcode='d'
              elseif(elond.le.xzr) then
                mcode='z'
              elseif(elond.lt.xhl) then
                mcode='d'
              elseif(elond.le.xhr) then
                mcode='h'
              else
                mcode='d'
              endif
            elseif(elatd.ge.37.19d0) then
              if(elond.lt.xzl) then
                mcode='d'
              elseif(elond.le.xzr) then
                mcode='z'
              elseif(elond.lt.xhl) then
                mcode='d'
              elseif(elond.le.xhr) then
                mcode='h'
              else
                mcode='d'
              endif
            elseif(elatd.ge.37.10d0) then
              if(elond.lt.xzl) then
                mcode='d'
              elseif(elond.le.xzr) then
                mcode='z'
              elseif(elond.lt.xhl) then
                mcode='d'
              elseif(elond.le.xhb) then
                mcode='h'
              else
                mcode='d'
              endif
            elseif(elatd.ge.36.87d0) then
              if(elond.lt.xzl) then
                mcode='d'
              elseif(elond.le.xzr) then
                mcode='z'
              else
                mcode='d'
              endif
            else
              if(elond.lt.xzl) then
                mcode='d'
              elseif(elond.le.xzb) then
                mcode='z'
              else
                mcode='d'
              endif
            endif
          elseif(elatd.ge.36.61d0) then
            mcode='d'
          elseif(elatd.ge.35.82d0) then
            xr=-119.76d0-1.200d0*(elatd-36.06d0)
            xl=-120.17d0-1.035d0*(elatd-35.82d0)
            xt=-120.76d0+1.545d0*(elatd-36.39d0)
            xb=-120.17d0+1.708d0*(elatd-35.82d0)
            if(elatd.ge.36.39d0) then
              if(elond.lt.xt) then
                mcode='d'
              elseif(elond.gt.xr) then
                mcode='d'
              else
                mcode='k'
              endif
            elseif(elatd.ge.36.06d0) then
              if(elond.lt.xl) then
                mcode='d'
              elseif(elond.gt.xr) then
                mcode='d'
              else
                mcode='k'
              endif
            else
              if(elond.lt.xl) then
                mcode='d'
              elseif(elond.gt.xb) then
                mcode='d'
              else
                mcode='k'
              endif
            endif
            imc=ichar(mcode)
            imk=ichar('k')
            if(imc.eq.imk) then
              vpa = 5.20d0
c              vpb = 0.120d0
              vpb = 0.088d0
              vpc = 7.50d0
              vsa = 2.20d0
              vsb = 0.042d0
              vsc = 4.50d0
              vh  = 18.d0
            endif
          else
            mcode='d'
          endif
        else
          vpa = 6.00d0
          vpb = 0.024d0
          vpc = 7.85d0
          vsa = 3.46d0
          vsb = 0.014d0
          vsc = 4.53d0
          vh  = 45.d0
          if(elond.lt.-119.2d0) then
            mcode='e'
          elseif(elond.lt.-118.7d0) then
            if(elatd.gt.37.75d0) then
              mcode='e'
            elseif(elatd.ge.37.35d0) then
              mcode='m'
            else
              mcode='e'
            endif
          elseif(elond.lt.-118.2d0) then
            if(elatd.gt.37.65d0) then
              mcode='e'
            elseif(elatd.ge.37.20d0) then
              mcode='s'
            else
              mcode='e'
            endif
          else
            mcode='e'
          endif
        endif
      endif
      return
c
c ...... extract station adjustments
c
  200 continue
      if(istna.eq.1) then
        safile=stn_adj_file
        call saextr(safile,stn,psa,ssa)
        return
      elseif(istna.eq.2) then
        psa=0.d0
        ssa=0.d0
        return
      elseif(mcode.eq.'x') then
        write(*,*) 'Must specify "-s" option when using "-v" option'
        stop
      endif
c
c ...... modified for file parameter input
c
	do i=1,lnblnk(models)
	  if(mcode.eq.models(i:i)) then
	    safile=model_dir(1:lnblnk(model_dir)) // '/' //
     1        'st_relp_stn_adj.' // models(i:i)
	  endif
	end do
	if(safile.eq.'') then
	  write(*,*) 'Illegal model code = ',mcode
	  stop
	endif
	call saextr(safile,stn,psa,ssa)
c
	return
	end

      subroutine saextr(safile,stn,psa,ssa)
c
c ...... extract station adjustments
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      character*80 safile
      character*4 stn,astn
      real*8 psa,ssa,psat,ssat
c
	lnb=lnblnk(safile)
	if(lnb.gt.0) then
	  open(21,file=safile(1:lnblnk(safile)),status='old')
	  psa=0.
	  ssa=0.
	else
	  psa=0.
	  ssa=0.
	  return
	endif
    1 read(21,'(1x,a4,2f6.2)',end=2) astn,psat,ssat
      iastn=16777216*ichar(astn(1:1))+65536*ichar(astn(2:2))+
     1           256*ichar(astn(3:3))+      ichar(astn(4:4))
      istn=16777216*ichar(stn(1:1))+65536*ichar(stn(2:2))+
     1          256*ichar(stn(3:3))+      ichar(stn(4:4))
      if(iastn.eq.istn) then
        psa=psat
        ssa=ssat
        go to 2
      endif
      go to 1
    2 continue
      close(21)
c
      return
      end
