      subroutine me(x,xm,xe,n)
      implicit real*8(a-h,o-z) 
      dimension x(n)
      xm=0.
      xe=0.
      do 900 ii=1,n
        xm=xm+x(ii)
c       print*,'xm x',xm,x(ii)
900   continue
      xm=xm/n
      do 800 ii=1,n
        xe=xe+(x(ii)-xm)**2
800   continue
      xe=xe/float(n)
      xe=sqrt(xe)
      return
      end
c
      subroutine ma22b(a,ia,n,w,e)
      implicit real*8(a-h,o-z)
      dimension a(1),w(1)
      print*,'entering ma22b',n
      call a22m(a,a,ia,n,w,e)
      return
      end
c
      subroutine a22m(a,b,ia,n,w,e)
      implicit real*8(a-h,o-z)
      dimension a(ia,1),w(1),b(1)
      l=0
      eps= 0.5d-7
401   format('a22m')
      print*,'entering a22m',n
      do 1 j=1,n
        do 1 i=1,j
          l=l+1
          b(l)=a(i,j)
1     continue
402   format('sinv')
      call sinv(a,n,eps,ierr)
      e=ierr
      ii=n
      jj=n
20    ii=jj
24    a(ii,jj)=b(l)
      l=l-1
      ii=ii-1
      if(ii)22,22,24
22    jj=jj-1
      if(jj)30,30,20
30    do 40  j=1,n
        do 40 i=1,j
          a(j,i)=a(i,j)
40    continue
      return
      end
c
      subroutine sinv(a,n,eps,ier)
      implicit real*8(a-h,o-z)
      real*8 din,work
      dimension a(1)
      print*,'entering sinv',n
      call mfsd(a,n,eps,ier)
      if (ier) 9,1,1
1     ipiv=n*(n+1)/2
      ind=ipiv
      do 6 i=1,n
      din=1.d0/dble(a(ipiv))
      a(ipiv)=din
      min=n
      kend=i-1
      lanf=n-kend
      if(kend)5,5,2
2     j=ind
      do 4 k=1,kend
      work=0.d0
      min=min-1
      lhor=ipiv
      lver=j
      do 3 l=lanf,min
      lver=lver+1
      lhor=lhor+l
3     work=work+dble(a(lver)*a(lhor))
      a(j)=-work*din
4     j=j-min
5     ipiv=ipiv-min
6     ind=ind-1
      do 8 i=1,n
      ipiv=ipiv+i
      j=ipiv
      do 8 k=i,n
      work=0.d0
      lhor=j
      do 7 l=k,n
      lver=lhor+k-i
      work=work+dble(a(lhor)*a(lver))
7     lhor=lhor+l
      a(j)=work
8     j=j+k
9     return
      end
c
      subroutine mfsd(a,n,eps,ier)
      implicit real*8(a-h,o-z)
      real*8 dpiv,dsum
      dimension a(1)
      print*,'entering mfsd',n
      if(n-1) 12,1,1
1     continue
      print*,'passing instruction1, n=',n
      ier=0
      kpiv=0
      do 11 k=1,n
      kpiv=kpiv+k
      ind=kpiv
      lend=k-1
      tol=abs(eps*a(kpiv))
      do 11 i=k,n
      dsum=0.d0
      if(lend) 2,4,2
2     do 3 l=1,lend
      lanf=kpiv-l
      lind=ind-l
3     dsum=dsum+dble(a(lanf)*a(lind))
4     dsum=dble(a(ind))-dsum
      if(i-k) 10,5,10
5     if(sngl(dsum)-tol) 6,6,9
6     if(dsum) 12,12,7
7     if(ier) 8,8,9
8     ier=k-1
9     dpiv=sqrt(dsum)
      a(kpiv)=dpiv
      dpiv=1.d0/dpiv
      goto 11
10    a(ind)=dsum*dpiv
11    ind=ind+i
      if(ier.ne.0) print 21,ier
      return
12    ier=-1
      print 22,n
21    format (' ret',i4)
22    format (' break, n=',i5)
      return
      end
      function de(t1,p1,t2,p2)
       implicit real*8(a-h,o-z)
      rad=0.017453293
      t=(t1+t2)*rad/2.
      de=((t2-t1)**2)+((p2-p1)**2)*(sin(t)**2)
      de=sqrt(de)
      return
      end
c
      subroutine pole1(eplat,eplong,stlat,stlong,xlatp,xlongp,phmp
     *  ,delta)
      implicit real*8(a-h,o-z)
      data radian/57.2957795/
c
c    returns latitude (xlatp) and longitude (xlonp) of pole in degrees
c
      th1=(90.-eplat)/radian
      th2=(90.-stlat)/radian
      ph1=eplong/radian
      ph2=stlong/radian
c
      cth1=cos(th1)
      cth2=cos(th2)
      sth1=sin(th1)
      sth2=sin(th2)
      cph1=cos(ph1)
      cph2=cos(ph2)
      sph1=sin(ph1)
      sph2=sin(ph2)
c
      cph21=cph1*cph2+sph1*sph2
      sph21=sph2*cph1-sph1*cph2
c
      cdel=sth1*sth2*cph21+cth1*cth2
      sdel=sqrt(1.-cdel**2)
      delta=atan2(sdel,cdel)*radian
