anchor.diff:            end do
           end do
   
  -        !print*, minval(alt_table(:,3)),maxval(alt_table(:,3))
  -        !Output data to a file
  -        !filename = "topo2.dat"
  -        !open(11, file=filename)
  -        !  do i=1, nlat_alt*nlon_alt
  -        !    write(11,*) alt_table(i,1), alt_table(i,2), alt_table(i,3)
  -        !  end do
  -        !close(11)
   
           deallocate(aux)
           deallocate(aux2)
  
positive.code: 
      end subroutine earth_elevation
  
  
      !--------------------------------------------------------------------------
      ! Moving average
      !--------------------------------------------------------------------------
      subroutine moving_average1(alt_table, gammas)
          real (kind=8), dimension(nlat_alt*nlon_alt,3), intent(inout) :: alt_table
          real (kind=8), intent(in) :: gammas
          integer :: i, j, k, m, n
          integer :: l, l2, l3, l4,l2aux,l3aux
          real (kind=8), dimension(:,:), allocatable  :: aux
          real (kind=8), dimension(:,:), allocatable  :: aux2
          real (kind=8) :: temp, alpha
  
          allocate (aux(nlat_alt, nlon_alt))
          allocate (aux2(nlat_alt, nlon_alt))
     
          k = 1
          do i = 1, nlat_alt
              do j = 1, nlon_alt
           
                  aux(i,j) = alt_table(k,3)
                  aux2(i,j) = alt_table(k,3)
                  k = k + 1
              end do
          end do
  
          l = 50
          l2 = 90  !90
          l3 = 90 !90
          l4 = 20
  
          do i = l, nlat_alt-l
              do j = l, nlon_alt-l
  
                  if(aux(i,j) == 500._r8/6320._r8)then
  
                      temp =  0._r8
                      do m = 0,l2-1
                          do n = 0,l2-1
                              temp = temp + aux(i-(l2-1)/2+m, j-(l2-1)/2+n)
                          end do
                      end do
                      temp = temp - aux(i,j)
                      aux2(i,j) = temp/(l2**2-1)
  	
                  end if
   
                  alpha = 5500._r8/6320._r8
                  !alpha = alpha**(0.25)
                  if(aux(i,j) == alpha*1._r8/gammas**4)then
                      ! if(aux(i,j) == 1._r8/gammas**4)then
  
                      temp =  0._r8
                      do m = 0,l3-1
                          do n = 0,l3-1
                              temp = temp + aux(i-(l3-1)/2+m, j-(l3-1)/2+n)
                          end do
                      end do
                      temp = temp - aux(i,j)
                      aux2(i,j) = temp/(l3**2-1)
  
                  end if
  
  
              !if(aux(i,j) > 500._r8/6320._r8)then
  
              !         temp =  0._r8
              !          do m = 0,l4-1
              !          do n = 0,l4-1
              !           temp = temp + aux(i-(l4-1)/2+m, j-(l4-1)/2+n)
              !        end do
              !     end do
              !    temp = temp - aux(i,j)
               !   aux2(i,j) = temp/(l4**2-1)
   
              !end if
  
  
  
  
              end do
          end do
  
          k = 1
          do i = 1, nlat_alt
              do j = 1, nlon_alt
                  alt_table(k,3) = aux2(i,j)
                  k = k + 1
              end do
          end do
  
          deallocate(aux)
          deallocate(aux2)
  
      end subroutine moving_average1
  
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
      subroutine moving_average2(alt_table, gammas)
          real (kind=8), dimension(nlat_alt*nlon_alt,3), intent(inout) :: alt_table
          real (kind=8), intent(in) :: gammas
          integer :: i, j, k, m, n
          integer :: l, l2
          real (kind=8), dimension(:,:), allocatable  :: aux
          real (kind=8), dimension(:,:), allocatable  :: aux2
          real (kind=8) :: temp
  
  
          allocate (aux(nlat_alt, nlon_alt))
          allocate (aux2(nlat_alt, nlon_alt))
     
  
          k = 1
          do i = 1, nlat_alt
              do j = 1, nlon_alt
           
                  aux(i,j) = alt_table(k,3)
                  aux2(i,j) = alt_table(k,3)
                  k = k + 1
              end do
          end do
  
          l = 17
          do i = l, nlat_alt-l
              do j = l, nlon_alt-l
  
  
                  if(aux(i,j) > 500._r8/6320._r8)then
  
                      temp =  0._r8
                      do m = 0,l-1
                          do n = 0,l-1
                              temp = temp + aux(i-(l-1)/2+m, j-(l-1)/2+n)
                          end do
                      end do
                      temp = temp - aux(i,j)
                      aux2(i,j) = temp/(l**2-1)
  
                  end if
  
              end do
          end do
  
          k = 1
          do i = 1, nlat_alt
              do j = 1, nlon_alt
                  alt_table(k,3) = aux2(i,j)
                  k = k + 1
              end do
          end do
  
          deallocate(aux)
          deallocate(aux2)
  
      end subroutine moving_average2
  
  
      !--------------------------------------------------------------------------
      ! Andes density data for interpolation
      !--------------------------------------------------------------------------
      subroutine andes_density_table3(alt_table)
          real (kind=8), dimension(nlat_alt*nlon_alt,3), intent(inout) :: alt_table
          integer :: i, j
          real (kind=8) :: maxx, minn, lat, lon, latc, lonc
          real (kind=8) :: dists, maxdist, gammas, epsilons
          real (kind=8), dimension(1:3) :: p
          real (kind=8), dimension(1:3) :: pc
          real (kind=8) :: coslat, alpha
  
          latc=-20.0*pi/180._r8
          lonc=-60.0*pi/180._r8
          maxdist = pi/6._r8
          epsilons = pi/12._r8
          alpha = 5500._r8/6320._r8
          !alpha = alpha**(0.25)
          minn = 0._r8
          maxx = 0._r8
  
          latc=-20.0*pi/180._r8
          lonc=-60.0*pi/180._r8
          coslat = dcos(latc)
          pc(1) = coslat*dcos(lonc)
          pc(2) = coslat*dsin(lonc)
          pc(3) = dsin(latc)
  
  
  
          do i = 1, nlat_alt*nlon_alt
              lat = alt_table(i,1)
              lon = alt_table(i,2)
              lat = lat*deg2rad
              lon = lon*deg2rad
  
              coslat = dcos(lat)
              p(1) = coslat*dcos(lon)
              p(2) = coslat*dsin(lon)
              p(3) = dsin (lat)
  
              !Distance to center
              dists = dsqrt((p(1)-pc(1))**2+(p(2)-pc(2))**2+(p(3)-pc(3))**2)
              ! dists = dsqrt((lat-latc)**2+(lon-lonc)**2)
        
              if(dists<=maxdist)then
                  if(alt_table(i,3) < minn) then
                      minn = alt_table(i,3)
  
                  else if(alt_table(i,3) > maxx) then
                      maxx = alt_table(i,3)
                  end if
              end if
          end do
  
  
          gammas = 3._r8
  
          do i = 1, nlat_alt*nlon_alt
              lat = alt_table(i,1)
              lon = alt_table(i,2)
              lat = lat*deg2rad
              lon = lon*deg2rad
  
              coslat = dcos(lat)
              p(1) = coslat*dcos(lon)
              p(2) = coslat*dsin(lon)
              p(3) = dsin(lat)
  
              !Distance to center
              dists = dsqrt((p(1)-pc(1))**2+(p(2)-pc(2))**2+(p(3)-pc(3))**2)
              !dists = dsqrt((lat-latc)**2+(lon-lonc)**2)
  
              if(dists<=maxdist)then
                  if(alt_table(i,3) <= 500._r8 )then
                      alt_table(i,3) = 500._r8
                  end if
   
                  alt_table(i,3) = alt_table(i,3)/maxx
  
              else
                  alt_table(i,3) = alpha*1._r8/gammas**4
                 !alt_table(i,3) = 1._r8/gammas**4
              end if
  
          end do
  
  
          call moving_average1(alt_table, gammas)
          call moving_average2(alt_table, gammas)
  
      end subroutine andes_density_table3
  
  
      !--------------------------------------------------------------------------
      !--------------------------------------------------------------------------
      !--------------------------------------------------------------------------
      !--------------------------------------------------------------------------
  
  
      subroutine andes_density_table2(alt_table)
          real (kind=8), dimension(nlat_alt*nlon_alt,3), intent(inout) :: alt_table
          integer :: i, j
          real (kind=8) :: maxx, minn, lat, lon, latc, lonc
          real (kind=8) :: dists, maxdist, gammas, epsilons
          real (kind=8), dimension(1:3) :: p
          real (kind=8), dimension(1:3) :: pc
          real (kind=8) :: coslat
  
          latc=-20.0*pi/180._r8
          lonc=-60.0*pi/180._r8
          maxdist = pi/6._r8
          epsilons = pi/12._r8
  
          coslat = dcos(latc)
          pc(1) = coslat*dcos(lonc)
          pc(2) = coslat*dsin(lonc)
          pc(3) = dsin(latc)
  
          minn = 0._r8
          maxx = 0._r8
  
  
          do i = 1, nlat_alt*nlon_alt
              lat = alt_table(i,1)
              lon = alt_table(i,2)
              lat = lat*deg2rad
              lon = lon*deg2rad
  
              coslat = dcos(lat)
              p(1) = coslat*dcos(lon)
              p(2) = coslat*dsin(lon)
              p(3) = dsin(lat)
        
              dists = dsqrt((p(1)-pc(1))**2+(p(2)-pc(2))**2+(p(3)-pc(3))**2)
  
              if(dists<=maxdist)then
                  if(alt_table(i,3) > maxx) then
                      maxx = alt_table(i,3)
                  end if
              end if
          end do
  
  
          do i = 1, nlat_alt*nlon_alt
              if(alt_table(i,3) <= 200._r8) then
                  alt_table(i,3) = 200._r8
              end if
              alt_table(i,3) = alt_table(i,3)/maxx
          end do
      end subroutine andes_density_table2
  
  
      subroutine andes_density_table(alt_table, iunit)
          real (kind=8), dimension(nlat_alt*nlon_alt,3), intent(inout) :: alt_table
          integer, intent(in) :: iunit 
          integer :: i, j, k, n, m, l,l2,l3
  
negative.code:         densf = interpol_densf([lat,lon], alt_table, -90._r8, 90._r8, -180._r8, 180._r8, nlat_alt, nlon_alt)
      end function smooth_andean_mountain
  
  
  end module refinter
  