Implement Intel OneMKL in /src/dynamics/operators_mod.F90: calc_ph()

[GrassBlock2016]

  subroutine calc_ph(block, dstate)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate
    real(r8), allocatable :: rd_o_cpd_vec(:,:,:)

    integer i, j, k

    call perf_start('calc_ph')

    associate (mesh    => block%mesh          , &
               mg_lev  => dstate%mg_lev       , & ! in
               dmg     => dstate%dmg          , & ! in
               qm      => tracers(block%id)%qm, & ! in
               ph_lev  => dstate%ph_lev       , & ! out
               pkh_lev => block%aux%pkh_lev   , & ! out
               ph      => dstate%ph           , & ! out
               phs     => dstate%phs          , & ! pointer
               ps      => dstate%ps           )   ! out
    k = mesh%half_kds

    allocate(rd_o_cpd_vec(mesh%full_ids:mesh%full_ide+1, mesh%full_jds:mesh%full_jde + merge(0, 1, mesh%has_north_pole()), mesh%half_kds + 1:mesh%half_kde))
    rd_o_cpd_vec(:,:,:) = rd_o_cpd

    ph_lev%d(:,:,k) = mg_lev%d(:,:,k)
    pkh_lev%d(:,:,k) = ph_lev%d(:,:,k)**rd_o_cpd
    do k = mesh%half_kds + 1, mesh%half_kde
      do j = mesh%full_jds, mesh%full_jde + merge(0, 1, mesh%has_north_pole())
        do i = mesh%full_ids, mesh%full_ide - 1, 2
          ph_lev%d(i,j,k) = ph_lev%d(i,j,k-1) + dmg%d(i,j,k-1) * (1 + qm%d(i,j,k-1))
          ! pkh_lev%d(i,j,k) = ph_lev%d(i,j,k)**rd_o_cpd
          ph_lev%d(i+1,j,k) = ph_lev%d(i+1,j,k-1) + dmg%d(i+1,j,k-1) * (1 + qm%d(i+1,j,k-1))
          ! pkh_lev%d(i+1,j,k) = ph_lev%d(i+1,j,k)**rd_o_cpd
        end do
        do i = mesh%full_ide - mod(mesh%full_ide - mesh%full_ids, 2), mesh%full_ide + 1
          ph_lev%d(i,j,k) = ph_lev%d(i,j,k-1) + dmg%d(i,j,k-1) * (1 + qm%d(i,j,k-1))
          ! pkh_lev%d(i,j,k) = ph_lev%d(i,j,k)**rd_o_cpd
        end do
      end do
    end do
    call vdpow((mesh%full_ide-mesh%full_ids+2)*(mesh%full_jde+merge(0, 1, mesh%has_north_pole())-mesh%full_jds+1)*(mesh%half_kde-mesh%half_kds), ph_lev%d, rd_o_cpd_vec, pkh_lev%d)
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds, mesh%full_jde + merge(0, 1, mesh%has_north_pole())
        do i = mesh%full_ids, mesh%full_ide + 1
          ph%d(i,j,k) = 0.5_r8 * (ph_lev%d(i,j,k) + ph_lev%d(i,j,k+1))
        end do
      end do
    end do
    ! NOTE: Move this to other place?
    if (hydrostatic) ps%d = phs%d
    end associate

    call perf_stop('calc_ph')

  end subroutine calc_ph

However, it will report errors while running bw.180x90 (No error while building).