c
      ccth=sth1*sth2*sph21/sdel
      scth=sqrt(1.-ccth**2)
      cth=atan2(scth,ccth)
c
      coc=1./(sdel*scth)
      scph=coco*(cth1*sth2*cph2-cth2*sth1*cph1)
      ccph=coco*(sth1*cth2*sph1-sth2*cth1*sph2)
c
      cph=atan2(scph,ccph)
      xlatp=90.-cth*radian
      xlongp=cph*radian
c
      v1=sth1*cph1+sth2*cph2
      v2=sth1*sph1+sth2*sph2
      v3=cth1+cth2
      vn=sqrt(v1**2+v2**2+v3**2)
      v1=v1/vn
      v2=v2/vn
      v3=v3/vn
      cthm=v3
      sthm=sqrt(1.-v3**2)
      cphm=v1/sthm
      sphm=v2/sthm
c
      cphmp=ccth*sthm*(cphm*ccph+sphm*scph)-scth*cthm
      sphmp=sthm*(sphm*ccph-cphm*scph)
      phmp=atan2(sphmp,cphmp)
      azmp=180.-phmp*radian
c
      cpsi=sth1*sph1*ccph-scph*cph1*sth1
      spsi=cth1/scth
      psi=atan2(spsi,cpsi)
      psi=psi*radian+90.
c
100   format(' phi,lam,phmp,psi,delta',5f8.2)
      psi1=180.-azmp-psi
      psi1=2.*psi1
      return
      end
c
      subroutine trapez(npoint,f,ds,sint)
      implicit real*8(a-h,o-z) 
      dimension f(1)
c    integration using the trapezoid method
      sint=0.
      nmax=npoint-1
      do 1 ns=2,nmax
1     sint=sint+2.*f(ns)
      sint=sint+f(1)+f(npoint)
      sint=sint*ds*0.5
      return
      end
c
      function ynterp(x,y,xi,n,nd)
c x(i) must be all different
c hiroo kanamori aug.4 1978
      implicit real*8 (a-h,o-z)
      dimension x(1), y(1)
      i=1
      j=n
      k=1
      if ( x(1) .gt. x(n) ) go to 11
1     if( xi-x(k)) 2,2,3
2     j=k
      go to 4
3     i=k
4     k=(i+j)/2
      if(k.ne.i) go to 1
      if( xi.gt.x(n)) j=n+1
      go to 21
11    if( x(k)-xi) 12,12,13
12    j=k
      go to 14
13    i=k
14    k=(i+j)/2
      if(k.ne.i) go to 11
      if( xi.lt. x(n)) j=n+1
21    if ( j.eq. 1) go to 23
      if(j.eq.n+1) go to 24
      ynterp=y(j-1)+(y(j)-y(j-1))*(xi-x(j-1))/(x(j)-x(j-1))
      return
23    ynterp=y(1)
      return
24    ynterp=y(n)
      return
      end
c
      subroutine pole(eplat,eplong,stlat,stlong,xlatp,xlongp,azmp
     *  ,delta)
      implicit real*8(a-h,o-z)
      data radian/57.2957795/
c
      th1=(90.-eplat)/radian
      th2=(90.-stlat)/radian
      ph1=eplong/radian
      ph2=stlong/radian
c
      cth1=cos(th1)
      cth2=cos(th2)
      sth1=sin(th1)
      sth2=sin(th2)
      cph1=cos(ph1)
      cph2=cos(ph2)
      sph1=sin(ph1)
      sph2=sin(ph2)
c
      cph21=cph1*cph2+sph1*sph2
      sph21=sph2*cph1-sph1*cph2
c
      cdel=sth1*sth2*cph21+cth1*cth2
      sdel=sqrt(1.-cdel**2)
      delta=atan2(sdel,cdel)*radian
c
      ccth=sth1*sth2*sph21/sdel
      scth=sqrt(1.-ccth**2)
c
      cth=atan2(scth,ccth)
c
      coco=1./(sdel*scth)
      scph=coco*(cth1*sth2*cph2-cth2*sth1*cph1)
      ccph=coco*(sth1*cth2*sph1-sth2*cth1*sph2)
c
      cph=atan2(scph,ccph)
c
      xlatp=90.-cth*radian
      xlongp=cph*radian
c
      v1=sth1*cph1+sth2*cph2
      v2=sth1*sph1+sth2*sph2
      v3=cth1+cth2
      vn=sqrt(v1**2+v2**2+v3**2)
      v1=v1/vn
      v2=v2/vn
      v3=v3/vn
      cthm=v3
      sthm=sqrt(1.-v3**2)
      cphm=v1/sthm
      sphm=v2/sthm
c
      cphmp=ccth*sthm*(cphm*ccph+sphm*scph)-scth*cthm
      sphmp=sthm*(sphm*ccph-cphm*scph)
c
      phmp=atan2(sphmp,cphmp)
c
      azmp=180.-phmp*radian
c
      cpsi=sth1*sph1*ccph-scph*cph1*sth1
      spsi=cth1/scth
      psi=atan2(spsi,cpsi)
      psi=psi*radian+90.
c
100   format(' phi,lam,phmp,psi,delta',5f8.2)
      psi1=180.-azmp-psi
      psi1=2.*psi1
      return
      end
