Merge branch 'dev/gfdl' into sync_clock_flag

config_src/infra/FMS1/MOM_domain_infra.F90

@@ -1716,13 +1716,13 @@ subroutine clone_MD_to_d2D(MD_in, mpp_domain, min_halo, halo_size, symmetric, &
                 symmetry=symmetric_dom, xextent=xextent, yextent=yextent, name=dom_name)
   endif
 
-  if ((MD_in%io_layout(1) + MD_in%io_layout(2) > 0) .and. &
-      (MD_in%layout(1)*MD_in%layout(2) > 1)) then
-    call mpp_define_io_domain(mpp_domain, MD_in%io_layout)
-  else
-    call mpp_define_io_domain(mpp_domain, (/ 1, 1 /) )
+  if (MD_in%layout(1) * MD_in%layout(2) > 1) then
+    if ((MD_in%io_layout(1) + MD_in%io_layout(2) > 0)) then
+      call mpp_define_io_domain(mpp_domain, MD_in%io_layout)
+    else
+      call mpp_define_io_domain(mpp_domain, [1, 1] )
+    endif
   endif
-
 end subroutine clone_MD_to_d2D
 
 !> Returns the index ranges that have been stored in a MOM_domain_type

config_src/infra/FMS1/MOM_io_infra.F90

@@ -137,15 +137,12 @@ logical function MOM_file_exists(filename, MOM_Domain)
 end function MOM_file_exists
 
 !> Returns true if the named file or its domain-decomposed variant exists.
-logical function FMS_file_exists(filename, domain, no_domain)
+logical function FMS_file_exists(filename)
   character(len=*),         intent(in) :: filename  !< The name of the file being inquired about
-  type(domain2d), optional, intent(in) :: domain    !< The mpp domain2d that describes the decomposition
-  logical,        optional, intent(in) :: no_domain !< This file does not use domain decomposition
-! This function uses the fms_io function file_exist to determine whether
-! a named file (or its decomposed variant) exists.
-
-  FMS_file_exists = file_exist(filename, domain, no_domain)
+  ! This function uses the fms_io function file_exist to determine whether
+  ! a named file (or its decomposed variant) exists.
 
+  FMS_file_exists = file_exist(filename)
 end function FMS_file_exists
 
 !> indicates whether an I/O handle is attached to an open file

config_src/infra/FMS2/MOM_io_infra.F90

@@ -17,6 +17,7 @@ module MOM_io_infra
 use fms2_io_mod,          only : get_variable_dimension_names, is_dimension_registered, get_dimension_size
 use fms2_io_mod,          only : is_dimension_unlimited, register_axis, unlimited
 use fms2_io_mod,          only : get_global_io_domain_indices
+use fms_io_utils_mod,     only : fms2_file_exist => file_exists
 
 use fms_mod,              only : write_version_number, open_namelist_file, check_nml_error
 use fms_io_mod,           only : file_exist, field_exist, field_size, read_data
@@ -170,15 +171,12 @@ logical function MOM_file_exists(filename, MOM_Domain)
 end function MOM_file_exists
 
 !> Returns true if the named file or its domain-decomposed variant exists.
-logical function FMS_file_exists(filename, domain, no_domain)
+logical function FMS_file_exists(filename)
   character(len=*),         intent(in) :: filename  !< The name of the file being inquired about
-  type(domain2d), optional, intent(in) :: domain    !< The mpp domain2d that describes the decomposition
-  logical,        optional, intent(in) :: no_domain !< This file does not use domain decomposition
-! This function uses the fms_io function file_exist to determine whether
-! a named file (or its decomposed variant) exists.
-
-  FMS_file_exists = file_exist(filename, domain, no_domain)
+  ! This function uses the fms_io function file_exist to determine whether
+  ! a named file (or its decomposed variant) exists.
 
+  FMS_file_exists = fms2_file_exist(filename)
 end function FMS_file_exists
 
 !> indicates whether an I/O handle is attached to an open file

src/diagnostics/MOM_sum_output.F90

