Merge branch 'dev/gfdl' of github.com:NOAA-GFDL/MOM6 into dev/gfdl

* 'dev/gfdl' of github.com:NOAA-GFDL/MOM6:
  Fix to out-of-bounds error
  Slight improvement in writeMessageAndDesc
  Corrected recent bug in writeMessageAndDesc
  Included H_subroundoff in a denominator
  +Added code to handle tab lengths in documentation
  Added explicit form of cpu_clock_id
  (*)Improve ePBL when EPBL_ORIGINAL_PE_CALC=False
  Alternate fix for divide-by-zero.
  Clean up some divide-by-zero issues in OBCs.
  Added a space to a description.
  Codecov: Set base as parent (target) branch
  Fixes an integer-kind mismatch in MOM_random, seed_from_time() (mom-ocean#1113)

.codecov.yml

@@ -3,9 +3,11 @@ coverage:
         project:
             default:
                 threshold: 100%
+                base: parent
         patch:
             default:
                 threshold: 100%
+                base: parent
 comment:
     # This must be set to the number of test cases (TCs)
     after_n_builds: 8

src/core/MOM_open_boundary.F90

@@ -4766,40 +4766,48 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
     segment=>OBC%segment(n)
     if (.not. associated(segment%tr_Reg)) cycle
     if (segment%is_E_or_W) then
+      I = segment%HI%IsdB
       do j=segment%HI%jsd,segment%HI%jed
-        I = segment%HI%IsdB
         ! ishift+I corresponds to the nearest interior tracer cell index
         ! idir switches the sign of the flow so that positive is into the reservoir
         if (segment%direction == OBC_DIRECTION_W) then
           ishift = 1 ; idir = -1
         else
           ishift = 0 ; idir = 1
         endif
+        ! Can keep this or take it out, either way
+        if (G%mask2dT(I+ishift,j) == 0.0) cycle
         ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep
         do m=1,ntr ; if (associated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz
-          u_L_out = max(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_out / (h(i+ishift,j,k)*G%dyCu(I,j)))
-          u_L_in  = min(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_in  / (h(i+ishift,j,k)*G%dyCu(I,j)))
+          u_L_out = max(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_out / &
+                    ((h(i+ishift,j,k) + GV%H_subroundoff)*G%dyCu(I,j)))
+          u_L_in  = min(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_in  / &
+                    ((h(i+ishift,j,k) + GV%H_subroundoff)*G%dyCu(I,j)))
           fac1 = 1.0 + (u_L_out-u_L_in)
           segment%tr_Reg%Tr(m)%tres(I,j,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(I,j,k) + &
                             (u_L_out*Reg%Tr(m)%t(I+ishift,j,k) - &
                              u_L_in*segment%tr_Reg%Tr(m)%t(I,j,k)))
           if (associated(OBC%tres_x)) OBC%tres_x(I,j,k,m) = segment%tr_Reg%Tr(m)%tres(I,j,k)
         enddo ; endif ; enddo
       enddo
-    else
+    elseif (segment%is_N_or_S) then
+      J = segment%HI%JsdB
       do i=segment%HI%isd,segment%HI%ied
-        J = segment%HI%JsdB
         ! jshift+J corresponds to the nearest interior tracer cell index
         ! jdir switches the sign of the flow so that positive is into the reservoir
         if (segment%direction == OBC_DIRECTION_S) then
           jshift = 1 ; jdir = -1
         else
           jshift = 0 ; jdir = 1
         endif
+        ! Can keep this or take it out, either way
+        if (G%mask2dT(i,j+jshift) == 0.0) cycle
         ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep
         do m=1,ntr ; if (associated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz
-          v_L_out = max(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_out / (h(i,j+jshift,k)*G%dxCv(i,J)))
-          v_L_in  = min(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_in  / (h(i,j+jshift,k)*G%dxCv(i,J)))
+          v_L_out = max(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_out / &
+                    ((h(i,j+jshift,k) + GV%H_subroundoff)*G%dxCv(i,J)))
+          v_L_in  = min(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_in  / &
+                    ((h(i,j+jshift,k) + GV%H_subroundoff)*G%dxCv(i,J)))
           fac1 = 1.0 + (v_L_out-v_L_in)
           segment%tr_Reg%Tr(m)%tres(i,J,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,J,k) + &
                             (v_L_out*Reg%Tr(m)%t(i,J+jshift,k) - &

src/framework/MOM_cpu_clock.F90

@@ -3,8 +3,9 @@ module MOM_cpu_clock
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
+use fms_mod, only : clock_flag_default
 use mpp_mod, only : cpu_clock_begin => mpp_clock_begin
-use mpp_mod, only : cpu_clock_end => mpp_clock_end, cpu_clock_id => mpp_clock_id
+use mpp_mod, only : cpu_clock_end => mpp_clock_end, mpp_clock_id
 use mpp_mod, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT, CLOCK_MODULE_DRIVER
 use mpp_mod, only : CLOCK_MODULE, CLOCK_ROUTINE, CLOCK_LOOP, CLOCK_INFRA
 use mpp_mod, only : CLOCK_SYNC => MPP_CLOCK_SYNC
@@ -15,4 +16,27 @@ module MOM_cpu_clock
 public :: CLOCK_COMPONENT, CLOCK_SUBCOMPONENT, CLOCK_MODULE_DRIVER, CLOCK_MODULE
 public :: CLOCK_ROUTINE, CLOCK_LOOP, CLOCK_INFRA, CLOCK_SYNC
 
+contains
+
+!> cpu_clock_id returns the integer handle for a named CPU clock.
+function cpu_clock_id( name, synchro_flag, grain )
+  character(len=*),  intent(in) :: name  !< The unique name of the CPU clock
+  integer, intent(in), optional :: synchro_flag !< An integer flag that controls whether the PEs
+                                       !! are synchronized before the cpu clocks start counting.
+                                       !! Synchronization occurs before the start of a clock if this
+                                       !! is odd, while additional (expensive) statistics can set
+                                       !! for other values. If absent, the default is taken from the
+                                       !! settings for FMS.
+  integer, intent(in), optional :: grain !< The timing granularity for this clock, usually set to
+                                       !! the values of CLOCK_COMPONENT, CLOCK_ROUTINE, CLOCK_LOOP, etc.
+  integer                       :: cpu_clock_id !< The integer CPU clock handle.
+
+  if (present(synchro_flag)) then
+    cpu_clock_id = mpp_clock_id(name, flags=synchro_flag, grain=grain)
+  else
+    cpu_clock_id = mpp_clock_id(name, flags=clock_flag_default, grain=grain)
+  endif
+
+end function cpu_clock_id
+
 end module MOM_cpu_clock

src/framework/MOM_document.F90

@@ -468,6 +468,9 @@ subroutine writeMessageAndDesc(doc, vmesg, desc, valueWasDefault, indent, &
   integer :: start_ind = 1             ! The starting index in the description for the next line.
   integer :: nl_ind, tab_ind, end_ind  ! The indices of new-lines, tabs, and the end of a line.
   integer :: len_text, len_tab, len_nl ! The lengths of the text string, tabs and new-lines.
+  integer :: len_cor                   ! The permitted length corrected for tab sizes in a line.
+  integer :: len_desc                  ! The non-whitespace length of the description.
+  integer :: substr_start              ! The starting index of a substring to search for tabs.
   integer :: indnt, msg_pad            ! Space counts used to format a message.
   logical :: msg_done, reset_msg_pad   ! Logicals used to format messages.
   logical :: all, short, layout, debug ! Flags indicating which files to write into.
@@ -494,16 +497,27 @@ subroutine writeMessageAndDesc(doc, vmesg, desc, valueWasDefault, indent, &
   do
     if (len_trim(desc(start_ind:)) < 1) exit
 
-    nl_ind = index(desc(start_ind:), "\n")
+    len_cor = len_text - msg_pad
 
+    substr_start = start_ind
+    len_desc = len_trim(desc)
+    do ! Adjust the available line length for anomalies in the size of tabs, counting \t as 2 spaces.
+      if (substr_start >= start_ind+len_cor) exit
+      tab_ind = index(desc(substr_start:min(len_desc,start_ind+len_cor)), "\t")
+      if (tab_ind == 0) exit
+      substr_start = substr_start + tab_ind
+      len_cor = len_cor + (len_tab - 2)
+    enddo
+
+    nl_ind = index(desc(start_ind:), "\n")
     end_ind = 0
-    if ((nl_ind > 0) .and. (len_trim(desc(start_ind:start_ind+nl_ind-2)) > len_text-msg_pad)) then
+    if ((nl_ind > 0) .and. (len_trim(desc(start_ind:start_ind+nl_ind-2)) > len_cor)) then
       ! This line is too long despite the new-line character.  Look for an earlier space to break.
-      end_ind = scan(desc(start_ind:start_ind+(len_text-msg_pad)), " ", back=.true.) - 1
+      end_ind = scan(desc(start_ind:start_ind+len_cor), " ", back=.true.) - 1
       if (end_ind > 0) nl_ind = 0
-    elseif ((nl_ind == 0) .and. (len_trim(desc(start_ind:)) > len_text-msg_pad)) then
+    elseif ((nl_ind == 0) .and. (len_trim(desc(start_ind:)) > len_cor)) then
       ! This line is too long and does not have a new-line character.  Look for a space to break.
-      end_ind = scan(desc(start_ind:start_ind+(len_text-msg_pad)), " ", back=.true.) - 1
+      end_ind = scan(desc(start_ind:start_ind+len_cor), " ", back=.true.) - 1
     endif
 
     reset_msg_pad = .false.

src/framework/MOM_random.F90

@@ -143,8 +143,9 @@ integer function seed_from_time(Time)
   call get_date(Time,yr,mo,dy,hr,mn,sc)
   s1 = sc + 61*(mn + 61*hr) + 379 ! Range 379 .. 89620
   ! Fun fact: 2147483647 is the eighth Mersenne prime.
-  ! This is not the reason for using 2147483647+1 here.
-  s2 = mod(dy + 32*(mo + 13*yr), 2147483648) ! Range 0 .. 2147483647
+  ! This is not the reason for using 2147483647 here. It is the
+  ! largest integer of kind=4.
+  s2 = modulo(dy + 32*(mo + 13*yr), 2147483647_4) ! Range 0 .. 2147483646
   seed_from_time = ieor(s1*4111, s2)
 
 end function seed_from_time

src/parameterizations/lateral/MOM_MEKE.F90

@@ -1063,7 +1063,7 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS)
                  "If true, use an alternative formula for computing the (equilibrium)"//&
                  "initial value of MEKE.", default=.false.)
   call get_param(param_file, mdl, "MEKE_EQUILIBRIUM_RESTORING", CS%MEKE_equilibrium_restoring, &
-                 "If true, restore MEKE back to its equilibrium value, which is calculated at"//&
+                 "If true, restore MEKE back to its equilibrium value, which is calculated at "//&
                  "each time step.", default=.false.)
   if (CS%MEKE_equilibrium_restoring) then
     call get_param(param_file, mdl, "MEKE_RESTORING_TIMESCALE", MEKE_restoring_timescale, &

src/parameterizations/vertical/MOM_bkgnd_mixing.F90

@@ -445,7 +445,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay, Kv, j, G, GV, US, CS)
 
   elseif ((.not. CS%Bryan_Lewis_diffusivity) .and. (.not.CS%bulkmixedlayer) .and. &
           (.not. CS%horiz_varying_background) .and. (CS%Kd /= CS%Kdml)) then
-    I_Hmix = 1.0 / CS%Hmix
+    I_Hmix = 1.0 / (CS%Hmix + GV%H_subroundoff*GV%H_to_Z)
     do i=is,ie ; depth(i) = 0.0 ; enddo
     do k=1,nz ; do i=is,ie
       depth_c = depth(i) + 0.5*GV%H_to_Z*h(i,j,k)

src/parameterizations/vertical/MOM_energetic_PBL.F90

@@ -772,7 +772,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
   integer :: OBL_it        ! Iteration counter
 
   real :: Surface_Scale ! Surface decay scale for vstar
-
+  logical :: calc_dT_expect ! If true calculate the expected changes in temperature and salinity.
+  logical :: calc_Te        ! If true calculate the expected final temperature and salinity values.
   logical :: debug=.false.  ! Change this hard-coded value for debugging.
 
   !  The following arrays are used only for debugging purposes.
@@ -788,7 +789,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
   if (.not. associated(CS)) call MOM_error(FATAL, "energetic_PBL: "//&
          "Module must be initialized before it is used.")
 
-  debug = .false. ; if (allocated(eCD%dT_expect) .or. allocated(eCD%dS_expect)) debug = .true.
+  calc_dT_expect = debug ; if (allocated(eCD%dT_expect) .or. allocated(eCD%dS_expect)) calc_dT_expect = .true.
+  calc_Te = (calc_dT_expect .or. (.not.CS%orig_PE_calc))
 
   h_neglect = GV%H_subroundoff
 
@@ -1285,7 +1287,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
                          dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE(k), dS_to_dPE(k), &
                          pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), &
                          dT_to_dColHt(k), dS_to_dColHt(k), &
-                         PE_chg=dPE_conv)
+                         PE_chg=dPE_conv, dPEc_dKd=dPEc_dKd)
               endif
               MKE_src = dMKE_max * (1.0 - exp(-MKE2_Hharm * Kddt_h_guess))
               dMKE_src_dK = dMKE_max * MKE2_Hharm * exp(-MKE2_Hharm * Kddt_h_guess)
@@ -1381,7 +1383,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
           htot  = htot + h(k)
         endif
 
-        if (debug) then
+        if (calc_Te) then
           if (k==2) then
             Te(1) = b1*(h(1)*T0(1))
             Se(1) = b1*(h(1)*S0(1))
@@ -1393,7 +1395,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
       enddo
       Kd(nz+1) = 0.0
 
-      if (debug) then
+      if (calc_dT_expect) then
         ! Complete the tridiagonal solve for Te.
         b1 = 1.0 / hp_a
         Te(nz) = b1 * (h(nz) * T0(nz) + Kddt_h(nz) * Te(nz-1))
@@ -1404,7 +1406,9 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs
           Se(k) = Se(k) + c1(K+1)*Se(k+1)
           eCD%dT_expect(k) = Te(k) - T0(k) ; eCD%dS_expect(k) = Se(k) - S0(k)
         enddo
+      endif
 
+      if (debug) then
         dPE_debug = 0.0
         do k=1,nz
           dPE_debug = dPE_debug + (dT_to_dPE(k) * (Te(k) - T0(k)) + &