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: 
          deallocate(aux)
          deallocate(aux2)
      end subroutine andes_density_table
  
  
      !--------------------------------------------------------------------------
      ! Andes density function
      !--------------------------------------------------------------------------
      real (kind=8) function dens_andes2(x, dens_table) result(densf)
          real (kind=8), dimension(1:2), intent(in) :: x
          real (kind=8), dimension(nlat_alt*nlon_alt,3), intent(in) :: dens_table
          real (kind=8) :: lat, lon
      
          lat = x(1)
          lon = x(2)
          lat = lat*rad2deg
          lon = lon*rad2deg
  
          densf = interpol_densf([lat,lon], dens_table, -90._r8, 90._r8, -180._r8, 180._r8, nlat_alt, nlon_alt)
  
      end function dens_andes2
  
  
  
      real (kind=8) function dens_andes1(x, dens_table) result(densf)
          real (kind=8), dimension(1:2), intent(in) :: x
          real (kind=8), dimension(nlat_alt*nlon_alt,3), intent(in) :: dens_table
          real (kind=8) :: lat, lon, latc, lonc
          real (kind=8) :: dists, maxdist, gammas, epsilons, coslat
          real (kind=8), dimension(1:3) :: p
          real (kind=8), dimension(1:3) :: pc
  
          latc=-20.0*pi/180._r8
          lonc=-60.0*pi/180._r8
          maxdist = pi/6._r8
          epsilons = pi/12._r8
          gammas = 3._r8
  
          lat = x(1)
          lon = x(2)
          coslat = dcos(lat)
          p(1) = coslat*dcos(lon)
          p(2) = coslat*dsin(lon)
          p(3) = dsin(lat)
  
          coslat = dcos(latc)
          pc(1) = coslat*dcos(lonc)
          pc(2) = coslat*dsin(lonc)
          pc(3) = dsin(latc)
  
          dists = dsqrt((p(1)-pc(1))**2+(p(2)-pc(2))**2+(p(3)-pc(3))**2)
  
          if(dists<=maxdist+epsilons)then
              lat = lat*rad2deg
              lon = lon*rad2deg
              densf = interpol_densf([lat,lon], dens_table, -90._r8, 90._r8, -180._r8, 180._r8, nlat_alt, nlon_alt)
          else
              densf = 1._r8/gammas**4
          end if
  
      end function dens_andes1
  
      !--------------------------------------------------------------------------
      ! Andes density function 2
      !--------------------------------------------------------------------------
  
  
      real (kind=8) function dens_andes(x, dens_table) result(dens_f)
          real (kind=8), dimension(1:2), intent(in) :: x
          real (kind=8), dimension(1:3) :: p, pc
          real (kind=8), dimension(nlat_alt*nlon_alt,3), intent(in) :: dens_table
  
