diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 6aaf8e0355..0e30dfc250 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -24,7 +24,6 @@ module MOM_ALE use MOM_string_functions, only : uppercase, extractWord use MOM_verticalGrid, only : verticalGrid_type -use regrid_ppoly_class, only : ppoly_t, ppoly_init, ppoly_destroy use regrid_edge_values, only : edgeValueArrays use regrid_edge_values, only : edge_values_implicit_h4 use regrid_edge_values, only : triDiagEdgeWorkAllocate, triDiagEdgeWorkDeallocate @@ -67,14 +66,6 @@ module MOM_ALE type, public :: ALE_CS private - ! Generic linear piecewise polynomial used for various purposes throughout - ! the code (the same ppoly is used to avoid having dynamical memory allocation) - type(ppoly_t) :: ppoly_linear - - ! Generic parabolic piecewise polynomial used for various purposes - ! throughout the code (the same ppoly is used to avoid having dynamical - ! memory allocation) - type(ppoly_t) :: ppoly_parab ! Indicate whether high-order boundary extrapolation should be used within ! boundary cells @@ -97,6 +88,10 @@ module MOM_ALE ! Used only for queries, not directly by this module integer :: nk + integer :: degree_linear=1 ! Degree of linear piecewise polynomial + integer :: degree_parab=2 ! Degree of parabolic piecewise polynomial + + ! Work space for communicating between regridding and remapping real ALLOCABLE_, dimension(NIMEM_,NJMEM_,NK_INTERFACE_) :: dzRegrid @@ -325,12 +320,18 @@ subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, tv, h ) ! Local variables integer :: i, j, k real :: hTmp(G%ke) + real, dimension(CS%nk,2) :: & + ppoly_linear_E !Edge value of polynomial + real, dimension(CS%nk,CS%degree_linear+1) :: & + ppoly_linear_coefficients !Coefficients of polynomial + ! NOTE: the variables 'CS%grid_generic' and 'CS%ppoly_linear' are declared at ! the module level. Memory is allocated once at the beginning of the run ! in 'ALE_memory_allocation'. ! Determine reconstruction within each column +!$OMP parallel do default(shared) private(hTmp,ppoly_linear_E,ppoly_linear_coefficients) do i = G%isc,G%iec+1 do j = G%jsc,G%jec+1 @@ -338,27 +339,27 @@ subroutine pressure_gradient_plm( CS, S_t, S_b, T_t, T_b, G, tv, h ) hTmp(:) = h(i,j,:) ! Reconstruct salinity profile - CS%ppoly_linear%E = 0.0 - CS%ppoly_linear%coefficients = 0.0 - call PLM_reconstruction( G%ke, hTmp, tv%S(i,j,:), CS%ppoly_linear ) + ppoly_linear_E = 0.0 + ppoly_linear_coefficients = 0.0 + call PLM_reconstruction( G%ke, hTmp, tv%S(i,j,:), ppoly_linear_E, ppoly_linear_coefficients ) if (CS%boundary_extrapolation_for_pressure) call & - PLM_boundary_extrapolation( G%ke, hTmp, tv%S(i,j,:), CS%ppoly_linear ) + PLM_boundary_extrapolation( G%ke, hTmp, tv%S(i,j,:), ppoly_linear_E, ppoly_linear_coefficients ) do k = 1,G%ke - S_t(i,j,k) = CS%ppoly_linear%E(k,1) - S_b(i,j,k) = CS%ppoly_linear%E(k,2) + S_t(i,j,k) = ppoly_linear_E(k,1) + S_b(i,j,k) = ppoly_linear_E(k,2) end do ! Reconstruct temperature profile - CS%ppoly_linear%E = 0.0 - CS%ppoly_linear%coefficients = 0.0 - call PLM_reconstruction( G%ke, hTmp, tv%T(i,j,:), CS%ppoly_linear ) + ppoly_linear_E = 0.0 + ppoly_linear_coefficients = 0.0 + call PLM_reconstruction( G%ke, hTmp, tv%T(i,j,:), ppoly_linear_E, ppoly_linear_coefficients ) if (CS%boundary_extrapolation_for_pressure) call & - PLM_boundary_extrapolation( G%ke, hTmp, tv%T(i,j,:), CS%ppoly_linear ) + PLM_boundary_extrapolation( G%ke, hTmp, tv%T(i,j,:), ppoly_linear_E, ppoly_linear_coefficients ) do k = 1,G%ke - T_t(i,j,k) = CS%ppoly_linear%E(k,1) - T_b(i,j,k) = CS%ppoly_linear%E(k,2) + T_t(i,j,k) = ppoly_linear_E(k,1) + T_b(i,j,k) = ppoly_linear_E(k,2) end do end do @@ -392,12 +393,18 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, tv, h ) ! Local variables integer :: i, j, k real :: hTmp(G%ke) + real, dimension(CS%nk,2) :: & + ppoly_parab_E !Edge value of polynomial + real, dimension(CS%nk,CS%degree_parab+1) :: & + ppoly_parab_coefficients !Coefficients of polynomial + ! NOTE: the variables 'CS%grid_generic' and 'CS%ppoly_parab' are declared at ! the module level. Memory is allocated once at the beginning of the run ! in 'ALE_memory_allocation'. ! Determine reconstruction within each column +!$OMP parallel do default(shared) private(hTmp,ppoly_parab_E,ppoly_parab_coefficients) do i = G%isc,G%iec+1 do j = G%jsc,G%jec+1 @@ -405,29 +412,29 @@ subroutine pressure_gradient_ppm( CS, S_t, S_b, T_t, T_b, G, tv, h ) hTmp(:) = h(i,j,:) ! Reconstruct salinity profile - CS%ppoly_parab%E = 0.0 - CS%ppoly_parab%coefficients = 0.0 - call edge_values_implicit_h4( G%ke, hTmp, tv%S(i,j,:), CS%edgeValueWrk, CS%ppoly_parab%E ) - call PPM_reconstruction( G%ke, hTmp, tv%S(i,j,:), CS%ppoly_parab ) + ppoly_parab_E = 0.0 + ppoly_parab_coefficients = 0.0 + call edge_values_implicit_h4( G%ke, hTmp, tv%S(i,j,:), CS%edgeValueWrk, ppoly_parab_E ) + call PPM_reconstruction( G%ke, hTmp, tv%S(i,j,:), ppoly_parab_E, ppoly_parab_coefficients ) if (CS%boundary_extrapolation_for_pressure) call & - PPM_boundary_extrapolation( G%ke, hTmp, tv%S(i,j,:), CS%ppoly_parab ) + PPM_boundary_extrapolation( G%ke, hTmp, tv%S(i,j,:), ppoly_parab_E, ppoly_parab_coefficients ) do k = 1,G%ke - S_t(i,j,k) = CS%ppoly_parab%E(k,1) - S_b(i,j,k) = CS%ppoly_parab%E(k,2) + S_t(i,j,k) = ppoly_parab_E(k,1) + S_b(i,j,k) = ppoly_parab_E(k,2) end do ! Reconstruct temperature profile - CS%ppoly_parab%E = 0.0 - CS%ppoly_parab%coefficients = 0.0 - call edge_values_implicit_h4( G%ke, hTmp, tv%T(i,j,:), CS%edgeValueWrk, CS%ppoly_parab%E ) - call PPM_reconstruction( G%ke, hTmp, tv%T(i,j,:), CS%ppoly_parab ) + ppoly_parab_E = 0.0 + ppoly_parab_coefficients = 0.0 + call edge_values_implicit_h4( G%ke, hTmp, tv%T(i,j,:), CS%edgeValueWrk, ppoly_parab_E ) + call PPM_reconstruction( G%ke, hTmp, tv%T(i,j,:), ppoly_parab_E, ppoly_parab_coefficients ) if (CS%boundary_extrapolation_for_pressure) call & - PPM_boundary_extrapolation( G%ke, hTmp, tv%T(i,j,:), CS%ppoly_parab ) + PPM_boundary_extrapolation( G%ke, hTmp, tv%T(i,j,:), ppoly_parab_E, ppoly_parab_coefficients ) do k = 1,G%ke - T_t(i,j,k) = CS%ppoly_parab%E(k,1) - T_b(i,j,k) = CS%ppoly_parab%E(k,2) + T_t(i,j,k) = ppoly_parab_E(k,1) + T_b(i,j,k) = ppoly_parab_E(k,2) end do end do @@ -462,12 +469,6 @@ subroutine ALE_memory_allocation( G, CS ) call triDiagEdgeWorkAllocate( nz, CS%edgeValueWrk ) call triDiagSlopeWorkAllocate( nz, CS%edgeSlopeWrk ) - ! Generic linear piecewise polynomial - call ppoly_init( CS%ppoly_linear, nz, 1 ) - - ! Generic parabolic piecewise polynomial - call ppoly_init( CS%ppoly_parab, nz, 2 ) - ! Work space ALLOC_(CS%dzRegrid(G%isd:G%ied,G%jsd:G%jed,nz+1)); CS%dzRegrid(:,:,:) = 0. @@ -488,10 +489,6 @@ subroutine ALE_memory_deallocation( CS ) call triDiagEdgeWorkDeallocate( CS%edgeValueWrk ) call triDiagSlopeWorkDeallocate( CS%edgeSlopeWrk ) - ! Piecewise polynomials - call ppoly_destroy( CS%ppoly_linear ) - call ppoly_destroy( CS%ppoly_parab ) - ! Work space DEALLOC_(CS%dzRegrid) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 16dd9705fa..873bf15f86 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -19,7 +19,6 @@ module MOM_regridding use MOM_EOS, only : calculate_density use MOM_string_functions,only : uppercase -use regrid_ppoly_class, only : ppoly_t, ppoly_init, ppoly_destroy use regrid_edge_values, only : edgeValueArrays use regrid_edge_values, only : edge_values_explicit_h2, edge_values_explicit_h4 use regrid_edge_values, only : edge_values_implicit_h4, edge_values_implicit_h6 @@ -51,17 +50,6 @@ module MOM_regridding type, public :: regridding_CS private - ! Generic linear piecewise polynomial used for various purposes throughout - ! the code (the same ppoly is used to avoid having dynamical memory allocation) - type(ppoly_t) :: ppoly_linear - - ! Generic parabolic piecewise polynomial used for various purposes - ! throughout the code (the same ppoly is used to avoid having dynamical - ! memory allocation) - type(ppoly_t) :: ppoly_parab - - type(ppoly_t) :: ppoly_i ! interpolation ppoly - type(ppoly_t) :: ppoly_r ! reconstruction ppoly ! This array is set by function setCoordinateResolution ! It contains the "resolution" or delta coordinate of the target @@ -76,6 +64,8 @@ module MOM_regridding integer :: nk ! Number of layers/levels + integer :: degree_i !Degree of interpolation polynomial + ! Indicates which grid to use in the vertical (z*, sigma, target interface ! densities) integer :: regridding_scheme @@ -398,7 +388,7 @@ subroutine buildGridZstar( CS, G, h, dzInterface ) nz = G%ke - +!$OMP parallel do default(private) shared(G,dzInterface,CS,nz,h) do j = G%jsc-1,G%jec+1 do i = G%isc-1,G%iec+1 @@ -595,7 +585,9 @@ subroutine buildGridRho( G, h, tv, dzInterface, remapCS, CS ) real :: threshold real :: max_thickness real :: correction - + real, dimension(CS%nk,2) :: ppoly_i_E !Edge value of polynomial + real, dimension(CS%nk,2) :: ppoly_i_S !Edge slope of polynomial + real, dimension(CS%nk,CS%degree_i+1) :: ppoly_i_coefficients !Coefficients of polynomial real, dimension(SZK_(G)) :: p_column, densities, T_column, S_column, Tmp_column integer, dimension(SZK_(G)) :: mapping real :: nominalDepth, totalThickness, dh @@ -692,7 +684,8 @@ subroutine buildGridRho( G, h, tv, dzInterface, remapCS, CS ) ! One regridding iteration call regridding_iteration( densities, CS%coordinateInterfaces, CS,& count_nonzero_layers, hTmp, xTmp, & - CS%ppoly_i, nz, h1, x1 ) + ppoly_i_E, ppoly_i_S, ppoly_i_coefficients, & + nz, h1, x1 ) ! Remap T and S from previous grid to new grid do k = 1,nz @@ -893,7 +886,7 @@ end subroutine build_grid_arbitrary ! Regridding iterations !------------------------------------------------------------------------------ subroutine regridding_iteration( densities, target_values, CS, & - n0, h0, x0, ppoly0, n1, h1, x1 ) + n0, h0, x0, ppoly0_E, ppoly0_S, ppoly0_coefficients, n1, h1, x1 ) ! ------------------------------------------------------------------------------ ! This routine performs one single iteration of the regridding based on target ! interface densities. Given the set of target values and cell densities, this @@ -910,7 +903,9 @@ subroutine regridding_iteration( densities, target_values, CS, & integer, intent(in) :: n0 ! Number of cells on source grid real, dimension(:), intent(in) :: h0 ! cell widths on source grid real, dimension(:), intent(in) :: x0 ! interface positions on source grid - type(ppoly_t), intent(inout) :: ppoly0 ! Piecewise polynomial for density interpolation + real, dimension(:,:),intent(inout) :: ppoly0_E !Edge value of polynomial + real, dimension(:,:),intent(inout) :: ppoly0_S !Edge slope of polynomial + real, dimension(:,:),intent(inout) :: ppoly0_coefficients !Coefficients of polynomial integer, intent(in) :: n1 ! Number of cells on target grid real, dimension(:), intent(out) :: h1 ! cell widths on target grid real, dimension(:), intent(out) :: x1 ! interface positions on target grid @@ -920,66 +915,66 @@ subroutine regridding_iteration( densities, target_values, CS, & integer :: degree, k ! Reset piecewise polynomials - ppoly0%E = 0.0 - ppoly0%S = 0.0 - ppoly0%coefficients = 0.0 + ppoly0_E = 0.0 + ppoly0_S = 0.0 + ppoly0_coefficients = 0.0 ! 1. Compute the interpolated profile of the density field and build grid select case ( CS%interpolation_scheme ) case ( INTERPOLATION_P1M_H2 ) degree = DEGREE_1 - call edge_values_explicit_h2( n0, h0, densities, ppoly0%E ) - call P1M_interpolation( n0, h0, densities, ppoly0 ) + call edge_values_explicit_h2( n0, h0, densities, ppoly0_E ) + call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call P1M_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) end if case ( INTERPOLATION_P1M_H4 ) degree = DEGREE_1 if ( n0 .ge. 4 ) then - call edge_values_explicit_h4( n0, h0, densities, ppoly0%E ) + call edge_values_explicit_h4( n0, h0, densities, ppoly0_E ) else - call edge_values_explicit_h2( n0, h0, densities, ppoly0%E ) + call edge_values_explicit_h2( n0, h0, densities, ppoly0_E ) end if - call P1M_interpolation( n0, h0, densities, ppoly0 ) + call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call P1M_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) end if case ( INTERPOLATION_P1M_IH4 ) degree = DEGREE_1 if ( n0 .ge. 4 ) then - call edge_values_implicit_h4( n0, h0, densities, CS%edgeValueWrk, ppoly0%E ) + call edge_values_implicit_h4( n0, h0, densities, CS%edgeValueWrk, ppoly0_E ) else - call edge_values_explicit_h2( n0, h0, densities, ppoly0%E ) + call edge_values_explicit_h2( n0, h0, densities, ppoly0_E ) end if - call P1M_interpolation( n0, h0, densities, ppoly0 ) + call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call P1M_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) end if case ( INTERPOLATION_PLM ) degree = DEGREE_1 - call PLM_reconstruction( n0, h0, densities, ppoly0 ) + call PLM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call PLM_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call PLM_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) end if case ( INTERPOLATION_PPM_H4 ) if ( n0 .ge. 4 ) then degree = DEGREE_2 - call edge_values_explicit_h4( n0, h0, densities, ppoly0%E ) - call PPM_reconstruction( n0, h0, densities, ppoly0 ) + call edge_values_explicit_h4( n0, h0, densities, ppoly0_E ) + call PPM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call PPM_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call PPM_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) end if else degree = DEGREE_1 - call edge_values_explicit_h2( n0, h0, densities, ppoly0%E ) - call P1M_interpolation( n0, h0, densities, ppoly0 ) + call edge_values_explicit_h2( n0, h0, densities, ppoly0_E ) + call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call P1M_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) end if end if @@ -987,17 +982,17 @@ subroutine regridding_iteration( densities, target_values, CS, & if ( n0 .ge. 4 ) then degree = DEGREE_2 - call edge_values_implicit_h4( n0, h0, densities, CS%edgeValueWrk, ppoly0%E ) - call PPM_reconstruction( n0, h0, densities, ppoly0 ) + call edge_values_implicit_h4( n0, h0, densities, CS%edgeValueWrk, ppoly0_E ) + call PPM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call PPM_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call PPM_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) end if else degree = DEGREE_1 - call edge_values_explicit_h2( n0, h0, densities, ppoly0%E ) - call P1M_interpolation( n0, h0, densities, ppoly0 ) + call edge_values_explicit_h2( n0, h0, densities, ppoly0_E ) + call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call P1M_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) end if end if @@ -1005,36 +1000,36 @@ subroutine regridding_iteration( densities, target_values, CS, & if ( n0 .ge. 4 ) then degree = DEGREE_3 - call edge_values_implicit_h4( n0, h0, densities, CS%edgeValueWrk, ppoly0%E ) - call edge_slopes_implicit_h3( n0, h0, densities, CS%edgeSlopeWrk, ppoly0%S ) - call P3M_interpolation( n0, h0, densities, ppoly0 ) + call edge_values_implicit_h4( n0, h0, densities, CS%edgeValueWrk, ppoly0_E ) + call edge_slopes_implicit_h3( n0, h0, densities, CS%edgeSlopeWrk, ppoly0_S ) + call P3M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_S, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call P3M_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call P3M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_S, ppoly0_coefficients ) end if else degree = DEGREE_1 - call edge_values_explicit_h2( n0, h0, densities, ppoly0%E ) - call P1M_interpolation( n0, h0, densities, ppoly0 ) + call edge_values_explicit_h2( n0, h0, densities, ppoly0_E ) + call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call P1M_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) end if end if case ( INTERPOLATION_P3M_IH6IH5 ) if ( n0 .ge. 6 ) then degree = DEGREE_3 - call edge_values_implicit_h6( n0, h0, densities, CS%edgeValueWrk, ppoly0%E ) - call edge_slopes_implicit_h5( n0, h0, densities, CS%edgeSlopeWrk, ppoly0%S ) - call P3M_interpolation( n0, h0, densities, ppoly0 ) + call edge_values_implicit_h6( n0, h0, densities, CS%edgeValueWrk, ppoly0_E ) + call edge_slopes_implicit_h5( n0, h0, densities, CS%edgeSlopeWrk, ppoly0_S ) + call P3M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_S, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call P3M_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call P3M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_S, ppoly0_coefficients ) end if else degree = DEGREE_1 - call edge_values_explicit_h2( n0, h0, densities, ppoly0%E ) - call P1M_interpolation( n0, h0, densities, ppoly0 ) + call edge_values_explicit_h2( n0, h0, densities, ppoly0_E ) + call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call P1M_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) end if end if @@ -1042,36 +1037,36 @@ subroutine regridding_iteration( densities, target_values, CS, & if ( n0 .ge. 4 ) then degree = DEGREE_4 - call edge_values_implicit_h4( n0, h0, densities, CS%edgeValueWrk, ppoly0%E ) - call edge_slopes_implicit_h3( n0, h0, densities, CS%edgeSlopeWrk, ppoly0%S ) - call PQM_reconstruction( n0, h0, densities, ppoly0 ) + call edge_values_implicit_h4( n0, h0, densities, CS%edgeValueWrk, ppoly0_E ) + call edge_slopes_implicit_h3( n0, h0, densities, CS%edgeSlopeWrk, ppoly0_S ) + call PQM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_S, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call PQM_boundary_extrapolation_v1( n0, h0, densities, ppoly0 ) + call PQM_boundary_extrapolation_v1( n0, h0, densities, ppoly0_E, ppoly0_S, ppoly0_coefficients ) end if else degree = DEGREE_1 - call edge_values_explicit_h2( n0, h0, densities, ppoly0%E ) - call P1M_interpolation( n0, h0, densities, ppoly0 ) + call edge_values_explicit_h2( n0, h0, densities, ppoly0_E ) + call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call P1M_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) end if end if case ( INTERPOLATION_PQM_IH6IH5 ) if ( n0 .ge. 6 ) then degree = DEGREE_4 - call edge_values_implicit_h6( n0, h0, densities, CS%edgeValueWrk, ppoly0%E ) - call edge_slopes_implicit_h5( n0, h0, densities, CS%edgeSlopeWrk, ppoly0%S ) - call PQM_reconstruction( n0, h0, densities, ppoly0 ) + call edge_values_implicit_h6( n0, h0, densities, CS%edgeValueWrk, ppoly0_E ) + call edge_slopes_implicit_h5( n0, h0, densities, CS%edgeSlopeWrk, ppoly0_S ) + call PQM_reconstruction( n0, h0, densities, ppoly0_E, ppoly0_S, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call PQM_boundary_extrapolation_v1( n0, h0, densities, ppoly0 ) + call PQM_boundary_extrapolation_v1( n0, h0, densities, ppoly0_E, ppoly0_S, ppoly0_coefficients ) end if else degree = DEGREE_1 - call edge_values_explicit_h2( n0, h0, densities, ppoly0%E ) - call P1M_interpolation( n0, h0, densities, ppoly0 ) + call edge_values_explicit_h2( n0, h0, densities, ppoly0_E ) + call P1M_interpolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) if ( CS%boundary_extrapolation) then - call P1M_boundary_extrapolation( n0, h0, densities, ppoly0 ) + call P1M_boundary_extrapolation( n0, h0, densities, ppoly0_E, ppoly0_coefficients ) end if end if @@ -1079,7 +1074,7 @@ subroutine regridding_iteration( densities, target_values, CS, & ! Based on global density profile, interpolate new grid and ! inflate vanished layers - call interpolate_grid( n0, h0, x0, ppoly0, target_values, degree, n1, h1, x1 ) + call interpolate_grid( n0, h0, x0, ppoly0_E, ppoly0_coefficients, target_values, degree, n1, h1, x1 ) call inflate_vanished_layers( CS%min_thickness, n1, h1 ) x1(1) = 0.0 do k = 1,n1 @@ -1092,7 +1087,7 @@ end subroutine regridding_iteration !------------------------------------------------------------------------------ ! Given target values (e.g., density), build new grid based on polynomial !------------------------------------------------------------------------------ -subroutine interpolate_grid( n0, h0, x0, ppoly0, target_values, degree, n1, h1, x1 ) +subroutine interpolate_grid( n0, h0, x0, ppoly0_E, ppoly0_coefficients, target_values, degree, n1, h1, x1 ) ! ------------------------------------------------------------------------------ ! Given the grid 'grid0' and the piecewise polynomial interpolant ! 'ppoly0' (possibly discontinuous), the coordinates of the new grid 'grid1' @@ -1103,7 +1098,8 @@ subroutine interpolate_grid( n0, h0, x0, ppoly0, target_values, degree, n1, h1, integer, intent(in) :: n0 real, dimension(:), intent(in) :: h0 real, dimension(:), intent(in) :: x0 - type(ppoly_t), intent(in) :: ppoly0 + real, dimension(:,:), intent(in) :: ppoly0_E !Edge value of polynomial + real, dimension(:,:), intent(in) :: ppoly0_coefficients !Coefficients of polynomial real, dimension(:), intent(in) :: target_values integer, intent(in) :: degree integer, intent(in) :: n1 @@ -1122,7 +1118,7 @@ subroutine interpolate_grid( n0, h0, x0, ppoly0, target_values, degree, n1, h1, ! Find coordinates for interior target values do k = 2,n1 t = target_values(k) - x1(k) = get_polynomial_coordinate ( n0, h0, ppoly0, t, degree ) + x1(k) = get_polynomial_coordinate ( n0, h0, ppoly0_E, ppoly0_coefficients, t, degree ) h1(k-1) = x1(k) - x1(k-1) end do h1(n1) = x1(n1+1) - x1(n1) @@ -1133,7 +1129,7 @@ end subroutine interpolate_grid !------------------------------------------------------------------------------ ! Given target value, find corresponding coordinate for given polynomial !------------------------------------------------------------------------------ -real function get_polynomial_coordinate ( N, h, ppoly, target_value, degree ) +real function get_polynomial_coordinate ( N, h, ppoly_E, ppoly_coefficients, target_value, degree ) ! ------------------------------------------------------------------------------ ! Here, 'ppoly' is assumed to be a piecewise discontinuous polynomial of degree ! 'degree' throughout the domain defined by 'grid'. A target value is given @@ -1154,7 +1150,8 @@ real function get_polynomial_coordinate ( N, h, ppoly, target_value, degree ) ! Arguments integer, intent(in) :: N real, intent(in) :: h(:) ! size n - type(ppoly_t), intent(in) :: ppoly + real, intent(in) :: ppoly_E(:,:) !Edge value of polynomial + real, intent(in) :: ppoly_coefficients(:,:) !Coefficients of polynomial real, intent(in) :: target_value integer, intent(in) :: degree @@ -1181,7 +1178,7 @@ real function get_polynomial_coordinate ( N, h, ppoly, target_value, degree ) ! If the target value is outside the range of all values, we ! force the target coordinate to be equal to the lowest or ! largest value, depending on which bound is overtaken - if ( target_value .LE. ppoly%E(1,1) ) then + if ( target_value .LE. ppoly_E(1,1) ) then x = 0. ! Left boundary is at x=0 get_polynomial_coordinate = x return ! return because there is no need to look further @@ -1192,8 +1189,8 @@ real function get_polynomial_coordinate ( N, h, ppoly, target_value, degree ) x = 0. ! Left boundary is at x=0 do k = 2,N x = x + h(k-1) ! Position of interface k - if ( ( target_value .GE. ppoly%E(k-1,2) ) .AND. & - ( target_value .LE. ppoly%E(k,1) ) ) then + if ( ( target_value .GE. ppoly_E(k-1,2) ) .AND. & + ( target_value .LE. ppoly_E(k,1) ) ) then get_polynomial_coordinate = x return ! return because there is no need to look further exit @@ -1204,7 +1201,7 @@ real function get_polynomial_coordinate ( N, h, ppoly, target_value, degree ) ! force the target coordinate to be equal to the lowest or ! largest value, depending on which bound is overtaken x = x + h(N) ! Position of right boundary - if ( target_value .GE. ppoly%E(N,2) ) then + if ( target_value .GE. ppoly_E(N,2) ) then get_polynomial_coordinate = x return ! return because there is no need to look further end if @@ -1218,8 +1215,8 @@ real function get_polynomial_coordinate ( N, h, ppoly, target_value, degree ) do k = 1,N x_l = x_r x_r = x_l + h(k) - if ( ( target_value .GT. ppoly%E(k,1) ) .AND. & - ( target_value .LT. ppoly%E(k,2) ) ) then + if ( ( target_value .GT. ppoly_E(k,1) ) .AND. & + ( target_value .LT. ppoly_E(k,2) ) ) then k_found = k exit end if @@ -1231,7 +1228,7 @@ real function get_polynomial_coordinate ( N, h, ppoly, target_value, degree ) ! means there is a major problem with the interpolant. This needs to be ! reported. if ( k_found .EQ. -1 ) then - write(*,*) target_value, ppoly%E(1,1), ppoly%E(N,2) + write(*,*) target_value, ppoly_E(1,1), ppoly_E(N,2) write(*,*) 'Could not find target coordinate in ' //& '"get_polynomial_coordinate". This is caused by an '//& 'inconsistent interpolant (perhaps not monotonically '//& @@ -1243,7 +1240,7 @@ real function get_polynomial_coordinate ( N, h, ppoly, target_value, degree ) ! the found cell a(:) = 0.0 do i = 1,degree+1 - a(i) = ppoly%coefficients(k_found,i) + a(i) = ppoly_coefficients(k_found,i) end do ! Guess value to start Newton-Raphson iterations (middle of cell) @@ -1656,18 +1653,6 @@ subroutine allocate_regridding( CS ) ! Target values allocate( CS%coordinateInterfaces(CS%nk+1) ) - ! Piecewise polynomials used for remapping - call ppoly_init( CS%ppoly_r, CS%nk, 4 ) - - ! Piecewise polynomials used for interpolation - call ppoly_init( CS%ppoly_i, CS%nk, 4 ) - - ! Generic linear piecewise polynomial - call ppoly_init( CS%ppoly_linear, CS%nk, 1 ) - - ! Generic parabolic piecewise polynomial - call ppoly_init( CS%ppoly_parab, CS%nk, 2 ) - ! Target resolution (for fixed coordinates) allocate( CS%coordinateResolution(CS%nk) ); CS%coordinateResolution(:) = -1.E30 @@ -1691,12 +1676,6 @@ subroutine regridding_memory_deallocation( CS ) ! Target values deallocate( CS%coordinateInterfaces ) - ! Piecewise polynomials - call ppoly_destroy( CS%ppoly_i ) - call ppoly_destroy( CS%ppoly_r ) - call ppoly_destroy( CS%ppoly_linear ) - call ppoly_destroy( CS%ppoly_parab ) - end subroutine regridding_memory_deallocation end module MOM_regridding diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 25bf00134b..c316fb8f91 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -12,7 +12,6 @@ module MOM_remapping use MOM_error_handler, only : MOM_error, FATAL use MOM_string_functions, only : uppercase use MOM_variables, only : ocean_grid_type, thermo_var_ptrs -use regrid_ppoly_class, only : ppoly_t, ppoly_init, ppoly_destroy use polynomial_functions, only : evaluation_polynomial, integration_polynomial use regrid_edge_values, only : edgeValueArrays use regrid_edge_values, only : edge_values_explicit_h4, edge_values_implicit_h4 @@ -36,13 +35,12 @@ module MOM_remapping type, public :: remapping_CS private ! Work arrays - type(ppoly_t) :: ppoly_r ! reconstruction ppoly type(edgeValueArrays) :: edgeValueWrk ! Work space for edge values type(edgeSlopeArrays) :: edgeSlopeWrk ! Work space for edge slopes ! Parameters integer :: nk = 0 ! Number of layers/levels in vertical integer :: remapping_scheme = -911 ! Determines which reconstruction to use - integer :: degree ! Degree of polynomial reconstruction + integer :: degree=0 ! Degree of polynomial reconstruction logical :: boundary_extrapolation = .true. ! If true, extrapolate boundaries end type @@ -121,6 +119,7 @@ subroutine remapping_main( CS, G, h, dxInterface, tv, u, v ) nz = G%ke ! Remap tracer +!$OMP parallel do default(shared) private(h1,dx,u_column) do j = G%jsc,G%jec do i = G%isc,G%iec if (G%mask2dT(i,j)>0.) then @@ -137,6 +136,7 @@ subroutine remapping_main( CS, G, h, dxInterface, tv, u, v ) ! Remap u velocity component if ( present(u) ) then +!$OMP parallel do default(shared) private(h1,dx,u_column) do j = G%jsc,G%jec do i = G%iscB,G%iecB if (G%mask2dCu(i,j)>0.) then @@ -152,6 +152,7 @@ subroutine remapping_main( CS, G, h, dxInterface, tv, u, v ) ! Remap v velocity component if ( present(v) ) then +!$OMP parallel do default(shared) private(h1,dx,u_column) do j = G%jscB,G%jecB do i = G%isc,G%iec if (G%mask2dCv(i,j)>0.) then @@ -374,6 +375,9 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 ) ! Local variables integer :: iMethod + real, dimension(CS%nk,2) :: ppoly_r_E !Edge value of polynomial + real, dimension(CS%nk,2) :: ppoly_r_S !Edge slope of polynomial + real, dimension(CS%nk,CS%degree+1) :: ppoly_r_coefficients !Coefficients of polynomial #ifdef __DO_SAFTEY_CHECKS__ integer :: k, remapping_scheme @@ -417,9 +421,9 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 ) iMethod = -999 ! Reset polynomial - CS%ppoly_r%E(:,:) = 0.0 - CS%ppoly_r%S(:,:) = 0.0 - CS%ppoly_r%coefficients(:,:) = 0.0 + ppoly_r_E(:,:) = 0.0 + ppoly_r_S(:,:) = 0.0 + ppoly_r_coefficients(:,:) = 0.0 remapping_scheme = CS%remapping_scheme if (n0<=1) then @@ -431,42 +435,42 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 ) endif select case ( remapping_scheme ) case ( REMAPPING_PCM ) - call PCM_reconstruction( n0, u0, CS%ppoly_r ) + call PCM_reconstruction( n0, u0, ppoly_r_E, ppoly_r_coefficients) iMethod = INTEGRATION_PCM case ( REMAPPING_PLM ) - call PLM_reconstruction( n0, h0, u0, CS%ppoly_r ) + call PLM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients ) if ( CS%boundary_extrapolation) then - call PLM_boundary_extrapolation( n0, h0, u0, CS%ppoly_r ) + call PLM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients) end if iMethod = INTEGRATION_PLM case ( REMAPPING_PPM_H4 ) - call edge_values_explicit_h4( n0, h0, u0, CS%ppoly_r%E ) - call PPM_reconstruction( n0, h0, u0, CS%ppoly_r ) + call edge_values_explicit_h4( n0, h0, u0, ppoly_r_E ) + call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients ) if ( CS%boundary_extrapolation) then - call PPM_boundary_extrapolation( n0, h0, u0, CS%ppoly_r ) + call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients ) end if iMethod = INTEGRATION_PPM case ( REMAPPING_PPM_IH4 ) - call edge_values_implicit_h4( n0, h0, u0, CS%edgeValueWrk, CS%ppoly_r%E ) - call PPM_reconstruction( n0, h0, u0, CS%ppoly_r ) + call edge_values_implicit_h4( n0, h0, u0, CS%edgeValueWrk, ppoly_r_E ) + call PPM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients ) if ( CS%boundary_extrapolation) then - call PPM_boundary_extrapolation( n0, h0, u0, CS%ppoly_r ) + call PPM_boundary_extrapolation( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients ) end if iMethod = INTEGRATION_PPM case ( REMAPPING_PQM_IH4IH3 ) - call edge_values_implicit_h4( n0, h0, u0, CS%edgeValueWrk, CS%ppoly_r%E ) - call edge_slopes_implicit_h3( n0, h0, u0, CS%edgeSlopeWrk, CS%ppoly_r%S ) - call PQM_reconstruction( n0, h0, u0, CS%ppoly_r ) + call edge_values_implicit_h4( n0, h0, u0, CS%edgeValueWrk, ppoly_r_E ) + call edge_slopes_implicit_h3( n0, h0, u0, CS%edgeSlopeWrk, ppoly_r_S ) + call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefficients ) if ( CS%boundary_extrapolation) then - call PQM_boundary_extrapolation_v1( n0, h0, u0, CS%ppoly_r ) + call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefficients ) end if iMethod = INTEGRATION_PQM case ( REMAPPING_PQM_IH6IH5 ) - call edge_values_implicit_h6( n0, h0, u0, CS%edgeValueWrk, CS%ppoly_r%E ) - call edge_slopes_implicit_h5( n0, h0, u0, CS%edgeSlopeWrk, CS%ppoly_r%S ) - call PQM_reconstruction( n0, h0, u0, CS%ppoly_r ) + call edge_values_implicit_h6( n0, h0, u0, CS%edgeValueWrk, ppoly_r_E ) + call edge_slopes_implicit_h5( n0, h0, u0, CS%edgeSlopeWrk, ppoly_r_S ) + call PQM_reconstruction( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefficients ) if ( CS%boundary_extrapolation) then - call PQM_boundary_extrapolation_v1( n0, h0, u0, CS%ppoly_r ) + call PQM_boundary_extrapolation_v1( n0, h0, u0, ppoly_r_E, ppoly_r_S, ppoly_r_coefficients ) end if iMethod = INTEGRATION_PQM case default @@ -474,7 +478,7 @@ subroutine remapping_core( CS, n0, h0, u0, n1, dx, u1 ) 'The selected remapping method is invalid' ) end select - call remapByDeltaZ( n0, h0, u0, CS%ppoly_r, n1, dx, iMethod, u1 ) + call remapByDeltaZ( n0, h0, u0, ppoly_r_E, ppoly_r_coefficients, n1, dx, iMethod, u1 ) ! call remapByProjection( n0, h0, u0, CS%ppoly_r, n1, h1, iMethod, u1 ) #ifdef __DO_SAFTEY_CHECKS__ @@ -527,12 +531,13 @@ end subroutine remapping_core ! ----------------------------------------------------------------------------- ! remapByProjection (integration of reconstructed profile) ! ----------------------------------------------------------------------------- -subroutine remapByProjection( n0, h0, u0, ppoly0, n1, h1, method, u1 ) +subroutine remapByProjection( n0, h0, u0, ppoly0_E, ppoly0_coefficients, n1, h1, method, u1 ) ! Arguments integer, intent(in) :: n0 ! number of cells in source grid real, intent(in) :: h0(:) ! source grid widths (size n0) real, intent(in) :: u0(:) ! source cell averages (size n0) - type(ppoly_t), intent(in) :: ppoly0 ! source piecewise polynomial + real, intent(in) :: ppoly0_E(:,:) !Edge value of polynomial + real, intent(in) :: ppoly0_coefficients(:,:) !Coefficients of polynomial integer, intent(in) :: n1 ! number of cells in target grid real, intent(in) :: h1(:) ! target grid widths (size n1) integer, intent(in) :: method ! remapping scheme to use @@ -551,7 +556,7 @@ subroutine remapByProjection( n0, h0, u0, ppoly0, n1, h1, method, u1 ) xL = xR xR = xL + h1(iTarget) - call integrateReconOnInterval( n0, h0, u0, ppoly0, method, & + call integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefficients, method, & xL, xR, h1(iTarget), u1(iTarget) ) end do ! end iTarget loop on target grid cells @@ -562,7 +567,7 @@ end subroutine remapByProjection ! ----------------------------------------------------------------------------- ! Remap using change in interface positions ! ----------------------------------------------------------------------------- -subroutine remapByDeltaZ( n0, h0, u0, ppoly0, n1, dx1, method, u1, h1 ) +subroutine remapByDeltaZ( n0, h0, u0, ppoly0_E, ppoly0_coefficients, n1, dx1, method, u1, h1 ) ! The new grid is defined relative to the original grid by change ! dx1(:) = xNew(:) - xOld(:) ! and the remapping calculated so that @@ -574,7 +579,8 @@ subroutine remapByDeltaZ( n0, h0, u0, ppoly0, n1, dx1, method, u1, h1 ) integer, intent(in) :: n0 ! number of cells in source grid real, intent(in) :: h0(:) ! source grid widths (size n0) real, intent(in) :: u0(:) ! source cell averages (size n0) - type(ppoly_t), intent(in) :: ppoly0 ! source piecewise polynomial + real, intent(in) :: ppoly0_E(:,:) !Edge value of polynomial + real, intent(in) :: ppoly0_coefficients(:,:) !Coefficients of polynomial integer, intent(in) :: n1 ! number of cells in target grid real, intent(in) :: dx1(:) ! target grid edge positions (size n1+1) integer :: method ! remapping scheme to use @@ -638,7 +644,7 @@ subroutine remapByDeltaZ( n0, h0, u0, ppoly0, n1, dx1, method, u1, h1 ) ! hFlux is the positive width of the remapped volume hFlux = abs(dx1(iTarget+1)) - call integrateReconOnInterval( n0, h0, u0, ppoly0, method, & + call integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefficients, method, & xL, xR, hFlux, uAve ) ! uAve is the average value of u, independent of sign of dx1 fluxR = dx1(iTarget+1)*uAve ! Includes sign of dx1 @@ -677,17 +683,18 @@ end subroutine remapByDeltaZ ! ----------------------------------------------------------------------------- ! integrate the reconstructed profile over a single cell ! ----------------------------------------------------------------------------- -subroutine integrateReconOnInterval( n0, h0, u0, ppoly0, method, & +subroutine integrateReconOnInterval( n0, h0, u0, ppoly0_E, ppoly0_coefficients, method, & xL, xR, hC, uAve ) ! Arguments - integer, intent(in) :: n0 ! number of cells in source grid - real, intent(in) :: h0(:) ! source grid sizes (size n0) - real, dimension(:), intent(in) :: u0 ! source cell averages - type(ppoly_t), intent(in) :: ppoly0 ! source piecewise polynomial - integer, intent(in) :: method ! remapping scheme to use - real, intent(in) :: xL, xR ! left/right edges of target cell - real, intent(in) :: hC ! cell width hC = xR - xL - real, intent(out) :: uAve ! average value on target cell + integer, intent(in) :: n0 ! number of cells in source grid + real, intent(in) :: h0(:) ! source grid sizes (size n0) + real, intent(in) :: u0(:) ! source cell averages + real, intent(in) :: ppoly0_E(:,:) !Edge value of polynomial + real, intent(in) :: ppoly0_coefficients(:,:) !Coefficients of polynomial + integer, intent(in) :: method ! remapping scheme to use + real, intent(in) :: xL, xR ! left/right edges of target cell + real, intent(in) :: hC ! cell width hC = xR - xL + real, intent(out) :: uAve ! average value on target cell ! Local variables integer :: j, k @@ -761,20 +768,20 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0, method, & ! value is set to be mean of the edge values (which should be the same). ! If it isn't, we simply interpolate. if ( h0(jL) == 0.0 ) then - uAve = 0.5 * ( ppoly0%E(jL,1) + ppoly0%E(jL,2) ) + uAve = 0.5 * ( ppoly0_E(jL,1) + ppoly0_E(jL,2) ) else ! WHY IS THIS NOT WRITTEN AS xi0 = ( xL - x0jLl ) / h0(jL) ---AJA xi0 = xL / h0(jL) - x0jLl / h0(jL) select case ( method ) case ( INTEGRATION_PCM ) - uAve = ppoly0%coefficients(jL,1) + uAve = ppoly0_coefficients(jL,1) case ( INTEGRATION_PLM ) - uAve = evaluation_polynomial( ppoly0%coefficients(jL,:), 2, xi0 ) + uAve = evaluation_polynomial( ppoly0_coefficients(jL,:), 2, xi0 ) case ( INTEGRATION_PPM ) - uAve = evaluation_polynomial( ppoly0%coefficients(jL,:), 3, xi0 ) + uAve = evaluation_polynomial( ppoly0_coefficients(jL,:), 3, xi0 ) case ( INTEGRATION_PQM ) - uAve = evaluation_polynomial( ppoly0%coefficients(jL,:), 5, xi0 ) + uAve = evaluation_polynomial( ppoly0_coefficients(jL,:), 5, xi0 ) case default call MOM_error( FATAL,'The selected integration method is invalid' ) end select @@ -829,16 +836,16 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0, method, & ! coordinates, hence: \int_xL^xR p(x) dx = h \int_xi0^xi1 p(xi) dxi select case ( method ) case ( INTEGRATION_PCM ) - q = ppoly0%coefficients(jL,1) * ( xR - xL ) + q = ppoly0_coefficients(jL,1) * ( xR - xL ) case ( INTEGRATION_PLM ) q = h0(jL) * & - integration_polynomial( xi0, xi1, ppoly0%coefficients(jL,:), 1 ) + integration_polynomial( xi0, xi1, ppoly0_coefficients(jL,:), 1 ) case ( INTEGRATION_PPM ) q = h0(jL) * & - integration_polynomial( xi0, xi1, ppoly0%coefficients(jL,:), 2 ) + integration_polynomial( xi0, xi1, ppoly0_coefficients(jL,:), 2 ) case ( INTEGRATION_PQM ) q = h0(jL) * & - integration_polynomial( xi0, xi1, ppoly0%coefficients(jL,:), 4 ) + integration_polynomial( xi0, xi1, ppoly0_coefficients(jL,:), 4 ) case default call MOM_error( FATAL,'The selected integration method is invalid' ) end select @@ -864,16 +871,16 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0, method, & xi1 = 1.0 select case ( method ) case ( INTEGRATION_PCM ) - q = q + ppoly0%coefficients(jL,1) * ( x0jLr - xL ) + q = q + ppoly0_coefficients(jL,1) * ( x0jLr - xL ) case ( INTEGRATION_PLM ) q = q + h0(jL) * & - integration_polynomial( xi0, xi1, ppoly0%coefficients(jL,:), 1 ) + integration_polynomial( xi0, xi1, ppoly0_coefficients(jL,:), 1 ) case ( INTEGRATION_PPM ) q = q + h0(jL) * & - integration_polynomial( xi0, xi1, ppoly0%coefficients(jL,:), 2 ) + integration_polynomial( xi0, xi1, ppoly0_coefficients(jL,:), 2 ) case ( INTEGRATION_PQM ) q = q + h0(jL) * & - integration_polynomial( xi0, xi1, ppoly0%coefficients(jL,:), 4 ) + integration_polynomial( xi0, xi1, ppoly0_coefficients(jL,:), 4 ) case default call MOM_error( FATAL, 'The selected integration method is invalid' ) end select @@ -892,16 +899,16 @@ subroutine integrateReconOnInterval( n0, h0, u0, ppoly0, method, & select case ( method ) case ( INTEGRATION_PCM ) - q = q + ppoly0%coefficients(jR,1) * ( xR - x0jRl ) + q = q + ppoly0_coefficients(jR,1) * ( xR - x0jRl ) case ( INTEGRATION_PLM ) q = q + h0(jR) * & - integration_polynomial( xi0, xi1, ppoly0%coefficients(jR,:), 1 ) + integration_polynomial( xi0, xi1, ppoly0_coefficients(jR,:), 1 ) case ( INTEGRATION_PPM ) q = q + h0(jR) * & - integration_polynomial( xi0, xi1, ppoly0%coefficients(jR,:), 2 ) + integration_polynomial( xi0, xi1, ppoly0_coefficients(jR,:), 2 ) case ( INTEGRATION_PQM ) q = q + h0(jR) * & - integration_polynomial( xi0, xi1, ppoly0%coefficients(jR,:), 4 ) + integration_polynomial( xi0, xi1, ppoly0_coefficients(jR,:), 4 ) case default call MOM_error( FATAL,'The selected integration method is invalid' ) end select @@ -1019,15 +1026,16 @@ subroutine setReconstructionType(string,CS) "Unrecognized choice for REMAPPING_SCHEME ("//trim(string)//").") end select - if (allocated(CS%ppoly_r%E) .and. degree/=CS%degree) then + if (CS%degree > 0 .and. degree/=CS%degree) then ! If the degree has changed then deallocate to force a re-allocation call end_remapping(CS) endif - CS%degree = degree - if (.not. allocated(CS%ppoly_r%E)) then + if ( CS%degree == 0 ) then call allocate_remapping( CS ) endif + CS%degree = degree + end subroutine setReconstructionType !------------------------------------------------------------------------------ @@ -1054,8 +1062,7 @@ end subroutine remapDisableBoundaryExtrapolation subroutine allocate_remapping( CS ) ! Arguments type(remapping_CS), intent(inout) :: CS - - call ppoly_init( CS%ppoly_r, CS%nk, CS%degree ) + call triDiagEdgeWorkAllocate( CS%nk, CS%edgeValueWrk ) call triDiagSlopeWorkAllocate( CS%nk, CS%edgeSlopeWrk ) @@ -1070,9 +1077,9 @@ subroutine end_remapping(CS) type(remapping_CS), intent(inout) :: CS ! Deallocate memory for grid - call ppoly_destroy( CS%ppoly_r ) call triDiagEdgeWorkDeallocate( CS%edgeValueWrk ) call triDiagSlopeWorkDeallocate( CS%edgeSlopeWrk ) + CS%degree = 0 end subroutine end_remapping @@ -1088,7 +1095,7 @@ logical function remappingUnitTests() data h1 /3*1./ ! 3 uniform layers with total depth of 3 data h2 /6*0.5/ ! 6 uniform layers with total depth of 3 type(remapping_CS) :: CS - type(ppoly_t) :: ppoly0 + real, allocatable, dimension(:,:) :: ppoly0_E, ppoly0_S, ppoly0_coefficients integer :: i real :: err @@ -1119,16 +1126,19 @@ logical function remappingUnitTests() write(*,*) 'h1 (by projection)' call dumpGrid(n1,h1,x1,u1) - call ppoly_init( ppoly0, n0, CS%degree ) - ppoly0%E(:,:) = 0.0 - ppoly0%S(:,:) = 0.0 - ppoly0%coefficients(:,:) = 0.0 + allocate(ppoly0_E(n0,2)) + allocate(ppoly0_S(n0,2)) + allocate(ppoly0_coefficients(n0,CS%degree+1)) + + ppoly0_E(:,:) = 0.0 + ppoly0_S(:,:) = 0.0 + ppoly0_coefficients(:,:) = 0.0 - call edge_values_explicit_h4( n0, h0, u0, ppoly0%E ) - call PPM_reconstruction( n0, h0, u0, ppoly0 ) - call PPM_boundary_extrapolation( n0, h0, u0, ppoly0 ) + call edge_values_explicit_h4( n0, h0, u0, ppoly0_E ) + call PPM_reconstruction( n0, h0, u0, ppoly0_E, ppoly0_coefficients ) + call PPM_boundary_extrapolation( n0, h0, u0, ppoly0_E, ppoly0_coefficients ) u1(:) = 0. - call remapByProjection( n0, h0, u0, ppoly0, & + call remapByProjection( n0, h0, u0, ppoly0_E, ppoly0_coefficients, & n1, h1, INTEGRATION_PPM, u1 ) do i=1,n1 err=u1(i)-8./3.*(0.5*real(1+n1)-real(i)) @@ -1136,7 +1146,7 @@ logical function remappingUnitTests() enddo u1(:) = 0. - call remapByDeltaZ( n0, h0, u0, ppoly0, & + call remapByDeltaZ( n0, h0, u0, ppoly0_E, ppoly0_coefficients, & n1, x1-x0(1:n1+1), & INTEGRATION_PPM, u1, hn1 ) write(*,*) 'h1 (by delta)' @@ -1150,7 +1160,7 @@ logical function remappingUnitTests() call buildGridFromH(n2, h2, x2) dx2(1:n0+1) = x2(1:n0+1) - x0 dx2(n0+2:n2+1) = x2(n0+2:n2+1) - x0(n0+1) - call remapByDeltaZ( n0, h0, u0, ppoly0, & + call remapByDeltaZ( n0, h0, u0, ppoly0_E, ppoly0_coefficients, & n2, dx2, & INTEGRATION_PPM, u2, hn2 ) write(*,*) 'h2' @@ -1163,6 +1173,8 @@ logical function remappingUnitTests() if (abs(err)>2.*epsilon(err)) remappingUnitTests = .true. enddo + deallocate(ppoly0_E, ppoly0_S, ppoly0_coefficients) + write(*,*) '==========================================================' contains diff --git a/src/ALE/P1M_functions.F90 b/src/ALE/P1M_functions.F90 index a15939b7d5..6450cb0cf6 100644 --- a/src/ALE/P1M_functions.F90 +++ b/src/ALE/P1M_functions.F90 @@ -25,7 +25,6 @@ module P1M_functions ! P1M_boundary_extrapolation (public) ! !============================================================================== -use regrid_ppoly_class, only : ppoly_t use regrid_edge_values, only : bound_edge_values, average_discontinuous_edge_values implicit none ; private @@ -41,7 +40,7 @@ module P1M_functions !------------------------------------------------------------------------------ ! p1m interpolation !------------------------------------------------------------------------------ -subroutine P1M_interpolation( N, h, u, ppoly ) +subroutine P1M_interpolation( N, h, u, ppoly_E, ppoly_coefficients ) ! ------------------------------------------------------------------------------ ! Linearly interpolate between edge values. ! The resulting piecewise interpolant is stored in 'ppoly'. @@ -60,27 +59,29 @@ subroutine P1M_interpolation( N, h, u, ppoly ) integer, intent(in) :: N ! Number of cells real, dimension(:), intent(in) :: h ! cell widths (size N) real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly + real, dimension(:,:), intent(inout) :: ppoly_E + real, dimension(:,:), intent(inout) :: ppoly_coefficients + ! Local variables integer :: k ! loop index real :: u0_l, u0_r ! edge values (left and right) ! Bound edge values (routine found in 'edge_values.F90') - call bound_edge_values( N, h, u, ppoly%E ) + call bound_edge_values( N, h, u, ppoly_E ) ! Systematically average discontinuous edge values (routine found in ! 'edge_values.F90') - call average_discontinuous_edge_values( N, ppoly%E ) + call average_discontinuous_edge_values( N, ppoly_E ) ! Loop on interior cells to build interpolants do k = 1,N - u0_l = ppoly%E(k,1) - u0_r = ppoly%E(k,2) + u0_l = ppoly_E(k,1) + u0_r = ppoly_E(k,2) - ppoly%coefficients(k,1) = u0_l - ppoly%coefficients(k,2) = u0_r - u0_l + ppoly_coefficients(k,1) = u0_l + ppoly_coefficients(k,2) = u0_r - u0_l end do ! end loop on interior cells @@ -90,7 +91,7 @@ end subroutine P1M_interpolation !------------------------------------------------------------------------------ ! p1m boundary extrapolation ! ----------------------------------------------------------------------------- -subroutine P1M_boundary_extrapolation( N, h, u, ppoly ) +subroutine P1M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) !------------------------------------------------------------------------------ ! Interpolation by linear polynomials within boundary cells. ! The left and right edge values in the left and right boundary cells, @@ -108,7 +109,8 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly ) integer, intent(in) :: N ! Number of cells real, dimension(:), intent(in) :: h ! cell widths (size N) real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly + real, dimension(:,:), intent(inout) :: ppoly_E + real, dimension(:,:), intent(inout) :: ppoly_coefficients ! Local variables real :: u0, u1 ! cell averages @@ -136,20 +138,20 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly ) ! by using the edge value in the neighboring cell. u0_r = u0 + 0.5 * slope - if ( (u1 - u0) * (ppoly%E(2,1) - u0_r) .LT. 0.0 ) then - slope = 2.0 * ( ppoly%E(2,1) - u0 ) + if ( (u1 - u0) * (ppoly_E(2,1) - u0_r) .LT. 0.0 ) then + slope = 2.0 * ( ppoly_E(2,1) - u0 ) end if ! Using the limited slope, the left edge value is reevaluated and ! the interpolant coefficients recomputed if ( h0 .NE. 0.0 ) then - ppoly%E(1,1) = u0 - 0.5 * slope + ppoly_E(1,1) = u0 - 0.5 * slope else - ppoly%E(1,1) = u0 + ppoly_E(1,1) = u0 end if - ppoly%coefficients(1,1) = ppoly%E(1,1) - ppoly%coefficients(1,2) = ppoly%E(1,2) - ppoly%E(1,1) + ppoly_coefficients(1,1) = ppoly_E(1,1) + ppoly_coefficients(1,2) = ppoly_E(1,2) - ppoly_E(1,1) ! ------------------------------------------ ! Right edge value in the left boundary cell @@ -164,18 +166,18 @@ subroutine P1M_boundary_extrapolation( N, h, u, ppoly ) u0_l = u1 - 0.5 * slope - if ( (u1 - u0) * (u0_l - ppoly%E(N-1,2)) .LT. 0.0 ) then - slope = 2.0 * ( u1 - ppoly%E(N-1,2) ) + if ( (u1 - u0) * (u0_l - ppoly_E(N-1,2)) .LT. 0.0 ) then + slope = 2.0 * ( u1 - ppoly_E(N-1,2) ) end if if ( h1 .NE. 0.0 ) then - ppoly%E(N,2) = u1 + 0.5 * slope + ppoly_E(N,2) = u1 + 0.5 * slope else - ppoly%E(N,2) = u1 + ppoly_E(N,2) = u1 end if - ppoly%coefficients(N,1) = ppoly%E(N,1) - ppoly%coefficients(N,2) = ppoly%E(N,2) - ppoly%E(N,1) + ppoly_coefficients(N,1) = ppoly_E(N,1) + ppoly_coefficients(N,2) = ppoly_E(N,2) - ppoly_E(N,1) end subroutine P1M_boundary_extrapolation diff --git a/src/ALE/P3M_functions.F90 b/src/ALE/P3M_functions.F90 index 74120ec65f..4b9d21871a 100644 --- a/src/ALE/P3M_functions.F90 +++ b/src/ALE/P3M_functions.F90 @@ -14,7 +14,6 @@ module P3M_functions ! cubic curve. ! !============================================================================== -use regrid_ppoly_class, only : ppoly_t use regrid_edge_values, only : bound_edge_values, average_discontinuous_edge_values implicit none ; private @@ -27,7 +26,7 @@ module P3M_functions !------------------------------------------------------------------------------ ! p3m interpolation ! ----------------------------------------------------------------------------- -subroutine P3M_interpolation( N, h, u, ppoly ) +subroutine P3M_interpolation( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients ) !------------------------------------------------------------------------------ ! Cubic interpolation between edges. ! @@ -39,10 +38,13 @@ subroutine P3M_interpolation( N, h, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: h ! cell widths (size N) - real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: h ! cell widths (size N) + real, dimension(:), intent(in) :: u ! cell averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E !Edge value of polynomial + real, dimension(:,:), intent(inout) :: ppoly_S !Edge slope of polynomial + real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial + ! Call the limiter for p3m, which takes care of everything from ! computing the coefficients of the cubic to monotonizing it. @@ -50,7 +52,7 @@ subroutine P3M_interpolation( N, h, u, ppoly ) ! 'P3M_interpolation' first but we do that to provide an homogeneous ! interface. - call P3M_limiter( N, h, u, ppoly ) + call P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients ) end subroutine P3M_interpolation @@ -58,7 +60,7 @@ end subroutine P3M_interpolation !------------------------------------------------------------------------------ ! p3m limiter ! ----------------------------------------------------------------------------- -subroutine P3M_limiter( N, h, u, ppoly ) +subroutine P3M_limiter( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients ) !------------------------------------------------------------------------------ ! The p3m limiter operates as follows: ! @@ -73,10 +75,14 @@ subroutine P3M_limiter( N, h, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: h ! cell widths (size N) - real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: h ! cell widths (size N) + real, dimension(:), intent(in) :: u ! cell averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E !Edge value of polynomial + real, dimension(:,:), intent(inout) :: ppoly_S !Edge slope of polynomial + real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial + +! real, dimension(:,:), intent(inout) :: ppoly_coefficients ! Local variables integer :: k ! loop index @@ -93,10 +99,10 @@ subroutine P3M_limiter( N, h, u, ppoly ) eps = 1e-10 ! 1. Bound edge values (boundary cells are assumed to be local extrema) - call bound_edge_values( N, h, u, ppoly%E ) + call bound_edge_values( N, h, u, ppoly_E ) ! 2. Systematically average discontinuous edge values - call average_discontinuous_edge_values( N, ppoly%E ) + call average_discontinuous_edge_values( N, ppoly_E ) ! 3. Loop on cells and do the following @@ -106,10 +112,10 @@ subroutine P3M_limiter( N, h, u, ppoly ) do k = 1,N ! Get edge values, edge slopes and cell width - u0_l = ppoly%E(k,1) - u0_r = ppoly%E(k,2) - u1_l = ppoly%S(k,1) - u1_r = ppoly%S(k,2) + u0_l = ppoly_E(k,1) + u0_r = ppoly_E(k,2) + u1_l = ppoly_S(k,1) + u1_r = ppoly_S(k,2) ! Get cell widths and cell averages (boundary cells are assumed to ! be local extrema for the sake of slopes) @@ -165,10 +171,10 @@ subroutine P3M_limiter( N, h, u, ppoly ) end if ! Build cubic interpolant (compute the coefficients) - call build_cubic_interpolant( h, k, ppoly ) + call build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coefficients ) ! Check whether cubic is monotonic - monotonic = is_cubic_monotonic( ppoly, k ) + monotonic = is_cubic_monotonic( ppoly_coefficients, k ) ! If cubic is not monotonic, monotonize it by modifiying the ! edge slopes, store the new edge slopes and recompute the @@ -178,11 +184,11 @@ subroutine P3M_limiter( N, h, u, ppoly ) end if ! Store edge slopes - ppoly%S(k,1) = u1_l - ppoly%S(k,2) = u1_r + ppoly_S(k,1) = u1_l + ppoly_S(k,2) = u1_r ! Recompute coefficients of cubic - call build_cubic_interpolant( h, k, ppoly ) + call build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coefficients ) end do ! loop on cells @@ -192,7 +198,7 @@ end subroutine P3M_limiter !------------------------------------------------------------------------------ ! p3m boundary extrapolation ! ----------------------------------------------------------------------------- -subroutine P3M_boundary_extrapolation( N, h, u, ppoly ) +subroutine P3M_boundary_extrapolation( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients ) !------------------------------------------------------------------------------ ! The following explanations apply to the left boundary cell. The same ! reasoning holds for the right boundary cell. @@ -207,10 +213,12 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: h ! cell widths (size N) - real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: h ! cell widths (size N) + real, dimension(:), intent(in) :: u ! cell averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E !Edge value of polynomial + real, dimension(:,:), intent(inout) :: ppoly_S !Edge slope of polynomial + real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial ! Local variables integer :: i0, i1 @@ -235,7 +243,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly ) ! Compute the left edge slope in neighboring cell and express it in ! the global coordinate system - b = ppoly%coefficients(i1,2) + b = ppoly_coefficients(i1,2) u1_r = b / h1 ! derivative evaluated at xi = 0.0, expressed w.r.t. x ! Limit the right slope by the PLM limited slope @@ -246,7 +254,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly ) ! The right edge value in the boundary cell is taken to be the left ! edge value in the neighboring cell - u0_r = ppoly%E(i1,1) + u0_r = ppoly_E(i1,1) ! Given the right edge value and slope, we determine the left ! edge value and slope by computing the parabola as determined by @@ -264,22 +272,22 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly ) end if ! Store edge values and slope, build cubic and check monotonicity - ppoly%E(i0,1) = u0_l - ppoly%E(i0,2) = u0_r - ppoly%S(i0,1) = u1_l - ppoly%S(i0,2) = u1_r + ppoly_E(i0,1) = u0_l + ppoly_E(i0,2) = u0_r + ppoly_S(i0,1) = u1_l + ppoly_S(i0,2) = u1_r ! Store edge values and slope, build cubic and check monotonicity - call build_cubic_interpolant( h, i0, ppoly ) - monotonic = is_cubic_monotonic( ppoly, i0 ) + call build_cubic_interpolant( h, i0, ppoly_E, ppoly_S, ppoly_coefficients ) + monotonic = is_cubic_monotonic( ppoly_coefficients, i0 ) if ( monotonic .EQ. 0 ) then call monotonize_cubic( h0, u0_l, u0_r, 0.0, slope, slope, u1_l, u1_r ) ! Rebuild cubic after monotonization - ppoly%S(i0,1) = u1_l - ppoly%S(i0,2) = u1_r - call build_cubic_interpolant( h, i0, ppoly ) + ppoly_S(i0,1) = u1_l + ppoly_S(i0,2) = u1_r + call build_cubic_interpolant( h, i0, ppoly_E, ppoly_S, ppoly_coefficients ) end if @@ -293,9 +301,9 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly ) ! Compute the right edge slope in neighboring cell and express it in ! the global coordinate system - b = ppoly%coefficients(i0,2) - c = ppoly%coefficients(i0,3) - d = ppoly%coefficients(i0,4) + b = ppoly_coefficients(i0,2) + c = ppoly_coefficients(i0,3) + d = ppoly_coefficients(i0,4) u1_l = (b + 2*c + 3*d) / h0 ! derivative evaluated at xi = 1.0 ! Limit the left slope by the PLM limited slope @@ -306,7 +314,7 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly ) ! The left edge value in the boundary cell is taken to be the right ! edge value in the neighboring cell - u0_l = ppoly%E(i0,2) + u0_l = ppoly_E(i0,2) ! Given the left edge value and slope, we determine the right ! edge value and slope by computing the parabola as determined by @@ -324,21 +332,21 @@ subroutine P3M_boundary_extrapolation( N, h, u, ppoly ) end if ! Store edge values and slope, build cubic and check monotonicity - ppoly%E(i1,1) = u0_l - ppoly%E(i1,2) = u0_r - ppoly%S(i1,1) = u1_l - ppoly%S(i1,2) = u1_r + ppoly_E(i1,1) = u0_l + ppoly_E(i1,2) = u0_r + ppoly_S(i1,1) = u1_l + ppoly_S(i1,2) = u1_r - call build_cubic_interpolant( h, i1, ppoly ) - monotonic = is_cubic_monotonic( ppoly, i1 ) + call build_cubic_interpolant( h, i1, ppoly_E, ppoly_S, ppoly_coefficients ) + monotonic = is_cubic_monotonic( ppoly_coefficients, i1 ) if ( monotonic .EQ. 0 ) then call monotonize_cubic( h1, u0_l, u0_r, slope, 0.0, slope, u1_l, u1_r ) ! Rebuild cubic after monotonization - ppoly%S(i1,1) = u1_l - ppoly%S(i1,2) = u1_r - call build_cubic_interpolant( h, i1, ppoly ) + ppoly_S(i1,1) = u1_l + ppoly_S(i1,2) = u1_r + call build_cubic_interpolant( h, i1, ppoly_E, ppoly_S, ppoly_coefficients ) end if @@ -348,7 +356,7 @@ end subroutine P3M_boundary_extrapolation !------------------------------------------------------------------------------ ! Build cubic interpolant in cell k ! ----------------------------------------------------------------------------- -subroutine build_cubic_interpolant( h, k, ppoly ) +subroutine build_cubic_interpolant( h, k, ppoly_E, ppoly_S, ppoly_coefficients ) !------------------------------------------------------------------------------ ! Given edge values and edge slopes, compute coefficients of cubic in cell k. ! @@ -357,9 +365,11 @@ subroutine build_cubic_interpolant( h, k, ppoly ) !------------------------------------------------------------------------------ ! Arguments - real, dimension(:), intent(in) :: h ! cell widths (size N) - integer, intent(in) :: k - type(ppoly_t), intent(inout) :: ppoly + real, dimension(:), intent(in) :: h ! cell widths (size N) + integer, intent(in) :: k + real, dimension(:,:), intent(in) :: ppoly_E !Edge value of polynomial + real, dimension(:,:), intent(in) :: ppoly_S !Edge slope of polynomial + real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial ! Local variables real :: u0_l, u0_r ! edge values @@ -369,21 +379,21 @@ subroutine build_cubic_interpolant( h, k, ppoly ) h_c = h(k) - u0_l = ppoly%E(k,1) - u0_r = ppoly%E(k,2) + u0_l = ppoly_E(k,1) + u0_r = ppoly_E(k,2) - u1_l = ppoly%S(k,1) * h_c - u1_r = ppoly%S(k,2) * h_c + u1_l = ppoly_S(k,1) * h_c + u1_r = ppoly_S(k,2) * h_c a0 = u0_l a1 = u1_l a2 = 3.0 * ( u0_r - u0_l ) - u1_r - 2.0 * u1_l a3 = u1_r + u1_l + 2.0 * ( u0_l - u0_r ) - ppoly%coefficients(k,1) = a0 - ppoly%coefficients(k,2) = a1 - ppoly%coefficients(k,3) = a2 - ppoly%coefficients(k,4) = a3 + ppoly_coefficients(k,1) = a0 + ppoly_coefficients(k,2) = a1 + ppoly_coefficients(k,3) = a2 + ppoly_coefficients(k,4) = a3 end subroutine build_cubic_interpolant @@ -391,7 +401,7 @@ end subroutine build_cubic_interpolant !------------------------------------------------------------------------------ ! Check whether cubic is monotonic ! ----------------------------------------------------------------------------- -integer function is_cubic_monotonic( ppoly, k ) +integer function is_cubic_monotonic( ppoly_coefficients, k ) !------------------------------------------------------------------------------ ! This function checks whether the cubic curve in cell k is monotonic. ! If so, returns 1. Otherwise, returns 0. @@ -402,8 +412,8 @@ integer function is_cubic_monotonic( ppoly, k ) !------------------------------------------------------------------------------ ! Arguments - type(ppoly_t), intent(in) :: ppoly - integer, intent(in) :: k + real, dimension(:,:), intent(in) :: ppoly_coefficients + integer, intent(in) :: k ! Local variables integer :: monotonic ! boolean indicating if monotonic or not @@ -417,10 +427,10 @@ integer function is_cubic_monotonic( ppoly, k ) ! to be equal to 0 or 1, respectively eps = 1e-14 - a0 = ppoly%coefficients(k,1) - a1 = ppoly%coefficients(k,2) - a2 = ppoly%coefficients(k,3) - a3 = ppoly%coefficients(k,4) + a0 = ppoly_coefficients(k,1) + a1 = ppoly_coefficients(k,2) + a2 = ppoly_coefficients(k,3) + a3 = ppoly_coefficients(k,4) a = a1 b = 2.0 * a2 diff --git a/src/ALE/PCM_functions.F90 b/src/ALE/PCM_functions.F90 index f0056f234f..ebade668f9 100644 --- a/src/ALE/PCM_functions.F90 +++ b/src/ALE/PCM_functions.F90 @@ -10,7 +10,6 @@ module PCM_functions ! reconstruction using the piecewise constant method (PCM). ! !============================================================================== -use regrid_ppoly_class, only : ppoly_t implicit none ; private @@ -21,7 +20,7 @@ module PCM_functions !------------------------------------------------------------------------------ ! pcm_reconstruction !------------------------------------------------------------------------------ -subroutine PCM_reconstruction( N, u, ppoly ) +subroutine PCM_reconstruction( N, u, ppoly_E, ppoly_coefficients ) !------------------------------------------------------------------------------ ! Reconstruction by constant polynomials within each cell. There is nothing to ! do but this routine is provided to ensure a homogeneous interface @@ -37,20 +36,21 @@ subroutine PCM_reconstruction( N, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: u ! cell averages - type(ppoly_t), intent(inout) :: ppoly + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: u ! cell averages + real, dimension(:,:), intent(inout) :: ppoly_E !Edge value of polynomial + real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial ! Local variables integer :: k ! The coefficients of the piecewise constant polynomial are simply ! the cell averages. - ppoly%coefficients(:,1) = u(:) + ppoly_coefficients(:,1) = u(:) ! The edge values are equal to the cell average do k = 1,N - ppoly%E(k,:) = u(k) + ppoly_E(k,:) = u(k) end do end subroutine PCM_reconstruction diff --git a/src/ALE/PLM_functions.F90 b/src/ALE/PLM_functions.F90 index 30a91a7985..1f686cac89 100644 --- a/src/ALE/PLM_functions.F90 +++ b/src/ALE/PLM_functions.F90 @@ -10,7 +10,6 @@ module PLM_functions ! reconstruction using the piecewise linear method (PLM). ! !============================================================================== -use regrid_ppoly_class, only : ppoly_t implicit none ; private @@ -21,7 +20,7 @@ module PLM_functions !------------------------------------------------------------------------------ ! PLM_reconstruction ! ----------------------------------------------------------------------------- -subroutine PLM_reconstruction( N, h, u, ppoly ) +subroutine PLM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients ) !------------------------------------------------------------------------------ ! Reconstruction by linear polynomials within each cell. ! @@ -34,10 +33,11 @@ subroutine PLM_reconstruction( N, h, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: h ! cell widths (size N) - real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: h ! cell widths (size N) + real, dimension(:), intent(in) :: u ! cell averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E + real, dimension(:,:), intent(inout) :: ppoly_coefficients ! Local variables integer :: k ! loop index @@ -82,10 +82,10 @@ subroutine PLM_reconstruction( N, h, u, ppoly ) a = u_c - 0.5 * slope b = slope - ppoly%coefficients(k,1) = a - ppoly%coefficients(k,2) = b - ppoly%E(k,1) = a - ppoly%E(k,2) = a+b + ppoly_coefficients(k,1) = a + ppoly_coefficients(k,2) = b + ppoly_E(k,1) = a + ppoly_E(k,2) = a+b end do ! end loop on interior cells @@ -98,31 +98,31 @@ subroutine PLM_reconstruction( N, h, u, ppoly ) slope = 0.0 a = u(1) - 0.5 * slope b = slope - ppoly%coefficients(1,1) = a - ppoly%coefficients(1,2) = b - ppoly%E(1,1) = a - ppoly%E(1,2) = a+b + ppoly_coefficients(1,1) = a + ppoly_coefficients(1,2) = b + ppoly_E(1,1) = a + ppoly_E(1,2) = a+b ! Right (bottom) boundary cell slope = u(N) - u(N-1) slope = 0.0 a = u(N) - 0.5 * slope b = slope - ppoly%coefficients(N,1) = a - ppoly%coefficients(N,2) = b - ppoly%E(N,1) = a - ppoly%E(N,2) = a+b + ppoly_coefficients(N,1) = a + ppoly_coefficients(N,2) = b + ppoly_E(N,1) = a + ppoly_E(N,2) = a+b ! Second pass: we need to check for nonmonotonic discontinuous edge values. ! When this occurs, the PLM slope is redefined so as to ensure monotonic edge ! values across edges. allocate( E_old(N,2) ) - E_old = ppoly%E + E_old = ppoly_E do k = 2,N-1 ! By default, the right and left slopes are both equal to the original slope - slope = ppoly%coefficients(k,2) + slope = ppoly_coefficients(k,2) sigma_l = slope sigma_r = slope @@ -147,10 +147,10 @@ subroutine PLM_reconstruction( N, h, u, ppoly ) a = u(k) - 0.5 * slope b = slope - ppoly%coefficients(k,1) = a - ppoly%coefficients(k,2) = b - ppoly%E(k,1) = a - ppoly%E(k,2) = a+b + ppoly_coefficients(k,1) = a + ppoly_coefficients(k,2) = b + ppoly_E(k,1) = a + ppoly_E(k,2) = a+b end do ! end loop on interior cells @@ -162,7 +162,7 @@ end subroutine PLM_reconstruction !------------------------------------------------------------------------------ ! plm boundary extrapolation ! ----------------------------------------------------------------------------- -subroutine PLM_boundary_extrapolation( N, h, u, ppoly ) +subroutine PLM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) !------------------------------------------------------------------------------ ! Reconstruction by linear polynomials within boundary cells. ! The left and right edge values in the left and right boundary cells, @@ -179,10 +179,11 @@ subroutine PLM_boundary_extrapolation( N, h, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: h ! cell widths (size N) - real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: h ! cell widths (size N) + real, dimension(:), intent(in) :: u ! cell averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E + real, dimension(:,:), intent(inout) :: ppoly_coefficients ! Local variables real :: u0, u1 ! cell averages @@ -199,17 +200,17 @@ subroutine PLM_boundary_extrapolation( N, h, u, ppoly ) u1 = u(2) ! The h2 scheme is used to compute the right edge value - ppoly%E(1,2) = (u0*h1 + u1*h0) / (h0 + h1) + ppoly_E(1,2) = (u0*h1 + u1*h0) / (h0 + h1) ! The standard PLM slope is computed as a first estimate for the ! reconstruction within the cell - slope = 2.0 * ( ppoly%E(1,2) - u0 ) + slope = 2.0 * ( ppoly_E(1,2) - u0 ) - ppoly%E(1,1) = u0 - 0.5 * slope - ppoly%E(1,2) = u0 + 0.5 * slope + ppoly_E(1,1) = u0 - 0.5 * slope + ppoly_E(1,2) = u0 + 0.5 * slope - ppoly%coefficients(1,1) = ppoly%E(1,1) - ppoly%coefficients(1,2) = ppoly%E(1,2) - ppoly%E(1,1) + ppoly_coefficients(1,1) = ppoly_E(1,1) + ppoly_coefficients(1,2) = ppoly_E(1,2) - ppoly_E(1,1) ! ------------------------------------------ ! Right edge value in the left boundary cell @@ -221,17 +222,17 @@ subroutine PLM_boundary_extrapolation( N, h, u, ppoly ) u1 = u(N) ! The h2 scheme is used to compute the right edge value - ppoly%E(N,1) = (u0*h1 + u1*h0) / (h0 + h1) + ppoly_E(N,1) = (u0*h1 + u1*h0) / (h0 + h1) ! The standard PLM slope is computed as a first estimate for the ! reconstruction within the cell - slope = 2.0 * ( u1 - ppoly%E(N,1) ) + slope = 2.0 * ( u1 - ppoly_E(N,1) ) - ppoly%E(N,1) = u1 - 0.5 * slope - ppoly%E(N,2) = u1 + 0.5 * slope + ppoly_E(N,1) = u1 - 0.5 * slope + ppoly_E(N,2) = u1 + 0.5 * slope - ppoly%coefficients(N,1) = ppoly%E(N,1) - ppoly%coefficients(N,2) = ppoly%E(N,2) - ppoly%E(N,1) + ppoly_coefficients(N,1) = ppoly_E(N,1) + ppoly_coefficients(N,2) = ppoly_E(N,2) - ppoly_E(N,1) end subroutine PLM_boundary_extrapolation diff --git a/src/ALE/PPM_functions.F90 b/src/ALE/PPM_functions.F90 index 90338503a9..31ec80a787 100644 --- a/src/ALE/PPM_functions.F90 +++ b/src/ALE/PPM_functions.F90 @@ -10,7 +10,6 @@ module PPM_functions ! reconstruction using the piecewise parabolic method (PPM). ! !============================================================================== -use regrid_ppoly_class, only : ppoly_t use regrid_edge_values, only : bound_edge_values, check_discontinuous_edge_values implicit none ; private @@ -22,7 +21,7 @@ module PPM_functions !------------------------------------------------------------------------------ ! PPM_reconstruction ! ----------------------------------------------------------------------------- -subroutine PPM_reconstruction( N, h, u, ppoly ) +subroutine PPM_reconstruction( N, h, u, ppoly_E, ppoly_coefficients ) !------------------------------------------------------------------------------ ! Reconstruction by quadratic polynomials within each cell. ! @@ -37,10 +36,11 @@ subroutine PPM_reconstruction( N, h, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: h ! cell widths (size N) - real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: h ! cell widths (size N) + real, dimension(:), intent(in) :: u ! cell averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E + real, dimension(:,:), intent(inout) :: ppoly_coefficients ! Local variables integer :: k ! loop index @@ -48,22 +48,22 @@ subroutine PPM_reconstruction( N, h, u, ppoly ) real :: a, b, c ! parabola coefficients ! PPM limiter - call PPM_limiter_standard( N, h, u, ppoly ) + call PPM_limiter_standard( N, h, u, ppoly_E ) ! Loop on cells to construct the parabola within each cell do k = 1,N - u0_l = ppoly%E(k,1) - u0_r = ppoly%E(k,2) + u0_l = ppoly_E(k,1) + u0_r = ppoly_E(k,2) a = u0_l b = 6.0 * u(k) - 4.0 * u0_l - 2.0 * u0_r c = 3.0 * ( u0_r + u0_l - 2.0 * u(k) ) ! Store coefficients - ppoly%coefficients(k,1) = a - ppoly%coefficients(k,2) = b - ppoly%coefficients(k,3) = c + ppoly_coefficients(k,1) = a + ppoly_coefficients(k,2) = b + ppoly_coefficients(k,3) = c end do ! end loop on interior cells @@ -73,7 +73,7 @@ end subroutine PPM_reconstruction !------------------------------------------------------------------------------ ! Limit ppm ! ----------------------------------------------------------------------------- -subroutine PPM_limiter_standard( N, h, u, ppoly ) +subroutine PPM_limiter_standard( N, h, u, ppoly_E ) !------------------------------------------------------------------------------ ! Standard PPM limiter (Colella & Woodward, JCP 1984). ! @@ -86,10 +86,10 @@ subroutine PPM_limiter_standard( N, h, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: h ! cell widths (size N) - real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: h ! cell widths (size N) + real, dimension(:), intent(in) :: u ! cell averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E ! Local variables integer :: k ! loop index @@ -98,10 +98,10 @@ subroutine PPM_limiter_standard( N, h, u, ppoly ) real :: expr1, expr2 ! Bound edge values - call bound_edge_values( N, h, u, ppoly%E ) + call bound_edge_values( N, h, u, ppoly_E ) ! Make discontinuous edge values monotonic - call check_discontinuous_edge_values( N, u, ppoly%E ) + call check_discontinuous_edge_values( N, u, ppoly_E ) ! Loop on interior cells to apply the standard ! PPM limiter (Colella & Woodward, JCP 84) @@ -112,8 +112,8 @@ subroutine PPM_limiter_standard( N, h, u, ppoly ) u_c = u(k) u_r = u(k+1) - u0_l = ppoly%E(k,1) - u0_r = ppoly%E(k,2) + u0_l = ppoly_E(k,1) + u0_r = ppoly_E(k,2) ! Auxiliary variables expr1 = (u0_r - u0_l) * (u_c - 0.5*(u0_l+u0_r)) @@ -133,14 +133,14 @@ subroutine PPM_limiter_standard( N, h, u, ppoly ) u0_r = 3.0 * u_c - 2.0 * u0_l end if - ppoly%E(k,1) = u0_l - ppoly%E(k,2) = u0_r + ppoly_E(k,1) = u0_l + ppoly_E(k,2) = u0_r end do ! end loop on interior cells ! PCM within boundary cells - ppoly%E(1,:) = u(1) - ppoly%E(N,:) = u(N) + ppoly_E(1,:) = u(1) + ppoly_E(N,:) = u(N) end subroutine PPM_limiter_standard @@ -148,7 +148,7 @@ end subroutine PPM_limiter_standard !------------------------------------------------------------------------------ ! ppm boundary extrapolation ! ----------------------------------------------------------------------------- -subroutine PPM_boundary_extrapolation( N, h, u, ppoly ) +subroutine PPM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) !------------------------------------------------------------------------------ ! Reconstruction by parabolas within boundary cells. ! @@ -172,10 +172,11 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: h ! cell widths (size N) - real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: h ! cell widths (size N) + real, dimension(:), intent(in) :: u ! cell averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E + real, dimension(:,:), intent(inout) :: ppoly_coefficients ! Local variables integer :: i0, i1 @@ -197,7 +198,7 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly ) ! Compute the left edge slope in neighboring cell and express it in ! the global coordinate system - b = ppoly%coefficients(i1,2) + b = ppoly_coefficients(i1,2) u1_r = b *(h0/h1) ! derivative evaluated at xi = 0.0, ! expressed w.r.t. xi (local coord. system) @@ -209,7 +210,7 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly ) ! The right edge value in the boundary cell is taken to be the left ! edge value in the neighboring cell - u0_r = ppoly%E(i1,1) + u0_r = ppoly_E(i1,1) ! Given the right edge value and slope, we determine the left ! edge value and slope by computing the parabola as determined by @@ -228,16 +229,16 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly ) u0_r = 3.0 * u0 - 2.0 * u0_l end if - ppoly%E(i0,1) = u0_l - ppoly%E(i0,2) = u0_r + ppoly_E(i0,1) = u0_l + ppoly_E(i0,2) = u0_r a = u0_l b = 6.0 * u0 - 4.0 * u0_l - 2.0 * u0_r c = 3.0 * ( u0_r + u0_l - 2.0 * u0 ) - ppoly%coefficients(i0,1) = a - ppoly%coefficients(i0,2) = b - ppoly%coefficients(i0,3) = c + ppoly_coefficients(i0,1) = a + ppoly_coefficients(i0,2) = b + ppoly_coefficients(i0,3) = c ! ----- Right boundary ----- i0 = N-1 @@ -249,8 +250,8 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly ) ! Compute the right edge slope in neighboring cell and express it in ! the global coordinate system - b = ppoly%coefficients(i0,2) - c = ppoly%coefficients(i0,3) + b = ppoly_coefficients(i0,2) + c = ppoly_coefficients(i0,3) u1_l = (b + 2*c) ! derivative evaluated at xi = 1.0 u1_l = u1_l * (h1/h0) @@ -262,7 +263,7 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly ) ! The left edge value in the boundary cell is taken to be the right ! edge value in the neighboring cell - u0_l = ppoly%E(i0,2) + u0_l = ppoly_E(i0,2) ! Given the left edge value and slope, we determine the right ! edge value and slope by computing the parabola as determined by @@ -281,16 +282,16 @@ subroutine PPM_boundary_extrapolation( N, h, u, ppoly ) u0_r = 3.0 * u1 - 2.0 * u0_l end if - ppoly%E(i1,1) = u0_l - ppoly%E(i1,2) = u0_r + ppoly_E(i1,1) = u0_l + ppoly_E(i1,2) = u0_r a = u0_l b = 6.0 * u1 - 4.0 * u0_l - 2.0 * u0_r c = 3.0 * ( u0_r + u0_l - 2.0 * u1 ) - ppoly%coefficients(i1,1) = a - ppoly%coefficients(i1,2) = b - ppoly%coefficients(i1,3) = c + ppoly_coefficients(i1,1) = a + ppoly_coefficients(i1,2) = b + ppoly_coefficients(i1,3) = c end subroutine PPM_boundary_extrapolation diff --git a/src/ALE/PQM_functions.F90 b/src/ALE/PQM_functions.F90 index b76fe81bb7..48e75f4fb2 100644 --- a/src/ALE/PQM_functions.F90 +++ b/src/ALE/PQM_functions.F90 @@ -10,7 +10,6 @@ module PQM_functions ! reconstruction using the piecewise quartic method (PQM). ! !============================================================================== -use regrid_ppoly_class, only : ppoly_t use regrid_edge_values, only : bound_edge_values, check_discontinuous_edge_values implicit none ; private @@ -22,7 +21,7 @@ module PQM_functions !------------------------------------------------------------------------------ ! PQM_reconstruction ! ----------------------------------------------------------------------------- -subroutine PQM_reconstruction( N, h, u, ppoly ) +subroutine PQM_reconstruction( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients ) !------------------------------------------------------------------------------ ! Reconstruction by quartic polynomials within each cell. ! @@ -35,11 +34,13 @@ subroutine PQM_reconstruction( N, h, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: h ! cell widths (size N) - real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly - + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: h ! cell widths (size N) + real, dimension(:), intent(in) :: u ! cell averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E !Edge value of polynomial + real, dimension(:,:), intent(inout) :: ppoly_S !Edge slope of polynomial + real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial + ! Local variables integer :: k ! loop index real :: h_c ! cell width @@ -48,16 +49,16 @@ subroutine PQM_reconstruction( N, h, u, ppoly ) real :: a, b, c, d, e ! parabola coefficients ! PQM limiter - call PQM_limiter( N, h, u, ppoly ) + call PQM_limiter( N, h, u, ppoly_E, ppoly_S ) ! Loop on cells to construct the cubic within each cell do k = 1,N - u0_l = ppoly%E(k,1) - u0_r = ppoly%E(k,2) + u0_l = ppoly_E(k,1) + u0_r = ppoly_E(k,2) - u1_l = ppoly%S(k,1) - u1_r = ppoly%S(k,2) + u1_l = ppoly_S(k,1) + u1_r = ppoly_S(k,2) h_c = h(k) @@ -68,11 +69,11 @@ subroutine PQM_reconstruction( N, h, u, ppoly ) e = 30.0 * u(k) + 2.5*h_c*(u1_r - u1_l) - 15.0*(u0_l + u0_r) ! Store coefficients - ppoly%coefficients(k,1) = a - ppoly%coefficients(k,2) = b - ppoly%coefficients(k,3) = c - ppoly%coefficients(k,4) = d - ppoly%coefficients(k,5) = e + ppoly_coefficients(k,1) = a + ppoly_coefficients(k,2) = b + ppoly_coefficients(k,3) = c + ppoly_coefficients(k,4) = d + ppoly_coefficients(k,5) = e end do ! end loop on cells @@ -82,7 +83,7 @@ end subroutine PQM_reconstruction !------------------------------------------------------------------------------ ! Limit pqm ! ----------------------------------------------------------------------------- -subroutine PQM_limiter( N, h, u, ppoly ) +subroutine PQM_limiter( N, h, u, ppoly_E, ppoly_S ) !------------------------------------------------------------------------------ ! Standard PQM limiter (White & Adcroft, JCP 2008). ! @@ -95,10 +96,11 @@ subroutine PQM_limiter( N, h, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: h ! cell widths (size N) - real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: h ! cell widths (size N) + real, dimension(:), intent(in) :: u ! cell averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E !Edge value of polynomial + real, dimension(:,:), intent(inout) :: ppoly_S !Edge slope of polynomial ! Local variables integer :: k ! loop index @@ -118,10 +120,10 @@ subroutine PQM_limiter( N, h, u, ppoly ) real :: x1, x2 ! Bound edge values - call bound_edge_values( N, h, u, ppoly%E ) + call bound_edge_values( N, h, u, ppoly_E ) ! Make discontinuous edge values monotonic (thru averaging) - call check_discontinuous_edge_values( N, u, ppoly%E ) + call check_discontinuous_edge_values( N, u, ppoly_E ) ! Loop on interior cells to apply the PQM limiter do k = 2,N-1 @@ -132,10 +134,10 @@ subroutine PQM_limiter( N, h, u, ppoly ) inflexion_r = 0 ! Get edge values, edge slopes and cell width - u0_l = ppoly%E(k,1) - u0_r = ppoly%E(k,2) - u1_l = ppoly%S(k,1) - u1_r = ppoly%S(k,2) + u0_l = ppoly_E(k,1) + u0_r = ppoly_E(k,2) + u1_l = ppoly_S(k,1) + u1_r = ppoly_S(k,2) ! Get cell widths and cell averages (boundary cells are assumed to ! be local extrema for the sake of slopes) @@ -336,19 +338,19 @@ subroutine PQM_limiter( N, h, u, ppoly ) end if ! clause to check where to collapse inflexion points ! Save edge values and edge slopes for reconstruction - ppoly%E(k,1) = u0_l - ppoly%E(k,2) = u0_r - ppoly%S(k,1) = u1_l - ppoly%S(k,2) = u1_r + ppoly_E(k,1) = u0_l + ppoly_E(k,2) = u0_r + ppoly_S(k,1) = u1_l + ppoly_S(k,2) = u1_r end do ! end loop on interior cells ! Constant reconstruction within boundary cells - ppoly%E(1,:) = u(1) - ppoly%S(1,:) = 0.0 + ppoly_E(1,:) = u(1) + ppoly_S(1,:) = 0.0 - ppoly%E(N,:) = u(N) - ppoly%S(N,:) = 0.0 + ppoly_E(N,:) = u(N) + ppoly_S(N,:) = 0.0 end subroutine PQM_limiter @@ -356,7 +358,7 @@ end subroutine PQM_limiter !------------------------------------------------------------------------------ ! pqm boundary extrapolation ! ----------------------------------------------------------------------------- -subroutine PQM_boundary_extrapolation( N, h, u, ppoly ) +subroutine PQM_boundary_extrapolation( N, h, u, ppoly_E, ppoly_coefficients ) !------------------------------------------------------------------------------ ! Reconstruction by parabolas within boundary cells. ! @@ -380,11 +382,12 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: h ! cell widths (size N) - real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly - + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: h ! cell widths (size N) + real, dimension(:), intent(in) :: u ! cell averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E !Edge value of polynomial + real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial + ! Local variables integer :: i0, i1 real :: u0, u1 @@ -405,7 +408,7 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly ) ! Compute the left edge slope in neighboring cell and express it in ! the global coordinate system - b = ppoly%coefficients(i1,2) + b = ppoly_coefficients(i1,2) u1_r = b *(h0/h1) ! derivative evaluated at xi = 0.0, ! expressed w.r.t. xi (local coord. system) @@ -417,7 +420,7 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly ) ! The right edge value in the boundary cell is taken to be the left ! edge value in the neighboring cell - u0_r = ppoly%E(i1,1) + u0_r = ppoly_E(i1,1) ! Given the right edge value and slope, we determine the left ! edge value and slope by computing the parabola as determined by @@ -436,19 +439,19 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly ) u0_r = 3.0 * u0 - 2.0 * u0_l end if - ppoly%E(i0,1) = u0_l - ppoly%E(i0,2) = u0_r + ppoly_E(i0,1) = u0_l + ppoly_E(i0,2) = u0_r a = u0_l b = 6.0 * u0 - 4.0 * u0_l - 2.0 * u0_r c = 3.0 * ( u0_r + u0_l - 2.0 * u0 ) ! The quartic is reduced to a parabola in the boundary cell - ppoly%coefficients(i0,1) = a - ppoly%coefficients(i0,2) = b - ppoly%coefficients(i0,3) = c - ppoly%coefficients(i0,4) = 0.0 - ppoly%coefficients(i0,5) = 0.0 + ppoly_coefficients(i0,1) = a + ppoly_coefficients(i0,2) = b + ppoly_coefficients(i0,3) = c + ppoly_coefficients(i0,4) = 0.0 + ppoly_coefficients(i0,5) = 0.0 ! ----- Right boundary ----- i0 = N-1 @@ -460,10 +463,10 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly ) ! Compute the right edge slope in neighboring cell and express it in ! the global coordinate system - b = ppoly%coefficients(i0,2) - c = ppoly%coefficients(i0,3) - d = ppoly%coefficients(i0,4) - e = ppoly%coefficients(i0,5) + b = ppoly_coefficients(i0,2) + c = ppoly_coefficients(i0,3) + d = ppoly_coefficients(i0,4) + e = ppoly_coefficients(i0,5) u1_l = (b + 2*c + 3*d + 4*e) ! derivative evaluated at xi = 1.0 u1_l = u1_l * (h1/h0) @@ -475,7 +478,7 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly ) ! The left edge value in the boundary cell is taken to be the right ! edge value in the neighboring cell - u0_l = ppoly%E(i0,2) + u0_l = ppoly_E(i0,2) ! Given the left edge value and slope, we determine the right ! edge value and slope by computing the parabola as determined by @@ -494,19 +497,19 @@ subroutine PQM_boundary_extrapolation( N, h, u, ppoly ) u0_r = 3.0 * u1 - 2.0 * u0_l end if - ppoly%E(i1,1) = u0_l - ppoly%E(i1,2) = u0_r + ppoly_E(i1,1) = u0_l + ppoly_E(i1,2) = u0_r a = u0_l b = 6.0 * u1 - 4.0 * u0_l - 2.0 * u0_r c = 3.0 * ( u0_r + u0_l - 2.0 * u1 ) ! The quartic is reduced to a parabola in the boundary cell - ppoly%coefficients(i1,1) = a - ppoly%coefficients(i1,2) = b - ppoly%coefficients(i1,3) = c - ppoly%coefficients(i1,4) = 0.0 - ppoly%coefficients(i1,5) = 0.0 + ppoly_coefficients(i1,1) = a + ppoly_coefficients(i1,2) = b + ppoly_coefficients(i1,3) = c + ppoly_coefficients(i1,4) = 0.0 + ppoly_coefficients(i1,5) = 0.0 end subroutine PQM_boundary_extrapolation @@ -514,7 +517,7 @@ end subroutine PQM_boundary_extrapolation !------------------------------------------------------------------------------ ! pqm boundary extrapolation using rational function ! ----------------------------------------------------------------------------- -subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly ) +subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly_E, ppoly_S, ppoly_coefficients ) !------------------------------------------------------------------------------ ! Reconstruction by parabolas within boundary cells. ! @@ -538,11 +541,13 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly ) !------------------------------------------------------------------------------ ! Arguments - integer, intent(in) :: N ! Number of cells - real, dimension(:), intent(in) :: h ! cell widths (size N) - real, dimension(:), intent(in) :: u ! cell averages (size N) - type(ppoly_t), intent(inout) :: ppoly - + integer, intent(in) :: N ! Number of cells + real, dimension(:), intent(in) :: h ! cell widths (size N) + real, dimension(:), intent(in) :: u ! cell averages (size N) + real, dimension(:,:), intent(inout) :: ppoly_E !Edge value of polynomial + real, dimension(:,:), intent(inout) :: ppoly_S !Edge slope of polynomial + real, dimension(:,:), intent(inout) :: ppoly_coefficients !Coefficients of polynomial + ! Local variables integer :: i0, i1 integer :: inflexion_l @@ -576,8 +581,8 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly ) ! The right edge value and slope of the boundary cell are taken to be the ! left edge value and slope of the adjacent cell - a = ppoly%coefficients(i1,1) - b = ppoly%coefficients(i1,2) + a = ppoly_coefficients(i1,1) + b = ppoly_coefficients(i1,2) u0_r = a ! edge value u1_r = b / h1 ! edge slope (w.r.t. global coord.) @@ -689,10 +694,10 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly ) end if ! Store edge values, edge slopes and coefficients - ppoly%E(i0,1) = u0_l - ppoly%E(i0,2) = u0_r - ppoly%S(i0,1) = u1_l - ppoly%S(i0,2) = u1_r + ppoly_E(i0,1) = u0_l + ppoly_E(i0,2) = u0_r + ppoly_S(i0,1) = u1_l + ppoly_S(i0,2) = u1_r a = u0_l b = h0 * u1_l @@ -701,11 +706,11 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly ) e = 30.0 * um + 2.5*h0*(u1_r - u1_l) - 15.0*(u0_l + u0_r) ! Store coefficients - ppoly%coefficients(i0,1) = a - ppoly%coefficients(i0,2) = b - ppoly%coefficients(i0,3) = c - ppoly%coefficients(i0,4) = d - ppoly%coefficients(i0,5) = e + ppoly_coefficients(i0,1) = a + ppoly_coefficients(i0,2) = b + ppoly_coefficients(i0,3) = c + ppoly_coefficients(i0,4) = d + ppoly_coefficients(i0,5) = e ! ----- Right boundary (BOTTOM) ----- i0 = N-1 @@ -723,11 +728,11 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly ) ! The left edge value and slope of the boundary cell are taken to be the ! right edge value and slope of the adjacent cell - a = ppoly%coefficients(i0,1) - b = ppoly%coefficients(i0,2) - c = ppoly%coefficients(i0,3) - d = ppoly%coefficients(i0,4) - e = ppoly%coefficients(i0,5) + a = ppoly_coefficients(i0,1) + b = ppoly_coefficients(i0,2) + c = ppoly_coefficients(i0,3) + d = ppoly_coefficients(i0,4) + e = ppoly_coefficients(i0,5) u0_l = a + b + c + d + e ! edge value u1_l = (b + 2*c + 3*d + 4*e) / h0 ! edge slope (w.r.t. global coord.) @@ -842,10 +847,10 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly ) end if ! Store edge values, edge slopes and coefficients - ppoly%E(i1,1) = u0_l - ppoly%E(i1,2) = u0_r - ppoly%S(i1,1) = u1_l - ppoly%S(i1,2) = u1_r + ppoly_E(i1,1) = u0_l + ppoly_E(i1,2) = u0_r + ppoly_S(i1,1) = u1_l + ppoly_S(i1,2) = u1_r a = u0_l b = h1 * u1_l @@ -853,11 +858,11 @@ subroutine PQM_boundary_extrapolation_v1( N, h, u, ppoly ) d = -60.0 * um + h1 *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l e = 30.0 * um + 2.5*h1*(u1_r - u1_l) - 15.0*(u0_l + u0_r) - ppoly%coefficients(i1,1) = a - ppoly%coefficients(i1,2) = b - ppoly%coefficients(i1,3) = c - ppoly%coefficients(i1,4) = d - ppoly%coefficients(i1,5) = e + ppoly_coefficients(i1,1) = a + ppoly_coefficients(i1,2) = b + ppoly_coefficients(i1,3) = c + ppoly_coefficients(i1,4) = d + ppoly_coefficients(i1,5) = e end subroutine PQM_boundary_extrapolation_v1