@@ -18,8 +18,6 @@ module MOM_sum_output
 use MOM_io,            only : axis_info, set_axis_info, delete_axis_info, get_filename_appendix
 use MOM_io,            only : attribute_info, set_attribute_info, delete_attribute_info
 use MOM_io,            only : APPEND_FILE, SINGLE_FILE, WRITEONLY_FILE
-use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type
-use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S
 use MOM_time_manager,  only : time_type, get_time, get_date, set_time, operator(>)
 use MOM_time_manager,  only : operator(+), operator(-), operator(*), operator(/)
 use MOM_time_manager,  only : operator(/=), operator(<=), operator(>=), operator(<)
@@ -297,7 +295,7 @@ end subroutine MOM_sum_output_end
 
 !>  This subroutine calculates and writes the total model energy, the energy and
 !! mass of each layer, and other globally integrated  physical quantities.
-subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_forcing)
+subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forcing)
   type(ocean_grid_type),   intent(in)    :: G   !< The ocean's grid structure.
   type(verticalGrid_type), intent(in)    :: GV  !< The ocean's vertical grid structure.
   type(unit_scale_type),   intent(in)    :: US  !< A dimensional unit scaling type
@@ -314,11 +312,10 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_
                                                 !! current execution.
   type(Sum_output_CS),     pointer       :: CS  !< The control structure returned by a
                                                 !! previous call to MOM_sum_output_init.
-  type(tracer_flow_control_CS), &
-                    optional, pointer    :: tracer_CSp !< tracer control structure.
-  type(ocean_OBC_type),         &
-                    optional, pointer    :: OBC !< Open boundaries control structure.
+  type(tracer_flow_control_CS), pointer  :: tracer_CSp !< Control structure with the tree of
+                                                !! all registered tracer packages
   type(time_type),  optional, intent(in) :: dt_forcing !< The forcing time step
+
   ! Local variables
   real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! The height of interfaces [Z ~> m].
   real :: areaTm(SZI_(G),SZJ_(G)) ! A masked version of areaT [L2 ~> m2].
@@ -409,17 +406,22 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_
   character(len=32)  :: mesg_intro, time_units, day_str, n_str, date_str
   logical :: date_stamped
   type(time_type) :: dt_force ! A time_type version of the forcing timestep.
-  real :: Tr_stocks(MAX_FIELDS_)
-  real :: Tr_min(MAX_FIELDS_), Tr_max(MAX_FIELDS_)
-  real :: Tr_min_x(MAX_FIELDS_), Tr_min_y(MAX_FIELDS_), Tr_min_z(MAX_FIELDS_)
-  real :: Tr_max_x(MAX_FIELDS_), Tr_max_y(MAX_FIELDS_), Tr_max_z(MAX_FIELDS_)
-  logical :: Tr_minmax_got(MAX_FIELDS_) = .false.
+  real :: Tr_stocks(MAX_FIELDS_) ! The total amounts of each of the registered tracers
+  real :: Tr_min(MAX_FIELDS_)   ! The global minimum unmasked value of the tracers
+  real :: Tr_max(MAX_FIELDS_)   ! The global maximum unmasked value of the tracers
+  real :: Tr_min_x(MAX_FIELDS_) ! The x-positions of the global tracer minima
+  real :: Tr_min_y(MAX_FIELDS_) ! The y-positions of the global tracer minima
+  real :: Tr_min_z(MAX_FIELDS_) ! The z-positions of the global tracer minima
+  real :: Tr_max_x(MAX_FIELDS_) ! The x-positions of the global tracer maxima
+  real :: Tr_max_y(MAX_FIELDS_) ! The y-positions of the global tracer maxima
+  real :: Tr_max_z(MAX_FIELDS_) ! The z-positions of the global tracer maxima
+  logical :: Tr_minmax_avail(MAX_FIELDS_) ! A flag indicating whether the global minimum and
+                                ! maximum information are available for each of the tracers
   character(len=40), dimension(MAX_FIELDS_) :: &
-    Tr_names, Tr_units
-  integer :: nTr_stocks
+    Tr_names, &          ! The short names for each of the tracers
+    Tr_units             ! The units for each of the tracers
+  integer :: nTr_stocks  ! The total number of tracers in all registered tracer packages
   integer :: iyear, imonth, iday, ihour, iminute, isecond, itick ! For call to get_date()
