Merge branch 'dev/gfdl' into unscale_chksum_arg

.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)
@@ -649,7 +651,7 @@ $(WORK)/%/restart/ocean.stats: $(BUILD)/symmetric/MOM6 | preproc
 # Not a true rule; only call this after `make test` to summarize test results.
 .PHONY: test.summary
 test.summary:
-	@./tools/report_test_results.sh $(WORK)/results
+	./tools/report_test_results.sh $(WORK)/results
 
 
 #---

.testing/tools/parse_perf.py

@@ -102,7 +102,7 @@ def parse_perf_report(perf_data_path):
 
                     if tok == '>':
                         bracks -= 1
-                    if tok == '(':
+                    if tok == ')':
                         parens -= 1
 
                 # Strip any whitespace tokens

ac/configure.ac

@@ -81,8 +81,8 @@ AS_IF([test "x$with_driver" != "x"],
 
 
 # Explicitly assume free-form Fortran
-AC_LANG(Fortran)
-AC_FC_SRCEXT(f90)
+AC_LANG([Fortran])
+AC_FC_SRCEXT([f90])
 
 
 # Determine MPI compiler wrappers

ac/deps/configure.fms.ac

@@ -10,7 +10,16 @@ AC_INIT(
 AC_CONFIG_SRCDIR([fms/fms.F90])
 AC_CONFIG_MACRO_DIR([m4])
 
+
 # C configuration
+
+# Autoconf assumes that LDFLAGS can be passed to CFLAGS, even though this is
+# not valid in some compilers.  This can cause basic CC tests to fail.
+# Since we do not link with CC, we can safely disable LDFLAGS for AC_PROG_CC.
+FC_LDFLAGS="$LDFLAGS"
+LDFLAGS=""
+
+# C compiler verification
 AC_PROG_CC
 AX_MPI
 CC=$MPICC
@@ -55,10 +64,13 @@ AC_CHECK_FUNCS([gettid], [], [
 # FMS 2019.01.03 uses __APPLE__ to disable Linux CPU affinity calls.
 AC_CHECK_FUNCS([sched_getaffinity], [], [AC_DEFINE([__APPLE__])])
 
+# Restore LDFLAGS
+LDFLAGS="$FC_LDFLAGS"
+
 
 # Standard Fortran configuration
-AC_LANG(Fortran)
-AC_FC_SRCEXT(f90)
+AC_LANG([Fortran])
+AC_FC_SRCEXT([f90])
 AC_PROG_FC
 
 

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/GFDL_ocean_BGC/generic_tracer_utils.F90

@@ -99,8 +99,9 @@ module g_tracer_utils
 contains
 
   !> Unknown
-  subroutine g_tracer_flux_init(g_tracer)
+  subroutine g_tracer_flux_init(g_tracer, verbosity)
     type(g_tracer_type), pointer :: g_tracer !< Pointer to this tracer node
+    integer, optional, intent(in) :: verbosity !< A 0-9 integer indicating a level of verbosity
   end subroutine g_tracer_flux_init
 
   !> Unknown

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

src/ALE/MOM_regridding.F90

@@ -211,6 +211,8 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
   real :: tmpReal  ! A temporary variable used in setting other variables [various]
   real :: P_Ref    ! The coordinate variable reference pression [R L2 T-2 ~> Pa]
   real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z).
+  real :: dz_extra      ! The thickness of an added layer to append to the woa09_dz profile when
+                        ! maximum_depth is large [m] (not in Z).
   real :: adaptTimeRatio, adaptZoomCoeff ! Temporary variables for input parameters [nondim]
   real :: adaptBuoyCoeff, adaptAlpha     ! Temporary variables for input parameters [nondim]
   real :: adaptZoom  ! The thickness of the near-surface zooming region with the adaptive coordinate [H ~> m or kg m-2]
@@ -311,7 +313,7 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
     param_name = create_coord_param(param_prefix, "DEF", param_suffix)
     coord_res_param = create_coord_param(param_prefix, "RES", param_suffix)
     string2 = 'UNIFORM'
-    if (maximum_depth>3000.) string2='WOA09' ! For convenience
+    if ((maximum_depth>3000.) .and. (maximum_depth<9250.)) string2='WOA09' ! For convenience
   endif
   call get_param(param_file, mdl, param_name, string, &
                  "Determines how to specify the coordinate "//&
@@ -458,20 +460,27 @@ subroutine initialize_regridding(CS, GV, US, max_depth, param_file, mdl, coord_m
     endif
   elseif (index(trim(string),'WOA09')==1) then
     if (len_trim(string)==5) then
-      tmpReal = 0. ; ke = 0
+      tmpReal = 0. ; ke = 0 ; dz_extra = 0.
       do while (tmpReal<maximum_depth)
         ke = ke + 1
+        if (ke > size(woa09_dz)) then
+          dz_extra = maximum_depth - tmpReal
+          exit
+        endif
         tmpReal = tmpReal + woa09_dz(ke)
       enddo
     elseif (index(trim(string),'WOA09:')==1) then
       if (len_trim(string)==6) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
                  'Expected string of form "WOA09:N" but got "'//trim(string)//'".')
       ke = extract_integer(string(7:len_trim(string)),'',1)
+      if (ke>40 .or. ke<1) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
+                   'For "WOA05:N" N must 0<N<41 but got "'//trim(string)//'".')
     endif
-    if (ke>40 .or. ke<1) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
-                 'For "WOA05:N" N must 0<N<41 but got "'//trim(string)//'".')
     allocate(dz(ke))
-    dz(1:ke) = woa09_dz(1:ke)
+    do k=1,min(ke, size(woa09_dz))
+      dz(k) = woa09_dz(k)
+    enddo
+    if (ke > size(woa09_dz)) dz(ke) = dz_extra
   else
     call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
       "Unrecognized coordinate configuration"//trim(string))