MODULE advect_mod
   USE icosa
   IMPLICIT NONE

   PRIVATE

   PUBLIC :: init_advect

CONTAINS

   !==========================================================================

   SUBROUTINE init_advect(normal, tangent, sqrt_leng)
      IMPLICIT NONE
      REAL(rstd), INTENT(OUT) :: normal(3*iim*jjm, 3)
      REAL(rstd), INTENT(OUT) :: tangent(3*iim*jjm, 3)
      REAL(rstd), INTENT(OUT) :: sqrt_leng(iim*jjm)

      INTEGER :: ij

!DIR$ SIMD
      DO ij = ij_begin, ij_end

         CALL cross_product2(xyz_v(ij + z_rdown, :), xyz_v(ij + z_rup, :), normal(ij + u_right, :))
         normal(ij + u_right, :) = normal(ij + u_right, :)/sqrt(sum(normal(ij + u_right, :)**2) + 1e-50)*ne(ij, right)

         CALL cross_product2(xyz_v(ij + z_up, :), xyz_v(ij + z_lup, :), normal(ij + u_lup, :))
         normal(ij + u_lup, :) = normal(ij + u_lup, :)/sqrt(sum(normal(ij + u_lup, :)**2) + 1e-50)*ne(ij, lup)

         CALL cross_product2(xyz_v(ij + z_ldown, :), xyz_v(ij + z_down, :), normal(ij + u_ldown, :))
         normal(ij + u_ldown, :) = normal(ij + u_ldown, :)/sqrt(sum(normal(ij + u_ldown, :)**2) + 1e-50)*ne(ij, ldown)

         tangent(ij + u_right, :) = xyz_v(ij + z_rup, :) - xyz_v(ij + z_rdown, :)
         tangent(ij + u_right, :) = tangent(ij + u_right, :)/sqrt(sum(tangent(ij + u_right, :)**2) + 1e-50)

         tangent(ij + u_lup, :) = xyz_v(ij + z_lup, :) - xyz_v(ij + z_up, :)
         tangent(ij + u_lup, :) = tangent(ij + u_lup, :)/sqrt(sum(tangent(ij + u_lup, :)**2) + 1e-50)

         tangent(ij + u_ldown, :) = xyz_v(ij + z_down, :) - xyz_v(ij + z_ldown, :)
         tangent(ij + u_ldown, :) = tangent(ij + u_ldown, :)/sqrt(sum(tangent(ij + u_ldown, :)**2) + 1e-50)

         sqrt_leng(ij) = sqrt(max(sum((xyz_v(ij + z_up, :) - xyz_i(ij, :))**2), sum((xyz_v(ij + z_down, :) - xyz_i(ij, :))**2), &
                                  sum((xyz_v(ij + z_rup, :) - xyz_i(ij, :))**2), sum((xyz_v(ij + z_rdown, :) - xyz_i(ij, :))**2), &
                                  sum((xyz_v(ij + z_lup, :) - xyz_i(ij, :))**2), sum((xyz_v(ij + z_ldown, :) - xyz_i(ij, :))**2)))
      END DO

   END SUBROUTINE init_advect

END MODULE advect_mod
