anchor.diff:         if(.not.allocated(mesh%densf_table))then
            allocate (mesh%densf_table(nlat_alt*nlon_alt, 3))
            call getunit(iunit4)
  -         call earth_elevation(iunit4, mesh%densf_table)
  -         call andes_density_table3(mesh%densf_table)
  +         call andes_density_table(mesh%densf_table,iunit4)
          end if
      end if
   
  
positive.code:        if(.not.allocated(mesh%densf_table))then
           allocate (mesh%densf_table(n_lat*n_lon, 3))
           call getunit(iunit3)
           filename2 = trim(datadir)//"densf_table.dat"
           call densftable(filename2, iunit3)
           open(iunit3, file=filename2)
           read(iunit3,*) ((mesh%densf_table(i,j), j=1,3), i=1,n_lat*n_lon)
           close(iunit3) 
         end if
      end if
  
negative.code: 
    !Intepolation for local refinement
    use refinter, only: &
     andes_density_table, &
     dens_andes, &
     densf, &
     densftable, &
     interpol_densf
    implicit none
  
  contains 
  