Line data    Source code

       1             : module gw_utils
       2             : 
       3             : !
       4             : ! This module contains utility code for the gravity wave modules.
       5             : !
       6             : 
       7             : implicit none
       8             : private
       9             : save
      10             : 
      11             : ! Real kind for gravity wave parameterization.
      12             : integer, public, parameter :: r8 = selected_real_kind(12)
      13             : 
      14             : ! Public interface
      15             : public :: get_unit_vector
      16             : public :: dot_2d
      17             : public :: midpoint_interp
      18             : 
      19             : contains
      20             : 
      21             : ! Take two components of a vector, and find the unit vector components and
      22             : ! total magnitude.
      23     1489176 : subroutine get_unit_vector(u, v, u_n, v_n, mag)
      24             :   real(r8), intent(in) :: u(:)
      25             :   real(r8), intent(in) :: v(:)
      26             :   real(r8), intent(out) :: u_n(:)
      27             :   real(r8), intent(out) :: v_n(:)
      28             :   real(r8), intent(out) :: mag(:)
      29             : 
      30             :   integer :: i
      31             : 
      32    24865776 :   mag = sqrt(u*u + v*v)
      33             : 
      34             :   ! Has to be a loop/if instead of a where, because floating point
      35             :   ! exceptions can trigger even on a masked divide-by-zero operation
      36             :   ! (especially on Intel).
      37    24865776 :   do i = 1, size(mag)
      38    24865776 :      if (mag(i) > 0._r8) then
      39    23376600 :         u_n(i) = u(i)/mag(i)
      40    23376600 :         v_n(i) = v(i)/mag(i)
      41             :      else
      42           0 :         u_n(i) = 0._r8
      43           0 :         v_n(i) = 0._r8
      44             :      end if
      45             :   end do
      46             : 
      47     1489176 : end subroutine get_unit_vector
      48             : 
      49             : ! Vectorized version of a 2D dot product (since the intrinsic dot_product
      50             : ! is more suitable for arrays of contiguous vectors).
      51    38718576 : function dot_2d(u1, v1, u2, v2)
      52             :   real(r8), intent(in) :: u1(:), v1(:)
      53             :   real(r8), intent(in) :: u2(:), v2(:)
      54             : 
      55             :   real(r8) :: dot_2d(size(u1))
      56             : 
      57   685228752 :   dot_2d = u1*u2 + v1*v2
      58             : 
      59             : end function dot_2d
      60             : 
      61             : ! Pure function that interpolates the values of the input array along
      62             : ! dimension 2. This is obviously not a very generic routine, unlike, say,
      63             : ! CAM's lininterp. But it's used often enough that it seems worth providing
      64             : ! here.
      65     5956704 : pure function midpoint_interp(arr) result(interp)
      66             :   real(r8), intent(in) :: arr(:,:)
      67             :   real(r8) :: interp(size(arr,1),size(arr,2)-1)
      68             : 
      69             :   integer :: i
      70             : 
      71   156363480 :   do i = 1, size(interp,2)
      72  2517400080 :      interp(:,i) = 0.5_r8 * (arr(:,i)+arr(:,i+1))
      73             :   end do
      74             : 
      75     5956704 : end function midpoint_interp
      76             : 
      77             : end module gw_utils