//------------------------------ Trace C++ ------------------------------


int Trace(float h0,int& k,float x[],float y[])
{ int i,j;
  int i1,j1,i2,j2,i3,j3,i4,j4;

  k=0;

  for(i=1;i<Mr;i++)                 // 1
  { if((h0-Ht[i][1])*(h0-Ht[i+1][1]) < 0.0 && Fg1[i][1]==0)
    { k=1;
      x[1]=1.0;
      y[1]=ConAv(h0,Ht[i][1],Ht[i+1][1],i,i+1);
      i1=i;j1=1;
      i2=i+1;j2=1;
      i3=i2;j3=2;
      i4=i1;j4=2;
      Fg1[i][1]=1;
    
      k=2;
      if((h0-Ht[i2][j2])*(h0-Ht[i3][j3]) < 0.0)
      { x[2]=ConAv(h0,Ht[i2][j2],Ht[i3][j3],j2,j3);
        y[2]=i2;
        i1=i3;j1=j3;
      }
      if((h0-Ht[i4][j4])*(h0-Ht[i3][j3]) < 0.0)
      { y[2]=ConAv(h0,Ht[i4][j4],Ht[i3][j3],i4,i3);
        x[2]=2;
        i1=i3;j1=j3;
        i2=i4;j2=j4;
      }
      if((h0-Ht[i4][j4])*(h0-Ht[i1][j1]) < 0.0)
      { x[2]=ConAv(h0,Ht[i4][j4],Ht[i1][j1],j4,j1);
        y[2]=i1;
        i2=i4;j2=j4;
      }
      break; 
    }
  }

  if(k==0)                 // 2
  { for(j=1;j<Nc;j++)
    { if((h0-Ht[1][j])*(h0-Ht[1][j+1]) < 0.0 && Fg2[j]==0)
      { k=1;
        y[1]=1.0;
        x[1]=ConAv(h0,Ht[1][j],Ht[1][j+1],j,j+1);
        i1=1;j1=j;
        i2=1;j2=j+1;
        i3=2;j3=j2;
        i4=2;j4=j1;
        Fg2[j]=1;
    
        k=2;
        if((h0-Ht[i2][j2])*(h0-Ht[i3][j3]) < 0.0)
        { y[2]=ConAv(h0,Ht[i2][j2],Ht[i3][j3],i2,i3);
          x[2]=j2;
          i1=i3;j1=j3;
        }
        if((h0-Ht[i4][j4])*(h0-Ht[i3][j3]) < 0.0)
        { x[2]=ConAv(h0,Ht[i4][j4],Ht[i3][j3],j4,j3);
          y[2]=2;
          i1=i3;j1=j3;
          i2=i4;j2=j4;
        }
        if((h0-Ht[i4][j4])*(h0-Ht[i1][j1]) < 0.0)
        { y[2]=ConAv(h0,Ht[i4][j4],Ht[i1][j1],i4,i1);
          x[2]=j1;
          Fg1[i1][j1]=1;
          i2=i4;j2=j4;
        }
        break; 
      }
    }
  }

  if(k==0)                 // 3
  { for(i=1;i<Mr;i++)
    { if((h0-Ht[i][Nc])*(h0-Ht[i+1][Nc) < 0.0 && Fg1[i][Nc]==0)
      { k=1;
        x[1]=Nc;
        y[1]=ConAv(h0,Ht[i][Nc],Ht[i+1][Nc],i,i+1);
        i1=i;j1=Nc;
        i2=i+1;j2=Nc;
        i3=i2;j3=Nc-1;
        i4=i1;j4=Nc-1;
        Fg1[i][Nc]=1;
    
        k=2;
        if((h0-Ht[i2][j2])*(h0-Ht[i3][j3]) < 0.0)
        { x[2]=ConAv(h0,Ht[i2][j2],Ht[i3][j3],j2,j3);
          y[2]=i2;
          i1=i3;j1=j3;
        }
        if((h0-Ht[i4][j4])*(h0-Ht[i3][j3]) < 0.0)
        { y[2]=ConAv(h0,Ht[i4][j4],Ht[i3][j3],i4,i3);
          x[2]=Nc-1;
          Fg1[i4][Nc-1]=1;
          i1=i3;j1=j3;
          i2=i4;j2=j4;
        }
        if((h0-Ht[i4][j4])*(h0-Ht[i1][j1]) < 0.0)
        { x[2]=ConAv(h0,Ht[i4][j4],Ht[i1][j1],j4,j1);
          y[2]=i1;
          i2=i4;j2=j4;
        }
        break; 
      }
    }
  }

  if(k==0)                 // 4
  { for(j=1;j<Nc;j++)
    { if((h0-Ht[Mr][j])*(h0-Ht[Mr][j+1]) < 0.0 && Fg1[Mr][j]==0)
      { k=1;
        y[1]=Mr;
        x[1]=ConAv(h0,Ht[Mr][j],Ht[Mr][j+1],j,j+1);
        i1=Mr;j1=j;
        i2=Mr;j2=j+1;
        i3=Mr-1;j3=j2;
        i4=Mr-1;j4=j1;
        Fg1[Mr][j]=1;
    
        k=2;
        if((h0-Ht[i2][j2])*(h0-Ht[i3][j3]) < 0.0)
        { y[2]=ConAv(h0,Ht[i2][j2],Ht[i3][j3],i2,i3);
          x[2]=j2;
          Fg1[i3][j3]=1;
          i1=i3;j1=j3;
        }
        if((h0-Ht[i4][j4])*(h0-Ht[i3][j3]) < 0.0)
        { x[2]=ConAv(h0,Ht[i4][j4],Ht[i3][j3],j4,j3);
          y[2]=Mr-1;
          i1=i3;j1=j3;
          i2=i4;j2=j4;
        }
        if((h0-Ht[i4][j4])*(h0-Ht[i1][j1]) < 0.0)
        { y[2]=ConAv(h0,Ht[i4][j4],Ht[i1][j1],i4,i1);
          x[2]=j1;
          Fg1[i4][j4]=1;
          i2=i4;j2=j4;
        }
        break; 
      }
    }
  }

  for(;;)
  { if(k!=0)
    { if(x[k]==1.0 || x[k]==(float)Nc) return(0);
      if(y[k]==1.0)
      { if(j1<j2) Fg2[j1]=1;
        else      Fg2[j2]=1;
        return(0);
      }
      if(y[k]==(float)Mr)
      { if(j1<j2) Fg1[Mr][j1]=1;
        else      Fg1[Mr][j2]=1;
        return(0);
      }

      if(j1==j2)
      { if(x[k-1]<x[k])
        { i3=i2;j3=j2+1;
          i4=i1;j4=j3;
        }
        else
        { i3=i2;j3=j2-1;
          i4=i1;j4=j3;
        }
      }
      else
      { if(y[k-1]<y[k])
        { i3=i2+1;j3=j2;
          i4=i3;j4=j1;
        }
        else
        { i3=i2-1;j3=j2;
          i4=i3;j4=j1;
        }
      }

      if(j1==j2)
      { k=k+1;
        if((h0-Ht[i2][j2])*(h0-Ht[i3][j3]) < 0.0)
        { x[k]=ConAv(h0,Ht[i2][j2],Ht[i3][j3],j2,j3);
          y[k]=i2;
          i1=i3;j1=j3;
        }
        if((h0-Ht[i4][j4])*(h0-Ht[i3][j3]) < 0.0)
        { y[k]=ConAv(h0,Ht[i4][j4],Ht[i3][j3],i4,i3);
          x[k]=j3;
          if(i3<i4) Fg1[i3][j3]=1;
          else      Fg1[i4][j4]=1;
          i1=i3;j1=j3;
          i2=i4;j2=j4;
        }
        if((h0-Ht[i4][j4])*(h0-Ht[i1][j1]) < 0.0)
        { x[k]=ConAv(h0,Ht[i4][j4],Ht[i1][j1],j4,j1);
          y[k]=i1;
          i2=i4;j2=j4;
        }
      }
      else
      { k=k+1;
        if((h0-Ht[i2][j2])*(h0-Ht[i3][j3]) < 0.0)
        { y[k]=ConAv(h0,Ht[i2][j2],Ht[i3][j3],i2,i3);
          x[k]=j2;
          if(i2<i3) Fg1[i2][j2]=1;
          else      Fg1[i3][j3]=1;
          i1=i3;j1=j3;
        }
        if((h0-Ht[i4][j4])*(h0-Ht[i3][j3]) < 0.0)
        { x[k]=ConAv(h0,Ht[i4][j4],Ht[i3][j3],j4,j3);
          y[k]=i3;
          i1=i3;j1=j3;
          i2=i4;j2=j4;
        }
        if((h0-Ht[i4][j4])*(h0-Ht[i1][j1]) < 0.0)
        { y[k]=ConAv(h0,Ht[i4][j4],Ht[i1][j1],i4,i1);
          x[k]=j1;
          if(i1<i4) Fg1[i1][j1]=1;
          else      Fg1[i4][j4]=1;
          i2=i4;j2=j4;
        }
      }
    }
  }


  for(j=2;j<=Nc-1;j++)
  { for(i=1;i<=Mr-1;i++)
    { if((h0-Ht[i][j])*(h0-Ht[i+1][j) < 0.0 && Fg1[i][j]==0)
      { k=1;
        x[1]=j;
        y[1]=ConAv(h0,Ht[i][j],Ht[i+1][j],i,i+1);
        i1=i;j1=j;
        i2=i+1;j2=j;
        i3=i2;j3=j+1;
        i4=i1;j4=j+1;
        Fg1[i][j]=1;
    
        k=2;
        if((h0-Ht[i2][j2])*(h0-Ht[i3][j3]) < 0.0)
        { x[2]=ConAv(h0,Ht[i2][j2],Ht[i3][j3],j2,j3);
          y[2]=i2;
          i1=i3;j1=j3;
          goto FM10;
        }
        if((h0-Ht[i4][j4])*(h0-Ht[i3][j3]) < 0.0)
        { y[2]=ConAv(h0,Ht[i4][j4],Ht[i3][j3],i4,i3);
          x[2]=j3;
          Fg1[i4][j4]=1;
          i1=i3;j1=j3;
          i2=i4;j2=j4;
          goto FM10;
        }
        if((h0-Ht[i4][j4])*(h0-Ht[i1][j1]) < 0.0)
        { x[2]=ConAv(h0,Ht[i4][j4],Ht[i1][j1],j4,j1);
          y[2]=i1;
          i2=i4;j2=j4;
          goto FM10;
        }
      }
    }
  }
  return(-1);



FM10:   ;

  for(;;)
  { if(x[1]==x[k] && y[1]==y[k]) return(0);

    if(j1==j2)
    { if(x[k-1]<x[k])
      { i3=i2;j3=j2+1;
        i4=i1;j4=j3;
      }
      else
      { i3=i2;j3=j2-1;
        i4=i1;j4=j3;
      }
    }
    else
    { if(y[k-1]<y[k])
      { i3=i2+1;j3=j2;
        i4=i3;j4=j1;
      }
      else
      { i3=i2-1;j3=j2;
        i4=i3;j4=j1;
      }
    }

    if(j1==j2)
    { k=k+1;
      if((h0-Ht[i2][j2])*(h0-Ht[i3][j3]) < 0.0)
      { x[k]=ConAv(h0,Ht[i2][j2],Ht[i3][j3],j2,j3);
        y[k]=i2;
        i1=i3;j1=j3;
      }
      if((h0-Ht[i4][j4])*(h0-Ht[i3][j3]) < 0.0)
      { y[k]=ConAv(h0,Ht[i4][j4],Ht[i3][j3],i4,i3);
        x[k]=j3;
        if(i3<i4) Fg1[i3][j3]=1;
        else      Fg1[i4][j4]=1;
        i1=i3;j1=j3;
        i2=i4;j2=j4;
      }
      if((h0-Ht[i4][j4])*(h0-Ht[i1][j1]) < 0.0)
      { x[k]=ConAv(h0,Ht[i4][j4],Ht[i1][j1],j4,j1);
        y[k]=i1;
        i2=i4;j2=j4;
      }
    }
    else
    { k=k+1;
      if((h0-Ht[i2][j2])*(h0-Ht[i3][j3]) < 0.0)
      { y[k]=ConAv(h0,Ht[i2][j2],Ht[i3][j3],i2,i3);
        x[k]=j2;
        if(i2<i3) Fg1[i2][j2]=1;
        else      Fg1[i3][j3]=1;
        i1=i3;j1=j3;
      }
      if((h0-Ht[i4][j4])*(h0-Ht[i3][j3]) < 0.0)
      { x[k]=ConAv(h0,Ht[i4][j4],Ht[i3][j3],j4,j3);
        y[k]=i3;
        i1=i3;j1=j3;
        i2=i4;j2=j4;
      }
      if((h0-Ht[i4][j4])*(h0-Ht[i1][j1]) < 0.0)
      { y[k]=ConAv(h0,Ht[i4][j4],Ht[i1][j1],i4,i1);
        x[k]=j1;
        if(i1<i4) Fg1[i1][j1]=1;
        else      Fg1[i4][j4]=1;
        i2=i4;j2=j4;
      }
    }    
  }

}


//------------------------------ ConAv C++ ------------------------------

float ConAv(float h0,float h1,float h2,int ij1,int ij2)
{ float xy;

  if(ij1<ij2) xy=ij1+(h0-h1)/(h2-h1);
  else        xy=ij2+(h0-h2)/(h1-h2);
  return(xy);
}