module gms_math_upwind3
  use datatype
  use mem_manager
  use gms_math_function
  implicit none
  interface uw3_x
    module procedure uw3_x_x
    module procedure uw3_x_xy
    module procedure uw3_x_xz
    module procedure uw3_x_xyz
  end interface

  interface uw3_y
    module procedure uw3_y_y
    module procedure uw3_y_xy
    module procedure uw3_y_yz
    module procedure uw3_y_xyz
  end interface

  interface uw3_z
    module procedure uw3_z_z
    module procedure uw3_z_xz
    module procedure uw3_z_yz
    module procedure uw3_z_xyz
  end interface

contains
  function uw3_x_x(input1, input2) result(output)
    type(var_x), intent(in) :: input1, input2
    type(var_x) :: output
    integer :: new_id

    call get_new_id_x(new_id)

    output%id = new_id

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)

	work_x(lb_axis1+2 : ub_axis1-2,:,:,new_id) & 
        = ( - 		work_x(lb_axis1+4:ub_axis1,:,:,input1%id) 	* work_x(lb_axis1+4:ub_axis1,:,:,input2%id) &
			+ 8.0D0*work_x(lb_axis1+3:ub_axis1-1,:,:,input1%id) * work_x(lb_axis1+3:ub_axis1-1,:,:,input2%id) &
			- 8.0D0*work_x(lb_axis1+1:ub_axis1-3,:,:,input1%id) * work_x(lb_axis1+1:ub_axis1-3,:,:,input2%id) &
            +		work_x(lb_axis1:ub_axis1-4,:,:,input1%id) 	* work_x(lb_axis1:ub_axis1-4,:,:,input2%id) ) / (12.0D0*dx) &
		  + sign(1.0D0, work_x(lb_axis1+2:ub_axis1-2,:,:,input1%id))*(dx**3)/12.0D0 &
		  * ( 	  work_x(lb_axis1+4:ub_axis1,:,:,input1%id)		* work_x(lb_axis1+4:ub_axis1,:,:,input2%id) &
		  - 4.0D0*work_x(lb_axis1+3:ub_axis1-1,:,:,input1%id)	* work_x(lb_axis1+3:ub_axis1-1,:,:,input2%id) &
		  + 6.0D0*work_x(lb_axis1+2:ub_axis1-2,:,:,input1%id)	* work_x(lb_axis1+2:ub_axis1-2,:,:,input2%id) &
		  - 4.0D0*work_x(lb_axis1+1:ub_axis1-3,:,:,input1%id)	* work_x(lb_axis1+1:ub_axis1-3,:,:,input2%id) &
          + 	  work_x(lb_axis1:ub_axis1-4,:,:,input1%id)		* work_x(lb_axis1:ub_axis1-4,:,:,input2%id) )/(dx**4)
    
  end function uw3_x_x

  function uw3_x_xy(input1, input2) result(output)
    type(var_xy), intent(in) :: input1, input2
    type(var_xy) :: output
    integer :: new_id

    call get_new_id_xy(new_id)

    output%id = new_id

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)

	work_xy(lb_axis1+2 : ub_axis1-2,:,:,new_id) & 
        = ( - 		work_xy(lb_axis1+4:ub_axis1,:,:,input1%id)	* work_xy(lb_axis1+4:ub_axis1,:,:,input2%id) &
			+ 8.0D0*work_xy(lb_axis1+3:ub_axis1-1,:,:,input1%id)* work_xy(lb_axis1+3:ub_axis1-1,:,:,input2%id) &
			- 8.0D0*work_xy(lb_axis1+1:ub_axis1-3,:,:,input1%id)* work_xy(lb_axis1+1:ub_axis1-3,:,:,input2%id) &
            +		work_xy(lb_axis1:ub_axis1-4,:,:,input1%id) 	* work_xy(lb_axis1:ub_axis1-4,:,:,input2%id) ) / (12.0D0*dx) &
		  + sign(1.0D0, work_xy(lb_axis1+2:ub_axis1-2,:,:,input1%id))*(dx**3)/12.0D0 &
		  * ( 	  work_xy(lb_axis1+4:ub_axis1,:,:,input1%id)	* work_xy(lb_axis1+4:ub_axis1,:,:,input2%id) &
		  - 4.0D0*work_xy(lb_axis1+3:ub_axis1-1,:,:,input1%id)	* work_xy(lb_axis1+3:ub_axis1-1,:,:,input2%id) &
		  + 6.0D0*work_xy(lb_axis1+2:ub_axis1-2,:,:,input1%id)	* work_xy(lb_axis1+2:ub_axis1-2,:,:,input2%id) &
		  - 4.0D0*work_xy(lb_axis1+1:ub_axis1-3,:,:,input1%id)	* work_xy(lb_axis1+1:ub_axis1-3,:,:,input2%id) &
          + 	  work_xy(lb_axis1:ub_axis1-4,:,:,input1%id)	* work_xy(lb_axis1:ub_axis1-4,:,:,input2%id) )/(dx**4)

  end function uw3_x_xy

  function uw3_x_xz(input1, input2) result(output)
    type(var_xz), intent(in) :: input1, input2
    type(var_xz) :: output
    integer :: new_id

    call get_new_id_xz(new_id)

    output%id = new_id

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)

	work_xz(lb_axis1+2 : ub_axis1-2,:,:,new_id) & 
        = ( - 		work_xz(lb_axis1+4:ub_axis1,:,:,input1%id)	* work_xz(lb_axis1+4:ub_axis1,:,:,input2%id) &
			+ 8.0D0*work_xz(lb_axis1+3:ub_axis1-1,:,:,input1%id)* work_xz(lb_axis1+3:ub_axis1-1,:,:,input2%id) &
			- 8.0D0*work_xz(lb_axis1+1:ub_axis1-3,:,:,input1%id)* work_xz(lb_axis1+1:ub_axis1-3,:,:,input2%id) &
            +		work_xz(lb_axis1:ub_axis1-4,:,:,input1%id) 	* work_xz(lb_axis1:ub_axis1-4,:,:,input2%id) ) / (12.0D0*dx) &
		  + sign(1.0D0, work_xz(lb_axis1+2:ub_axis1-2,:,:,input1%id))*(dx**3)/12.0D0 &
		  * ( 	  work_xz(lb_axis1+4:ub_axis1,:,:,input1%id)	* work_xz(lb_axis1+4:ub_axis1,:,:,input2%id) &
		  - 4.0D0*work_xz(lb_axis1+3:ub_axis1-1,:,:,input1%id)	* work_xz(lb_axis1+3:ub_axis1-1,:,:,input2%id) &
		  + 6.0D0*work_xz(lb_axis1+2:ub_axis1-2,:,:,input1%id)	* work_xz(lb_axis1+2:ub_axis1-2,:,:,input2%id) &
		  - 4.0D0*work_xz(lb_axis1+1:ub_axis1-3,:,:,input1%id)	* work_xz(lb_axis1+1:ub_axis1-3,:,:,input2%id) &
          + 	  work_xz(lb_axis1:ub_axis1-4,:,:,input1%id)	* work_xz(lb_axis1:ub_axis1-4,:,:,input2%id) )/(dx**4)

  end function uw3_x_xz

  function uw3_x_xyz(input1, input2) result(output)
    type(var_xyz), intent(in) :: input1, input2
    type(var_xyz) :: output
    integer :: new_id

    call get_new_id_xyz(new_id)

    output%id = new_id

	output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)

	work_xyz(lb_axis1+2 : ub_axis1-2,:,:,new_id) & 
        = ( - 		work_xyz(lb_axis1+4:ub_axis1,:,:,input1%id)		* work_xyz(lb_axis1+4:ub_axis1,:,:,input2%id) &
			+ 8.0D0*work_xyz(lb_axis1+3:ub_axis1-1,:,:,input1%id)	* work_xyz(lb_axis1+3:ub_axis1-1,:,:,input2%id) &
			- 8.0D0*work_xyz(lb_axis1+1:ub_axis1-3,:,:,input1%id)	* work_xyz(lb_axis1+1:ub_axis1-3,:,:,input2%id) &
            +		work_xyz(lb_axis1:ub_axis1-4,:,:,input1%id) 	* work_xyz(lb_axis1:ub_axis1-4,:,:,input2%id) ) / (12.0D0*dx) &
		  + sign(1.0D0, work_xyz(lb_axis1+2:ub_axis1-2,:,:,input1%id))*(dx**3)/12.0D0 &
		  * ( 	  work_xyz(lb_axis1+4:ub_axis1,:,:,input1%id)		* work_xyz(lb_axis1+4:ub_axis1,:,:,input2%id) &
		  - 4.0D0*work_xyz(lb_axis1+3:ub_axis1-1,:,:,input1%id)		* work_xyz(lb_axis1+3:ub_axis1-1,:,:,input2%id) &
		  + 6.0D0*work_xyz(lb_axis1+2:ub_axis1-2,:,:,input1%id)		* work_xyz(lb_axis1+2:ub_axis1-2,:,:,input2%id) &
		  - 4.0D0*work_xyz(lb_axis1+1:ub_axis1-3,:,:,input1%id)		* work_xyz(lb_axis1+1:ub_axis1-3,:,:,input2%id) &
          + 	  work_xyz(lb_axis1:ub_axis1-4,:,:,input1%id)		* work_xyz(lb_axis1:ub_axis1-4,:,:,input2%id) )/(dx**4)

  end function uw3_x_xyz

  function uw3_y_y(input1, input2) result(output)
    type(var_y), intent(in) :: input1, input2
    type(var_y) :: output
    integer :: new_id

    call get_new_id_y(new_id)

    output%id = new_id

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)

	work_y(:,lb_axis2+2 : ub_axis2-2,:,new_id) & 
        = ( - 		work_y(:,lb_axis2+4:ub_axis2,:,input1%id) 	* work_y(:,lb_axis2+4:ub_axis2,:,input2%id) &
			+ 8.0D0*work_y(:,lb_axis2+3:ub_axis2-1,:,input1%id) * work_y(:,lb_axis2+3:ub_axis2-1,:,input2%id) &
			- 8.0D0*work_y(:,lb_axis2+1:ub_axis2-3,:,input1%id) * work_y(:,lb_axis2+1:ub_axis2-3,:,input2%id) &
            +		work_y(:,lb_axis2:ub_axis2-4,:,input1%id) 	* work_y(:,lb_axis2:ub_axis2-4,:,input2%id) ) / (12.0D0*dy) &
		  + sign(1.0D0, work_y(:,lb_axis2+2:ub_axis2-2,:,input1%id))*(dy**3)/12.0D0 &
		  * ( 	  work_y(:,lb_axis2+4:ub_axis2,:,input1%id)		* work_y(:,lb_axis2+4:ub_axis2,:,input2%id) &
		  - 4.0D0*work_y(:,lb_axis2+3:ub_axis2-1,:,input1%id)	* work_y(:,lb_axis2+3:ub_axis2-1,:,input2%id) &
		  + 6.0D0*work_y(:,lb_axis2+2:ub_axis2-2,:,input1%id)	* work_y(:,lb_axis2+2:ub_axis2-2,:,input2%id) &
		  - 4.0D0*work_y(:,lb_axis2+1:ub_axis2-3,:,input1%id)	* work_y(:,lb_axis2+1:ub_axis2-3,:,input2%id) &
          + 	  work_y(:,lb_axis2:ub_axis2-4,:,input1%id)		* work_y(:,lb_axis2:ub_axis2-4,:,input2%id) )/(dy**4)

  end function uw3_y_y

  function uw3_y_xy(input1, input2) result(output)
    type(var_xy), intent(in) :: input1, input2
    type(var_xy) :: output
    integer :: new_id

    call get_new_id_xy(new_id)

    output%id = new_id

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)

	work_xy(:,lb_axis2+2 : ub_axis2-2,:,new_id) & 
        = ( - 		work_xy(:,lb_axis2+4:ub_axis2,:,input1%id)	* work_xy(:,lb_axis2+4:ub_axis2,:,input2%id) &
			+ 8.0D0*work_xy(:,lb_axis2+3:ub_axis2-1,:,input1%id)* work_xy(:,lb_axis2+3:ub_axis2-1,:,input2%id) &
			- 8.0D0*work_xy(:,lb_axis2+1:ub_axis2-3,:,input1%id)* work_xy(:,lb_axis2+1:ub_axis2-3,:,input2%id) &
            +		work_xy(:,lb_axis2:ub_axis2-4,:,input1%id)	* work_xy(:,lb_axis2:ub_axis2-4,:,input2%id) ) / (12.0D0*dy) &
		  + sign(1.0D0, work_xy(:,lb_axis2+2:ub_axis2-2,:,input1%id))*(dy**3)/12.0D0 &
		  * ( 	  work_xy(:,lb_axis2+4:ub_axis2,:,input1%id)	* work_xy(:,lb_axis2+4:ub_axis2,:,input2%id) &
		  - 4.0D0*work_xy(:,lb_axis2+3:ub_axis2-1,:,input1%id)	* work_xy(:,lb_axis2+3:ub_axis2-1,:,input2%id) &
		  + 6.0D0*work_xy(:,lb_axis2+2:ub_axis2-2,:,input1%id)	* work_xy(:,lb_axis2+2:ub_axis2-2,:,input2%id) &
		  - 4.0D0*work_xy(:,lb_axis2+1:ub_axis2-3,:,input1%id)	* work_xy(:,lb_axis2+1:ub_axis2-3,:,input2%id) &
          + 	  work_xy(:,lb_axis2:ub_axis2-4,:,input1%id)	* work_xy(:,lb_axis2:ub_axis2-4,:,input2%id) )/(dy**4)

  end function uw3_y_xy

  function uw3_y_yz(input1, input2) result(output)
    type(var_yz), intent(in) :: input1, input2
    type(var_yz) :: output
    integer :: new_id

    call get_new_id_yz(new_id)

    output%id = new_id

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)

	work_yz(:,lb_axis2+2 : ub_axis2-2,:,new_id) & 
        = ( - 		work_yz(:,lb_axis2+4:ub_axis2,:,input1%id)	* work_yz(:,lb_axis2+4:ub_axis2,:,input2%id) &
			+ 8.0D0*work_yz(:,lb_axis2+3:ub_axis2-1,:,input1%id)* work_yz(:,lb_axis2+3:ub_axis2-1,:,input2%id) &
			- 8.0D0*work_yz(:,lb_axis2+1:ub_axis2-3,:,input1%id)* work_yz(:,lb_axis2+1:ub_axis2-3,:,input2%id) &
            +		work_yz(:,lb_axis2:ub_axis2-4,:,input1%id) 	* work_yz(:,lb_axis2:ub_axis2-4,:,input2%id) ) / (12.0D0*dy) &
		  + sign(1.0D0, work_yz(:,lb_axis2+2:ub_axis2-2,:,input1%id))*(dy**3)/12.0D0 &
		  * ( 	  work_yz(:,lb_axis2+4:ub_axis2,:,input1%id)	* work_yz(:,lb_axis2+4:ub_axis2,:,input2%id) &
		  - 4.0D0*work_yz(:,lb_axis2+3:ub_axis2-1,:,input1%id)	* work_yz(:,lb_axis2+3:ub_axis2-1,:,input2%id) &
		  + 6.0D0*work_yz(:,lb_axis2+2:ub_axis2-2,:,input1%id)	* work_yz(:,lb_axis2+2:ub_axis2-2,:,input2%id) &
		  - 4.0D0*work_yz(:,lb_axis2+1:ub_axis2-3,:,input1%id)	* work_yz(:,lb_axis2+1:ub_axis2-3,:,input2%id) &
          + 	  work_yz(:,lb_axis2:ub_axis2-4,:,input1%id)	* work_yz(:,lb_axis2:ub_axis2-4,:,input2%id) )/(dy**4)

  end function uw3_y_yz

  function uw3_y_xyz(input1, input2) result(output)
    type(var_xyz), intent(in) :: input1, input2
    type(var_xyz) :: output
    integer :: new_id

    call get_new_id_xyz(new_id)

    output%id = new_id

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)

	work_xyz(:,lb_axis2+2 : ub_axis2-2,:,new_id) & 
        = ( - 		work_xyz(:,lb_axis2+4:ub_axis2,:,input1%id)		* work_xyz(:,lb_axis2+4:ub_axis2,:,input2%id) &
			+ 8.0D0*work_xyz(:,lb_axis2+3:ub_axis2-1,:,input1%id)	* work_xyz(:,lb_axis2+3:ub_axis2-1,:,input2%id) &
			- 8.0D0*work_xyz(:,lb_axis2+1:ub_axis2-3,:,input1%id)	* work_xyz(:,lb_axis2+1:ub_axis2-3,:,input2%id) &
            +		work_xyz(:,lb_axis2:ub_axis2-4,:,input1%id) 	* work_xyz(:,lb_axis2:ub_axis2-4,:,input2%id) ) / (12.0D0*dy) &
		  + sign(1.0D0, work_xyz(:,lb_axis2+2:ub_axis2-2,:,input1%id))*(dy**3)/12.0D0 &
		  * ( 	  work_xyz(:,lb_axis2+4:ub_axis2,:,input1%id)		* work_xyz(:,lb_axis2+4:ub_axis2,:,input2%id) &
		  - 4.0D0*work_xyz(:,lb_axis2+3:ub_axis2-1,:,input1%id)		* work_xyz(:,lb_axis2+3:ub_axis2-1,:,input2%id) &
		  + 6.0D0*work_xyz(:,lb_axis2+2:ub_axis2-2,:,input1%id)		* work_xyz(:,lb_axis2+2:ub_axis2-2,:,input2%id) &
		  - 4.0D0*work_xyz(:,lb_axis2+1:ub_axis2-3,:,input1%id)		* work_xyz(:,lb_axis2+1:ub_axis2-3,:,input2%id) &
          + 	  work_xyz(:,lb_axis2:ub_axis2-4,:,input1%id)		* work_xyz(:,lb_axis2:ub_axis2-4,:,input2%id) )/(dy**4)

  end function uw3_y_xyz

  function uw3_z_z(input1, input2) result(output)
    type(var_z), intent(in) :: input1, input2
    type(var_z) :: output
    integer :: new_id

    call get_new_id_z(new_id)

    output%id = new_id

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)

	work_z(:,:,lb_axis3+2 : ub_axis3-2,new_id) & 
        = ( - 		work_z(:,:,lb_axis3+4:ub_axis3,input1%id) 	* work_z(:,:,lb_axis3+4:ub_axis3,input2%id) &
			+ 8.0D0*work_z(:,:,lb_axis3+3:ub_axis3-1,input1%id) * work_z(:,:,lb_axis3+3:ub_axis3-1,input2%id) &
			- 8.0D0*work_z(:,:,lb_axis3+1:ub_axis3-3,input1%id) * work_z(:,:,lb_axis3+1:ub_axis3-3,input2%id) &
            +		work_z(:,:,lb_axis3:ub_axis3-4,input1%id) 	* work_z(:,:,lb_axis3:ub_axis3-4,input2%id) ) / (12.0D0*dz) &
		  + sign(1.0D0, work_z(:,:,lb_axis3+2:ub_axis3-2,input1%id))*(dz**3)/12.0D0 &
		  * ( 	  work_z(:,:,lb_axis3+4:ub_axis3,input1%id)		* work_z(:,:,lb_axis3+4:ub_axis3,input2%id) &
		  - 4.0D0*work_z(:,:,lb_axis3+3:ub_axis3-1,input1%id)	* work_z(:,:,lb_axis3+3:ub_axis3-1,input2%id) &
		  + 6.0D0*work_z(:,:,lb_axis3+2:ub_axis3-2,input1%id)	* work_z(:,:,lb_axis3+2:ub_axis3-2,input2%id) &
		  - 4.0D0*work_z(:,:,lb_axis3+1:ub_axis3-3,input1%id)	* work_z(:,:,lb_axis3+1:ub_axis3-3,input2%id) &
          + 	  work_z(:,:,lb_axis3:ub_axis3-4,input1%id)		* work_z(:,:,lb_axis3:ub_axis3-4,input2%id) )/(dz**4)

  end function uw3_z_z

  function uw3_z_xz(input1, input2) result(output)
    type(var_xz), intent(in) :: input1, input2
    type(var_xz) :: output
    integer :: new_id

    call get_new_id_xz(new_id)

    output%id = new_id

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)

	work_xz(:,:,lb_axis3+2 : ub_axis3-2,new_id) & 
        = ( - 		work_xz(:,:,lb_axis3+4:ub_axis3,input1%id) 	* work_xz(:,:,lb_axis3+4:ub_axis3,input2%id) &
			+ 8.0D0*work_xz(:,:,lb_axis3+3:ub_axis3-1,input1%id)* work_xz(:,:,lb_axis3+3:ub_axis3-1,input2%id) &
			- 8.0D0*work_xz(:,:,lb_axis3+1:ub_axis3-3,input1%id)* work_xz(:,:,lb_axis3+1:ub_axis3-3,input2%id) &
            +		work_xz(:,:,lb_axis3:ub_axis3-4,input1%id) 	* work_xz(:,:,lb_axis3:ub_axis3-4,input2%id) ) / (12.0D0*dz) &
		  + sign(1.0D0, work_xz(:,:,lb_axis3+2:ub_axis3-2,input1%id))*(dz**3)/12.0D0 &
		  * ( 	  work_xz(:,:,lb_axis3+4:ub_axis3,input1%id)	* work_xz(:,:,lb_axis3+4:ub_axis3,input2%id) &
		  - 4.0D0*work_xz(:,:,lb_axis3+3:ub_axis3-1,input1%id)	* work_xz(:,:,lb_axis3+3:ub_axis3-1,input2%id) &
		  + 6.0D0*work_xz(:,:,lb_axis3+2:ub_axis3-2,input1%id)	* work_xz(:,:,lb_axis3+2:ub_axis3-2,input2%id) &
		  - 4.0D0*work_xz(:,:,lb_axis3+1:ub_axis3-3,input1%id)	* work_xz(:,:,lb_axis3+1:ub_axis3-3,input2%id) &
          + 	  work_xz(:,:,lb_axis3:ub_axis3-4,input1%id)	* work_xz(:,:,lb_axis3:ub_axis3-4,input2%id) )/(dz**4)

  end function uw3_z_xz

  function uw3_z_yz(input1, input2) result(output)
    type(var_yz), intent(in) :: input1, input2
    type(var_yz) :: output
    integer :: new_id

    call get_new_id_yz(new_id)

    output%id = new_id

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)

	work_yz(:,:,lb_axis3+2 : ub_axis3-2,new_id) & 
        = ( - 		work_yz(:,:,lb_axis3+4:ub_axis3,input1%id) 	* work_yz(:,:,lb_axis3+4:ub_axis3,input2%id) &
			+ 8.0D0*work_yz(:,:,lb_axis3+3:ub_axis3-1,input1%id)* work_yz(:,:,lb_axis3+3:ub_axis3-1,input2%id) &
			- 8.0D0*work_yz(:,:,lb_axis3+1:ub_axis3-3,input1%id)* work_yz(:,:,lb_axis3+1:ub_axis3-3,input2%id) &
            +		work_yz(:,:,lb_axis3:ub_axis3-4,input1%id) 	* work_yz(:,:,lb_axis3:ub_axis3-4,input2%id) ) / (12.0D0*dz) &
		  + sign(1.0D0, work_yz(:,:,lb_axis3+2:ub_axis3-2,input1%id))*(dz**3)/12.0D0 &
		  * ( 	  work_yz(:,:,lb_axis3+4:ub_axis3,input1%id)	* work_yz(:,:,lb_axis3+4:ub_axis3,input2%id) &
		  - 4.0D0*work_yz(:,:,lb_axis3+3:ub_axis3-1,input1%id)	* work_yz(:,:,lb_axis3+3:ub_axis3-1,input2%id) &
		  + 6.0D0*work_yz(:,:,lb_axis3+2:ub_axis3-2,input1%id)	* work_yz(:,:,lb_axis3+2:ub_axis3-2,input2%id) &
		  - 4.0D0*work_yz(:,:,lb_axis3+1:ub_axis3-3,input1%id)	* work_yz(:,:,lb_axis3+1:ub_axis3-3,input2%id) &
          + 	  work_yz(:,:,lb_axis3:ub_axis3-4,input1%id)	* work_yz(:,:,lb_axis3:ub_axis3-4,input2%id) )/(dz**4)

  end function uw3_z_yz

  function uw3_z_xyz(input1, input2) result(output)
    type(var_xyz), intent(in) :: input1, input2
    type(var_xyz) :: output
    integer :: new_id

    call get_new_id_xyz(new_id)

    output%id = new_id

    output%grid(1) = input2%grid(1)
    output%grid(2) = input2%grid(2)
    output%grid(3) = input2%grid(3)

	work_xyz(:,:,lb_axis3+2 : ub_axis3-2,new_id) & 
        = ( - 		work_xyz(:,:,lb_axis3+4:ub_axis3,input1%id) 	* work_xyz(:,:,lb_axis3+4:ub_axis3,input2%id) &
			+ 8.0D0*work_xyz(:,:,lb_axis3+3:ub_axis3-1,input1%id)	* work_xyz(:,:,lb_axis3+3:ub_axis3-1,input2%id) &
			- 8.0D0*work_xyz(:,:,lb_axis3+1:ub_axis3-3,input1%id)	* work_xyz(:,:,lb_axis3+1:ub_axis3-3,input2%id) &
            +		work_xyz(:,:,lb_axis3:ub_axis3-4,input1%id) 	* work_xyz(:,:,lb_axis3:ub_axis3-4,input2%id) ) / (12.0D0*dz) &
		  + sign(1.0D0, work_xyz(:,:,lb_axis3+2:ub_axis3-2,input1%id))*(dz**3)/12.0D0 &
		  * ( 	  work_xyz(:,:,lb_axis3+4:ub_axis3,input1%id)		* work_xyz(:,:,lb_axis3+4:ub_axis3,input2%id) &
		  - 4.0D0*work_xyz(:,:,lb_axis3+3:ub_axis3-1,input1%id)		* work_xyz(:,:,lb_axis3+3:ub_axis3-1,input2%id) &
		  + 6.0D0*work_xyz(:,:,lb_axis3+2:ub_axis3-2,input1%id)		* work_xyz(:,:,lb_axis3+2:ub_axis3-2,input2%id) &
		  - 4.0D0*work_xyz(:,:,lb_axis3+1:ub_axis3-3,input1%id)		* work_xyz(:,:,lb_axis3+1:ub_axis3-3,input2%id) &
          + 	  work_xyz(:,:,lb_axis3:ub_axis3-4,input1%id)		* work_xyz(:,:,lb_axis3:ub_axis3-4,input2%id) )/(dz**4)

  end function uw3_z_xyz

end module gms_math_upwind3