! this is the w-node version
      subroutine interpolag_Sdep(jt)
! This subroutine computes the values of F_LM and F_MM 
! at positions (x-u*dt) for use in the subroutine lagrangian
use types,only:rprec
use param,only:ld,nx,ny,nz,dx,dy,dz,dt
use sgsmodule
implicit none
integer :: jx, jy, jz ,jt
integer :: jjx, jjy, jjz
real(kind=rprec), dimension(ld,ny,nz) :: xp, yp, zp
real(kind=rprec) :: frac_x,frac_y,frac_z
real(kind=rprec) :: comp_x,comp_y,comp_z
real(kind=rprec), dimension(nx+2,ny+2,nz+2) :: FF_LM,FF_MM
real(kind=rprec), dimension(nx+2,ny+2,nz+2) :: FF_QN,FF_NN

integer :: addx, addy, addz
integer :: jxaddx, jyaddy, jzaddz
real(kind=rprec), dimension(ld,ny,nz) :: u_temp, v_temp
real(kind=rprec)::test_dyn_count
      
!     creates dummy arrays FF_LM and FF_MM to use in the subroutine
 
do jx=1,nx
do jy=1,ny
do jz=1,nz
   FF_LM(jx+1,jy+1,jz+1) = F_LM(jx,jy,jz)
   FF_MM(jx+1,jy+1,jz+1) = F_MM(jx,jy,jz)
   FF_QN(jx+1,jy+1,jz+1) = F_QN(jx,jy,jz)
   FF_NN(jx+1,jy+1,jz+1) = F_NN(jx,jy,jz)
end do
end do
end do

!This is a bit like witch craft but the following lines do take care
!of all the edges including the corners

FF_LM(1,:,:) = FF_LM(nx+1,:,:)
FF_LM(nx+2,:,:) = FF_LM(2,:,:)

FF_LM(:,1,:) = FF_LM(:,ny+1,:) 
FF_LM(:,ny+2,:) = FF_LM(:,2,:) 

FF_LM(:,:,1) = FF_LM(:,:,2) 
FF_LM(:,:,nz+2) = FF_LM(:,:,nz+1) 

FF_MM(1,:,:) = FF_MM(nx+1,:,:)
FF_MM(nx+2,:,:) = FF_MM(2,:,:)

FF_MM(:,1,:) = FF_MM(:,ny+1,:) 
FF_MM(:,ny+2,:) = FF_MM(:,2,:) 

FF_MM(:,:,1) = FF_MM(:,:,2) 
FF_MM(:,:,nz+2) = FF_MM(:,:,nz+1) 

FF_QN(1,:,:) = FF_QN(nx+1,:,:)
FF_QN(nx+2,:,:) = FF_QN(2,:,:)

FF_QN(:,1,:) = FF_QN(:,ny+1,:) 
FF_QN(:,ny+2,:) = FF_QN(:,2,:) 

FF_QN(:,:,1) = FF_QN(:,:,2) 
FF_QN(:,:,nz+2) = FF_QN(:,:,nz+1) 

FF_NN(1,:,:) = FF_NN(nx+1,:,:)
FF_NN(nx+2,:,:) = FF_NN(2,:,:)

FF_NN(:,1,:) = FF_NN(:,ny+1,:) 
FF_NN(:,ny+2,:) = FF_NN(:,2,:) 

FF_NN(:,:,1) = FF_NN(:,:,2) 
FF_NN(:,:,nz+2) = FF_NN(:,:,nz+1) 

! end of witch craft
!	puts u_lag and and v_lag on w nodes
u_temp = u_lag/real(dyn_count,kind=rprec)
v_temp = v_lag/real(dyn_count,kind=rprec)

do jz= 2,nz
   u_lag (:,:,jz) = 0.5_rprec*(u_temp(:,:,jz) + u_temp(:,:,jz-1))
   v_lag (:,:,jz) = 0.5_rprec*(v_temp(:,:,jz) + v_temp(:,:,jz-1))
end do
u_lag (:,:,1) = u_temp(:,:,1)
v_lag (:,:,1) = v_temp(:,:,1)
w_lag = w_lag/real(dyn_count,kind=rprec)
w_lag (:,:,1) = .25_rprec*w_lag (:,:,2)
!     computes the 3-D inverse displacement arrays that describe
!     the location where the point was at the previous step
xp = -u_lag*dt*real(dyn_count,kind=rprec)/dx    !! is u already normalized
yp = -v_lag*dt*real(dyn_count,kind=rprec)/dy   
zp = -w_lag*dt*real(dyn_count,kind=rprec)/dz

!	because first plane is on u,v,p nodes
!	this corrects for the fact that the first cell
!	in the z direction has height dz/2
!	it doubles the zp fraction if this fraction relates to the ecll beneath it
zp(:,:,2)=zp(:,:,2)+min(real(zp(jx,jy,2),kind=rprec),0._rprec)
zp(:,:,1)=zp(:,:,1)+max(real(zp(jx,jy,1),kind=rprec),0._rprec)

