Merge remote-tracking branch 'gfdl/dev/gfdl' into merge_gfdl

.testing/Makefile

@@ -265,10 +265,12 @@ $(BUILD)/timing/Makefile: MOM_ACFLAGS += --with-driver=timing_tests
 
 
 # Build executables
+.NOTPARALLEL:$(foreach e,$(UNIT_EXECS),$(BUILD)/unit/$(e))
 $(BUILD)/unit/test_%: $(BUILD)/unit/Makefile FORCE
 	cd $(@D) && $(TIME) $(MAKE) $(@F) -j
 $(BUILD)/unit/Makefile: $(foreach e,$(UNIT_EXECS),../config_src/drivers/unit_tests/$(e).F90)
 
+.NOTPARALLEL:$(foreach e,$(TIMING_EXECS),$(BUILD)/timing/$(e))
 $(BUILD)/timing/time_%: $(BUILD)/timing/Makefile FORCE
 	cd $(@D) && $(TIME) $(MAKE) $(@F) -j
 $(BUILD)/timing/Makefile: $(foreach e,$(TIMING_EXECS),../config_src/drivers/timing_tests/$(e).F90)

config_src/drivers/timing_tests/time_MOM_remapping.F90

@@ -0,0 +1,100 @@
+program time_MOM_remapping
+
+! This file is part of MOM6. See LICENSE.md for the license.
+
+use MOM_remapping, only : remapping_CS
+use MOM_remapping, only : initialize_remapping
+use MOM_remapping, only : remapping_core_h
+
+implicit none
+
+type(remapping_CS) :: CS
+integer, parameter :: nk=75, nij=20*20, nits=10, nsamp=100, nschemes = 2
+character(len=10) :: scheme_labels(nschemes)
+real, dimension(nschemes) :: timings ! Time for nits of nij calls for each scheme [s]
+real, dimension(nschemes) :: tmean ! Mean time for a call [s]
+real, dimension(nschemes) :: tstd ! Standard deviation of time for a call [s]
+real, dimension(nschemes) :: tmin ! Shortest time for a call [s]
+real, dimension(nschemes) :: tmax ! Longest time for a call [s]
+real, dimension(:,:), allocatable :: u0, u1 ! Source/target values [arbitrary but same units as each other]
+real, dimension(:,:), allocatable :: h0, h1 ! Source target thicknesses [0..1]
+real :: start, finish ! Times [s]
+real :: h0sum, h1sum ! Totals of h0 and h1 [nondim]
+integer :: ij, k, isamp, iter, ischeme ! Indices and counters
+integer :: seed_size ! Number of integers used by seed
+integer, allocatable :: seed(:) ! Random number seed
+
+! Set seed for random numbers
+call random_seed(size=seed_size)
+allocate( seed(seed_Size) )
+seed(:) = 102030405
+call random_seed(put=seed)
+
+scheme_labels(1) = 'PCM'
+scheme_labels(2) = 'PLM'
+
+! Set up some test data (note: using k,i indexing rather than i,k)
+allocate( u0(nk,nij), h0(nk,nij), u1(nk,nij), h1(nk,nij) )
+call random_number(u0) ! In range 0-1
+call random_number(h0) ! In range 0-1
+call random_number(h1) ! In range 0-1
+do ij = 1, nij
+  h0(:,ij) = max(0., h0(:,ij) - 0.05) ! Make 5% of values equal to zero
+  h1(:,ij) = max(0., h1(:,ij) - 0.05) ! Make 5% of values equal to zero
+  h0sum = h0(1,ij)
+  h1sum = h1(1,ij)
+  do k = 2, nk
+    h0sum = h0sum + h0(k,ij)
+    h1sum = h1sum + h1(k,ij)
+  enddo
+  h0(:,ij) = h0(:,ij) / h0sum
+  h1(:,ij) = h1(:,ij) / h1sum
+enddo
+
+! Loop over many samples of timing loop to collect statistics
+tmean(:) = 0.
+tstd(:) = 0.
+tmin(:) = 1.e9
+tmax(:) = 0.
+do isamp = 1, nsamp
+  ! Time reconstruction + remapping
+  do ischeme = 1, nschemes
+    call initialize_remapping(CS, remapping_scheme=trim(scheme_labels(ischeme)))
+    call cpu_time(start)
+    do iter = 1, nits ! Make many passes to reduce sampling error
+      do ij = 1, nij ! Calling nij times to make similar to cost in MOM_ALE()
+        call remapping_core_h(CS, nk, h0(:,ij), u0(:,ij), nk, h1(:,ij), u1(:,ij))
+      enddo
+    enddo
+    call cpu_time(finish)
+    timings(ischeme) = (finish-start)/real(nits*nij) ! Average time per call
+  enddo
+  tmean(:) = tmean(:) + timings(:)
+  tstd(:) = tstd(:) + timings(:)**2 ! tstd contains sum of squares here
+  tmin(:) = min( tmin(:), timings(:) )
+  tmax(:) = max( tmax(:), timings(:) )
+enddo
+tmean(:) = tmean(:) / real(nsamp) ! convert to mean
+tstd(:) = tstd(:) / real(nsamp) ! convert to mean of squares
+tstd(:) = tstd(:) - tmean(:)**2  ! convert to variance
+tstd(:) = sqrt( tstd(:) * real(nsamp) / real(nsamp-1) ) ! convert to standard deviation
+
+
+! Display results in YAML
+write(*,'(a)') "{"
+do ischeme = 1, nschemes
+  write(*,"(2x,5a)") '"MOM_remapping remapping_core_h(remapping_scheme=', &
+                     trim(scheme_labels(ischeme)), ')": {'
+  write(*,"(4x,a,1pe11.4,',')") '"min": ',tmin(ischeme)
+  write(*,"(4x,a,1pe11.4,',')") '"mean":',tmean(ischeme)
+  write(*,"(4x,a,1pe11.4,',')") '"std": ',tstd(ischeme)
+  write(*,"(4x,a,i7,',')") '"n_samples": ',nsamp
+  if (ischeme.ne.nschemes) then
+    write(*,"(4x,a,1pe11.4,'},')") '"max": ',tmax(ischeme)
+  else
+    write(*,"(4x,a,1pe11.4,'}')") '"max": ',tmax(ischeme)
+  endif
+enddo
+write(*,'(a)') "}"
+
+end program time_MOM_remapping