-  logical :: local_open_BC
-  type(OBC_segment_type), pointer :: segment => NULL()
 
  ! A description for output of each of the fields.
   type(vardesc) :: vars(NUM_FIELDS+MAX_FIELDS_)
@@ -479,16 +481,10 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_
     vars(17) = var_desc("Heat_anom","Joules","Anomalous Total Heat Change",'1','1')
   endif
 
-  local_open_BC = .false.
-  if (present(OBC)) then ; if (associated(OBC)) then
-    local_open_BC = (OBC%open_u_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally)
-  endif ; endif
-
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
   isr = is - (G%isd-1) ; ier = ie - (G%isd-1) ; jsr = js - (G%jsd-1) ; jer = je - (G%jsd-1)
 
-
   HL2_to_kg = GV%H_to_kg_m2*US%L_to_m**2
 
   if (.not.associated(CS)) call MOM_error(FATAL, &
@@ -504,34 +500,6 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_
       tmp1(i,j,k) = h(i,j,k) * (HL2_to_kg*areaTm(i,j))
     enddo ; enddo ; enddo
 
-    ! This block avoids using the points beyond an open boundary condition
-    ! in the accumulation of mass, but perhaps it would be unnecessary if there
-    ! were a more judicious use of masks in the loops 4 or 7 lines above.
-    if (local_open_BC) then
-      do ns=1, OBC%number_of_segments
-        segment => OBC%segment(ns)
-        if (.not. segment%on_pe .or. segment%specified) cycle
-        I=segment%HI%IsdB ; J=segment%HI%JsdB
-        if (segment%direction == OBC_DIRECTION_E) then
-          do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
-            tmp1(i+1,j,k) = 0.0
-          enddo ; enddo
-        elseif (segment%direction == OBC_DIRECTION_W) then
-          do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
-            tmp1(i,j,k) = 0.0
-          enddo ; enddo
-        elseif (segment%direction == OBC_DIRECTION_N) then
-          do k=1,nz ; do i=segment%HI%isd,segment%HI%ied
-            tmp1(i,j+1,k) = 0.0
-          enddo ; enddo
-        elseif (segment%direction == OBC_DIRECTION_S) then
-          do k=1,nz ; do i=segment%HI%isd,segment%HI%ied
-            tmp1(i,j,k) = 0.0
-          enddo ; enddo
-        endif
-      enddo
-    endif
-
     mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, EFP_sum=mass_EFP)
     do k=1,nz ; vol_lay(k) = (US%m_to_L**2*GV%H_to_Z/GV%H_to_kg_m2)*mass_lay(k) ; enddo
   else
@@ -558,19 +526,18 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_
   endif ! Boussinesq
 
   nTr_stocks = 0
-  if (present(tracer_CSp)) then
-    call call_tracer_stocks(h, Tr_stocks, G, GV, tracer_CSp, stock_names=Tr_names, &
-                            stock_units=Tr_units, num_stocks=nTr_stocks,&
-                            got_min_max=Tr_minmax_got, global_min=Tr_min, global_max=Tr_max, &
-                            xgmin=Tr_min_x, ygmin=Tr_min_y, zgmin=Tr_min_z,&
-                            xgmax=Tr_max_x, ygmax=Tr_max_y, zgmax=Tr_max_z)
-    if (nTr_stocks > 0) then
-      do m=1,nTr_stocks
-        vars(num_nc_fields+m) = var_desc(Tr_names(m), units=Tr_units(m), &
-                      longname=Tr_names(m), hor_grid='1', z_grid='1')
-      enddo
-      num_nc_fields = num_nc_fields + nTr_stocks
-    endif
+  Tr_minmax_avail(:) = .false.
+  call call_tracer_stocks(h, Tr_stocks, G, GV, tracer_CSp, stock_names=Tr_names, &
+                          stock_units=Tr_units, num_stocks=nTr_stocks,&
+                          got_min_max=Tr_minmax_avail, global_min=Tr_min, global_max=Tr_max, &
+                          xgmin=Tr_min_x, ygmin=Tr_min_y, zgmin=Tr_min_z,&
+                          xgmax=Tr_max_x, ygmax=Tr_max_y, zgmax=Tr_max_z)
+  if (nTr_stocks > 0) then
+    do m=1,nTr_stocks
+      vars(num_nc_fields+m) = var_desc(Tr_names(m), units=Tr_units(m), &
+                    longname=Tr_names(m), hor_grid='1', z_grid='1')
+    enddo
+    num_nc_fields = num_nc_fields + nTr_stocks
   endif
 
   if (CS%previous_calls == 0) then
