! this is the w-node version
subroutine interpolag_Ssim(jt)
! This subroutine computes the values of F_LM and F_MM 
! at positions (x-u*dt) for use in the subroutine lagrangian
!use sgsmodule,only:u_lag,v_lag,w_lag
use types,only:rprec
use param,only:ld,nx,ny,nz,dx,dy,dz,dt,cs_count
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) :: lagran_dt, consta,test_cs_count
integer :: addx, addy, addz
integer :: jxaddx, jyaddy, jzaddz
real(kind=rprec), dimension(ld,ny,nz) :: u_temp, v_temp
!     creates dummy arrays FF_LM and FF_MM to use in the subroutine
!$omp parallel do default(shared) private(jx,jy,jz)	
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)
end do
end do
end do
!$omp end parallel 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) 

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

!$omp parallel do default(shared) private(jz)	
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
!$omp end parallel do
u_lag (:,:,1) = u_temp(:,:,1)
v_lag (:,:,1) = v_temp(:,:,1)
w_lag = w_lag/real(cs_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(cs_count,kind=rprec)/dx    !! is u already normalized
yp = -v_lag*dt*real(cs_count,kind=rprec)/dy   
zp = -w_lag*dt*real(cs_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 cell beneath it
do jx= 1,nx
do jy=1,ny
   zp(jx,jy,2) = zp(jx,jy,2) + dmin1(zp(jx,jy,2),0._rprec)
   zp(jx,jy,1) = zp(jy,jy,1) + dmax1(zp(jx,jy,1),0._rprec)
end do
end do

test_cs_count=maxval(abs(xp))
print*, 'Lagrangian CFL condition'
print*, test_cs_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(sign(1._rprec,xp(jx,jy,jz)))
     addy = int(sign(1._rprec,yp(jx,jy,jz)))
     addz = int(sign(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,1000) == 0) .AND. (jt>55000)) then    
!TS        write(80,5180)jt*dt,float(jx),float(jy),float(jz),real(frac_x),real*frac_y),real(frac_z)
        write(80,*)jt*dt,jx,jy,jz,real(frac_x),real(frac_y),real(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 interpoated 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 interpoated 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)
  end do
 end do
end do

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