source: branches/GRISLIv3/SOURCES/celltest_tracer.f90 @ 455

    subroutine CellTest(time_tra, ic, jc, kc, v_x, v_y, v_E, xtracer, &
      ytracer, E_k, im, jm, km, fx, fy, fz)



      !!  Find the cell origin index for the previous time_tra step
      use module3d_phy, only: e,uxbar,uybar,flotmx,flotmy,ux,uy,gzmx,gzmy
      use geography, only: nx,ny,nz,dx,dy
      use tracer_vars, only: xgrid,ygrid

      implicit none

      integer, intent(IN) :: time_tra
      integer, intent(IN) :: ic, jc, kc
      real, intent(IN) :: v_x, v_y, v_E
      real, intent(IN) :: xtracer, ytracer
      real, intent(INOUT) :: E_k
      integer, intent(OUT) :: im, jm, km
      real, intent(OUT) :: fx, fy, fz

      integer :: kk   ! indice vertical pour definition de E, mais on veut conserver la valeur de k
!      real,dimension(nz) :: E   ! vertical coordinate in ice, scaled to H zeta

      E(1)=0.
      E(NZ)=1.
      do KK=1,NZ
         if ((KK.ne.1).and.(KK.ne.NZ)) E(KK)=(KK-1.)/(NZ-1.)
      end do


      if(v_x>0.0) then
        im = ic-1
        do
          if(xtracer>=xgrid(im) .and. xtracer<=xgrid(im+1)) exit
          !! PRINT *,"CellTest(x-): Multiple decrement required (",ic,",",jc,",",kc,")"
          im = im-1
          if(im<1) then
            print *,"CellTest(x-;",time_tra,"): Out-of_range decrement (",ic,",",jc,",",kc,")", v_x
!           STOP
           im =ic
           exit
         end if
        end do
      else
        im = ic
        do
          if(xtracer>=xgrid(im) .and. xtracer<=xgrid(im+1)) exit
          !! PRINT *,"CellTest(x+): Multiple increment required (",ic,",",jc,",",kc,")"
          im = im+1
          if(im==nx-1) then
            print *,"CellTest(x+;",time_tra,"): Out-of_range increment (",ic,",",jc,",",kc,")", v_x
!           STOP
             im = ic -1
             exit
          end if
        end do
      end if

      fx = abs(xtracer-xgrid(im))/dx
  if (fx>1.00)    fx = fx - aint(fx)

!     if((fx>1.0).or.(fx<0.0)) THEN
!       PRINT *,"Cell test(!!!)", ic,jc, kc,im, fx
!     END IF
!----------------------------------------------------------------------------------

      if(v_y>0.0) then
        jm = jc -1
        do
          if(ytracer>=ygrid(jm) .and. ytracer<=ygrid(jm+1)) exit
          !! PRINT *,"CellTest(y-): Multiple decrement required (",ic,",",jc,",",kc,")"
          jm = jm-1
          if(jm<1) then
            print *,"CellTest(y-;",time_tra,"): Out-of_range decrement (",ic,",",jc,",",kc,")", v_y
!           STOP
              jm =jc
              exit
          end if
        end do
      else
        jm = jc
        do
          if(ytracer>=ygrid(jm) .and. ytracer<=ygrid(jm+1)) exit
          !! PRINT *,"CellTest(y+): Multiple increment required (",ic,",",jc,",",kc,")"
          jm = jm+1
          if(jm==ny-1) then
            print *,"CellTest(y+;",time_tra,"): Out-of-range increment (",ic,",",jc,",",kc,")", v_y
!           print*, ytracer, ygrid(jm)
            print*, uybar(ic,jc), uybar(ic,jc+1), flotmy(ic,jc), flotmy(ic, jc+1)
            print *, uy(ic,jc,:)
            print *, uy(ic,jc+1,:)
            print*, gzmy(ic, jc), gzmy(ic,jc+1)
!           STOP
             jm = jc-1
             exit
          end if
        end do
      end if

      fy = abs(ytracer-ygrid(jm))/dy
  if (fy>1.00)    fy = fy - aint(fy)

!     if((fy>1.0).or.(fy<0.0)) THEN
!       PRINT *,"Cell test(!!!)", ic,jc,kc,im,jm, fy
!     END IF
 
!-------------------------------------------------------------------------------

      !! Note that v_E is "velocity" in E coordinate system which is
      !! positive DOWN
      ! also note that in Catritz code k=1 for surface

    
  if(E_k < E(1)) then
        E_k = E(1)
        km=1     
   
  else if(E_k > E(nz)) then
        E_k = E(nz)
        km=nz-1

  else   ! this used to be a loop on its own, now shortcuts if E_k is special 

      if(v_E>=0.0) then   !! Downward flow of ice in E system
        if (kc==1) then
          !! No action, pending tactical discussion with GKCC
          km = 1
        else
          km=kc - 1
          do
            if( (E_k >= E(km)) .and. (E_k <= E(km+1)) ) exit
            km = km - 1
            if(km==0) then
              print *,"CellTest(z+;",time_tra,"): Out-of-range increment (",ic,",",jc,",",kc,")", v_E
!             STOP
              km = kc
              exit
            end if
          end do
        end if
      else                             !! Upward flow of ice in E system
        if(kc==nz) then    
        ! this is incomplete: ice comes from below, cannot use above layer, fix in interpolate
          km = kc-1
        else
          km = kc
         do
          if( (E_k >= E(km)) .and. (E_k <= E(km+1)) ) exit
          !! PRINT *,"CellTest(z,",time_tra,"): Decrementing km=",km
          km = km + 1
          if(km==nz) then
            print *,"CellTest(z-;",time_tra,"): Out-of-range decrement (",ic,",",jc,",",kc,")", v_E
!           STOP
             km = kc-1
             exit
          end if
         end do
       end if

     end if  ! velocity direction

      fz = abs((E_k-E(km))/(E(km+1)-E(km)))
  if (fz>1.00)    fz = fz - aint(fz)

  end if 

    end subroutine CellTest
!===========================================================================