negative.code: 
      end subroutine earth_elevation
  
      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
        real (kind=8) :: maxx, lat, lon, latc, lonc
        real (kind=8) :: dists, maxdist, maxdist2, gammas, epsilons, Fo
        real (kind=8), dimension(1:3) :: p
        real (kind=8), dimension(1:3) :: pc
        real (kind=8) :: coslat, a, b, c, sx, corte, alpha, beta,corte2
        real (kind=8), dimension(:,:), allocatable  :: topo, topo2, aux, aux2
        real (kind=8), dimension(:,:), allocatable  :: grad_table, grad
        real (kind=8) :: temp, dx, dy, gx, gy, dlat, dlon, x, lambda
        character (len=100) :: filename
  
        call earth_elevation(iunit, alt_table)
        allocate (topo(nlat_alt, nlon_alt))
        allocate (topo2(nlat_alt, nlon_alt))
        allocate (aux(nlat_alt, nlon_alt))
        allocate (aux2(nlat_alt, nlon_alt))
        
        ! Parameters
        latc=-20.0*deg2rad
        lonc=-70.0*deg2rad
        coslat = dcos(latc)
        pc(1) = coslat*dcos(lonc)
        pc(2) = coslat*dsin(lonc)
        pc(3) = dsin(latc)
  
        maxdist = 45.d0*deg2rad
        maxdist2 = 28.d0*deg2rad
        epsilons = 15.d0*deg2rad
        gammas = 2._r8
  
        corte = 500.d0
        corte2 = 0.d0
  
        dlat = pi/(nlat_alt-1)
        dlon = pi2/(nlon_alt-1)
  
        !call bump(maxdist2*rad2deg,(maxdist2+epsilons)*rad2deg)
        !c = int_phi(1.d0, 0.d0, 1.d0)
  
        ! Saves the topography data
        k = 1
        do i = 1, nlat_alt
          do j = 1, nlon_alt
            topo(i,j) = alt_table(k,3)
            k = k + 1
          end do
        end do
  
        aux = 0.d0
        do i = 1, nlat_alt
          lat = -pio2 + (i-1)*dlat
          do j = 1, nlon_alt
            lon = -pi + (j-1)*dlon
            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)
            !aux3(i,j) = xsi((dists-maxdist2)/epsilons,c)
            if(dists<=maxdist2)then
              aux(i,j) = 1.d0
             else if(dists<=maxdist2+epsilons)then
               sx = (maxdist2+epsilons-dists)/epsilons
               aux(i,j) = sx
            endif
          enddo
        enddo
  
  
  
        do i = 1, nlat_alt
          lat = -pio2 + (i-1)*dlat
          do j = 1, nlon_alt
            lon = -pi + (j-1)*dlon
            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(topo(i,j)<=corte .or. lon*rad2deg>-62.d0 .or.lon*rad2deg<-82.d0 .or.lat*rad2deg>15.d0 .or.lat*rad2deg<-55.d0)then
                topo(i,j) = 0.d0
              ! Remove Mount Roraima
              else if(lat*rad2deg > -5.d0 .and. lat*rad2deg < 10.d0 .and. lon*rad2deg > -68.d0 .and. lon*rad2deg < -50.d0 )then
                topo(i,j) = 0.d0
              end if    
            else 
              topo(i,j) = 0.d0
            end if
  
            if(topo(i,j)>0.d0)then
              topo(i,j) = 1.d0
            endif 
          enddo
        enddo
  
        !topo = topo/maxval(topo)
        !filename = "topo.dat"
        !open(4, file=filename)
        !  do i=1, nlat_alt
        !      write(4,*) topo(i,:)
        !  end do
        !close(4)
  
  
        Fo = 0.25d0
        do l = 1, 500
          do i = 2, nlat_alt-1
            do j = 2, nlon_alt-1
              topo2(i,j) = (1-4.d0*Fo)*topo(i,j) + Fo*(topo(i+1,j)+topo(i-1,j)+topo(i,j+1)+topo(i,j-1))
            enddo
          enddo
          topo = topo2
        enddo
  
        l = 55
        !do i = l, nlat_alt-l
        !  do j = l, nlon_alt-l
        !    temp =  0._r8
        !    do m = 0, l-1
        !      do n = 0, l-1
        !        temp = temp + topo(i-(l-1)/2+m, j-(l-1)/2+n)
        !      end do
        !    end do
        !    topo2(i,j) = temp/l**2
        !  end do
        !end do
  
  
        l = 55!15
        do i = l, nlat_alt-l
          do j = l, nlon_alt-l
            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
            aux2(i,j) = temp/l**2
          end do
        end do
  
        aux = aux2
        aux = aux/maxval(aux)
  
        topo = topo2
        topo = topo/maxval(topo)
  
        !aux = 1.d0/gammas**4 + (1-1.d0/gammas**4)*aux
        !topo = 1.d0/gammas**4 + (1-1.d0/gammas**4)*topo      
        lambda = 0.6d0 
        topo = 1.d0/gammas**4 + (1.d0-1.d0/gammas**4)*(lambda*topo+(1.d0-lambda)*aux)
  
        !filename = "topo2.dat"
        !open(4, file=filename)
        !  do i=1, nlat_alt
        !      write(4,*) topo(i,:)
        !     !write(4,*) aux(i,:)
        !  end do
        !close(4)
  
        !stop
  
        ! Update values in the table
        k = 1
        do i = 1, nlat_alt
          do j = 1, nlon_alt
              alt_table(k,3) = topo(i,j)
              !alt_table(k,3) = aux(i,j)
              k = k + 1
          end do
        end do
  
        deallocate(topo)
        deallocate(topo2)
        deallocate(aux)
        deallocate(aux2)
      end subroutine andes_density_table
  
      subroutine andes_density_table2(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
  