WIP remapping class

src/ALE/MOM_remapping.F90

@@ -18,6 +18,8 @@ module MOM_remapping
 use PQM_functions, only : PQM_reconstruction, PQM_boundary_extrapolation_v1
 use MOM_hybgen_remap, only : hybgen_plm_coefs, hybgen_ppm_coefs, hybgen_weno_coefs
 
+use Recon1d_PCM, only : PCM
+
 implicit none ; private
 
 !> Container for remapping parameters
@@ -1624,6 +1626,7 @@ logical function remapping_unit_tests(verbose)
   type(testing) :: test ! Unit testing convenience functions
   integer :: om4
   character(len=4) :: om4_tag
+  type(PCM) :: PCM
 
   call test%set( verbose=verbose ) ! Sets the verbosity flag in test
 
@@ -2257,6 +2260,9 @@ logical function remapping_unit_tests(verbose)
                      3, (/0.,0.,0./), (/0.,0.,0./), &
                      3, (/0.,0.,0./), (/0.,0.,0./) )
 
+  if (verbose) write(test%stdout,*) '- - - - - - - - - - Recon1d PCM tests  - - - - - - - - -'
+  call test%test( PCM%unit_tests(verbose, test%stdout, test%stderr), 'PCM unit test')
+
   remapping_unit_tests = test%summarize('remapping_unit_tests')
 
 end function remapping_unit_tests

src/ALE/Recon1d_PCM.F90