config_src/drivers/unit_tests/test_MOM_remapping.F90

@@ -0,0 +1,7 @@
+program test_MOM_remapping
+
+use MOM_remapping, only : remapping_unit_tests
+
+if (remapping_unit_tests(.true.)) stop 1
+
+end program test_MOM_remapping

config_src/external/drifters/MOM_particles_types.F90

@@ -3,6 +3,7 @@ module particles_types_mod
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+use, intrinsic :: iso_fortran_env, only : int64
 use MOM_grid, only : ocean_grid_type
 use MOM_domains, only: domain2D
 
@@ -75,7 +76,7 @@ module particles_types_mod
   real :: vvel_old !< Previous meridional velocity component (m/s)
   integer :: year  !< Year of this record
   integer :: particle_num  !< Current particle number
-  integer(kind=8) :: id = -1 !< Particle Identifier
+  integer(kind=int64) :: id = -1 !< Particle Identifier
   type(xyt), pointer :: next=>null()  !< Pointer to the next position in the list
 end type xyt
 
@@ -98,8 +99,8 @@ module particles_types_mod
   real :: start_day       !< origination position (degrees) and day
   integer :: start_year                         !< origination year
   real :: halo_part  !< equal to zero for particles on the computational domain, and 1 for particles on the halo
-  integer(kind=8) :: id                      !< particle identifier
-  integer(kind=8) :: drifter_num             !< particle identifier
+  integer(kind=int64) :: id                      !< particle identifier
+  integer(kind=int64) :: drifter_num             !< particle identifier
   integer :: ine                           !< nearest i-index in NE direction (for convenience)
   integer :: jne                           !< nearest j-index in NE direction (for convenience)
   real :: xi                               !< non-dimensional x-coordinate within current cell (0..1)
