diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index b62865598e..9950769f97 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -84,6 +84,27 @@ subroutine turbotmp_ppm_reconstruction_y_bridge(bx, h_in, h_S, h_N, mask2dT, & end subroutine turbotmp_ppm_reconstruction_y_bridge end interface + interface + !> Bridge for the PPM_reconstruction_x subroutine + subroutine turbotmp_ppm_reconstruction_x_bridge(bx, h_in, h_W, h_E, mask2dT, & + h_min, monotonic, simple_2nd, obc) bind(C) + use iso_c_binding, only : c_double, c_bool, c_ptr + use array_mod, only : RealArray_c + use box_mod, only : Box_c + implicit none + + type(Box_C), intent(in) :: bx !< Index space over which to iterate + type(RealArray_C), intent(in) :: h_in !< Layer thickness + type(RealArray_C), intent(inout) :: h_W !< West edge thickness + type(RealArray_C), intent(inout) :: h_E !< East edge thickness + type(RealArray_C), intent(in) :: mask2dT !< Mask (0 land, 1 ocean) + real(c_double), intent(in), value :: h_min !< Minimum thickness + logical(c_bool), intent(in), value :: monotonic !< Use CW84 limiter + logical(c_bool), intent(in), value :: simple_2nd !< Use 2nd order scheme + type(c_ptr), intent(in), value :: obc !< Pointer to OBC structure + end subroutine turbotmp_ppm_reconstruction_x_bridge + end interface + #include public continuity_PPM, continuity_PPM_init, continuity_PPM_stencil @@ -511,6 +532,7 @@ subroutine zonal_edge_thickness(h_in, h_W, h_E, G, GV, US, CS, OBC, LB_in) integer :: isl, iel, jsl, jel, stencil type(Box_t) :: bx, bxH character(len=256) :: mesg + type(RealArray_t) :: h_in_a, h_W_a, h_E_a, mask2dT_a if (present(LB_in)) then LB = LB_in @@ -553,8 +575,21 @@ subroutine zonal_edge_thickness(h_in, h_W, h_E, G, GV, US, CS, OBC, LB_in) h_W(i,j,k) = h_in(i,j,k) ; h_E(i,j,k) = h_in(i,j,k) enddo else - call PPM_reconstruction_x(bxH, h_in, h_W, h_E, G, GV, G%mask2dT, & + ! Duplicate the arrays + call h_in_a%dup(h_in) + call h_W_a%dup(h_W) + call h_E_a%dup(h_E) + call mask2dT_a%dup(G%mask2dT) + + ! Copy data into array containers + call h_in_a%copy2Array(h_in) + call mask2dT_a%copy2Array(G%mask2dT) + + ! Calculate the reconstruction + call PPM_reconstruction_x(bxH, h_in_a, h_W_a, h_E_a, mask2dT_a, & 2.0*GV%Angstrom_H, CS%monotonic, CS%simple_2nd, OBC) + call h_W_a%copy2F(h_W) + call h_E_a%copy2F(h_E) endif ! Free memory associated with index boxes @@ -563,6 +598,12 @@ subroutine zonal_edge_thickness(h_in, h_W, h_E, G, GV, US, CS, OBC, LB_in) call cpu_clock_end(id_clock_reconstruct) + ! Free up temporary containers + call h_in_a%free() + call h_W_a%free() + call h_E_a%free() + call mask2dT_a%free() + end subroutine zonal_edge_thickness @@ -2453,17 +2494,15 @@ subroutine set_merid_BT_cont(v, h_in, h_S, h_N, BT_cont, vh_tot_0, dvhdv_tot_0, end subroutine set_merid_BT_cont !> Calculates left/right edge values for PPM reconstruction. -subroutine PPM_reconstruction_x(bxH, h_in, h_W, h_E, G, GV, mask2dT, h_min, monotonic, simple_2nd, OBC) +subroutine PPM_reconstruction_x_fortran(bxH, h_in_a, h_W_a, h_E_a, mask2dT_a, h_min, monotonic, simple_2nd, OBC) type(Box_t), intent(in) :: bxH !< H-grid iteration Box - type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure. - type(verticalgrid_type), intent(in) :: GV !< Ocean's vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness [H ~> m or kg m-2]. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h_W !< West edge thickness in the reconstruction, - !! [H ~> m or kg m-2]. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h_E !< East edge thickness in the reconstruction, - !! [H ~> m or kg m-2]. - real, dimension(SZI_(g),SZJ_(G)), intent(in) :: mask2dT !< 0 for land points and 1 for ocean points - !! on the h-grid [nondim] + type(RealArray_t), intent(in) :: h_in_a !< Layer thickness [H ~> m or kg m-2]. + type(RealArray_t), intent(inout) :: h_W_a !< West edge thickness in the reconstruction + !! [H ~> m or kg m-2]. + type(RealArray_t), intent(inout) :: h_E_a !< East edge thickness in the reconstruction + !! [H ~> m or kg m-2]. + type(RealArray_t), intent(in) :: mask2dT_a !< 0 for land points and 1 for ocean points + !! on the h-grid [nondim] real, intent(in) :: h_min !< The minimum thickness !! that can be obtained by a concave parabolic fit [H ~> m or kg m-2] logical, intent(in) :: monotonic !< If true, use the @@ -2475,7 +2514,7 @@ subroutine PPM_reconstruction_x(bxH, h_in, h_W, h_E, G, GV, mask2dT, h_min, mono type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. ! Local variables with useful mnemonic names. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: slp ! The slopes per grid point [H ~> m or kg m-2] + real, dimension(:,:,:), allocatable :: slp ! The slopes per grid point [H ~> m or kg m-2] real, parameter :: oneSixth = 1./6. ! [nondim] real :: h_ip1, h_im1 ! Neighboring thicknesses or sensibly extrapolated values [H ~> m or kg m-2] real :: dMx, dMn ! The difference between the local thickness and the maximum (dMx) or @@ -2484,9 +2523,21 @@ subroutine PPM_reconstruction_x(bxH, h_in, h_W, h_E, G, GV, mask2dT, h_min, mono integer :: n logical :: local_open_BC type(OBC_segment_type), pointer :: segment => NULL() - integer, parameter :: ndims = 3 type(Box_t) :: bx, bxE - type(RealArray_t) :: h_in_a, h_W_a, h_E_a + + real, dimension(:,:), contiguous, pointer :: mask2dT + real, dimension(:,:,:), contiguous, pointer :: h_in + real, dimension(:,:,:), contiguous, pointer :: h_W + real, dimension(:,:,:), contiguous, pointer :: h_E + + ! Get the views for containers (subroutine arguments) + call h_W_a%view(h_W) + call h_E_a%view(h_E) + call h_in_a%view(h_in) + call mask2dT_a%view(mask2dT) + + ! Allocate local slope array using bounds from h_in_a + allocate(slp(h_in_a%lb(1):h_in_a%ub(1),h_in_a%lb(2):h_in_a%ub(2),h_in_a%lb(3):h_in_a%ub(3))) local_open_BC = .false. if (associated(OBC)) then @@ -2579,38 +2630,22 @@ subroutine PPM_reconstruction_x(bxH, h_in, h_W, h_E, G, GV, mask2dT, h_min, mono enddo endif - ! Duplicate the arrays - call h_in_a%dup(h_in) - call h_W_a%dup(h_W) - call h_E_a%dup(h_E) - - ! Copy data into array containers - call h_in_a%copy2Array(h_in) - call h_W_a%copy2Array(h_W) - call h_E_a%copy2Array(h_E) - if (monotonic) then call PPM_limit_cw84(bx, h_in_a, h_W_a, h_E_a) else call PPM_limit_pos(bx, h_in_a, h_W_a, h_E_a, h_min) endif - ! Copy data back to Fortran arrays - call h_W_a%copy2F(h_W) - call h_E_a%copy2F(h_E) - - ! Free up temporary containers - call h_in_a%free() - call h_W_a%free() - call h_E_a%free() + ! Deallocate local temporary array + if(allocated(slp)) deallocate(slp) - ! Free up the iteration space boxes + ! Deallocate local iteration boxes call bx%free() call bxE%free() return -end subroutine PPM_reconstruction_x +end subroutine PPM_reconstruction_x_fortran !> Calculates left/right edge values for PPM reconstruction. subroutine PPM_reconstruction_y_fortran(bxH, h_in_a, h_S_a, h_N_a, mask2dT_a, h_min, monotonic, simple_2nd, OBC) @@ -3244,6 +3279,90 @@ subroutine PPM_reconstruction_y(bxH, h_in_a, h_S_a, h_N_a, mask2dT_a, h_min, mon end subroutine PPM_reconstruction_y +!< shim for PPM_reconstruction_x +subroutine PPM_reconstruction_x(bxH, h_in_a, h_W_a, h_E_a, mask2dT_a, h_min, monotonic, simple_2nd, OBC) + implicit none + + type(Box_t), intent(in) :: bxH !< H-grid iteration Box + type(RealArray_t), intent(in) :: h_in_a !< Layer thickness + type(RealArray_t), intent(inout) :: h_W_a !< West edge thickness + type(RealArray_t), intent(inout) :: h_E_a !< East edge thickness + type(RealArray_t), intent(in) :: mask2dT_a !< Mask (0 land, 1 ocean) + real, intent(in) :: h_min !< Minimum thickness + logical, intent(in) :: monotonic !< Use CW84 limiter + logical, intent(in) :: simple_2nd !< Use simple 2nd order + type(ocean_OBC_type), pointer :: OBC !< Open boundary control + + ! local variables + integer :: mode + type(Box_C) :: bx_c + type(RealArray_C) :: h_in_c, h_W_c, h_E_c, mask2dT_c + type(c_ptr) :: OBC_c + logical(c_bool) :: monotonic_c, simple_2nd_c + integer :: rc + type (io_recorder) :: rec + logical :: capture + character(len=80) :: kernel + character(len=100) :: dir + character(len=256) :: binFile, metaFile + + kernel="ppm_reconstruction_x" + + mode = getenv_mode("PPM_RECONSTRUCTION_X_MODE", default=TIMH_runFORTRAN) + + select case (mode) + + case (TIMH_capture) + capture = .FALSE. + if((.not. already_recorded(TRIM(kernel))) .and. is_root_pe()) capture = .TRUE. + + if(capture) then + dir = "capture" + rc = mkdir_posix(TRIM(dir) // c_null_char, int(o'755', c_int)) + + binFile = TRIM(dir) // "/" // TRIM(kernel) // ".bin" + metaFile = TRIM(dir) // "/" // TRIM(kernel) // ".meta" + + call rec%open_write(binFile, metaFile) + + call rec%add("_bxH", bxH) + call rec%add("_h_in", h_in_a ) + call rec%add("_h_W_before", h_W_a ) + call rec%add("_h_E_before", h_E_a ) + call rec%add("_mask2d_t", mask2dT_a ) + call rec%add("_h_min", h_min) + call rec%add("_monotonic", monotonic) + call rec%add("_simple_2nd", simple_2nd) + !call rec%add("OBC", OBC) + endif + + call ppm_reconstruction_x_fortran(bxH, h_in_a, h_W_a, h_E_a, & + mask2dT_a, h_min, monotonic, simple_2nd, OBC) + if(capture) then + call rec%add("_h_W_after", h_W_a) + call rec%add("_h_E_after", h_E_a) + call rec%close() + call mark_recorded(TRIM(kernel)) + endif +#ifdef _TIM + case (TIMH_runAMREX) + + ! create C-compatible descriptors + bx_c = bxH%to_c(); h_in_c = h_in_a%to_c(); h_W_c = h_W_a%to_c(); + h_E_c = h_E_a%to_c(); mask2dT_c = mask2dT_a%to_c() + monotonic_c = monotonic; simple_2nd_c = simple_2nd + if(associated(OBC)) then; OBC_c = c_loc(OBC); else; OBC_c = c_null_ptr; endif + ! Call C++ bridge to execute AMReX code + call turbotmp_ppm_reconstruction_x_bridge(bx_c, h_in_c, h_W_c, h_E_c, mask2dT_c, & + h_min, monotonic_c, simple_2nd_c, OBC_c) +#endif + case default + call ppm_reconstruction_x_fortran(bxH, h_in_a, h_W_a, h_E_a, & + mask2dT_a, h_min, monotonic, simple_2nd, OBC) + end select + +end subroutine PPM_reconstruction_x + !> \namespace mom_continuity_ppm !! !! This module contains the subroutines that advect layer