program SHADMIT
      use SHTOOLS

      implicit none
      integer, parameter ::   degmax = 3600
      real*8 ::               c1ilm(2,degmax+1, degmax+1), admit(degmax+1), corr(degmax+1)
      real*8 ::               c2ilm(2,degmax+1, degmax+1), admit_err(degmax+1)
      real*8 ::               topo(2,degmax+1, degmax+1), grav(2,degmax+1, degmax+1)
      real*8 ::               param(2), r, r0_pot, gm
      integer ::              lmax, lmax1, lmax2, i, l
      character*80 ::         shfile1, shfile2
      
      write(*,*) 'Input file name of Spherical Harmonics (Topography)'
      read(*,*) shfile1
      write(*,*) 'Input file name of Spherical Harmonics (Gravity)'
      read(*,*) shfile2
      call SHRead(shfile1, c1ilm, lmax1)
      write(*,*) 'Max. Degree of Topography is ...', lmax1
      call SHRead(shfile2, c2ilm, lmax2, header=param(1:2))
      gm = param(1)
      r0_pot= param(2)
      write(*,*) 'Max. Degree of Gravity is ...', lmax2
      write(*,*) 'Input Max. Degree for Admittance Computation'
      read(*,*) lmax
      topo = 0.d0
      grav = 0.d0
      r = c1ilm(1,1,1)
      do l=0, lmax            
                grav(1:2,l+1,1:l+1) = c2ilm(1:2,l+1,1:l+1) * dble(l+1) * (r0_pot/r)**(l+2)        
      enddo 
      grav = grav * 1.0d5 * gm / (r0_pot**2)
      topo = c1ilm / 1.0d3
      call SHAdmitCorr ( grav, topo, lmax, admit, corr, admit_err )
      open (20,file='shadmit.out',status='unknown')
      do i=3, lmax+1
        write(20,*) i-1, admit(i), corr(i), admit_err(i)
      enddo
      close(20)
      stop
end program  SHADMIT
