Merge branch 'NOAA-GFDL:dev/gfdl' into evp_in_BTCS_tmp

Author: theresa-cordero

.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
 
 
 #---

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/FMS_cap/MOM_surface_forcing_gfdl.F90

@@ -635,16 +635,20 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
   endif
 
   ! Set the wind stresses and ustar.
-  if (associated(fluxes%ustar) .and. associated(fluxes%ustar_gustless) .and. associated(fluxes%tau_mag)) then
+  if (associated(fluxes%ustar) .and. associated(fluxes%ustar_gustless) .and. associated(fluxes%tau_mag) &
+      .and. associated(fluxes%tau_mag_gustless) ) then
     call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, ustar=fluxes%ustar, &
-                              mag_tau=fluxes%tau_mag, gustless_ustar=fluxes%ustar_gustless)
+                              mag_tau=fluxes%tau_mag, gustless_ustar=fluxes%ustar_gustless, &
+                              gustless_mag_tau=fluxes%tau_mag_gustless)
   else
     if (associated(fluxes%ustar)) &
       call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, ustar=fluxes%ustar)
     if (associated(fluxes%ustar_gustless)) &
       call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, gustless_ustar=fluxes%ustar_gustless)
     if (associated(fluxes%tau_mag)) &
       call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, mag_tau=fluxes%tau_mag)
+    if (associated(fluxes%tau_mag_gustless)) &
+      call extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, gustless_mag_tau=fluxes%tau_mag_gustless)
   endif
 
   if (coupler_type_initialized(fluxes%tr_fluxes) .and. &
@@ -908,7 +912,7 @@ end subroutine convert_IOB_to_forces
 !! Ice_ocean_boundary_type into optional argument arrays, including changes of units, sign
 !! conventions, and putting the fields into arrays with MOM-standard sized halos.
 subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy, ustar, &
-                                gustless_ustar, mag_tau, tau_halo)
+                                gustless_ustar, mag_tau, gustless_mag_tau, tau_halo)
   type(ice_ocean_boundary_type), &
                    target, intent(in)    :: IOB  !< An ice-ocean boundary type with fluxes to drive
                                                  !! the ocean in a coupled model
@@ -931,6 +935,9 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
   real, dimension(SZI_(G),SZJ_(G)), &
                  optional, intent(inout) :: mag_tau !< The magintude of the wind stress at tracer points
                                                  !! including subgridscale variability and gustiness [R Z L T-2 ~> Pa]
+  real, dimension(SZI_(G),SZJ_(G)), &
+                 optional, intent(out) :: gustless_mag_tau !< The magintude of the wind stress at tracer points
+                                                 !! without any contributions from gustiness [R Z L T-2 ~> Pa]
   integer,       optional, intent(in)    :: tau_halo !< The halo size of wind stresses to set, 0 by default.
 
   ! Local variables
@@ -947,7 +954,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
   real :: tau_mag       ! magnitude of the wind stress [R Z L T-2 ~> Pa]
   real :: stress_conversion ! A unit conversion factor from Pa times any stress multiplier [R Z L T-2 Pa-1 ~> 1]
 
-  logical :: do_ustar, do_gustless, do_tau_mag
+  logical :: do_ustar, do_gustless, do_tau_mag, do_gustless_tau_mag
   integer :: wind_stagger  ! AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains)
   integer :: i, j, is, ie, js, je, ish, ieh, jsh, jeh, Isqh, Ieqh, Jsqh, Jeqh, i0, j0, halo
 
@@ -960,7 +967,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
   IRho0 = US%L_to_Z / CS%Rho0
   stress_conversion = US%Pa_to_RLZ_T2 * CS%wind_stress_multiplier
 
-  do_ustar = present(ustar) ; do_gustless = present(gustless_ustar) ; do_tau_mag = present(mag_tau)
+  do_ustar = present(ustar) ; do_gustless = present(gustless_ustar)
+  do_tau_mag = present(mag_tau) ; do_gustless_tau_mag = present(gustless_mag_tau)
 
   wind_stagger = CS%wind_stagger
   if ((IOB%wind_stagger == AGRID) .or. (IOB%wind_stagger == BGRID_NE) .or. &
@@ -973,7 +981,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
 
   ! Set surface momentum stress related fields as a function of staggering.
   if (present(taux) .or. present(tauy) .or. &
-      ((do_ustar .or. do_tau_mag .or. do_gustless) .and. .not.associated(IOB%stress_mag)) ) then
+      ((do_ustar .or. do_tau_mag .or. do_gustless .or. do_gustless_tau_mag) &
+       .and. .not.associated(IOB%stress_mag)) ) then
 
     if (wind_stagger == BGRID_NE) then
       taux_in_B(:,:) = 0.0 ; tauy_in_B(:,:) = 0.0
@@ -1053,7 +1062,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
     endif   ! endif for extracting wind stress fields with various staggerings
   endif
 
-  if (do_ustar .or. do_tau_mag .or. do_gustless) then
+  if (do_ustar .or. do_tau_mag .or. do_gustless .or. do_gustless_tau_mag) then
     ! Set surface friction velocity directly or as a function of staggering.
     ! ustar is required for the bulk mixed layer formulation and other turbulent mixing
     ! parametizations. The background gustiness (for example with a relatively small value
@@ -1071,6 +1080,8 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
         endif
         if (do_tau_mag) &
           mag_tau(i,j) = gustiness + US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0)
+        if (do_gustless_tau_mag) &
+          gustless_mag_tau(i,j) = US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0)
         if (do_ustar) &
           ustar(i,j) = sqrt(gustiness*IRho0 + IRho0*US%Pa_to_RLZ_T2*IOB%stress_mag(i-i0,j-j0))
       enddo ; enddo ; endif
@@ -1097,6 +1108,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
         endif
         if (do_ustar) ustar(i,j) = sqrt(gustiness*IRho0 + IRho0 * tau_mag)
         if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag
+        if (do_gustless_tau_mag) gustless_mag_tau(i,j) = tau_mag
         if (CS%answer_date < 20190101) then
           if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0)
         else
@@ -1110,6 +1122,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
         if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0.0)) gustiness = CS%gust(i,j)
         if (do_ustar) ustar(i,j) = sqrt(gustiness*IRho0 + IRho0 * tau_mag)
         if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag
+        if (do_gustless_tau_mag) gustless_mag_tau(i,j) = tau_mag
         if (CS%answer_date < 20190101) then
           if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0)
         else
@@ -1132,6 +1145,7 @@ subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, US, CS, taux, tauy,
 
         if (do_ustar) ustar(i,j) = sqrt(gustiness*IRho0 + IRho0 * tau_mag)
         if (do_tau_mag) mag_tau(i,j) = gustiness + tau_mag
+        if (do_gustless_tau_mag) gustless_mag_tau(i,j) = tau_mag
         if (CS%answer_date < 20190101) then
           if (do_gustless) gustless_ustar(i,j) = sqrt(US%L_to_Z*tau_mag / CS%Rho0)
         else

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