@@ -0,0 +1,163 @@
+!> The equation of state using the Jackett and McDougall fits to the UNESCO EOS
+module Recon1d_PCM
+
+! This file is part of MOM6. See LICENSE.md for the license.
+
+use Recon1d_type, only : Recon1d, testing
+
+implicit none ; private
+
+public PCM
+
+!> The EOS_base implementation of the UNESCO equation of state
+type, extends (Recon1d) :: PCM
+
+  real, allocatable :: h(:) !< Grid spacing (thickness) [typically H]
+
+contains
+  !> Implementation of the PCM constructor
+  procedure :: construct => construct_PCM
+  !> Implementation of the PCM reconstruction
+  procedure :: reconstruct => PCM_reconstruct
+  !> Implementation of function returning the PCM edge values
+  procedure :: lr_edge => PCM_lr_edge
+  !> Implementation of the PCM average over an interval [A]
+  procedure :: average => PCM_average
+  !> Implementation of finding the PCM position of a value
+  procedure :: x_for_f => PCM_x_for_f
+  !> Implementation of unit tests for the PCM reconstruction
+  procedure :: unit_tests => PCM_unit_tests
+
+end type PCM
+
+contains
+
+!> Initialize a 1D PCM reconstruction for n cells
+subroutine construct_PCM(this, n)
+  class(PCM), intent(out) :: this !< This reconstruction
+  integer,    intent(in)  :: n    !< Number of cells in this column
+
+  this%n = n
+
+  allocate( this%u_mean(n) )
+  allocate( this%h(n) )
+
+  if (this%do_legacy) then
+    this%degree = 1
+    allocate( this%poly_coef(n,1) )
+  endif
+
+end subroutine construct_PCM
+
+!> Calculate a 1D PCM reconstructions based on h(:) and u(:)
+subroutine PCM_reconstruct(this, h, u)
+  class(PCM), intent(inout) :: this !< This reconstruction
+  real,       intent(in)    :: h(*) !< Grid spacing (thickness) [typically H]
+  real,       intent(in)    :: u(*) !< Cell mean values [A]
+  ! Local variables
+  integer :: k
+
+  do k = 1, this%n
+    this%u_mean(k) = u(k)
+    this%h(k) = h(k)
+  enddo
+
+  if (this%do_legacy) then
+    do k = 1, this%n
+      this%poly_coef(k,1) = u(k)
+    enddo
+  endif
+
+end subroutine PCM_reconstruct
+
+!> Returns the left/right edge values in cell k of a 1D PCM reconstruction
+subroutine PCM_lr_edge(this, k, u_left, u_right)
+  class(PCM), intent(in)  :: this    !< This reconstruction
+  integer,    intent(in)  :: k       !< Cell number
+  real,       intent(out) :: u_left  !< Left edge value [A]
+  real,       intent(out) :: u_right !< Right edge value [A]
+
+  u_left = this%u_mean(k)
+  u_right = this%u_mean(k)
+
+end subroutine PCM_lr_edge
+
+!> Position at which 1D PCM reconstruction has a value f in cell k [0..1]
+real function PCM_x_for_f(this, k, f)
+  class(PCM), intent(in) :: this !< This reconstruction
+  integer,    intent(in) :: k    !< Cell number
+  real,       intent(in) :: f    !< Value of reconstruction to solve for [A]
+
+  PCM_x_for_f = 0.5 ! For PCM, every value between x=0 and x=1 has the same value
+                    ! but x=0.5 is the average position
+  if ( f < this%u_mean(k) ) then
+    PCM_x_for_f = 0.
+  elseif ( f > this%u_mean(k) ) then
+    PCM_x_for_f = 1.
+  endif
+
+end function PCM_x_for_f
+
+!> Average between xa and xb for cell k of a 1D PCM reconstruction [A]
+real function PCM_average(this, k, xa, xb)
+  class(PCM), intent(in) :: this !< This reconstruction
+  integer,    intent(in) :: k    !< Cell number
+  real,       intent(in) :: xa   !< Start of averageing interval on element (0 to 1)
+  real,       intent(in) :: xb   !< End of averageing interval on element (0 to 1)
+
+  PCM_average = this%u_mean(k)
+
+end function PCM_average
+
+!> Runs PCM reconstruction unit tests and returns True for any fails, False otherwise
+logical function PCM_unit_tests(this, verbose, stdout, stderr)
+  class(PCM), intent(inout) :: this    !< This reconstruction
+  logical,    intent(in)    :: verbose !< True, if verbose
+  integer,    intent(in)    :: stdout  !< I/O channel for stdout
+  integer,    intent(in)    :: stderr  !< I/O channel for stderr
+  ! Local variables
+  real, allocatable :: ul(:), ur(:), um(:) ! test values [A]
+  type(testing) :: test ! convenience functions
+  integer :: k
+
+  call test%set( verbose=verbose ) ! Sets the verbosity flag in test
+
+  call this%construct(3)
+  call test%test( this%n /= 3, 'Setting number of levels')
+  allocate( um(3), ul(3), ur(3) )
+
+  call this%reconstruct( (/2.,2.,2./), (/1.,3.,5./) )
+  call test%real_arr(3, this%u_mean, (/1.,3.,5./), 'Setting cell values')
+  do k = 1, 3
+    um(k) = this%cell_mean(k)
+  enddo
+  call test%real_arr(3, um, (/1.,3.,5./), 'Return cell mean')
+
+  do k = 1, 3
+    call this%lr_edge(k, ul(k), ur(k))
+  enddo
+  call test%real_arr(3, ul, (/1.,3.,5./), 'Return left edge')
+  call test%real_arr(3, ur, (/1.,3.,5./), 'Return right edge')
+
+  do k = 1, 3
+    um(k) = this%average(k, 0.5, 0.75)
+  enddo
+  call test%real_arr(3, um, (/1.,3.,5./), 'Return interval average')
+
+  do k = 1, 3
+    ul(k) = this%x_for_f(k, real(2*k-2))
+    um(k) = this%x_for_f(k, real(2*k-1))
+    ur(k) = this%x_for_f(k, real(2*k-0))
+  enddo
+  call test%real_arr(3, ul, (/0.,0.,0./), 'Return position of f<')
+  call test%real_arr(3, um, (/0.5,0.5,0.5/), 'Return position of f=')
+  call test%real_arr(3, ur, (/1.,1.,1./), 'Return position of f>')
+
+  PCM_unit_tests = test%summarize('PCM_unit_tests')
+
+end function PCM_unit_tests
+
+!> \namespace recon1d_pcm
+!!
+
+end module Recon1d_PCM