From 8e5a1b758902fad7e33194361cd6da8c22944676 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 13 Oct 2021 16:11:44 -0400 Subject: [PATCH 1/5] FMS1: Don't create IO domains for single-PE domain FMS mpp domain creation is done in the `clone_MD_to_d2D` function, and currently an IO domain is always created within the MPP domain. This has caused problems with single-PE runs in FMS1, where the `write_field` logic was not able to reach the part which removes halo data if an IO domain was present, and halo data was being written to the restart files. This issue arose when `PARALLEL_RESTARTFILES` was set to `True` for single-PE tests. This does not appear to be a problem with FMS2, and no action is needed by the FMS team. --- config_src/infra/FMS1/MOM_domain_infra.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/config_src/infra/FMS1/MOM_domain_infra.F90 b/config_src/infra/FMS1/MOM_domain_infra.F90 index 590637158f..7eff4597f3 100644 --- a/config_src/infra/FMS1/MOM_domain_infra.F90 +++ b/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 From 6ff0caec66ecd3817a8cded0b82a53f165584942 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Fri, 22 Oct 2021 14:58:18 -0400 Subject: [PATCH 2/5] Use FMS2 `file_exists`, remove domain args This patch uses the FMS2 `file_exists` function when using the FMS2 infra. Previously, the FMS1 version of this function was being used. This patch also removes the `mpp_domain` and `no_domain` arguments from the direct `file_exist` calls, which were not used by any known MOM6 configurations. (Nor would we expect them to be, since explicit references to FMS should not exist outside of the infra layer.) Since FMS2 does not use these arguments, their removal also creates a more meaningful interface between the two frameworks. Motivation: An issue with the FMS1 `file_exists` under the FMS2 infra was discovered in the UFS model on Hera. It was only reproducible in submitted jobs, and not for interactive jobs, and only with the GCC 9.2 compiler. (Other GCC versions were not tested.) One potential explanation is that it is related to the `save` attribute of the domain pointer, `d_ptr`. In the case above, `d_ptr` pointed to the MOM input domain for the failed cases. For the other working cases, `d_ptr` pointed to a `NULL()` value and behavior was normal. It is possible that `d_ptr` is inconsisently updated when FMS1 and FMS2 IO operations are used together, which should probably be considered undefined behavior. --- config_src/infra/FMS1/MOM_io_infra.F90 | 11 ++++------- config_src/infra/FMS2/MOM_io_infra.F90 | 12 +++++------- 2 files changed, 9 insertions(+), 14 deletions(-) diff --git a/config_src/infra/FMS1/MOM_io_infra.F90 b/config_src/infra/FMS1/MOM_io_infra.F90 index 1501f3171b..f956f9fa51 100644 --- a/config_src/infra/FMS1/MOM_io_infra.F90 +++ b/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 diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index 0b8c19d836..62a43ab99b 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/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 From 7e8d6e2f4ac3c78d536b1c7362bae4fa293db4ce Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 19 Oct 2021 09:06:56 -0400 Subject: [PATCH 3/5] +Eliminate unused arguments in diagnostics code Eliminated the unused optional OBC argument to write_energy() and several unused optional arguments to wave_speed() and wave_speeds() that are set instead via arguments to wave_speed_init() that store these values in a wave_speed_CS type. Also made the optional row_scale argument to tridiag_det() and the tracer_CSp argument to write_energy() that were always present in calls into mandatory arguments. All answers are bitwise identical, and solutions do not change. --- src/diagnostics/MOM_sum_output.F90 | 95 ++++++++++-------------------- src/diagnostics/MOM_wave_speed.F90 | 30 ++-------- 2 files changed, 37 insertions(+), 88 deletions(-) diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index d190cee7a3..5f144af4d5 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/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,")" )') & diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 6a4d9660d7..a468f36658 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/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 From a2aaad6df3bacf81ab4cfc0530cd99b26396aa28 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 13 Oct 2021 16:13:07 -0400 Subject: [PATCH 4/5] Always create single-threaded files in 1 PE runs Modified create_file to always specify that the threading will be for a single PE to do the writing to a single file if there is only a single ocean PE, to avoid some incorrect behavior inside of the FMS IO modules. This can fix restart problems in some single-processor cases with an inconsistent setting of the optional threading argument, but in all cases that ran correctly before, all answers are bitwise identical. --- src/framework/MOM_io.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 00eeb4cf89..563f9f9f8a 100644 --- a/src/framework/MOM_io.F90 +++ b/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) From 6dd1d14b3cc0cf5e8abd39f61e781e54c3e39e22 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 20 Oct 2021 10:26:42 -0400 Subject: [PATCH 5/5] +Optional arg cleanup in horizontal param code Cleaned up 13 falsely optional or unused arguments in the horizontal parameterization code, and related changes. This includes: - Made the previously optional OBC pointer arguments that were always being used in calls to 3 routines in MOM_lateral_mixing_coeffs.F90 into mandatory arguments. Because these are pointers, the deciding factor of whether to use them is really whether they are associated. - Made an internal optional argument that was always being used mandatory in 2 routines in MOM_internal_tides.F90. - Made 2 internal optional arguments that were always being used mandatory in thickness_diffuse_full(). - Eliminated the unused deta_tidal_deta argument to calc_tidal_forcing() and made the m_to_Z argument to the same routine mandatory. The former value is instead obtained by a call to tidal_sensitivity. - Eliminated 3 unused arguments and made an optional argument that was always used mandatory for find_deficit_ratios() in MOM_regularize_layers.F90. --- .../lateral/MOM_internal_tides.F90 | 17 ++-- .../lateral/MOM_lateral_mixing_coeffs.F90 | 22 +++--- .../lateral/MOM_thickness_diffuse.F90 | 43 ++++------ .../lateral/MOM_tidal_forcing.F90 | 29 +++---- .../vertical/MOM_regularize_layers.F90 | 78 +++---------------- 5 files changed, 55 insertions(+), 134 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 860b5233c9..2f9bb1f653 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -1905,7 +1905,7 @@ subroutine PPM_reconstruction_x(h_in, h_l, h_r, G, LB, simple_2nd) real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_l !< Left edge value of reconstruction (2D). real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_r !< Right edge value of reconstruction (2D). type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds. - logical, optional, intent(in) :: simple_2nd !< If true, use the arithmetic mean + logical, intent(in) :: simple_2nd !< If true, use the arithmetic mean !! energy densities as default edge values !! for a simple 2nd order scheme. ! Local variables @@ -1913,15 +1913,14 @@ subroutine PPM_reconstruction_x(h_in, h_l, h_r, G, LB, simple_2nd) real, parameter :: oneSixth = 1./6. real :: h_ip1, h_im1 real :: dMx, dMn - logical :: use_CW84, use_2nd + logical :: use_CW84 character(len=256) :: mesg ! The text of an error message integer :: i, j, isl, iel, jsl, jel, stencil - use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd isl = LB%ish-1 ; iel = LB%ieh+1 ; jsl = LB%jsh ; jel = LB%jeh ! This is the stencil of the reconstruction, not the scheme overall. - stencil = 2 ; if (use_2nd) stencil = 1 + stencil = 2 ; if (simple_2nd) stencil = 1 if ((isl-stencil < G%isd) .or. (iel+stencil > G%ied)) then write(mesg,'("In MOM_internal_tides, PPM_reconstruction_x called with a ", & @@ -1936,7 +1935,7 @@ subroutine PPM_reconstruction_x(h_in, h_l, h_r, G, LB, simple_2nd) call MOM_error(FATAL,mesg) endif - if (use_2nd) then + if (simple_2nd) then do j=jsl,jel ; do i=isl,iel h_im1 = G%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-G%mask2dT(i-1,j)) * h_in(i,j) h_ip1 = G%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-G%mask2dT(i+1,j)) * h_in(i,j) @@ -1981,7 +1980,7 @@ subroutine PPM_reconstruction_y(h_in, h_l, h_r, G, LB, simple_2nd) real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_l !< Left edge value of reconstruction (2D). real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_r !< Right edge value of reconstruction (2D). type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds. - logical, optional, intent(in) :: simple_2nd !< If true, use the arithmetic mean + logical, intent(in) :: simple_2nd !< If true, use the arithmetic mean !! energy densities as default edge values !! for a simple 2nd order scheme. ! Local variables @@ -1989,15 +1988,13 @@ subroutine PPM_reconstruction_y(h_in, h_l, h_r, G, LB, simple_2nd) real, parameter :: oneSixth = 1./6. real :: h_jp1, h_jm1 real :: dMx, dMn - logical :: use_2nd character(len=256) :: mesg ! The text of an error message integer :: i, j, isl, iel, jsl, jel, stencil - use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd isl = LB%ish ; iel = LB%ieh ; jsl = LB%jsh-1 ; jel = LB%jeh+1 ! This is the stencil of the reconstruction, not the scheme overall. - stencil = 2 ; if (use_2nd) stencil = 1 + stencil = 2 ; if (simple_2nd) stencil = 1 if ((isl < G%isd) .or. (iel > G%ied)) then write(mesg,'("In MOM_internal_tides, PPM_reconstruction_y called with a ", & @@ -2012,7 +2009,7 @@ subroutine PPM_reconstruction_y(h_in, h_l, h_r, G, LB, simple_2nd) call MOM_error(FATAL,mesg) endif - if (use_2nd) then + if (simple_2nd) then do j=jsl,jel ; do i=isl,iel h_jm1 = G%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-G%mask2dT(i,j-1)) * h_in(i,j) h_jp1 = G%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-G%mask2dT(i,j+1)) * h_in(i,j) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index d0df4b81ba..9306233112 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -451,7 +451,7 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables real, intent(in) :: dt !< Time increment [T ~> s] type(VarMix_CS), pointer :: CS !< Variable mixing coefficients - type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. ! Local variables real, dimension(SZI_(G), SZJ_(G),SZK_(GV)+1) :: & e ! The interface heights relative to mean sea level [Z ~> m]. @@ -477,10 +477,10 @@ subroutine calc_slope_functions(h, tv, dt, G, GV, US, CS, OBC) if (CS%use_stored_slopes) then call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*CS%kappa_smooth, & CS%slope_x, CS%slope_y, N2_u=N2_u, N2_v=N2_v, halo=1, OBC=OBC) - call calc_Visbeck_coeffs_old(h, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, US, CS, OBC=OBC) + call calc_Visbeck_coeffs_old(h, CS%slope_x, CS%slope_y, N2_u, N2_v, G, GV, US, CS, OBC) else !call calc_isoneutral_slopes(G, GV, h, e, tv, dt*CS%kappa_smooth, CS%slope_x, CS%slope_y) - call calc_slope_functions_using_just_e(h, G, GV, US, CS, e, .true., OBC=OBC) + call calc_slope_functions_using_just_e(h, G, GV, US, CS, e, .true., OBC) endif endif endif @@ -515,7 +515,7 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C !! at v-points [L2 Z-2 T-2 ~> s-2] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(VarMix_CS), pointer :: CS !< Variable mixing coefficients - type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. ! Local variables real :: S2 ! Interface slope squared [nondim] @@ -543,10 +543,10 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C local_open_u_BC = .false. local_open_v_BC = .false. - if (present(OBC)) then ; if (associated(OBC)) then + if (associated(OBC)) then local_open_u_BC = OBC%open_u_BCs_exist_globally local_open_v_BC = OBC%open_v_BCs_exist_globally - endif ; endif + endif S2max = CS%Visbeck_S_max**2 @@ -673,8 +673,8 @@ subroutine calc_Eady_growth_rate_2D(CS, G, GV, US, OBC, h, e, dzu, dzv, dzSxN, d 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 - type(ocean_OBC_type), pointer, intent(in) :: OBC !< Open boundaries control structure. -real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Interface height [Z ~> m] + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Interface height [Z ~> m] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface height [Z ~> m] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: dzu !< dz at u-points [Z ~> m] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: dzv !< dz at v-points [Z ~> m] @@ -859,7 +859,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface position [Z ~> m] logical, intent(in) :: calculate_slopes !< If true, calculate slopes !! internally otherwise use slopes stored in CS - type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. ! Local variables real :: E_x(SZIB_(G), SZJ_(G)) ! X-slope of interface at u points [nondim] (for diagnostics) real :: E_y(SZI_(G), SZJB_(G)) ! Y-slope of interface at v points [nondim] (for diagnostics) @@ -890,10 +890,10 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop local_open_u_BC = .false. local_open_v_BC = .false. - if (present(OBC)) then ; if (associated(OBC)) then + if (associated(OBC)) then local_open_u_BC = OBC%open_u_BCs_exist_globally local_open_v_BC = OBC%open_v_BCs_exist_globally - endif ; endif + endif one_meter = 1.0 * GV%m_to_H h_neglect = GV%H_subroundoff diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index daeb64fab9..78425676b1 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -580,11 +580,11 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV real, intent(in) :: dt !< Time increment [T ~> s] type(MEKE_type), pointer :: MEKE !< MEKE control structure type(thickness_diffuse_CS), pointer :: CS !< Control structure for thickness diffusion - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), optional, intent(in) :: int_slope_u !< Ratio that determine how much of + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: int_slope_u !< Ratio that determine how much of !! the isopycnal slopes are taken directly from !! the interface slopes without consideration of !! density gradients [nondim]. - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), optional, intent(in) :: int_slope_v !< Ratio that determine how much of + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: int_slope_v !< Ratio that determine how much of !! the isopycnal slopes are taken directly from !! the interface slopes without consideration of !! density gradients [nondim]. @@ -694,7 +694,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV real, dimension(SZIB_(G), SZJ_(G),SZK_(GV)+1) :: diag_sfn_x, diag_sfn_unlim_x ! Diagnostics real, dimension(SZI_(G), SZJB_(G),SZK_(GV)+1) :: diag_sfn_y, diag_sfn_unlim_y ! Diagnostics - logical :: present_int_slope_u, present_int_slope_v logical :: present_slope_x, present_slope_y, calc_derivatives integer, dimension(2) :: EOSdom_u ! The shifted i-computational domain to use for equation of ! state calculations at u-points. @@ -715,8 +714,6 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV N2_floor = CS%N2_floor*US%Z_to_L**2 use_EOS = associated(tv%eqn_of_state) - present_int_slope_u = PRESENT(int_slope_u) - present_int_slope_v = PRESENT(int_slope_v) present_slope_x = PRESENT(slope_x) present_slope_y = PRESENT(slope_y) use_Stanley = CS%Stanley_det_coeff >= 0. @@ -818,12 +815,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV EOSdom_u(1) = (is-1) - (G%IsdB-1) ; EOSdom_u(2) = ie - (G%IsdB-1) !$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, & -!$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect, & -!$OMP I_slope_max2,h_neglect2,present_int_slope_u, & -!$OMP int_slope_u,KH_u,uhtot,h_frac,h_avail_rsum, & -!$OMP uhD,h_avail,G_scale,Work_u,CS,slope_x,cg1, & -!$OMP diag_sfn_x, diag_sfn_unlim_x,N2_floor,EOSdom_u, & -!$OMP use_stanley, Tsgs2, & +!$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,I_slope_max2, & +!$OMP h_neglect2,int_slope_u,KH_u,uhtot,h_frac,h_avail_rsum, & +!$OMP uhD,h_avail,G_scale,Work_u,CS,slope_x,cg1,diag_sfn_x, & +!$OMP diag_sfn_unlim_x,N2_floor,EOSdom_u,use_stanley, Tsgs2, & !$OMP present_slope_x,G_rho0,Slope_x_PE,hN2_x_PE) & !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, & !$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, & @@ -941,11 +936,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! Adjust real slope by weights that bias towards slope of interfaces ! that ignore density gradients along layers. - if (present_int_slope_u) then - Slope = (1.0 - int_slope_u(I,j,K)) * Slope + & - int_slope_u(I,j,K) * ((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)) - slope2_Ratio_u(I,K) = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio_u(I,K) - endif + Slope = (1.0 - int_slope_u(I,j,K)) * Slope + & + int_slope_u(I,j,K) * ((e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)) + slope2_Ratio_u(I,K) = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio_u(I,K) Slope_x_PE(I,j,k) = MIN(Slope,CS%slope_max) hN2_x_PE(I,j,k) = hN2_u(I,K) @@ -1090,12 +1083,10 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! Calculate the meridional fluxes and gradients. EOSdom_v(:) = EOS_domain(G%HI) !$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, & -!$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect, & -!$OMP I_slope_max2,h_neglect2,present_int_slope_v, & -!$OMP int_slope_v,KH_v,vhtot,h_frac,h_avail_rsum, & -!$OMP vhD,h_avail,G_scale,Work_v,CS,slope_y,cg1, & -!$OMP diag_sfn_y,diag_sfn_unlim_y,N2_floor,EOSdom_v,& -!$OMP use_stanley, Tsgs2, & +!$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect,I_slope_max2, & +!$OMP h_neglect2,int_slope_v,KH_v,vhtot,h_frac,h_avail_rsum, & +!$OMP vhD,h_avail,G_scale,Work_v,CS,slope_y,cg1,diag_sfn_y, & +!$OMP diag_sfn_unlim_y,N2_floor,EOSdom_v,use_stanley,Tsgs2, & !$OMP present_slope_y,G_rho0,Slope_y_PE,hN2_y_PE) & !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, & !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, & @@ -1211,11 +1202,9 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV ! Adjust real slope by weights that bias towards slope of interfaces ! that ignore density gradients along layers. - if (present_int_slope_v) then - Slope = (1.0 - int_slope_v(i,J,K)) * Slope + & - int_slope_v(i,J,K) * ((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)) - slope2_Ratio_v(i,K) = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio_v(i,K) - endif + Slope = (1.0 - int_slope_v(i,J,K)) * Slope + & + int_slope_v(i,J,K) * ((e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)) + slope2_Ratio_v(i,K) = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio_v(i,K) Slope_y_PE(i,J,k) = MIN(Slope,CS%slope_max) hN2_y_PE(i,J,k) = hN2_v(i,K) diff --git a/src/parameterizations/lateral/MOM_tidal_forcing.F90 b/src/parameterizations/lateral/MOM_tidal_forcing.F90 index 862b622d56..0d92a14d2a 100644 --- a/src/parameterizations/lateral/MOM_tidal_forcing.F90 +++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -579,7 +579,7 @@ end subroutine tidal_forcing_sensitivity !! height. For now, eta and eta_tidal are both geopotential heights in depth !! units, but probably the input for eta should really be replaced with the !! column mass anomalies. -subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, deta_tidal_deta, m_to_Z) +subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, m_to_Z) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(time_type), intent(in) :: Time !< The time for the caluculation. real, dimension(SZI_(G),SZJ_(G)), intent(in) :: eta !< The sea surface height anomaly from @@ -588,16 +588,12 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, deta_tidal_deta, m_to !! anomalies [Z ~> m]. type(tidal_forcing_CS), pointer :: CS !< The control structure returned by a !! previous call to tidal_forcing_init. - real, optional, intent(out) :: deta_tidal_deta !< The partial derivative of - !! eta_tidal with the local value of - !! eta [nondim]. - real, optional, intent(in) :: m_to_Z !< A scaling factor from m to the units of eta. + real, intent(in) :: m_to_Z !< A scaling factor from m to the units of eta. ! Local variables real :: now ! The relative time in seconds. real :: amp_cosomegat, amp_sinomegat real :: cosomegat, sinomegat - real :: m_Z ! A scaling factor from m to depth units. real :: eta_prop ! The nondimenional constant of proportionality beteen eta and eta_tidal. integer :: i, j, c, m, is, ie, js, je, Isq, Ieq, Jsq, Jeq is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -622,21 +618,14 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, deta_tidal_deta, m_to eta_prop = 0.0 endif - if (present(deta_tidal_deta)) then - deta_tidal_deta = eta_prop - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ; eta_tidal(i,j) = 0.0 ; enddo ; enddo - else - do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta_tidal(i,j) = eta_prop*eta(i,j) - enddo ; enddo - endif - - m_Z = 1.0 ; if (present(m_to_Z)) m_Z = m_to_Z + do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 + eta_tidal(i,j) = eta_prop*eta(i,j) + enddo ; enddo do c=1,CS%nc m = CS%struct(c) - amp_cosomegat = m_Z*CS%amp(c)*CS%love_no(c) * cos(CS%freq(c)*now + CS%phase0(c)) - amp_sinomegat = m_Z*CS%amp(c)*CS%love_no(c) * sin(CS%freq(c)*now + CS%phase0(c)) + amp_cosomegat = m_to_Z*CS%amp(c)*CS%love_no(c) * cos(CS%freq(c)*now + CS%phase0(c)) + amp_sinomegat = m_to_Z*CS%amp(c)*CS%love_no(c) * sin(CS%freq(c)*now + CS%phase0(c)) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 eta_tidal(i,j) = eta_tidal(i,j) + (amp_cosomegat*CS%cos_struct(i,j,m) + & amp_sinomegat*CS%sin_struct(i,j,m)) @@ -647,7 +636,7 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, deta_tidal_deta, m_to cosomegat = cos(CS%freq(c)*now) sinomegat = sin(CS%freq(c)*now) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta_tidal(i,j) = eta_tidal(i,j) + m_Z*CS%ampsal(i,j,c) * & + eta_tidal(i,j) = eta_tidal(i,j) + m_to_Z*CS%ampsal(i,j,c) * & (cosomegat*CS%cosphasesal(i,j,c) + sinomegat*CS%sinphasesal(i,j,c)) enddo ; enddo enddo ; endif @@ -656,7 +645,7 @@ subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, CS, deta_tidal_deta, m_to cosomegat = cos(CS%freq(c)*now) sinomegat = sin(CS%freq(c)*now) do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 - eta_tidal(i,j) = eta_tidal(i,j) - m_Z*CS%SAL_SCALAR*CS%amp_prev(i,j,c) * & + eta_tidal(i,j) = eta_tidal(i,j) - m_to_Z*CS%SAL_SCALAR*CS%amp_prev(i,j,c) * & (cosomegat*CS%cosphase_prev(i,j,c) + sinomegat*CS%sinphase_prev(i,j,c)) enddo ; enddo enddo ; endif diff --git a/src/parameterizations/vertical/MOM_regularize_layers.F90 b/src/parameterizations/vertical/MOM_regularize_layers.F90 index f67fb48fc7..af92e522a2 100644 --- a/src/parameterizations/vertical/MOM_regularize_layers.F90 +++ b/src/parameterizations/vertical/MOM_regularize_layers.F90 @@ -219,7 +219,7 @@ subroutine regularize_surface(h, tv, dt, ea, eb, G, GV, US, CS) e(i,j,K+1) = e(i,j,K) - h(i,j,k) enddo ; enddo ; enddo - call find_deficit_ratios(e, def_rat_u, def_rat_v, G, GV, CS, h=h) + call find_deficit_ratios(e, def_rat_u, def_rat_v, G, GV, CS, h) ! Determine which columns are problematic do j=js,je ; do_j(j) = .false. ; enddo @@ -612,8 +612,7 @@ end subroutine regularize_surface !! thickness at velocity points differ from the arithmetic means, relative to !! the the arithmetic means, after eliminating thickness variations that are !! solely due to topography and aggregating all interior layers into one. -subroutine find_deficit_ratios(e, def_rat_u, def_rat_v, G, GV, CS, & - def_rat_u_2lay, def_rat_v_2lay, halo, h) +subroutine find_deficit_ratios(e, def_rat_u, def_rat_v, G, GV, CS, h) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & @@ -626,30 +625,18 @@ subroutine find_deficit_ratios(e, def_rat_u, def_rat_v, G, GV, CS, & !! [nondim]. type(regularize_layers_CS), pointer :: CS !< The control structure returned by a !! previous call to regularize_layers_init. - real, dimension(SZIB_(G),SZJ_(G)), & - optional, intent(out) :: def_rat_u_2lay !< The thickness deficit ratio at u - !! points when the mixed and buffer layers - !! are aggregated into 1 layer [nondim]. - real, dimension(SZI_(G),SZJB_(G)), & - optional, intent(out) :: def_rat_v_2lay !< The thickness deficit ratio at v - !! pointswhen the mixed and buffer layers - !! are aggregated into 1 layer [nondim]. - integer, optional, intent(in) :: halo !< An extra-wide halo size, 0 by default. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - optional, intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. - !! If h is not present, vertical differences - !! in interface heights are used instead. + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + ! Local variables real, dimension(SZIB_(G),SZJ_(G)) :: & h_def_u, & ! The vertically summed thickness deficits at u-points [H ~> m or kg m-2]. - h_norm_u, & ! The vertically summed arithmetic mean thickness by which + h_norm_u ! The vertically summed arithmetic mean thickness by which ! h_def_u is normalized [H ~> m or kg m-2]. - h_def2_u real, dimension(SZI_(G),SZJB_(G)) :: & h_def_v, & ! The vertically summed thickness deficits at v-points [H ~> m or kg m-2]. - h_norm_v, & ! The vertically summed arithmetic mean thickness by which + h_norm_v ! The vertically summed arithmetic mean thickness by which ! h_def_v is normalized [H ~> m or kg m-2]. - h_def2_v real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: Hmix_min ! A local copy of CS%Hmix_min [H ~> m or kg m-2]. @@ -657,9 +644,6 @@ subroutine find_deficit_ratios(e, def_rat_u, def_rat_v, G, GV, CS, & integer :: i, j, k, is, ie, js, je, nz, nkmb is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke - if (present(halo)) then - is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo - endif nkmb = GV%nk_rho_varies h_neglect = GV%H_subroundoff Hmix_min = CS%Hmix_min @@ -677,22 +661,8 @@ subroutine find_deficit_ratios(e, def_rat_u, def_rat_v, G, GV, CS, & h_def_u(I,j) = 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect) h_norm_u(I,j) = 0.5*(h1+h2) enddo ; enddo - if (present(def_rat_u_2lay)) then ; do j=js,je ; do I=is-1,ie - ! This is a particular metric of the aggregation into two layers. - h1 = e(i,j,1)-e(i,j,nkmb+1) ; h2 = e(i+1,j,1)-e(i+1,j,nkmb+1) - if (e(i,j,nkmb+1) < e(i+1,j,nz+1)) then - if (h1 > h2) h1 = max(e(i,j,1)-e(i+1,j,nz+1), h2) - elseif (e(i+1,j,nkmb+1) < e(i,j,nz+1)) then - if (h2 > h1) h2 = max(e(i+1,j,1)-e(i,j,nz+1), h1) - endif - h_def2_u(I,j) = h_def_u(I,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect) - enddo ; enddo ; endif do k=1,nkmb ; do j=js,je ; do I=is-1,ie - if (present(h)) then - h1 = h(i,j,k) ; h2 = h(i+1,j,k) - else - h1 = e(i,j,K)-e(i,j,K+1) ; h2 = e(i+1,j,K)-e(i+1,j,K+1) - endif + h1 = h(i,j,k) ; h2 = h(i+1,j,k) ! Thickness deficits can not arise simply because a layer's bottom is bounded ! by the bathymetry. if (e(i,j,K+1) < e(i+1,j,nz+1)) then @@ -703,15 +673,10 @@ subroutine find_deficit_ratios(e, def_rat_u, def_rat_v, G, GV, CS, & h_def_u(I,j) = h_def_u(I,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect) h_norm_u(I,j) = h_norm_u(I,j) + 0.5*(h1+h2) enddo ; enddo ; enddo - if (present(def_rat_u_2lay)) then ; do j=js,je ; do I=is-1,ie - def_rat_u(I,j) = G%mask2dCu(I,j) * h_def_u(I,j) / & - (max(Hmix_min, h_norm_u(I,j)) + h_neglect) - def_rat_u_2lay(I,j) = G%mask2dCu(I,j) * h_def2_u(I,j) / & - (max(Hmix_min, h_norm_u(I,j)) + h_neglect) - enddo ; enddo ; else ; do j=js,je ; do I=is-1,ie + do j=js,je ; do I=is-1,ie def_rat_u(I,j) = G%mask2dCu(I,j) * h_def_u(I,j) / & (max(Hmix_min, h_norm_u(I,j)) + h_neglect) - enddo ; enddo ; endif + enddo ; enddo ! Determine which meridional faces are problematic. do J=js-1,je ; do i=is,ie @@ -726,22 +691,8 @@ subroutine find_deficit_ratios(e, def_rat_u, def_rat_v, G, GV, CS, & h_def_v(i,J) = 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect) h_norm_v(i,J) = 0.5*(h1+h2) enddo ; enddo - if (present(def_rat_v_2lay)) then ; do J=js-1,je ; do i=is,ie - ! This is a particular metric of the aggregation into two layers. - h1 = e(i,j,1)-e(i,j,nkmb+1) ; h2 = e(i,j+1,1)-e(i,j+1,nkmb+1) - if (e(i,j,nkmb+1) < e(i,j+1,nz+1)) then - if (h1 > h2) h1 = max(e(i,j,1)-e(i,j+1,nz+1), h2) - elseif (e(i,j+1,nkmb+1) < e(i,j,nz+1)) then - if (h2 > h1) h2 = max(e(i,j+1,1)-e(i,j,nz+1), h1) - endif - h_def2_v(i,J) = h_def_v(i,J) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect) - enddo ; enddo ; endif do k=1,nkmb ; do J=js-1,je ; do i=is,ie - if (present(h)) then - h1 = h(i,j,k) ; h2 = h(i,j+1,k) - else - h1 = e(i,j,K)-e(i,j,K+1) ; h2 = e(i,j+1,K)-e(i,j+1,K+1) - endif + h1 = h(i,j,k) ; h2 = h(i,j+1,k) ! Thickness deficits can not arise simply because a layer's bottom is bounded ! by the bathymetry. if (e(i,j,K+1) < e(i,j+1,nz+1)) then @@ -752,15 +703,10 @@ subroutine find_deficit_ratios(e, def_rat_u, def_rat_v, G, GV, CS, & h_def_v(i,J) = h_def_v(i,J) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect) h_norm_v(i,J) = h_norm_v(i,J) + 0.5*(h1+h2) enddo ; enddo ; enddo - if (present(def_rat_v_2lay)) then ; do J=js-1,je ; do i=is,ie - def_rat_v(i,J) = G%mask2dCv(i,J) * h_def_v(i,J) / & - (max(Hmix_min, h_norm_v(i,J)) + h_neglect) - def_rat_v_2lay(i,J) = G%mask2dCv(i,J) * h_def2_v(i,J) / & - (max(Hmix_min, h_norm_v(i,J)) + h_neglect) - enddo ; enddo ; else ; do J=js-1,je ; do i=is,ie + do J=js-1,je ; do i=is,ie def_rat_v(i,J) = G%mask2dCv(i,J) * h_def_v(i,J) / & (max(Hmix_min, h_norm_v(i,J)) + h_neglect) - enddo ; enddo ; endif + enddo ; enddo end subroutine find_deficit_ratios