source: trunk/SOURCES/interpolate_tracer.f90 @ 10

subroutine interpolate(ic, jc, kc, im, jm, km, fx,fy,fz, E_k)
  
  use module3d_phy
  use tracer_vars
  use climat_perturb_mois_mod ! on en a besoin, notamment pour NFT

  implicit none
  integer, intent(IN) :: ic, jc, kc, im, jm, km
  real, intent(INOUT) :: fx, fy, fz, E_k
  real :: fx_1, fy_1, fz_1, a_xy, a_xy_p, f_zt, a_tmp, ac , al
  real :: fra1, fra2, fra0, accucumul1, accucumul2, accucumul0, &
       thinning, thinning2, f_k, a_xy_m, fzz, fprime, fsec
  integer :: index1, index2, index0
 
  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 
  
  !! Note that input (im,jm,km) refer to points on (ij,jj) grid
  !! whereas CellTest search was performed on (i,j) grid


  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

  fx_1 = 1.0-fx
  fy_1 = 1.0-fy
  fz_1 = 1.0-fz

  ! start with bilinear interp in horizontal plane
  ! above
  a_xy   = fx_1*fy_1*tdepk(km,im,jm)+fx*fy_1*tdepk(km,im+1,jm)+ &
       fx_1*fy*tdepk(km,im,jm+1)+fx*fy*tdepk(km,im+1,jm+1)
  ! below
  a_xy_p = fx_1*fy_1*tdepk(km+1,im,jm)+fx*fy_1*tdepk(km+1,im+1,jm)+ &
       fx_1*fy*tdepk(km+1,im,jm+1)+fx*fy*tdepk(km+1,im+1,jm+1)


  if  ( ((a_xy - a_xy_p) < 200.).or.(H(ic+1,jc+1)<100.)) then
     tdep(kc,ic,jc) = fz_1*a_xy+fz*a_xy_p
     
  else if ((a_xy_p<time_max_accu +1000.) &
       .or.(fz<1.E-7).or.(fz_1<1.E-7) ) then   ! added 25/03/04: don t bother for < 1/50mm
     
     a_tmp = fz_1*a_xy+fz*a_xy_p
     
     if  ((km==1).and.(fz<0.5)) then
        tdep(kc,ic,jc) = a_tmp
        
     else  ! cubic interp from the closest neighbour
        if (fz>0.5) then  ! furthest to km, use km+2 unless bottom
           a_xy_m= a_xy
           a_xy = a_xy_p
           if (km+1==nz) then 
              a_xy_p = a_xy + 1.3*(a_xy - a_xy_m) ! like articially old 22nd layer
           else
              a_xy_p = fx_1*fy_1*tdepk(km+2,im,jm)+fx*fy_1*tdepk(km+2,im+1,jm)+ &
                   fx_1*fy*tdepk(km+2,im,jm+1)+fx*fy*tdepk(km+2,im+1,jm+1)
           end if
           fzz= - fz_1  ! because km+1 closer than km and E_D<E(km+1)
        else       
           a_xy_m = fx_1*fy_1*tdepk(km-1,im,jm)+fx*fy_1*tdepk(km-1,im+1,jm)+ &
                fx_1*fy*tdepk(km-1,im,jm+1)+fx*fy*tdepk(km-1,im+1,jm+1)
           fzz = fz
        end if
        fprime = - abs(a_xy_p - a_xy_m)/2.  ! /E(km+1)-E(km) included in fz
        fsec = - abs(a_xy_p - 2.*a_xy + a_xy_m)   ! idem
        ! there is a risk of overshooting, limit change to 1.3x linear change
        al  = abs(tdepk(kc,ic,jc) - a_tmp)
        ac  = abs(tdepk(kc,ic,jc) - (a_xy + fzz*fprime + fzz**2*fsec/2.))
        if (ac<(1.3*al)) then   
           tdep(kc,ic,jc) = a_xy + fzz*fprime + fzz**2*fsec/2.
        else
           tdep(kc,ic,jc) = a_tmp
        end if
     end if
     
     
  else
     
     index1  = floor(-a_xy/100.)           ! gives the youngest: ind1+1 < tdepk < ind1
     fra1 = real((-a_xy/100.) - index1)
     accucumul1 = real(  -fra1*(accucumul(index1)-accucumul(index1+1))  +  accucumul(index1) )
     index2  = floor(-a_xy_p/100.)
     fra2 = real((-a_xy_p/100.) - index2)
     accucumul2 = real( -fra2 * (accucumul(index2)-accucumul(index2+1)) + accucumul(index2) )
     
     !print*, 'ic,jc,kc :',ic,jc,kc
     !print*, 'a_xy et a_xy_p', a_xy, a_xy_p
     !print*, 'index1 et 2',index1,index2
     !print*, 'fra1 et 2', fra1,fra2
     !print*, 'accucumul 1 et 2 :', accucumul1, accucumul2,accucumul2-accucumul1  !aurel neem
     thinning =    ( real( -(E(km+1)-E(km)) / (accucumul2 - accucumul1)) )      ! checked
     !print*, 'thinning ', thinning  !aurel neem

     if (fz < 0.5) then
        !======= closer to km, probably upward motion, rather seldom NL thinks
        if (km>1) then
           a_xy_m = fx_1*fy_1*tdepk(km-1,im,jm)+fx*fy_1*tdepk(km-1,im+1,jm)+ &
                fx_1*fy*tdepk(km-1,im,jm+1)+fx*fy*tdepk(km-1,im+1,jm+1)
           ! should be a_xy_m - a_xy > 0
           if ((a_xy_m>time_max_accu+ 500.).and.(   (a_xy_m - a_xy) > 200.)) then
              index0  = floor(-a_xy_m/100.)
              fra0 = real((-a_xy_m/100.) - index0)
              accucumul0 = real( -fra0 * (accucumul(index0)-accucumul(index0+1)) &
                   + accucumul(index0) )
              ! acc0 > acc1

              !print*, 'ic,jc,kc :',ic,jc,kc
              !print*, 'accucumul 0 et 1 : ', accucumul0, accucumul1,accucumul1-accucumul0 !aurel neem
              thinning2 =    ( real( -(E(km)-E(km-1)) / (accucumul1 - accucumul0)) )            ! ok
              !print*, 'thinning2 ', thinning2  !aurel neem

              thinning =  max(1.e-8, (0.5-fz/2.0)*thinning2 + (0.5+fz/2.0)*thinning)
              
              ! fz/c2 says my theory, though for Greenland I find that /c3 or /c4 may be better...
           end if
        end if
 
        f_k = real( accucumul1 - (E_k - E(km) ) / thinning )           ! ok: f_k < acc1
        ! downward search
        index0=index1-1   ! no inversion in stratigraphy
        L4: do
           index0 = index0 + 1
           if ((index0>=NFT-1).or.(accucumul(index0)  < f_k )) exit L4  
           ! index < ind(t) < index-1,  t(i+1) < t < t(i)
        end do L4
        index0 = index0-1
        
        !=========================================
     else                               !  fz> 0.5, probably downwards
        
        if (km<nz-1) then
           a_xy_m = fx_1*fy_1*tdepk(km+2,im,jm)+fx*fy_1*tdepk(km+2,im+1,jm)+ &
                fx_1*fy*tdepk(km+2,im,jm+1)+fx*fy*tdepk(km+2,im+1,jm+1)
           if ((a_xy_m>time_max_accu+ 500.).and.(   (a_xy_m - a_xy_p) < -200.)) then
              index0  = floor(-a_xy_m/100.)
              fra0 = real((-a_xy_m/100.) - index0)
              accucumul0 = real( -fra0 * (accucumul(index0)-accucumul(index0+1)) &
                   + accucumul(index0) )
             ! print*, 'ic,jc,kc :',ic,jc,kc
             ! print*, 'accucumul 0 et 2 : ', accucumul0, accucumul2,accucumul0-accucumul2 !aurel neem
              thinning2 =     ( real( -(E(km+2)-E(km+1)) / (accucumul0 - accucumul2)) )      
             ! print*, 'thinning2 ', thinning2  !aurel neem
              ! E decreases upwards!
              thinning =  max(1.e-8, (0.5-fz_1/2.0)*thinning2 + (0.5+fz_1/2.0)*thinning)
           end if
        end if
        
        f_k = real( accucumul2 + (E(km+1) - E_k) / thinning )
        ! upward search
        index0=index2+2   ! no inversion in stratigraphy
        L5: do
           index0 = index0 - 1
           if ((index0<=1).or.(accucumul(index0)> f_k )) exit L5
           ! index < ind(t) < index+1,  t(i+1) < t < t(i)
        end do L5
     end if
     ! now we know that we have index-1 < tdep(k) < index   from accucumul condition
     
     if ((index0==1).or.(index0==NFT-1)) then
        fzz=2.
     else  
        fzz = real( (accucumul(index0) - f_k)/ (accucumul(index0)-accucumul(index0+1)) )
     end if

     if ((fzz>1.0).or.(fzz<0.0)) then
        tdep(kc,ic,jc) = fz_1*a_xy + fz*a_xy_p
     else
        tdep(kc,ic,jc) = ( fzz + real(index0) ) * real(-100)
     end if
     
  end if