@@ -147,7 +148,7 @@ module particles_types_mod
   logical :: ignore_traj=.False. !< If true, then model does not write trajectory data at all
   logical :: use_new_predictive_corrective =.False. !< Flag to use Bob's predictive corrective particle scheme
   !Added by Alon
-  integer(kind=8) :: debug_particle_with_id = -1 !< If positive, monitors a part with this id
+  integer(kind=int64) :: debug_particle_with_id = -1 !< If positive, monitors a part with this id
   type(buffer), pointer :: obuffer_n=>null() !< Buffer for outgoing parts to the north
   type(buffer), pointer :: ibuffer_n=>null() !< Buffer for incoming parts from the north
   type(buffer), pointer :: obuffer_s=>null() !< Buffer for outgoing parts to the south

config_src/infra/FMS1/MOM_diag_manager_infra.F90

@@ -8,6 +8,7 @@ module MOM_diag_manager_infra
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+use, intrinsic :: iso_fortran_env, only : real64
 use diag_axis_mod,    only : fms_axis_init=>diag_axis_init
 use diag_axis_mod,    only : fms_get_diag_axis_name => get_diag_axis_name
 use diag_axis_mod,    only : EAST, NORTH
@@ -359,7 +360,7 @@ end function send_data_infra_3d
 logical function send_data_infra_2d_r8(diag_field_id, field, is_in, ie_in, js_in, je_in, &
                                        time, mask, rmask, weight, err_msg)
   integer,                           intent(in) :: diag_field_id !< The diagnostic manager identifier for this field
-  real(kind=8), dimension(:,:),      intent(in) :: field !< A 2-d array of values being recorded
+  real(kind=real64), dimension(:,:), intent(in) :: field !< A 2-d array of values being recorded
   integer,                 optional, intent(in) :: is_in !< The starting i-index for the data being recorded
   integer,                 optional, intent(in) :: ie_in !< The end i-index for the data being recorded
   integer,                 optional, intent(in) :: js_in !< The starting j-index for the data being recorded
@@ -382,7 +383,7 @@ end function send_data_infra_2d_r8
 logical function send_data_infra_3d_r8(diag_field_id, field, is_in, ie_in, js_in, je_in, ks_in, ke_in, &
                                     time, mask, rmask, weight, err_msg)
   integer,                             intent(in) :: diag_field_id !< The diagnostic manager identifier for this field
-  real(kind=8), dimension(:,:,:),      intent(in) :: field !< A rank 1 array of floating point values being recorded
+  real(kind=real64), dimension(:,:,:), intent(in) :: field !< A rank 1 array of floating point values being recorded
   integer,                   optional, intent(in) :: is_in !< The starting i-index for the data being recorded
   integer,                   optional, intent(in) :: ie_in !< The end i-index for the data being recorded
   integer,                   optional, intent(in) :: js_in !< The starting j-index for the data being recorded

config_src/infra/FMS2/MOM_diag_manager_infra.F90

@@ -8,6 +8,7 @@ module MOM_diag_manager_infra
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+use, intrinsic :: iso_fortran_env, only : real64
 use diag_axis_mod,    only : fms_axis_init=>diag_axis_init
 use diag_axis_mod,    only : fms_get_diag_axis_name => get_diag_axis_name
 use diag_axis_mod,    only : EAST, NORTH
@@ -361,7 +362,7 @@ end function send_data_infra_3d
 logical function send_data_infra_2d_r8(diag_field_id, field, is_in, ie_in, js_in, je_in, &
                                        time, mask, rmask, weight, err_msg)
   integer,                           intent(in) :: diag_field_id !< The diagnostic manager identifier for this field
-  real(kind=8), dimension(:,:),      intent(in) :: field !< A 2-d array of values being recorded
+  real(kind=real64), dimension(:,:), intent(in) :: field !< A 2-d array of values being recorded
   integer,                 optional, intent(in) :: is_in !< The starting i-index for the data being recorded
   integer,                 optional, intent(in) :: ie_in !< The end i-index for the data being recorded
   integer,                 optional, intent(in) :: js_in !< The starting j-index for the data being recorded