test_dyn_count=maxval(abs(xp))
!print*, 'test_dyn_count'
!print*, test_dyn_count
!      if ((test_dyn_count>0.7) .OR. (test_dyn_count<0.5)) then
!      jt_count = 1
!      dyn_count = int (0.6 / test_dyn_count * float(dyn_count))
!      dyn_count = max(1,dyn_count)
!      dyn_count = min(10,dyn_count)
!      print*, 'dyn_count corrected to', dyn_count
!      end if
do jx=1,nx
jjx = jx+1
do jy=1,ny
jjy = jy+1
do jz=1,nz
   jjz = jz+1
!     the are the values to add to the indices jx, jy, and jz
!     the are +1 or -1 depending on what cube should be used for interpolation
   addx = int(dsign(1._rprec,xp(jx,jy,jz)))
   addy = int(dsign(1._rprec,yp(jx,jy,jz)))
   addz = int(dsign(1._rprec,zp(jx,jy,jz)))
   
   jxaddx = jjx + addx 
   jyaddy = jjy + addy
   jzaddz = jjz + addz
!     computes the relative weights given to F_** in the cube depending on point location
   comp_x = abs(xp(jx,jy,jz))
   comp_y = abs(yp(jx,jy,jz))
   comp_z = abs(zp(jx,jy,jz)) 
     
   if ((mod(jt,500) == 0) .AND. (jt>55000)) then    
      write(80,5180) jt*dt,float(jx),float(jy),float(jz),frac_x,frac_y,frac_z
   end if
    
5180 format(120(E14.5))
  
   frac_x = 1._rprec - comp_x
   frac_y = 1._rprec - comp_y
   frac_z = 1._rprec - comp_z

!     computes interpolated F_LM

   F_LM(jx,jy,jz)=frac_x*frac_y*&
        (FF_LM(jjx,jjy,jjz)*frac_z+FF_LM(jjx,jjy,jzaddz)*comp_z)&
        + frac_x*comp_y*&
        (FF_LM(jjx,jyaddy,jjz)*frac_z+FF_LM(jjx,jyaddy,jzaddz)*comp_z)&
        + comp_x*frac_y*&
        (FF_LM(jxaddx,jjy,jjz)*frac_z+FF_LM(jxaddx,jjy,jzaddz)*comp_z)&
        + comp_x*comp_y*&
        (FF_LM(jxaddx,jyaddy,jjz)*frac_z&
        +FF_LM(jxaddx,jyaddy,jzaddz)*comp_z)

!     computes interpolated F_MM
   F_MM(jx,jy,jz)=frac_x*frac_y*&
        (FF_MM(jjx,jjy,jjz)*frac_z+FF_MM(jjx,jjy,jzaddz)*comp_z)&
        + frac_x*comp_y*&
        (FF_MM(jjx,jyaddy,jjz)*frac_z+FF_MM(jjx,jyaddy,jzaddz)*comp_z)&
        + comp_x*frac_y*&
        (FF_MM(jxaddx,jjy,jjz)*frac_z+FF_MM(jxaddx,jjy,jzaddz)*comp_z)&
        + comp_x*comp_y*&
        (FF_MM(jxaddx,jyaddy,jjz)*frac_z&
        +FF_MM(jxaddx,jyaddy,jzaddz)*comp_z)

!     computes interpolated F_QN
   F_QN(jx,jy,jz)=frac_x*frac_y*&
        (FF_QN(jjx,jjy,jjz)*frac_z+FF_QN(jjx,jjy,jzaddz)*comp_z)&
        + frac_x*comp_y*&
        (FF_QN(jjx,jyaddy,jjz)*frac_z+FF_QN(jjx,jyaddy,jzaddz)*comp_z)&
        + comp_x*frac_y*&
        (FF_QN(jxaddx,jjy,jjz)*frac_z+FF_QN(jxaddx,jjy,jzaddz)*comp_z)&
        + comp_x*comp_y*&
        (FF_QN(jxaddx,jyaddy,jjz)*frac_z&
        +FF_QN(jxaddx,jyaddy,jzaddz)*comp_z)

!     computes interpolated F_NN
   F_NN(jx,jy,jz)=frac_x*frac_y*&
        (FF_NN(jjx,jjy,jjz)*frac_z+FF_NN(jjx,jjy,jzaddz)*comp_z)&
        + frac_x*comp_y*&
        (FF_NN(jjx,jyaddy,jjz)*frac_z+FF_NN(jjx,jyaddy,jzaddz)*comp_z)&
        + comp_x*frac_y*&
        (FF_NN(jxaddx,jjy,jjz)*frac_z+FF_NN(jxaddx,jjy,jzaddz)*comp_z)&
        + comp_x*comp_y*&
        (FF_NN(jxaddx,jyaddy,jjz)*frac_z&
        +FF_NN(jxaddx,jyaddy,jzaddz)*comp_z)
end do
end do
end do

u_lag = 0._rprec
v_lag = 0._rprec
w_lag = 0._rprec
end subroutine interpolag_Sdep
