module gms_math_upwind1 use datatype use mem_manager use gms_math_function implicit none interface uw1_x module procedure uw1_x_x module procedure uw1_x_xy module procedure uw1_x_xz module procedure uw1_x_xyz end interface interface uw1_y module procedure uw1_y_y module procedure uw1_y_xy module procedure uw1_y_yz module procedure uw1_y_xyz end interface interface uw1_z module procedure uw1_z_z module procedure uw1_z_xz module procedure uw1_z_yz module procedure uw1_z_xyz end interface contains function uw1_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+1 : ub_axis1-1,:,:,new_id) & = ( work_x(lb_axis1+2:ub_axis1,:,:,input1%id) * work_x(lb_axis1+2:ub_axis1,:,:,input2%id) & - work_x(lb_axis1:ub_axis1-2,:,:,input1%id) * work_x(lb_axis1:ub_axis1-2,:,:,input2%id) ) / (2.0D0*dx) & - sign(1.0D0, work_x(lb_axis1+1:ub_axis1-1,:,:,input1%id))*dx/2.0D0 & * ( work_x(lb_axis1+2:ub_axis1,:,:,input1%id)* work_x(lb_axis1+2:ub_axis1,:,:,input2%id) & - 2.0D0*work_x(lb_axis1+1:ub_axis1-1,:,:,input1%id)* work_x(lb_axis1+1:ub_axis1-1,:,:,input2%id) & + work_x(lb_axis1:ub_axis1-2,:,:,input1%id)* work_x(lb_axis1:ub_axis1-2,:,:,input2%id) )/(dx**2) end function uw1_x_x function uw1_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+1 : ub_axis1-1,:,:,new_id) & = ( work_xy(lb_axis1+2:ub_axis1,:,:,input1%id) * work_xy(lb_axis1+2:ub_axis1,:,:,input2%id) & - work_xy(lb_axis1:ub_axis1-2,:,:,input1%id) * work_xy(lb_axis1:ub_axis1-2,:,:,input2%id) ) / (2.0D0*dx) & - sign(1.0D0, work_xy(lb_axis1+1:ub_axis1-1,:,:,input1%id))*dx/2.0D0 & * ( work_xy(lb_axis1+2:ub_axis1,:,:,input1%id)* work_xy(lb_axis1+2:ub_axis1,:,:,input2%id) & - 2.0D0*work_xy(lb_axis1+1:ub_axis1-1,:,:,input1%id)* work_xy(lb_axis1+1:ub_axis1-1,:,:,input2%id) & + work_xy(lb_axis1:ub_axis1-2,:,:,input1%id)* work_xy(lb_axis1:ub_axis1-2,:,:,input2%id) )/(dx**2) end function uw1_x_xy function uw1_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+1 : ub_axis1-1,:,:,new_id) & = ( work_xz(lb_axis1+2:ub_axis1,:,:,input1%id) * work_xz(lb_axis1+2:ub_axis1,:,:,input2%id) & - work_xz(lb_axis1:ub_axis1-2,:,:,input1%id) * work_xz(lb_axis1:ub_axis1-2,:,:,input2%id) ) / (2.0D0*dx) & - sign(1.0D0, work_xz(lb_axis1+1:ub_axis1-1,:,:,input1%id))*dx/2.0D0 & * ( work_xz(lb_axis1+2:ub_axis1,:,:,input1%id)* work_xz(lb_axis1+2:ub_axis1,:,:,input2%id) & - 2.0D0*work_xz(lb_axis1+1:ub_axis1-1,:,:,input1%id)* work_xz(lb_axis1+1:ub_axis1-1,:,:,input2%id) & + work_xz(lb_axis1:ub_axis1-2,:,:,input1%id)* work_xz(lb_axis1:ub_axis1-2,:,:,input2%id) )/(dx**2) end function uw1_x_xz function uw1_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+1 : ub_axis1-1,:,:,new_id) & = ( work_xyz(lb_axis1+2:ub_axis1,:,:,input1%id) * work_xyz(lb_axis1+2:ub_axis1,:,:,input2%id) & - work_xyz(lb_axis1:ub_axis1-2,:,:,input1%id) * work_xyz(lb_axis1:ub_axis1-2,:,:,input2%id) ) / (2.0D0*dx) & - sign( 1.0D0, work_xyz(lb_axis1+1:ub_axis1-1,:,:,input1%id) )*dx/2.0D0 & * ( work_xyz(lb_axis1+2:ub_axis1,:,:,input1%id) * work_xyz(lb_axis1+2:ub_axis1,:,:,input2%id) & - 2.0D0*work_xyz(lb_axis1+1:ub_axis1-1,:,:,input1%id) * work_xyz(lb_axis1+1:ub_axis1-1,:,:,input2%id) & + work_xyz(lb_axis1:ub_axis1-2,:,:,input1%id) * work_xyz(lb_axis1:ub_axis1-2,:,:,input2%id) )/(dx**2) end function uw1_x_xyz function uw1_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+1 : ub_axis2-1,:,new_id) & = ( work_y(:,lb_axis2+2:ub_axis2,:,input1%id) * work_y(:,lb_axis2+2:ub_axis2,:,input2%id) & - work_y(:,lb_axis2:ub_axis2-2,:,input1%id) * work_y(:,lb_axis2:ub_axis2-2,:,input2%id) ) / (2.0D0*dy) & - sign(1.0D0, work_y(:,lb_axis2+1:ub_axis2-1,:,input1%id))*dy/2.0D0 & * ( work_y(:,lb_axis2+2:ub_axis2,:,input1%id)* work_y(:,lb_axis2+2:ub_axis2,:,input2%id) & - 2.0D0*work_y(:,lb_axis2+1:ub_axis2-1,:,input1%id)* work_y(:,lb_axis2+1:ub_axis2-1,:,input2%id) & + work_y(:,lb_axis2:ub_axis2-2,:,input1%id)* work_y(:,lb_axis2:ub_axis2-2,:,input2%id) )/(dy**2) end function uw1_y_y function uw1_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+1 : ub_axis2-1,:,new_id) & = ( work_xy(:,lb_axis2+2:ub_axis2,:,input1%id) * work_xy(:,lb_axis2+2:ub_axis2,:,input2%id) & - work_xy(:,lb_axis2:ub_axis2-2,:,input1%id) * work_xy(:,lb_axis2:ub_axis2-2,:,input2%id) ) / (2.0D0*dy) & - sign(1.0D0, work_xy(:,lb_axis2+1:ub_axis2-1,:,input1%id))*dy/2.0D0 & * ( work_xy(:,lb_axis2+2:ub_axis2,:,input1%id)* work_xy(:,lb_axis2+2:ub_axis2,:,input2%id) & - 2.0D0*work_xy(:,lb_axis2+1:ub_axis2-1,:,input1%id)* work_xy(:,lb_axis2+1:ub_axis2-1,:,input2%id) & + work_xy(:,lb_axis2:ub_axis2-2,:,input1%id)* work_xy(:,lb_axis2:ub_axis2-2,:,input2%id) )/(dy**2) end function uw1_y_xy function uw1_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+1 : ub_axis2-1,:,new_id) & = ( work_yz(:,lb_axis2+2:ub_axis2,:,input1%id) * work_yz(:,lb_axis2+2:ub_axis2,:,input2%id) & - work_yz(:,lb_axis2:ub_axis2-2,:,input1%id) * work_yz(:,lb_axis2:ub_axis2-2,:,input2%id) ) / (2.0D0*dy) & - sign(1.0D0, work_yz(:,lb_axis2+1:ub_axis2-1,:,input1%id))*dy/2.0D0 & * ( work_yz(:,lb_axis2+2:ub_axis2,:,input1%id)* work_yz(:,lb_axis2+2:ub_axis2,:,input2%id) & - 2.0D0*work_yz(:,lb_axis2+1:ub_axis2-1,:,input1%id)* work_yz(:,lb_axis2+1:ub_axis2-1,:,input2%id) & + work_yz(:,lb_axis2:ub_axis2-2,:,input1%id)* work_yz(:,lb_axis2:ub_axis2-2,:,input2%id) )/(dy**2) end function uw1_y_yz function uw1_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+1 : ub_axis2-1,:,new_id) & = ( work_xyz(:,lb_axis2+2:ub_axis2,:,input1%id) * work_xyz(:,lb_axis2+2:ub_axis2,:,input2%id) & - work_xyz(:,lb_axis2:ub_axis2-2,:,input1%id) * work_xyz(:,lb_axis2:ub_axis2-2,:,input2%id) ) / (2.0D0*dy) & - sign(1.0D0, work_xyz(:,lb_axis2+1:ub_axis2-1,:,input1%id))*dy/2.0D0 & * ( work_xyz(:,lb_axis2+2:ub_axis2,:,input1%id)* work_xyz(:,lb_axis2+2:ub_axis2,:,input2%id) & - 2.0D0*work_xyz(:,lb_axis2+1:ub_axis2-1,:,input1%id)*work_xyz(:,lb_axis2+1:ub_axis2-1,:,input2%id) & + work_xyz(:,lb_axis2:ub_axis2-2,:,input1%id)* work_xyz(:,lb_axis2:ub_axis2-2,:,input2%id) )/(dy**2) end function uw1_y_xyz function uw1_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+1 : ub_axis3-1,new_id) & = ( work_z(:,:,lb_axis3+2:ub_axis3,input1%id) * work_z(:,:,lb_axis3+2:ub_axis3,input2%id) & - work_z(:,:,lb_axis3:ub_axis3-2,input1%id) * work_z(:,:,lb_axis3:ub_axis3-2,input2%id) ) / (2.0D0*dz) & - sign(1.0D0, work_z(:,:,lb_axis3+1:ub_axis3-1,input1%id))*dz/2.0D0 & * ( work_z(:,:,lb_axis3+2:ub_axis3,input1%id)* work_z(:,:,lb_axis3+2:ub_axis3,input2%id) & - 2.0D0*work_z(:,:,lb_axis3+1:ub_axis3-1,input1%id)* work_z(:,:,lb_axis3+1:ub_axis3-1,input2%id) & + work_z(:,:,lb_axis3:ub_axis3-2,input1%id)* work_z(:,:,lb_axis3:ub_axis3-2,input2%id) )/(dz**2) end function uw1_z_z function uw1_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+1 : ub_axis3-1,new_id) & = ( work_xz(:,:,lb_axis3+2:ub_axis3,input1%id) * work_xz(:,:,lb_axis3+2:ub_axis3,input2%id) & - work_xz(:,:,lb_axis3:ub_axis3-2,input1%id) * work_xz(:,:,lb_axis3:ub_axis3-2,input2%id) ) / (2.0D0*dz) & - sign(1.0D0, work_xz(:,:,lb_axis3+1:ub_axis3-1,input1%id))*dz/2.0D0 & * ( work_xz(:,:,lb_axis3+2:ub_axis3,input1%id)* work_xz(:,:,lb_axis3+2:ub_axis3,input2%id) & - 2.0D0*work_xz(:,:,lb_axis3+1:ub_axis3-1,input1%id)* work_xz(:,:,lb_axis3+1:ub_axis3-1,input2%id) & + work_xz(:,:,lb_axis3:ub_axis3-2,input1%id)* work_xz(:,:,lb_axis3:ub_axis3-2,input2%id) )/(dz**2) end function uw1_z_xz function uw1_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+1 : ub_axis3-1,new_id) & = ( work_yz(:,:,lb_axis3+2:ub_axis3,input1%id) * work_yz(:,:,lb_axis3+2:ub_axis3,input2%id) & - work_yz(:,:,lb_axis3:ub_axis3-2,input1%id) * work_yz(:,:,lb_axis3:ub_axis3-2,input2%id) ) / (2.0D0*dz) & - sign(1.0D0, work_yz(:,:,lb_axis3+1:ub_axis3-1,input1%id))*dz/2.0D0 & * ( work_yz(:,:,lb_axis3+2:ub_axis3,input1%id)* work_yz(:,:,lb_axis3+2:ub_axis3,input2%id) & - 2.0D0*work_yz(:,:,lb_axis3+1:ub_axis3-1,input1%id)* work_yz(:,:,lb_axis3+1:ub_axis3-1,input2%id) & + work_yz(:,:,lb_axis3:ub_axis3-2,input1%id)* work_yz(:,:,lb_axis3:ub_axis3-2,input2%id) )/(dz**2) end function uw1_z_yz function uw1_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+1 : ub_axis3-1,new_id) & = ( work_xyz(:,:,lb_axis3+2:ub_axis3,input1%id) * work_xyz(:,:,lb_axis3+2:ub_axis3,input2%id) & - work_xyz(:,:,lb_axis3:ub_axis3-2,input1%id) * work_xyz(:,:,lb_axis3:ub_axis3-2,input2%id) ) / (2.0D0*dz) & - sign(1.0D0, work_xyz(:,:,lb_axis3+1:ub_axis3-1,input1%id))*dz/2.0D0 & * ( work_xyz(:,:,lb_axis3+2:ub_axis3,input1%id)* work_xyz(:,:,lb_axis3+2:ub_axis3,input2%id) & - 2.0D0*work_xyz(:,:,lb_axis3+1:ub_axis3-1,input1%id)*work_xyz(:,:,lb_axis3+1:ub_axis3-1,input2%id) & + work_xyz(:,:,lb_axis3:ub_axis3-2,input1%id)* work_xyz(:,:,lb_axis3:ub_axis3-2,input2%id) )/(dz**2) end function uw1_z_xyz end module gms_math_upwind1