@@ -384,7 +385,7 @@ end function send_data_infra_2d_r8
 logical function send_data_infra_3d_r8(diag_field_id, field, is_in, ie_in, js_in, je_in, ks_in, ke_in, &
                                     time, mask, rmask, weight, err_msg)
   integer,                             intent(in) :: diag_field_id !< The diagnostic manager identifier for this field
-  real(kind=8), dimension(:,:,:),      intent(in) :: field !< A rank 1 array of floating point values being recorded
+  real(kind=real64), dimension(:,:,:), intent(in) :: field !< A rank 1 array of floating point values being recorded
   integer,                   optional, intent(in) :: is_in !< The starting i-index for the data being recorded
   integer,                   optional, intent(in) :: ie_in !< The end i-index for the data being recorded
   integer,                   optional, intent(in) :: js_in !< The starting j-index for the data being recorded

src/ALE/MOM_ALE.F90

@@ -178,6 +178,7 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   logical           :: local_logical
   logical           :: remap_boundary_extrap
   logical           :: init_boundary_extrap
+  logical           :: om4_remap_via_sub_cells
   type(hybgen_regrid_CS), pointer :: hybgen_regridCS => NULL() ! Control structure for hybgen regridding
                                                          ! for sharing parameters.
 
@@ -234,6 +235,11 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
   call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
                  "This sets the default value for the various _ANSWER_DATE parameters.", &
                  default=99991231)
+  call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
+                 "This selects the remapping algorithm used in OM4 that does not use "//&
+                 "the full reconstruction for the top- and lower-most sub-layers, but instead "//&
+                 "assumes they are always vanished (untrue) and so just uses their edge values. "//&
+                 "We recommend setting this option to false.", default=.true.)
   call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", CS%answer_date, &
                  "The vintage of the expressions and order of arithmetic to use for remapping.  "//&
                  "Values below 20190101 result in the use of older, less accurate expressions "//&
@@ -247,12 +253,14 @@ subroutine ALE_init( param_file, GV, US, max_depth, CS)
                              check_reconstruction=check_reconstruction, &
                              check_remapping=check_remapping, &
                              force_bounds_in_subcell=force_bounds_in_subcell, &
+                             om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
                              answer_date=CS%answer_date)
   call initialize_remapping( CS%vel_remapCS, vel_string, &
                              boundary_extrapolation=init_boundary_extrap, &
                              check_reconstruction=check_reconstruction, &
                              check_remapping=check_remapping, &
                              force_bounds_in_subcell=force_bounds_in_subcell, &
+                             om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
                              answer_date=CS%answer_date)
 
   call get_param(param_file, mdl, "PARTIAL_CELL_VELOCITY_REMAP", CS%partial_cell_vel_remap, &
@@ -326,6 +334,21 @@ subroutine ALE_set_extrap_boundaries( param_file, CS)
   call remapping_set_param(CS%remapCS, boundary_extrapolation=remap_boundary_extrap)
 end subroutine ALE_set_extrap_boundaries
 
+!> Sets the remapping algorithm to that of OM4
+!!
+!! The remapping aglorithm used in OM4 made poor assumptions about the reconstructions
+!! in the top/bottom layers, namely that they were always vanished and could be
+!! represented solely by their upper/lower edge value respectively.
+!! Passing .false. here uses the full reconstruction of those top and bottom layers
+!! and properly sample those layers.
+subroutine ALE_set_OM4_remap_algorithm( CS, om4_remap_via_sub_cells )
+  type(ALE_CS), pointer :: CS !< Module control structure
+  logical, intent(in)   :: om4_remap_via_sub_cells !< If true, use OM4 remapping algorithm
+
+  call remapping_set_param(CS%remapCS, om4_remap_via_sub_cells =om4_remap_via_sub_cells )
+
+end subroutine ALE_set_OM4_remap_algorithm
+
 !> Initialize diagnostics for the ALE module.
 subroutine ALE_register_diags(Time, G, GV, US, diag, CS)
   type(time_type),target,     intent(in)  :: Time  !< Time structure