program SHLOCALIZEDADMIT
      use SHTOOLS

      implicit none
      integer, parameter ::   degmax = 3600
      real*8 ::               c1ilm(2,degmax+1, degmax+1), c2ilm(2,degmax+1, degmax+1)
      real*8 ::               topo(2,degmax+1, degmax+1)
      real*8, allocatable ::  grav(:,:,:), admit(:), corr(:), admit_error(:), corr_error(:)
      real*8, allocatable ::  tapers(:,:), eigenvalues(:)
      real*8 ::               param(2), r, r0_pot, gm, pi, theta0, alpha, lat, lon, sn
      integer ::              lmax, lmax1, lmax2, i, l, lwin, lmaxwin, astat(5), K, option1
      integer, allocatable :: taper_order(:)
      character*80 ::         shfile1, shfile2

      pi = acos(-1.d0)
      
      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
      
!============= Localized options

      print*, "Number of tapers to use > "
      read(*,*) K

      print*, "Latitude of feature of interest (degrees) > "  
      read(*,*) lat           
      print*, "Longitude (degrees) > "        
      read(*,*) lon                   
      print*, "Angular radius of localization window (degrees) > "    
      read(*,*) theta0        
      theta0 = theta0*pi/180.0d0
                        
      print*, "Create localization window using"      
      print*, "(1) Desired concentration factor (alpha^2)"    
      print*, "(2) Desired (approximate) Shannon number; (Lwin+1) Theta0 / pi)"       
      print*, "(3) Desired Spectral bandwidth (Lwin)" 
      read(*,*) option1               
        
      if (option1==1) then            
              print*, "Input desired concentration factor of the Kth taper > "                
              read(*,*) alpha         
              lwin = SHFindLWin(theta0, 0, alpha, K)          
              print*, "Corresponding spherical harmonic bandwidth = ", lwin           
              print*, "Corresponding approximate Shannon number = ", (lwin+1)*theta0/pi       
      elseif (option1==2) then                
              print*, "Input Shannon number > "               
              read(*,*) sn            
              lwin = nint(sn*pi/theta0) - 1           
              print*, "Corresponding spherical harmonic bandwidth = ", lwin   
      elseif (option1==3) then                
              print*, "Input Lwin > "         
              read(*,*) lwin          
              print*, "Corresponding approximate Shannon number = ", (lwin+1)*theta0/pi       
      else            
              stop    
      endif           
      
      allocate(tapers(lwin+1, (lwin+1)**2), stat = astat(1))
      allocate(taper_order((lwin+1)**2), stat = astat(2))
      allocate(eigenvalues((lwin+1)**2), stat = astat(3))
        
      if (astat(1) /= 0 .or. astat(2) /= 0 .or. astat(3) /= 0 ) then
              print*, "Problem allocatig arrays for tapers, taper_order, or eigenvalues", astat(1), astat(2), astat(3)
              stop
      endif
        
      call SHReturnTapers(theta0, lwin, tapers, eigenvalues, taper_order)
        
      print*, "Concentration factor of first taper = ", eigenvalues(1)
      print*, "Taper order = ", taper_order(1)
        
      lmaxwin = lmax+lwin             
      r = c1ilm(1,1,1)

!===============

     allocate(grav(2,lmax+1,lmax+1), stat=astat(1))  
     allocate(admit(lmax+lwin+1), stat=astat(2))     
     allocate(corr(lmax+lwin+1), stat=astat(3))      
     allocate(admit_error(lmax+lwin+1), stat=astat(4))
     allocate(corr_error(lmax+lwin+1), stat=astat(5))
        
     grav = 0.0d0
     admit = 0.0d0
     corr = 0.0d0
     admit_error = 0.0d0
     corr_error = 0.0d0              

     if (astat(1) /= 0 .or. astat(2) /= 0 .or. astat(3) /= 0 .or. astat(4) /= 0 &
             .or. astat(5) /=0) then         
             print*, "Problem allocating memory"             
             stop    
     endif  

!================

      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 SHLocalizedAdmitCorr(tapers, taper_order, lwin, lat, lon, grav, topo, lmax, admit, &               
              corr, K, admit_error=admit_error, corr_error=corr_error) 
      open (20,file='shlocalizedadmit.out',status='unknown')

      if (K==1) then  
              do l = lwin+1, lmax-lwin             
                      write(20,*) l, admit(l+1), corr(l+1)
              enddo   
      else
              do l = lwin+1, lmax - lwin           
                      write(20,*) l, admit(l+1), admit_error(l+1), corr(l+1), corr_error(l+1)
              enddo   
      endif

      close(20)
 
      deallocate(grav)        
      deallocate(admit)       
      deallocate(admit_error) 
      deallocate(corr)
      deallocate(corr_error)
      deallocate(tapers)
      deallocate(taper_order)
      deallocate(eigenvalues)


      stop
end program  SHLOCALIZEDADMIT