@@ -884,7 +851,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_
          write(stdout,'("      Total ",a,": ",ES24.16,X,a)') &
               trim(Tr_names(m)), Tr_stocks(m), trim(Tr_units(m))
 
-         if (Tr_minmax_got(m)) then
+         if (Tr_minmax_avail(m)) then
            write(stdout,'(64X,"Global Min:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")"  )') &
                 Tr_min(m),Tr_min_x(m),Tr_min_y(m),Tr_min_z(m)
            write(stdout,'(64X,"Global Max:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")"  )') &

src/diagnostics/MOM_wave_speed.F90

@@ -55,7 +55,7 @@ module MOM_wave_speed
 
 !> Calculates the wave speed of the first baroclinic mode.
 subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_N2_column_fraction, &
-                      mono_N2_depth, modal_structure, better_speed_est, min_speed, wave_speed_tol)
+                      mono_N2_depth, modal_structure)
   type(ocean_grid_type),            intent(in)  :: G  !< Ocean grid structure
   type(verticalGrid_type),          intent(in)  :: GV !< Vertical grid structure
   type(unit_scale_type),            intent(in)  :: US !< A dimensional unit scaling type
@@ -76,12 +76,6 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_
                                           !! modal structure [Z ~> m].
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           optional, intent(out) :: modal_structure !< Normalized model structure [nondim]
-  logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
-                                     !! mode speed as the starting point for iterations.
-  real,    optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
-                                     !! below which 0 is returned [L T-1 ~> m s-1].
-  real,    optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
-                                     !! wave speeds [nondim]
 
   ! Local variables
   real, dimension(SZK_(GV)+1) :: &
@@ -181,10 +175,10 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_
   Z_to_pres = GV%Z_to_H * (GV%H_to_RZ * GV%g_Earth)
   use_EOS = associated(tv%eqn_of_state)
 
-  better_est = CS%better_cg1_est ; if (present(better_speed_est)) better_est = better_speed_est
+  better_est = CS%better_cg1_est
 
   if (better_est) then
-    tol_solve = CS%wave_speed_tol ; if (present(wave_speed_tol)) tol_solve = wave_speed_tol
+    tol_solve = CS%wave_speed_tol
     tol_Hfrac  = 0.1*tol_solve ; tol_merge = tol_solve / real(nz)
   else
     tol_solve = 0.001 ; tol_Hfrac  = 0.0001 ; tol_merge = 0.001
@@ -197,7 +191,7 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_
   ! worst possible oceanic case of g'H < 0.5*10m/s2*1e4m = 5.e4 m2/s2 < 1024**2*c2_scale, suggesting
   ! that c2_scale can safely be set to 1/(16*1024**2), which would decrease the stable floor on
   ! min_speed to ~6.9e-8 m/s for 90 layers or 2.33e-7 m/s for 1000 layers.
-  cg1_min2 = CS%min_speed2 ; if (present(min_speed)) cg1_min2 = min_speed**2
+  cg1_min2 = CS%min_speed2
   rescale = 1024.0**4 ; I_rescale = 1.0/rescale
   c2_scale = US%m_s_to_L_T**2 / 4096.0**2 ! Other powers of 2 give identical results.
 
@@ -638,8 +632,7 @@ subroutine tdma6(n, a, c, lam, y)
 end subroutine tdma6
 
 !> Calculates the wave speeds for the first few barolinic modes.
-subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos, better_speed_est, &
-                       min_speed, wave_speed_tol)
+subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos)
   type(ocean_grid_type),                    intent(in)  :: G !< Ocean grid structure
   type(verticalGrid_type),                  intent(in)  :: GV !< Vertical grid structure
   type(unit_scale_type),                    intent(in)  :: US !< A dimensional unit scaling type
