anchor.diff:      !Offset between mass center and hexag centers
       offsethx%n=mesh%nv
       offsethx%pos=0
  -    allocate(offsethx%f(1:offsethx%n))
  -    do i=1,mesh%nv
  -       offsethx%f(i)=arclen(mesh%hx(i)%b%p,mesh%v(i)%p)
  -    end do
  +
  +    if(trim(mesh%pos)=="readref_andes")then
  +       if(.not.allocated(mesh%densf_table))then
  +         print*, "Reading Andes density table ..."
  +         allocate (mesh%densf_table(nlat_alt*nlon_alt, 3))
  +         call getunit(iunit2)
  +         call earth_elevation(iunit2, mesh%densf_table)
  +         !call andes_density_table(mesh%densf_table)
  +         !call andes_density_table2(mesh%densf_table)
  +         call andes_density_table3(mesh%densf_table)
  +       end if
  +
  +         allocate(offsethx%f(1:offsethx%n))
  +         do i=1,mesh%nv
  +           p = vorbarycenterdens(i,mesh) 
  +           offsethx%f(i)=arclen(p,mesh%v(i)%p)
  +         end do
  +       deallocate (mesh%densf_table)
  +    else    
  +        allocate(offsethx%f(1:offsethx%n))
  +        do i=1,mesh%nv
  +          offsethx%f(i)=arclen(mesh%hx(i)%b%p,mesh%v(i)%p)
  +        end do
  +    end if
  +
   
       minoffsethx=minval(offsethx%f(1:offsethx%n))
       maxoffsethx=maxval(offsethx%f(1:offsethx%n))
  
positive.code:     print*, "Set lat (degrees):"
      !read (*, *) tlat
      !Point
      !tlon=-64.5
      !tlat=-13.5
  
      !tlon=35.9121339133462_r8
  
negative.code:        pi, &
         r8, &
         rad2deg, &
         nlat_alt, &
         nlon_alt, &
         simulcase
  
    !Data structures
  