!=============== end tdep-interp ====================

!       f_zt = (tdep(kc,ic,jc) - a_xy) / (a_xy_p - a_xy)
! NL 21/11/02: use only fz, not f_zt, i don't think f_zt is right

      a_xy   = fx_1*fy_1*xdepk(km,im,jm)+fx*fy_1*xdepk(km,im+1,jm)+ &
               fx_1*fy*xdepk(km,im,jm+1)+fx*fy*xdepk(km,im+1,jm+1)
      a_xy_p = fx_1*fy_1*xdepk(km+1,im,jm)+fx*fy_1*xdepk(km+1,im+1,jm)+ &
               fx_1*fy*xdepk(km+1,im,jm+1)+fx*fy*xdepk(km+1,im+1,jm+1)

      xdep(kc,ic,jc) = fz_1*a_xy+fz*a_xy_p

!========= YDEP ===================

      a_xy   = fx_1*fy_1*ydepk(km,im,jm)+fx*fy_1*ydepk(km,im+1,jm)+ &
               fx_1*fy*ydepk(km,im,jm+1)+fx*fy*ydepk(km,im+1,jm+1)
      a_xy_p = fx_1*fy_1*ydepk(km+1,im,jm)+fx*fy_1*ydepk(km+1,im+1,jm)+ &
               fx_1*fy*ydepk(km+1,im,jm+1)+fx*fy*ydepk(km+1,im+1,jm+1)

      ydep(kc,ic,jc) = fz_1*a_xy+fz*a_xy_p


end subroutine interpolate
!==========================================================================