@@ -650,12 +643,6 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos, better_spee
   type(wave_speed_CS), optional,            pointer     :: CS !< Control structure for MOM_wave_speed
   logical,             optional,            intent(in)  :: full_halos !< If true, do the calculation
                                                                       !! over the entire computational domain.
-  logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
-                                     !! mode speed as the starting point for iterations.
-  real,    optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
-                                     !! below which 0 is returned [L T-1 ~> m s-1].
-  real,    optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
-                                     !! wave speeds [nondim]
 
   ! Local variables
   real, dimension(SZK_(GV)+1) :: &
@@ -757,16 +744,13 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos, better_spee
   c2_scale = US%m_s_to_L_T**2 / 4096.0**2 ! Other powers of 2 give identical results.
 
   better_est = .false. ; if (present(CS)) better_est = CS%better_cg1_est
-  if (present(better_speed_est)) better_est = better_speed_est
   if (better_est) then
     tol_solve = 0.001 ; if (present(CS)) tol_solve = CS%wave_speed_tol
-    if (present(wave_speed_tol)) tol_solve = wave_speed_tol
     tol_Hfrac  = 0.1*tol_solve ; tol_merge = tol_solve / real(nz)
   else
     tol_solve = 0.001 ; tol_Hfrac  = 0.0001 ; tol_merge = 0.001
   endif
   cg1_min2 = 0.0 ; if (present(CS)) cg1_min2 = CS%min_speed2
-  if (present(min_speed)) cg1_min2 = min_speed**2
 
   ! Zero out all wave speeds.  Values over land or for columns that are too weakly stratified
   ! are not changed from this zero value.
@@ -1151,18 +1135,16 @@ subroutine tridiag_det(a, c, ks, ke, lam, det, ddet, row_scale)
   real,               intent(in) :: lam   !< Value subtracted from b
   real,               intent(out):: det   !< Determinant
   real,               intent(out):: ddet  !< Derivative of determinant with lam
-  real,     optional, intent(in) :: row_scale !< A scaling factor of the rows of the
+  real,               intent(in) :: row_scale !< A scaling factor of the rows of the
                                       !! matrix to limit the growth of the determinant
   ! Local variables
   real :: detKm1, detKm2   ! Cumulative value of the determinant for the previous two layers.
   real :: ddetKm1, ddetKm2 ! Derivative of the cumulative determinant with lam for the previous two layers.
   real, parameter :: rescale = 1024.0**4 ! max value of determinant allowed before rescaling
-  real :: rscl      ! A rescaling factor that is applied succesively to each row.
   real :: I_rescale ! inverse of rescale
   integer :: k      ! row (layer interface) index
 
   I_rescale = 1.0 / rescale
-  rscl = 1.0 ; if (present(row_scale)) rscl = row_scale
 
   detKm1 = 1.0 ; ddetKm1 = 0.0
   det = (a(ks)+c(ks)) - lam ; ddet = -1.0

src/framework/MOM_io.F90

@@ -6,7 +6,7 @@ module MOM_io
 use MOM_array_transform,  only : allocate_rotated_array, rotate_array
 use MOM_array_transform,  only : rotate_array_pair, rotate_vector
 use MOM_domains,          only : MOM_domain_type, domain1D, broadcast, get_domain_components
-use MOM_domains,          only : rescale_comp_data, AGRID, BGRID_NE, CGRID_NE
+use MOM_domains,          only : rescale_comp_data, num_PEs, AGRID, BGRID_NE, CGRID_NE
 use MOM_dyn_horgrid,      only : dyn_horgrid_type
 use MOM_ensemble_manager, only : get_ensemble_id
 use MOM_error_handler,    only : MOM_error, NOTE, FATAL, WARNING, is_root_PE
@@ -236,6 +236,8 @@ subroutine create_file(IO_handle, filename, vars, novars, fields, threading, tim
     IsgB = dG%IsgB ; IegB = dG%IegB ; JsgB = dG%JsgB ; JegB = dG%JegB
   endif
 
+  if (domain_set .and. (num_PEs() == 1)) thread = SINGLE_FILE
+
   one_file = .true.
   if (domain_set) one_file = (thread == SINGLE_FILE)