From ed1bd5f4a50d8a086bfa8be55ae491b2c1bbf997 Mon Sep 17 00:00:00 2001 From: Hemant Khatri Date: Thu, 30 Apr 2020 18:34:23 -0400 Subject: [PATCH 01/40] "Barotropic velocity tendency" diagnostic added --- src/core/MOM_barotropic.F90 | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 5f97f5933a..6eef429f92 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -294,6 +294,7 @@ module MOM_barotropic integer :: id_visc_rem_u = -1, id_visc_rem_v = -1, id_eta_cor = -1 integer :: id_ubt = -1, id_vbt = -1, id_eta_bt = -1, id_ubtav = -1, id_vbtav = -1 integer :: id_ubt_st = -1, id_vbt_st = -1, id_eta_st = -1 + integer :: id_ubtdt = -1, id_vbtdt = -1 integer :: id_ubt_hifreq = -1, id_vbt_hifreq = -1, id_eta_hifreq = -1 integer :: id_uhbt_hifreq = -1, id_vhbt_hifreq = -1, id_eta_pred_hifreq = -1 integer :: id_gtotn = -1, id_gtots = -1, id_gtote = -1, id_gtotw = -1 @@ -479,6 +480,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, real :: q(SZIBW_(CS),SZJBW_(CS)) ! A pseudo potential vorticity [T-1 Z-1 ~> s-1 m-1]. real, dimension(SZIBW_(CS),SZJW_(CS)) :: & ubt, & ! The zonal barotropic velocity [L T-1 ~> m s-1]. + ubt_st, & ! The zonal barotropic velocity at the start of timestep [L T-1 ~> m s-1]. + ubt_dt, & ! The zonal barotropic velocity tendency [L T-2 ~> m s-2]. bt_rem_u, & ! The fraction of the barotropic zonal velocity that remains ! after a time step, the remainder being lost to bottom drag. ! bt_rem_u is a nondimensional number between 0 and 1. @@ -512,6 +515,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! spacing [H L ~> m2 or kg m-1]. real, dimension(SZIW_(CS),SZJBW_(CS)) :: & vbt, & ! The meridional barotropic velocity [L T-1 ~> m s-1]. + vbt_st, & ! The meridional barotropic velocity at the start of timestep [L T-1 ~> m s-1]. + vbt_dt, & ! The meridional barotropic velocity tendency [L T-2 ~> m s-2]. bt_rem_v, & ! The fraction of the barotropic meridional velocity that ! remains after a time step, the rest being lost to bottom ! drag. bt_rem_v is a nondimensional number between 0 and 1. @@ -1479,6 +1484,14 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, call uvchksum("BT visc_rem_[uv]", visc_rem_u, visc_rem_v, G%HI, haloshift=1) endif + if (CS%id_ubtdt > 0) then + ubt_st(:,:) = ubt(:,:) + endif + + if (CS%id_vbtdt > 0) then + vbt_st(:,:) = vbt(:,:) + endif + if (query_averaging_enabled(CS%diag)) then if (CS%id_eta_st > 0) call post_data(CS%id_eta_st, eta(isd:ied,jsd:jed), CS%diag) if (CS%id_ubt_st > 0) call post_data(CS%id_ubt_st, ubt(IsdB:IedB,jsd:jed), CS%diag) @@ -2201,6 +2214,19 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, enddo ; enddo call post_data(CS%id_Corv_bt, Corv_bt_sum(isd:ied,JsdB:JedB), CS%diag) endif + if (CS%id_ubtdt > 0) then + do j=js,je ; do I=is-1,ie + ubt_dt(I,j) = (ubt_wtd(I,j) - ubt_st(I,j))/dt + enddo ; enddo + call post_data(CS%id_ubtdt, ubt_dt(IsdB:IedB,jsd:jed), CS%diag) + endif + if (CS%id_vbtdt > 0) then + do J=js-1,je ; do i=is,ie + vbt_dt(i,J) = (vbt_wtd(i,J) - vbt_st(i,J))/dt + enddo ; enddo + call post_data(CS%id_vbtdt, vbt_dt(isd:ied,JsdB:JedB), CS%diag) + endif + if (CS%id_ubtforce > 0) call post_data(CS%id_ubtforce, BT_force_u(IsdB:IedB,jsd:jed), CS%diag) if (CS%id_vbtforce > 0) call post_data(CS%id_vbtforce, BT_force_v(isd:ied,JsdB:JedB), CS%diag) if (CS%id_uaccel > 0) call post_data(CS%id_uaccel, u_accel_bt(IsdB:IedB,jsd:jed), CS%diag) @@ -4203,6 +4229,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, 'Barotropic zonal acceleration from baroclinic terms', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_vbtforce = register_diag_field('ocean_model', 'vbtforce', diag%axesCv1, Time, & 'Barotropic meridional acceleration from baroclinic terms', 'm s-2', conversion=US%L_T2_to_m_s2) + CS%id_ubtdt = register_diag_field('ocean_model', 'ubt_dt', diag%axesCu1, Time, & + 'Barotropic zonal acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) + CS%id_vbtdt = register_diag_field('ocean_model', 'vbt_dt', diag%axesCv1, Time, & + 'Barotropic meridional acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) CS%id_eta_bt = register_diag_field('ocean_model', 'eta_bt', diag%axesT1, Time, & 'Barotropic end SSH', thickness_units, conversion=GV%H_to_m) From e05054ef10abd72772a4593f1fa14413613803f1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 7 May 2020 07:57:40 -0400 Subject: [PATCH 02/40] +Fix units in comments and parameter descriptions Corrected the reported units for several parameters and the units in a number of comments. This changes some comments in MOM_parameter_doc files. All answers are bitwise identical. --- src/core/MOM_isopycnal_slopes.F90 | 4 ++-- src/core/MOM_open_boundary.F90 | 24 ++++++++----------- .../MOM_state_initialization.F90 | 11 ++++----- .../vertical/MOM_tidal_mixing.F90 | 2 +- 4 files changed, 18 insertions(+), 23 deletions(-) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 53f0d59294..fa60fb821d 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -242,7 +242,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & slope_x(I,j,K) = 0.0 endif - if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of Brunt-Vaisala frequency [s-2] + if (present_N2_u) N2_u(I,j,k) = G_Rho0 * drdz * G%mask2dCu(I,j) ! Square of buoyancy frequency [T-2 ~> s-2] else ! With .not.use_EOS, the layers are constant density. slope_x(I,j,K) = (Z_to_L*(e(i,j,K)-e(i+1,j,K))) * G%IdxCu(I,j) @@ -328,7 +328,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & slope_y(i,J,K) = 0.0 endif - if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of Brunt-Vaisala frequency [s-2] + if (present_N2_v) N2_v(i,J,k) = G_Rho0 * drdz * G%mask2dCv(i,J) ! Square of buoyancy frequency [T-2 ~> s-2] else ! With .not.use_EOS, the layers are constant density. slope_y(i,J,K) = (Z_to_L*(e(i,j,K)-e(i,j+1,K))) * G%IdyCv(i,J) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index ffede1c0c2..5b6dc168f4 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -265,10 +265,8 @@ module MOM_open_boundary !! velocities (or speed of characteristics) at the !! new time level (1) or the running mean (0) for velocities. !! Valid values range from 0 to 1, with a default of 0.3. - real :: rx_max !< The maximum magnitude of the baroclinic radiation - !! velocity (or speed of characteristics) [m s-1]. The - !! default value is 10 m s-1. - !### The description above seems inconsistent with the code, and the units should be [nondim]. + real :: rx_max !< The maximum magnitude of the baroclinic radiation velocity (or speed of + !! characteristics) in units of grid points per timestep [nondim]. logical :: OBC_pe !< Is there an open boundary on this tile? type(remapping_CS), pointer :: remap_CS !< ALE remapping control structure for segments only type(OBC_registry_type), pointer :: OBC_Reg => NULL() !< Registry type for boundaries @@ -500,13 +498,11 @@ subroutine open_boundary_config(G, US, param_file, OBC) call initialize_segment_data(G, OBC, param_file) if (open_boundary_query(OBC, apply_open_OBC=.true.)) then - !### I think that OBC%rx_max as used is actually nondimensional, with effective - ! units of grid points per time step. call get_param(param_file, mdl, "OBC_RADIATION_MAX", OBC%rx_max, & - "The maximum magnitude of the baroclinic radiation "//& - "velocity (or speed of characteristics). This is only "//& + "The maximum magnitude of the baroclinic radiation velocity (or speed of "//& + "characteristics), in gridpoints per timestep. This is only "//& "used if one of the open boundary segments is using Orlanski.", & - units="m s-1", default=10.0) !### Should the units here be "nondim"? + units="nondim", default=10.0) !### Should the default be changed to 1.0? call get_param(param_file, mdl, "OBC_RAD_VEL_WT", OBC%gamma_uv, & "The relative weighting for the baroclinic radiation "//& "velocities (or speed of characteristics) at the new "//& @@ -1804,9 +1800,9 @@ subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv, US) I=segment%HI%IsdB do j=segment%HI%jsd,segment%HI%jed if (segment%direction == OBC_DIRECTION_E) then - areaCu(I,j) = G%areaT(i,j) ! Both of these are in [L2] + areaCu(I,j) = G%areaT(i,j) ! Both of these are in [L2 ~> m2] else ! West - areaCu(I,j) = G%areaT(i+1,j) ! Both of these are in [L2] + areaCu(I,j) = G%areaT(i+1,j) ! Both of these are in [L2 ~> m2] endif enddo else @@ -1814,9 +1810,9 @@ subroutine open_boundary_impose_land_mask(OBC, G, areaCu, areaCv, US) J=segment%HI%JsdB do i=segment%HI%isd,segment%HI%ied if (segment%direction == OBC_DIRECTION_S) then - areaCv(i,J) = G%areaT(i,j+1) ! Both of these are in [L2] + areaCv(i,J) = G%areaT(i,j+1) ! Both of these are in [L2 ~> m2] else ! North - areaCu(i,J) = G%areaT(i,j) ! Both of these are in [L2] + areaCu(i,J) = G%areaT(i,j) ! Both of these are in [L2 ~> m2] endif enddo endif @@ -1906,7 +1902,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, US, dt) real :: dhdt, dhdx, dhdy ! One-point differences in time or space [L T-1 ~> m s-1] real :: gamma_u, gamma_2 ! Fractional weightings of new values [nondim] real :: tau ! A local nudging timescale [T ~> s] - real :: rx_max, ry_max ! coefficients for radiation [nondim] or [L2 T-2 ~> m2 s-2] + real :: rx_max, ry_max ! coefficients for radiation [nondim] real :: rx_new, rx_avg ! coefficients for radiation [nondim] or [L2 T-2 ~> m2 s-2] real :: ry_new, ry_avg ! coefficients for radiation [nondim] or [L2 T-2 ~> m2 s-2] real :: cff_new, cff_avg ! denominator in oblique [L2 T-2 ~> m2 s-2] diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index beeaf6e46a..07d928d76b 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -1108,11 +1108,11 @@ subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read_params) just_read = .false. ; if (present(just_read_params)) just_read = just_read_params call get_param(PF, mdl, "SURFACE_PRESSURE_FILE", p_surf_file, & - "The initial condition file for the surface height.", & + "The initial condition file for the surface pressure exerted by ice.", & fail_if_missing=.not.just_read, do_not_log=just_read) call get_param(PF, mdl, "SURFACE_PRESSURE_VAR", p_surf_var, & - "The initial condition variable for the surface height.", & - units="kg m-2", default="", do_not_log=just_read) !### The units here should be Pa? + "The initial condition variable for the surface pressure exerted by ice.", & + units="Pa", default="", do_not_log=just_read) call get_param(PF, mdl, "INPUTDIR", inputdir, default=".", do_not_log=.true.) filename = trim(slasher(inputdir))//trim(p_surf_file) if (.not.just_read) call log_param(PF, mdl, "!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename) @@ -2177,14 +2177,13 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, G, GV, US, PF, just_read_param ! to the North/South Pole past the limits of the input data, they are extrapolated using the average ! value at the northernmost/southernmost latitude. - call horiz_interp_and_extrap_tracer(tfilename, potemp_var, 1.0, 1, & G, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, & - tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018,ongrid=pre_gridded) + tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018, ongrid=pre_gridded) call horiz_interp_and_extrap_tracer(sfilename, salin_var, 1.0, 1, & G, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, reentrant_x, & - tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018,ongrid=pre_gridded) + tripolar_n, homogenize, m_to_Z=US%m_to_Z, answers_2018=hor_regrid_answers_2018, ongrid=pre_gridded) kd = size(z_in,1) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 27b316e144..a788a964f6 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -590,7 +590,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, CS) if (CS%use_CVMix_tidal) then CS%id_N2_int = register_diag_field('ocean_model','N2_int',diag%axesTi,Time, & - 'Bouyancy frequency squared, at interfaces', 's-2') !###, conversion=US%s_to_T**2) + 'Bouyancy frequency squared, at interfaces', 's-2', conversion=US%s_to_T**2) !> TODO: add units CS%id_Simmons_coeff = register_diag_field('ocean_model','Simmons_coeff',diag%axesT1,Time, & 'time-invariant portion of the tidal mixing coefficient using the Simmons', '') From a041a95f1c9d33409c822802f4e3978603e70705 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 7 May 2020 08:04:44 -0400 Subject: [PATCH 03/40] (*)Corrected loop bounds in build_adapt_column Corrected the upper loop bound in a call to calculate_density_derivs in build_adapt_column. This code is not yet tested, so all answers are bitwise identical in the MOM6-examples test cases. --- src/ALE/coord_adapt.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ALE/coord_adapt.F90 b/src/ALE/coord_adapt.F90 index 42ae0ee245..8fa4b09fc5 100644 --- a/src/ALE/coord_adapt.F90 +++ b/src/ALE/coord_adapt.F90 @@ -209,7 +209,7 @@ subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNex ! a positive curvature means we're too light relative to adjacent columns, ! so del2sigma needs to be positive too (push the interface deeper) call calculate_density_derivs(tInt(i,j,:), sInt(i,j,:), zInt(i,j,:) * (GV%H_to_RZ * GV%g_Earth), & - alpha, beta, tv%eqn_of_state, (/1,nz/) ) !### This should be (/1,nz+1/) - see 25 lines below. + alpha, beta, tv%eqn_of_state, (/1,nz+1/) ) do K = 2, nz ! TODO make lower bound here configurable dh_d2s(K) = del2sigma(K) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / & From 3e28d61258712e66cdc3b3dd2ae20c007c392dfa Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 7 May 2020 08:05:26 -0400 Subject: [PATCH 04/40] Corrected rescaling with CONST_SEA_LEVEL Corrected parentheses in calculating the balancing_flux with an ice-shelf. This corrects the dimensional rescaling with the CONST_SEA_LEVEL option within the ice-shelf code, but does not change answers in any MOM6-examples test cases. --- src/ice_shelf/MOM_ice_shelf.F90 | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 891d6b3ea7..026745ee91 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1051,10 +1051,10 @@ subroutine add_shelf_flux(G, US, CS, sfc_state, fluxes) enddo ; enddo balancing_area = global_area_integral(bal_frac, G) - if (balancing_area > 0.0) then !### Examine whether the rescaling should be inside the parenthesis. - balancing_flux = US%kg_m2s_to_RZ_T*(global_area_integral(ISS%water_flux, G, scale=US%RZ_T_to_kg_m2s, & - area=ISS%area_shelf_h) + & - delta_mass_shelf ) / balancing_area + if (balancing_area > 0.0) then + balancing_flux = ( US%kg_m2s_to_RZ_T*global_area_integral(ISS%water_flux, G, scale=US%RZ_T_to_kg_m2s, & + area=ISS%area_shelf_h) + & + delta_mass_shelf ) / balancing_area else balancing_flux = 0.0 endif @@ -1166,7 +1166,6 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, forces, fl write(0,*) 'IG: ', G%isd, G%isc, G%iec, G%ied, G%jsd, G%jsc, G%jsd, G%jed endif - CS%Time = Time ! ### This might not be in the right place? CS%diag => diag ! Are we being called from the solo ice-sheet driver? When called by the ocean From b7e37761f588493af21b3f47f76cb0b1f25d624e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 7 May 2020 08:12:01 -0400 Subject: [PATCH 05/40] (*)Fix dimensional inconsistency in calc_CVMix_shear Corrected a dimensionally inconsistent expression in the thickness that determines when a layer has vanished in calculate_CVMix_shear. This would change answers for some ranges of dimensional rescaling or for non-Boussinesq configurations using the CVMix shear code. All answers in the MOM6-examples test cases are bitwise identical. --- src/parameterizations/vertical/MOM_CVMix_shear.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_CVMix_shear.F90 b/src/parameterizations/vertical/MOM_CVMix_shear.F90 index 6f1a629ab4..f099305f0c 100644 --- a/src/parameterizations/vertical/MOM_CVMix_shear.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_shear.F90 @@ -88,10 +88,11 @@ subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS ) real, dimension(G%ke+1) :: Ri_Grad !< Gradient Richardson number [nondim] real, dimension(G%ke+1) :: Kvisc !< Vertical viscosity at interfaces [m2 s-1] real, dimension(G%ke+1) :: Kdiff !< Diapycnal diffusivity at interfaces [m2 s-1] - real, parameter :: epsln = 1.e-10 !< Threshold to identify vanished layers [m] + real :: epsln !< Threshold to identify vanished layers [H ~> m or kg m-2] ! some constants GoRho = US%L_to_Z**2 * GV%g_Earth / GV%Rho0 + epsln = 1.e-10 * GV%m_to_H do j = G%jsc, G%jec do i = G%isc, G%iec @@ -148,7 +149,6 @@ subroutine calculate_CVMix_shear(u_H, v_H, h, tv, kd, kv, G, GV, US, CS ) if (CS%smooth_ri) then ! 1) fill Ri_grad in vanished layers with adjacent value - !### For dimensional consistency, epsln needs to be epsln*GV%m_to_H. do k = 2, G%ke if (h(i,j,k) <= epsln) Ri_grad(k) = Ri_grad(k-1) enddo From 450e954a7dd243500586d8a944fad5971c2fe019 Mon Sep 17 00:00:00 2001 From: Raphael Dussin Date: Thu, 7 May 2020 09:30:12 -0400 Subject: [PATCH 06/40] fix out of bounds when nrows == 1 --- src/diagnostics/MOM_wave_speed.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90 index 65a23e0fa2..9da2963c16 100644 --- a/src/diagnostics/MOM_wave_speed.F90 +++ b/src/diagnostics/MOM_wave_speed.F90 @@ -1246,7 +1246,7 @@ subroutine tridiag_det(a, b, c, nrows, lam, det_out, ddet_out, row_scale) rscl = 1.0 ; if (present(row_scale)) rscl = row_scale det(1) = 1.0 ; ddet(1) = 0.0 - det(2) = b(2)-lam ; ddet(2) = -1.0 + if (nrows > 1) then ; det(2) = b(2)-lam ; ddet(2) = -1.0 ; endif do n=3,nrows det(n) = rscl*(b(n)-lam)*det(n-1) - rscl*(a(n)*c(n-1))*det(n-2) ddet(n) = rscl*(b(n)-lam)*ddet(n-1) - rscl*(a(n)*c(n-1))*ddet(n-2) - det(n-1) From e16d7ece85be1a5a47789f4bfbd1d5a25653ec01 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 7 May 2020 12:54:46 -0400 Subject: [PATCH 07/40] Reduced the halo size in an update for eta Reduced the halo size for an update to eta to 1. There had been a note in the code suggesting that this should be possible for the past 3 years, but that for unknown reasons this changed the answers with the circle_OBCs test case. At some point whatever problems with the OBCs that caused this unexpected dependency on overly large halos has been corrected, and this halo update can now use its expected size of 1. All answers are bitwise identical. --- src/core/MOM_dynamics_split_RK2.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 2517846b4c..5a20e60b04 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -388,10 +388,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s !--- begin set up for group halo pass cont_stencil = continuity_stencil(CS%continuity_CSp) - !### Apart from circle_OBCs halo for eta could be 1, but halo>=3 is required - !### to match circle_OBCs solutions. Why? call cpu_clock_begin(id_clock_pass) - call create_group_pass(CS%pass_eta, eta, G%Domain) !### , halo=1) + call create_group_pass(CS%pass_eta, eta, G%Domain, halo=1) call create_group_pass(CS%pass_visc_rem, CS%visc_rem_u, CS%visc_rem_v, G%Domain, & To_All+SCALAR_PAIR, CGRID_NE, halo=max(1,cont_stencil)) call create_group_pass(CS%pass_uvp, up, vp, G%Domain, halo=max(1,cont_stencil)) From ce3df026e6c50741f5f67aa2401f2ba7535e9a92 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 7 May 2020 15:42:51 -0400 Subject: [PATCH 08/40] +Optionally correct timesteps in unsplit viscosity Added the option to correct the timesteps used for the first predictor step viscosity and the surface boundary layer viscosities with the unsplit time stepping. This includes adding the runtime parameter FIX_UNSPLIT_DT_VISC_BUG, which leads to changes in the MOM_parameter_doc files when SPLIT=False and USE_RK2=False. Several unneeded line continuations were removed. By default, all answers are bitwise identical. --- src/core/MOM_dynamics_unsplit.F90 | 56 +++++++++++++-------------- src/core/MOM_dynamics_unsplit_RK2.F90 | 3 +- 2 files changed, 30 insertions(+), 29 deletions(-) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 8c6e7d4299..25b25d5667 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -120,6 +120,10 @@ module MOM_dynamics_unsplit real, pointer, dimension(:,:) :: tauy_bot => NULL() !< frictional y-bottom stress from the ocean !! to the seafloor [R L Z T-2 ~> Pa] + logical :: use_correct_dt_visc ! If true, use the correct timestep in the viscous terms applied + ! in the first predictor step with the unsplit time stepping scheme, + ! and in the calculation of the turbulent mixed layer properties + ! for viscosity. The default should be true, but it is false. logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: module_is_initialized = .false. !< Record whether this mouled has been initialzed. @@ -228,6 +232,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp, vpp ! Predicted meridional velocities [L T-1 ~> m s-1] real, dimension(:,:), pointer :: p_surf => NULL() real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s]. + real :: dt_visc ! The time step for a part of the update due to viscosity [T ~> s]. logical :: dyn_p_surf integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -255,8 +260,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! diffu = horizontal viscosity terms (u,h) call enable_averages(dt, Time_local, CS%diag) call cpu_clock_begin(id_clock_horvisc) - call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, & - G, GV, US, CS%hor_visc_CSp) + call horizontal_viscosity(u, v, h, CS%diffu, CS%diffv, MEKE, Varmix, G, GV, US, CS%hor_visc_CSp) call cpu_clock_end(id_clock_horvisc) call disable_averaging(CS%diag) @@ -323,31 +327,29 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! up = u + dt_pred * (PFu + CAu) call cpu_clock_begin(id_clock_mom_update) do k=1,nz ; do j=js,je ; do I=Isq,Ieq - up(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt_pred * & - (CS%PFu(I,j,k) + CS%CAu(I,j,k))) + up(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt_pred * (CS%PFu(I,j,k) + CS%CAu(I,j,k))) enddo ; enddo ; enddo do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - vp(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt_pred * & - (CS%PFv(i,J,k) + CS%CAv(i,J,k))) + vp(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt_pred * (CS%PFv(i,J,k) + CS%CAv(i,J,k))) enddo ; enddo ; enddo call cpu_clock_end(id_clock_mom_update) if (CS%debug) then call MOM_state_chksum("Predictor 1", up, vp, h_av, uh, vh, G, GV, US) - call MOM_accel_chksum("Predictor 1 accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv,& + call MOM_accel_chksum("Predictor 1 accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv, & CS%diffu, CS%diffv, G, GV, US) endif ! up <- up + dt/2 d/dz visc d/dz up call cpu_clock_begin(id_clock_vertvisc) call enable_averages(dt, Time_local, CS%diag) - call set_viscous_ML(u, v, h_av, tv, forces, visc, dt*0.5, G, GV, US, & - CS%set_visc_CSp) + dt_visc = 0.5*dt ; if (CS%use_correct_dt_visc) dt_visc = dt + call set_viscous_ML(u, v, h_av, tv, forces, visc, dt_visc, G, GV, US, CS%set_visc_CSp) call disable_averaging(CS%diag) - !### I think that the time steps in the next two calls should be dt_pred. - call vertvisc_coef(up, vp, h_av, forces, visc, dt*0.5, G, GV, US, & - CS%vertvisc_CSp, CS%OBC) - call vertvisc(up, vp, h_av, forces, visc, dt*0.5, CS%OBC, CS%ADp, CS%CDp, & + + dt_visc = 0.5*dt ; if (CS%use_correct_dt_visc) dt_visc = dt_pred + call vertvisc_coef(up, vp, h_av, forces, visc, dt_visc, G, GV, US, CS%vertvisc_CSp, CS%OBC) + call vertvisc(up, vp, h_av, forces, visc, dt_visc, CS%OBC, CS%ADp, CS%CDp, & G, GV, US, CS%vertvisc_CSp, Waves=Waves) call cpu_clock_end(id_clock_vertvisc) call pass_vector(up, vp, G%Domain, clock=id_clock_pass) @@ -355,8 +357,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! uh = up * hp ! h_av = hp + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(up, vp, hp, h_av, uh, vh, (0.5*dt), G, GV, US, & - CS%continuity_CSp, OBC=CS%OBC) + call continuity(up, vp, hp, h_av, uh, vh, (0.5*dt), G, GV, US, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call pass_var(h_av, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -392,25 +393,22 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! upp = u + dt/2 * ( PFu + CAu ) call cpu_clock_begin(id_clock_mom_update) do k=1,nz ; do j=js,je ; do I=Isq,Ieq - upp(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * 0.5 * & - (CS%PFu(I,j,k) + CS%CAu(I,j,k))) + upp(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * 0.5 * (CS%PFu(I,j,k) + CS%CAu(I,j,k))) enddo ; enddo ; enddo do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - vpp(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt * 0.5 * & - (CS%PFv(i,J,k) + CS%CAv(i,J,k))) + vpp(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt * 0.5 * (CS%PFv(i,J,k) + CS%CAv(i,J,k))) enddo ; enddo ; enddo call cpu_clock_end(id_clock_mom_update) if (CS%debug) then call MOM_state_chksum("Predictor 2", upp, vpp, h_av, uh, vh, G, GV, US) - call MOM_accel_chksum("Predictor 2 accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv,& + call MOM_accel_chksum("Predictor 2 accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv, & CS%diffu, CS%diffv, G, GV, US) endif ! upp <- upp + dt/2 d/dz visc d/dz upp call cpu_clock_begin(id_clock_vertvisc) - call vertvisc_coef(upp, vpp, hp, forces, visc, dt*0.5, G, GV, US, & - CS%vertvisc_CSp, CS%OBC) + call vertvisc_coef(upp, vpp, hp, forces, visc, dt*0.5, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc(upp, vpp, hp, forces, visc, dt*0.5, CS%OBC, CS%ADp, CS%CDp, & G, GV, US, CS%vertvisc_CSp, Waves=Waves) call cpu_clock_end(id_clock_vertvisc) @@ -419,8 +417,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! uh = upp * hp ! h = hp + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(upp, vpp, hp, h, uh, vh, (dt*0.5), G, GV, US, & - CS%continuity_CSp, OBC=CS%OBC) + call continuity(upp, vpp, hp, h, uh, vh, (dt*0.5), G, GV, US, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call pass_var(h, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -470,12 +467,10 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & call open_boundary_zero_normal_flow(CS%OBC, G, CS%CAu, CS%CAv) endif do k=1,nz ; do j=js,je ; do I=Isq,Ieq - u(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * & - (CS%PFu(I,j,k) + CS%CAu(I,j,k))) + u(I,j,k) = G%mask2dCu(I,j) * (u(I,j,k) + dt * (CS%PFu(I,j,k) + CS%CAu(I,j,k))) enddo ; enddo ; enddo do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie - v(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt * & - (CS%PFv(i,J,k) + CS%CAv(i,J,k))) + v(i,J,k) = G%mask2dCv(i,J) * (v(i,J,k) + dt * (CS%PFv(i,J,k) + CS%CAv(i,J,k))) enddo ; enddo ; enddo ! u <- u + dt d/dz visc d/dz u @@ -639,6 +634,11 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS default=.false., debuggingParam=.true.) call get_param(param_file, mdl, "TIDES", use_tides, & "If true, apply tidal momentum forcing.", default=.false.) + call get_param(param_file, mdl, "FIX_UNSPLIT_DT_VISC_BUG", CS%use_correct_dt_visc, & + "If true, use the correct timestep in the viscous terms applied in the first "//& + "predictor step with the unsplit time stepping scheme, and in the calculation "//& + "of the turbulent mixed layer properties for viscosity. "//& + "The default should be true.", default=.false.) allocate(CS%taux_bot(IsdB:IedB,jsd:jed)) ; CS%taux_bot(:,:) = 0.0 allocate(CS%tauy_bot(isd:ied,JsdB:JedB)) ; CS%tauy_bot(:,:) = 0.0 diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index d3adfaa194..085611c51b 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -339,9 +339,10 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! up[n-1/2] <- up*[n-1/2] + dt/2 d/dz visc d/dz up[n-1/2] call cpu_clock_begin(id_clock_vertvisc) call enable_averages(dt, Time_local, CS%diag) - call set_viscous_ML(up, vp, h_av, tv, forces, visc, dt_pred, G, GV, US, & + call set_viscous_ML(up, vp, h_av, tv, forces, visc, dt, G, GV, US, & CS%set_visc_CSp) call disable_averaging(CS%diag) + call vertvisc_coef(up, vp, h_av, forces, visc, dt_pred, G, GV, US, & CS%vertvisc_CSp, CS%OBC) call vertvisc(up, vp, h_av, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, & From c39ab9e4c2eed6e1a5e68366e9196598def69b0e Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 7 May 2020 15:53:40 -0400 Subject: [PATCH 09/40] +Use time types to manage offline tracers Revised the offline tracer code to use time types to manage the offline tracer time stepping. There is a new optional argument to step_offline. The offline tracers are not adequately tested, but the code compiles and the logic should be correct. All answers in the MOM6-examples test cases are bitwise identical. --- src/core/MOM.F90 | 27 +++++++++++++++------------ src/tracer/MOM_offline_main.F90 | 18 +++++++++++------- 2 files changed, 26 insertions(+), 19 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index aff4860a21..1553ff7e1a 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -45,7 +45,7 @@ module MOM use MOM_spatial_means, only : global_mass_integral use MOM_time_manager, only : time_type, real_to_time, time_type_to_real, operator(+) use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/) -use MOM_time_manager, only : operator(>=), increment_date +use MOM_time_manager, only : operator(>=), operator(==), increment_date use MOM_unit_tests, only : unit_tests use coupler_types_mod, only : coupler_type_send_data, coupler_1d_bc_type, coupler_type_spawn @@ -1404,7 +1404,8 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS logical :: skip_diffusion integer :: id_eta_diff_end - integer, pointer :: accumulated_time => NULL() + type(time_type), pointer :: accumulated_time => NULL() + type(time_type), pointer :: vertical_time => NULL() integer :: i,j,k integer :: is, ie, js, je, isd, ied, jsd, jed @@ -1426,32 +1427,30 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS call cpu_clock_begin(id_clock_offline_tracer) call extract_offline_main(CS%offline_CSp, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, & - dt_offline, dt_offline_vertical, skip_diffusion) + vertical_time, dt_offline, dt_offline_vertical, skip_diffusion) Time_end = increment_date(Time_start, seconds=floor(time_interval+0.001)) call enable_averaging(time_interval, Time_end, CS%diag) ! Check to see if this is the first iteration of the offline interval - if (accumulated_time==0) then + if (accumulated_time == real_to_time(0.0)) then first_iter = .true. else ! This is probably unnecessary but is used to guard against unwanted behavior first_iter = .false. endif - ! Check to see if vertical tracer functions should be done - if ( mod(accumulated_time, floor(US%T_to_s*dt_offline_vertical + 1e-6)) == 0 ) then + ! Check to see if vertical tracer functions should be done + if (first_iter .or. (accumulated_time >= vertical_time)) then do_vertical = .true. + vertical_time = accumulated_time + real_to_time(US%T_to_s*dt_offline_vertical) else do_vertical = .false. endif ! Increment the amount of time elapsed since last read and check if it's time to roll around - accumulated_time = mod(accumulated_time + int(time_interval), floor(US%T_to_s*dt_offline+1e-6)) - if (accumulated_time==0) then - last_iter = .true. - else - last_iter = .false. - endif + accumulated_time = accumulated_time + real_to_time(time_interval) + + last_iter = (accumulated_time >= real_to_time(US%T_to_s*dt_offline)) if (CS%use_ALE_algorithm) then ! If this is the first iteration in the offline timestep, then we need to read in fields and @@ -1566,6 +1565,10 @@ subroutine step_offline(forces, fluxes, sfc_state, Time_start, time_interval, CS fluxes%fluxes_used = .true. + if (last_iter) then + accumulated_time = real_to_time(0.0) + endif + call cpu_clock_end(id_clock_offline_tracer) end subroutine step_offline diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index 65f83ecfea..c0ee07d434 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -28,7 +28,7 @@ module MOM_offline_main use MOM_offline_aux, only : distribute_residual_uh_upwards, distribute_residual_vh_upwards use MOM_opacity, only : opacity_CS, optics_type use MOM_open_boundary, only : ocean_OBC_type -use MOM_time_manager, only : time_type +use MOM_time_manager, only : time_type, real_to_time use MOM_tracer_advect, only : tracer_advect_CS, advect_tracer use MOM_tracer_diabatic, only : applyTracerBoundaryFluxesInOut use MOM_tracer_flow_control, only : tracer_flow_control_CS, call_tracer_column_fns, call_tracer_stocks @@ -79,7 +79,8 @@ module MOM_offline_main integer :: start_index !< Timelevel to start integer :: iter_no !< Timelevel to start integer :: numtime !< How many timelevels in the input fields - integer :: accumulated_time !< Length of time accumulated in the current offline interval + type(time_type) :: accumulated_time !< Length of time accumulated in the current offline interval + type(time_type) :: vertical_time !< The next value of accumulate_time at which to apply vertical processes ! Index of each of the variables to be read in with separate indices for each variable if they ! are set off from each other in time integer :: ridx_sum = -1 !< Read index offset of the summed variables @@ -1200,7 +1201,7 @@ end subroutine post_offline_convergence_diags !> Extracts members of the offline main control structure. All arguments are optional except !! the control structure itself -subroutine extract_offline_main(CS, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, & +subroutine extract_offline_main(CS, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, vertical_time, & dt_offline, dt_offline_vertical, skip_diffusion) type(offline_transport_CS), target, intent(in ) :: CS !< Offline control structure ! Returned optional arguments @@ -1212,9 +1213,10 @@ subroutine extract_offline_main(CS, uhtr, vhtr, eatr, ebtr, h_end, accumulated_t !! one time step [H ~> m or kg m-2] real, dimension(:,:,:), optional, pointer :: h_end !< Thicknesses at the end of offline timestep !! [H ~> m or kg m-2] - !### Why are the following variables integers? - integer, optional, pointer :: accumulated_time !< Length of time accumulated in the - !! current offline interval [s] + type(time_type), optional, pointer :: accumulated_time !< Length of time accumulated in the + !! current offline interval + type(time_type), optional, pointer :: vertical_time !< The next value of accumulate_time at which to + !! vertical processes real, optional, intent( out) :: dt_offline !< Timestep used for offline tracers [T ~> s] real, optional, intent( out) :: dt_offline_vertical !< Timestep used for calls to tracer !! vertical physics [T ~> s] @@ -1229,6 +1231,7 @@ subroutine extract_offline_main(CS, uhtr, vhtr, eatr, ebtr, h_end, accumulated_t ! Pointers to integer members which need to be modified if (present(accumulated_time)) accumulated_time => CS%accumulated_time + if (present(vertical_time)) vertical_time => CS%vertical_time ! Return value of non-modified integers if (present(dt_offline)) dt_offline = CS%dt_offline @@ -1414,7 +1417,8 @@ subroutine offline_transport_init(param_file, CS, diabatic_CSp, G, GV, US) end select ! Set the accumulated time to zero - CS%accumulated_time = 0 + CS%accumulated_time = real_to_time(0.0) + CS%vertical_time = CS%accumulated_time ! Set the starting read index for time-averaged and snapshotted fields CS%ridx_sum = CS%start_index if (CS%fields_are_offset) CS%ridx_snap = next_modulo_time(CS%start_index,CS%numtime) From ea17a977f7c7935c753a483b6f24c7bbdb6fcdd4 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 7 May 2020 19:08:30 -0400 Subject: [PATCH 10/40] Corrected units in comments in coord_rho Corrected or added the units of variables in the coord_rho module. All answers are bitwise identical. --- src/ALE/MOM_regridding.F90 | 2 +- src/ALE/coord_rho.F90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index 9bc71dd15f..ed6e66e0ae 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -1353,7 +1353,7 @@ subroutine build_rho_grid( G, GV, US, h, tv, dzInterface, remapCS, CS ) integer :: nz integer :: i, j, k real :: nominalDepth ! Depth of the bottom of the ocean, positive downward [H ~> m or kg m-2] - real, dimension(SZK_(GV)+1) :: zOld, zNew + real, dimension(SZK_(GV)+1) :: zOld, zNew ! Old and new interface heights [H ~> m or kg m-2] real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] #ifdef __DO_SAFETY_CHECKS__ real :: totalThickness diff --git a/src/ALE/coord_rho.F90 b/src/ALE/coord_rho.F90 index 5299c74b1b..dce802ff3c 100644 --- a/src/ALE/coord_rho.F90 +++ b/src/ALE/coord_rho.F90 @@ -91,7 +91,7 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & h_neglect, h_neglect_edge) type(rho_CS), intent(in) :: CS !< coord_rho control structure integer, intent(in) :: nz !< Number of levels on source grid (i.e. length of h, T, S) - real, intent(in) :: depth !< Depth of ocean bottom (positive in m) + real, intent(in) :: depth !< Depth of ocean bottom (positive downward) [H ~> m or kg m-2] real, dimension(nz), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] real, dimension(nz), intent(in) :: T !< Temperature for source column [degC] real, dimension(nz), intent(in) :: S !< Salinity for source column [ppt] @@ -111,7 +111,7 @@ subroutine build_rho_column(CS, nz, depth, h, T, S, eqn_of_state, z_interface, & real, dimension(nz) :: densities ! Layer density [R ~> kg m-3] real, dimension(nz+1) :: xTmp ! Temporary positions [H ~> m or kg m-2] real, dimension(CS%nk) :: h_new ! New thicknesses [H ~> m or kg m-2] - real, dimension(CS%nk+1) :: x1 + real, dimension(CS%nk+1) :: x1 ! Interface heights [H ~> m or kg m-2] ! Construct source column with vanished layers removed (stored in h_nv) call copy_finite_thicknesses(nz, h, CS%min_thickness, count_nonzero_layers, h_nv, mapping) From 0d983b57934faaaa1dc71f67cb5ac0184e203089 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 7 May 2020 19:11:32 -0400 Subject: [PATCH 11/40] (*)Fixed dimensionally inconsistent diagnostic remap Corrected a dimensionally inconsistent statement when remapping diagnostics to rho-space. The ocean bottom depth had been passed in [m] but should have been passed in [H], which could change diagnostics for non-Boussinesq models. Also added multiple comments documenting the units of variables in MOM_regridding. All solutions are bitwise identical, but some rho-space diagnostics could change, especially for non-Boussinesq models. --- src/framework/MOM_diag_remap.F90 | 101 +++++++++++++++---------------- 1 file changed, 50 insertions(+), 51 deletions(-) diff --git a/src/framework/MOM_diag_remap.F90 b/src/framework/MOM_diag_remap.F90 index 0fe937a173..4e12abaa5b 100644 --- a/src/framework/MOM_diag_remap.F90 +++ b/src/framework/MOM_diag_remap.F90 @@ -31,11 +31,11 @@ ! ! For interpolation between h and u grids, we use the following relations: ! -! h->u: f_u[ig] = 0.5 * (f_h[ ig ] + f_h[ig+1]) -! f_u[i1] = 0.5 * (f_h[i1-1] + f_h[ i1 ]) +! h->u: f_u(ig) = 0.5 * (f_h( ig ) + f_h(ig+1)) +! f_u(i1) = 0.5 * (f_h(i1-1) + f_h( i1 )) ! -! u->h: f_h[ig] = 0.5 * (f_u[ig-1] + f_u[ ig ]) -! f_h[i1] = 0.5 * (f_u[ i1 ] + f_u[i1+1]) +! u->h: f_h(ig) = 0.5 * (f_u(ig-1) + f_u( ig )) +! f_h(i1) = 0.5 * (f_u( i1 ) + f_u(i1+1)) ! ! where ig is the grid index and i1 is the 1-based index. That is, a 1-based ! u-point is ahead of its matching h-point in non-symmetric mode, but behind @@ -110,11 +110,11 @@ module MOM_diag_remap character(len=16) :: diag_coord_name = '' !< A name for the purpose of run-time parameters character(len=8) :: diag_module_suffix = '' !< The suffix for the module to appear in diag_table type(remapping_CS) :: remap_cs !< Remapping control structure use for this axes - type(regridding_CS) :: regrid_cs !< Regridding control structure that defines the coordiantes for this axes + type(regridding_CS) :: regrid_cs !< Regridding control structure that defines the coordinates for this axes integer :: nz = 0 !< Number of vertical levels used for remapping - real, dimension(:,:,:), allocatable :: h !< Remap grid thicknesses - real, dimension(:,:,:), allocatable :: h_extensive !< Remap grid thicknesses for extensive variables - real, dimension(:), allocatable :: dz !< Nominal layer thicknesses + real, dimension(:,:,:), allocatable :: h !< Remap grid thicknesses [H ~> m or kg m-2] + real, dimension(:,:,:), allocatable :: h_extensive !< Remap grid thicknesses for extensive + !! variables [H ~> m or kg m-2] integer :: interface_axes_id = 0 !< Vertical axes id for remapping at interfaces integer :: layer_axes_id = 0 !< Vertical axes id for remapping on layers logical :: answers_2018 !< If true, use the order of arithmetic and expressions for remapping @@ -151,7 +151,6 @@ subroutine diag_remap_end(remap_cs) type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diag remapping control structure if (allocated(remap_cs%h)) deallocate(remap_cs%h) - if (allocated(remap_cs%dz)) deallocate(remap_cs%dz) remap_cs%configured = .false. remap_cs%initialized = .false. remap_cs%used = .false. @@ -274,15 +273,15 @@ function diag_remap_axes_configured(remap_cs) !! coordinates then technically we should also regenerate the !! target grid whenever T/S change. subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_target) - type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diagnostic coordinate control structure - type(ocean_grid_type), pointer :: G !< The ocean's grid type - type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(:, :, :), intent(in) :: h !< New thickness [H ~> m or kg m-2] - real, dimension(:, :, :), intent(in) :: T !< New temperatures [degC] - real, dimension(:, :, :), intent(in) :: S !< New salinities [ppt] - type(EOS_type), pointer :: eqn_of_state !< A pointer to the equation of state - real, dimension(:, :, :), intent(inout) :: h_target !< Where to store the new diagnostic array + type(diag_remap_ctrl), intent(inout) :: remap_cs !< Diagnostic coordinate control structure + type(ocean_grid_type), pointer :: G !< The ocean's grid type + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(:,:,:), intent(in) :: h !< New thickness [H ~> m or kg m-2] + real, dimension(:,:,:), intent(in) :: T !< New temperatures [degC] + real, dimension(:,:,:), intent(in) :: S !< New salinities [ppt] + type(EOS_type), pointer :: eqn_of_state !< A pointer to the equation of state + real, dimension(:,:,:), intent(inout) :: h_target !< The new diagnostic thicknesses [H ~> m or kg m-2] ! Local variables real, dimension(remap_cs%nz + 1) :: zInterfaces ! Interface positions [H ~> m or kg m-2] @@ -328,9 +327,8 @@ subroutine diag_remap_update(remap_cs, G, GV, US, h, T, S, eqn_of_state, h_targe call build_sigma_column(get_sigma_CS(remap_cs%regrid_cs), & GV%Z_to_H*G%bathyT(i,j), sum(h(i,j,:)), zInterfaces) elseif (remap_cs%vertical_coord == coordinateMode('RHO')) then -!### I think that the conversion factor in the 2nd line should be GV%Z_to_H call build_rho_column(get_rho_CS(remap_cs%regrid_cs), G%ke, & - US%Z_to_m*G%bathyT(i,j), h(i,j,:), T(i,j,:), S(i,j,:), & + GV%Z_to_H*G%bathyT(i,j), h(i,j,:), T(i,j,:), S(i,j,:), & eqn_of_state, zInterfaces, h_neglect, h_neglect_edge) elseif (remap_cs%vertical_coord == coordinateMode('SLIGHT')) then ! call build_slight_column(remap_cs%regrid_cs,remap_cs%remap_cs, nz, & @@ -354,16 +352,16 @@ subroutine diag_remap_do_remap(remap_cs, G, GV, h, staggered_in_x, staggered_in_ type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - real, dimension(:,:,:), intent(in) :: h !< The current thicknesses + real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points - real, dimension(:,:,:), pointer :: mask !< A mask for the field - real, dimension(:,:,:), intent(in) :: field(:,:,:) !< The diagnostic field to be remapped - real, dimension(:,:,:), intent(inout) :: remapped_field !< Field remapped to new coordinate + real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim] + real, dimension(:,:,:), intent(in) :: field(:,:,:) !< The diagnostic field to be remapped [A] + real, dimension(:,:,:), intent(inout) :: remapped_field !< Field remapped to new coordinate [A] ! Local variables - real, dimension(remap_cs%nz) :: h_dest - real, dimension(size(h,3)) :: h_src - real :: h_neglect, h_neglect_edge + real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] + real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] + real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] integer :: nz_src, nz_dest integer :: i, j, k !< Grid index integer :: i1, j1 !< 1-based index @@ -446,14 +444,15 @@ end subroutine diag_remap_do_remap !> Calculate masks for target grid subroutine diag_remap_calc_hmask(remap_cs, G, mask) - type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(:,:,:), intent(out) :: mask !< h-point mask for target grid + type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + real, dimension(:,:,:), intent(out) :: mask !< h-point mask for target grid [nondim] ! Local variables - real, dimension(remap_cs%nz) :: h_dest + real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] integer :: i, j, k logical :: mask_vanished_layers - real :: h_tot, h_err + real :: h_tot ! Sum of all thicknesses [H ~> m or kg m-2] + real :: h_err ! An estimate of a negligible thickness [H ~> m or kg m-2] call assert(remap_cs%initialized, 'diag_remap_calc_hmask: remap_cs not initialized.') @@ -492,16 +491,16 @@ subroutine vertically_reintegrate_diag_field(remap_cs, G, h, h_target, staggered mask, field, reintegrated_field) type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(:,:,:), intent(in) :: h !< The thicknesses of the source grid - real, dimension(:,:,:), intent(in) :: h_target !< The thicknesses of the target grid + real, dimension(:,:,:), intent(in) :: h !< The thicknesses of the source grid [H ~> m or kg m-2] + real, dimension(:,:,:), intent(in) :: h_target !< The thicknesses of the target grid [H ~> m or kg m-2] logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points - real, dimension(:,:,:), pointer :: mask !< A mask for the field - real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped - real, dimension(:,:,:), intent(inout) :: reintegrated_field !< Field argument remapped to alternative coordinate + real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim] + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(:,:,:), intent(inout) :: reintegrated_field !< Field argument remapped to alternative coordinate [A] ! Local variables - real, dimension(remap_cs%nz) :: h_dest - real, dimension(size(h,3)) :: h_src + real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] + real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] integer :: nz_src, nz_dest integer :: i, j, k !< Grid index integer :: i1, j1 !< 1-based index @@ -572,16 +571,16 @@ end subroutine vertically_reintegrate_diag_field subroutine vertically_interpolate_diag_field(remap_cs, G, h, staggered_in_x, staggered_in_y, & mask, field, interpolated_field) type(diag_remap_ctrl), intent(in) :: remap_cs !< Diagnostic coodinate control structure - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure - real, dimension(:,:,:), intent(in) :: h !< The current thicknesses + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] logical, intent(in) :: staggered_in_x !< True is the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True is the y-axis location is at v or q points - real, dimension(:,:,:), pointer :: mask !< A mask for the field - real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped - real, dimension(:,:,:), intent(inout) :: interpolated_field !< Field argument remapped to alternative coordinate + real, dimension(:,:,:), pointer :: mask !< A mask for the field [nondim] + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(:,:,:), intent(inout) :: interpolated_field !< Field argument remapped to alternative coordinate [A] ! Local variables - real, dimension(remap_cs%nz) :: h_dest - real, dimension(size(h,3)) :: h_src + real, dimension(remap_cs%nz) :: h_dest ! Destination thicknesses [H ~> m or kg m-2] + real, dimension(size(h,3)) :: h_src ! A column of source thicknesses [H ~> m or kg m-2] integer :: nz_src, nz_dest integer :: i, j, k !< Grid index integer :: i1, j1 !< 1-based index @@ -656,20 +655,20 @@ subroutine horizontally_average_diag_field(G, GV, h, staggered_in_x, staggered_i averaged_mask) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean vertical grid structure - real, dimension(:,:,:), intent(in) :: h !< The current thicknesses + real, dimension(:,:,:), intent(in) :: h !< The current thicknesses [H ~> m or kg m-2] logical, intent(in) :: staggered_in_x !< True if the x-axis location is at u or q points logical, intent(in) :: staggered_in_y !< True if the y-axis location is at v or q points logical, intent(in) :: is_layer !< True if the z-axis location is at h points logical, intent(in) :: is_extensive !< True if the z-direction is spatially integrated (over layers) - real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped - real, dimension(:), intent(inout) :: averaged_field !< Field argument horizontally averaged - logical, dimension(:), intent(inout) :: averaged_mask !< Mask for horizontally averaged field + real, dimension(:,:,:), intent(in) :: field !< The diagnostic field to be remapped [A] + real, dimension(:), intent(inout) :: averaged_field !< Field argument horizontally averaged [A] + logical, dimension(:), intent(inout) :: averaged_mask !< Mask for horizontally averaged field [nondim] ! Local variables real, dimension(G%isc:G%iec, G%jsc:G%jec, size(field,3)) :: volume, stuff real, dimension(size(field, 3)) :: vol_sum, stuff_sum ! nz+1 is needed for interface averages type(EFP_type), dimension(2*size(field,3)) :: sums_EFP ! Sums of volume or stuff by layer - real :: height + real :: height ! An average thickness attributed to an velocity point [H ~> m or kg m-2] integer :: i, j, k, nz integer :: i1, j1 !< 1-based index From 949b04ff88d07ffb11377d5265a72aa2acaaa9d9 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 7 May 2020 19:51:33 -0400 Subject: [PATCH 12/40] dOxygenized a comment in MOM_dyn_unsplit_CS --- src/core/MOM_dynamics_unsplit.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 25b25d5667..99cdc932e2 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -120,10 +120,10 @@ module MOM_dynamics_unsplit real, pointer, dimension(:,:) :: tauy_bot => NULL() !< frictional y-bottom stress from the ocean !! to the seafloor [R L Z T-2 ~> Pa] - logical :: use_correct_dt_visc ! If true, use the correct timestep in the viscous terms applied - ! in the first predictor step with the unsplit time stepping scheme, - ! and in the calculation of the turbulent mixed layer properties - ! for viscosity. The default should be true, but it is false. + logical :: use_correct_dt_visc !< If true, use the correct timestep in the viscous terms applied + !! in the first predictor step with the unsplit time stepping scheme, + !! and in the calculation of the turbulent mixed layer properties + !! for viscosity. The default should be true, but it is false. logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: module_is_initialized = .false. !< Record whether this mouled has been initialzed. From 6ed79a43b2aef5af8138f336cc3ff72f1051a450 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 8 May 2020 07:15:47 -0400 Subject: [PATCH 13/40] +Added FIX_UNSPLIT_DT_VISC_BUG to unsplit_RK2 Add the runtime parameter FIX_UNSPLIT_DT_VISC_BUG to MOM_dynamics_unsplit_RK2 to preserve previous answers by default, analogously to what was done in a recent commit for MOM_dynamics_unsplit. Some unneeded line continuations were also removed for brevity. By default all answers are identical to those on the dev/gfdl branch. --- src/core/MOM_dynamics_unsplit.F90 | 10 ++++---- src/core/MOM_dynamics_unsplit_RK2.F90 | 36 ++++++++++++++------------- 2 files changed, 24 insertions(+), 22 deletions(-) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 99cdc932e2..d6a5186be3 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -629,16 +629,16 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS CS%diag => diag + call get_param(param_file, mdl, "FIX_UNSPLIT_DT_VISC_BUG", CS%use_correct_dt_visc, & + "If true, use the correct timestep in the viscous terms applied in the first "//& + "predictor step with the unsplit time stepping scheme, and in the calculation "//& + "of the turbulent mixed layer properties for viscosity with unsplit or "//& + "unsplit_RK2. The default should be true.", default=.false.) call get_param(param_file, mdl, "DEBUG", CS%debug, & "If true, write out verbose debugging data.", & default=.false., debuggingParam=.true.) call get_param(param_file, mdl, "TIDES", use_tides, & "If true, apply tidal momentum forcing.", default=.false.) - call get_param(param_file, mdl, "FIX_UNSPLIT_DT_VISC_BUG", CS%use_correct_dt_visc, & - "If true, use the correct timestep in the viscous terms applied in the first "//& - "predictor step with the unsplit time stepping scheme, and in the calculation "//& - "of the turbulent mixed layer properties for viscosity. "//& - "The default should be true.", default=.false.) allocate(CS%taux_bot(IsdB:IedB,jsd:jed)) ; CS%taux_bot(:,:) = 0.0 allocate(CS%tauy_bot(isd:ied,JsdB:JedB)) ; CS%tauy_bot(:,:) = 0.0 diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 085611c51b..e3ec48ff58 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -123,6 +123,9 @@ module MOM_dynamics_unsplit_RK2 !! the extent to which the treatment of gravity waves !! is forward-backward (0) or simulated backward !! Euler (1). 0 is almost always used. + logical :: use_correct_dt_visc !< If true, use the correct timestep in the calculation of the + !! turbulent mixed layer properties for viscosity. + !! The default should be true, but it is false. logical :: debug !< If true, write verbose checksums for debugging purposes. logical :: module_is_initialized = .false. !< Record whether this mouled has been initialzed. @@ -238,8 +241,8 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up ! Predicted zonal velocities [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp ! Predicted meridional velocities [L T-1 ~> m s-1] real, dimension(:,:), pointer :: p_surf => NULL() - real :: dt_pred ! The time step for the predictor part of the baroclinic - ! time stepping [T ~> s]. + real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s] + real :: dt_visc ! The time step for a part of the update due to viscosity [T ~> s] logical :: dyn_p_surf integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -280,17 +283,15 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, call cpu_clock_begin(id_clock_continuity) ! This is a duplicate calculation of the last continuity from the previous step ! and could/should be optimized out. -AJA - call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, GV, US, & - CS%continuity_CSp, OBC=CS%OBC) + call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, GV, US, CS%continuity_CSp, OBC=CS%OBC) call cpu_clock_end(id_clock_continuity) call pass_var(hp, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) ! h_av = (h + hp)/2 (used in PV denominator) call cpu_clock_begin(id_clock_mom_update) - do k=1,nz - do j=js-2,je+2 ; do i=is-2,ie+2 - h_av(i,j,k) = (h_in(i,j,k) + hp(i,j,k)) * 0.5 + do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2 + h_av(i,j,k) = (h_in(i,j,k) + hp(i,j,k)) * 0.5 enddo ; enddo ; enddo call cpu_clock_end(id_clock_mom_update) @@ -305,8 +306,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, if (dyn_p_surf) then ; do j=js-2,je+2 ; do i=is-2,ie+2 p_surf(i,j) = 0.5*p_surf_begin(i,j) + 0.5*p_surf_end(i,j) enddo ; enddo ; endif - call PressureForce(h_in, tv, CS%PFu, CS%PFv, G, GV, US, & - CS%PressureForce_CSp, CS%ALE_CSp, p_surf) + call PressureForce(h_in, tv, CS%PFu, CS%PFv, G, GV, US, CS%PressureForce_CSp, CS%ALE_CSp, p_surf) call cpu_clock_end(id_clock_pres) call pass_vector(CS%PFu, CS%PFv, G%Domain, clock=id_clock_pass) call pass_vector(CS%CAu, CS%CAv, G%Domain, clock=id_clock_pass) @@ -339,12 +339,11 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! up[n-1/2] <- up*[n-1/2] + dt/2 d/dz visc d/dz up[n-1/2] call cpu_clock_begin(id_clock_vertvisc) call enable_averages(dt, Time_local, CS%diag) - call set_viscous_ML(up, vp, h_av, tv, forces, visc, dt, G, GV, US, & - CS%set_visc_CSp) + dt_visc = dt_pred ; if (CS%use_correct_dt_visc) dt_visc = dt + call set_viscous_ML(up, vp, h_av, tv, forces, visc, dt_visc, G, GV, US, CS%set_visc_CSp) call disable_averaging(CS%diag) - call vertvisc_coef(up, vp, h_av, forces, visc, dt_pred, G, GV, US, & - CS%vertvisc_CSp, CS%OBC) + call vertvisc_coef(up, vp, h_av, forces, visc, dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc(up, vp, h_av, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, & G, GV, US, CS%vertvisc_CSp) call cpu_clock_end(id_clock_vertvisc) @@ -395,12 +394,10 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! up[n] <- up* + dt d/dz visc d/dz up ! u[n] <- u*[n] + dt d/dz visc d/dz u[n] call cpu_clock_begin(id_clock_vertvisc) - call vertvisc_coef(up, vp, h_av, forces, visc, dt, G, GV, US, & - CS%vertvisc_CSp, CS%OBC) + call vertvisc_coef(up, vp, h_av, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc(up, vp, h_av, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, & G, GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot) - call vertvisc_coef(u_in, v_in, h_av, forces, visc, dt, G, GV, US, & - CS%vertvisc_CSp, CS%OBC) + call vertvisc_coef(u_in, v_in, h_av, forces, visc, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC) call vertvisc(u_in, v_in, h_av, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp,& G, GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot) call cpu_clock_end(id_clock_vertvisc) @@ -594,6 +591,11 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag "If SPLIT is false and USE_RK2 is true, BEGW can be "//& "between 0 and 0.5 to damp gravity waves.", & units="nondim", default=0.0) + call get_param(param_file, mdl, "FIX_UNSPLIT_DT_VISC_BUG", CS%use_correct_dt_visc, & + "If true, use the correct timestep in the viscous terms applied in the first "//& + "predictor step with the unsplit time stepping scheme, and in the calculation "//& + "of the turbulent mixed layer properties for viscosity with unsplit or "//& + "unsplit_RK2. The default should be true.", default=.false.) call get_param(param_file, mdl, "DEBUG", CS%debug, & "If true, write out verbose debugging data.", & default=.false., debuggingParam=.true.) From b8c3570260e74940bf0af06e6502e61d5293ac0f Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sat, 9 May 2020 11:27:23 -0400 Subject: [PATCH 14/40] +Add false position iteration for MLD to ePBL Added the option to replace bisection with false position iteration when determining the mixed layer depth. The new option is selected with the new runtime parameter EPBL_MLD_BISECTION = False. When EPBL_MLD_TOLERANCE = 1 m, the new option reduces the average number of iterations from 9.53 to 4.23 in the OM4_05 test case. When EPBL_MLD_TOLERANCE = 0.001 m, the average is reduced from 19.4 iterations to 7.65 with the new option. By default all answers are bitwise identical, but there are changes to the entries in the MOM_parameter_doc files, both from the new option and from not reporting other ePBL options when they are not valid options. --- .../vertical/MOM_energetic_PBL.F90 | 99 ++++++++++++++----- 1 file changed, 77 insertions(+), 22 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index a9e68736e7..2a66c60c88 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -4,6 +4,7 @@ module MOM_energetic_PBL ! This file is part of MOM6. See LICENSE.md for the license. use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE +use MOM_coms, only : EFP_type, real_to_EFP, EFP_to_real, operator(+), assignment(=), EFP_sum_across_PEs use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc use MOM_diag_mediator, only : time_type, diag_ctrl use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type @@ -33,12 +34,12 @@ module MOM_energetic_PBL type, public :: energetic_PBL_CS ; private !/ Constants - real :: VonKar = 0.41 !< The von Karman coefficient. This should be runtime, but because - !! it is runtime in KPP and set to 0.4 it might change answers. - real :: omega !< The Earth's rotation rate [T-1 ~> s-1]. - real :: omega_frac !< When setting the decay scale for turbulence, use this fraction of - !! the absolute rotation rate blended with the local value of f, as - !! sqrt((1-of)*f^2 + of*4*omega^2) [nondim]. + real :: VonKar = 0.41 !< The von Karman coefficient. This should be a runtime parameter, + !! but because it is set to 0.4 at runtime in KPP it might change answers. + real :: omega !< The Earth's rotation rate [T-1 ~> s-1]. + real :: omega_frac !< When setting the decay scale for turbulence, use this fraction of + !! the absolute rotation rate blended with the local value of f, as + !! sqrt((1-omega_frac)*f^2 + omega_frac*4*omega^2) [nondim]. !/ Convection related terms real :: nstar !< The fraction of the TKE input to the mixed layer available to drive @@ -47,9 +48,14 @@ module MOM_energetic_PBL !! TKE produced by buoyancy. !/ Mixing Length terms - logical :: Use_MLD_iteration=.false. !< False to use old ePBL method. - logical :: MLD_iteration_guess=.false. !< False to default to guessing half the - !! ocean depth for the iteration. + logical :: Use_MLD_iteration !< If true, use the proximity to the bottom of the actively turbulent + !! surface boundary layer to constrain the mixing lengths. + logical :: MLD_iteration_guess !< False to default to guessing half the + !! ocean depth for the first iteration. + logical :: MLD_bisection !! If true, use bisection with the iterative determination of the + !! self-consistent mixed layer depth. Otherwise use the false position + !! after a maximum and minimum bound have been evaluated and the + !! returned value from the previous guess or bisection before this. integer :: max_MLD_its !< The maximum number of iterations that can be used to find a !! self-consistent mixed layer depth with Use_MLD_iteration. real :: MixLenExponent !< Exponent in the mixing length shape-function. @@ -179,6 +185,8 @@ module MOM_energetic_PBL LA, & !< Langmuir number [nondim] LA_MOD !< Modified Langmuir number [nondim] + type(EFP_type), dimension(2) :: sum_its !< The total number of iterations and columns worked on + real, allocatable, dimension(:,:,:) :: & Velocity_Scale, & !< The velocity scale used in getting Kd [Z T-1 ~> m s-1] Mixing_Length !< The length scale used in getting Kd [Z ~> m] @@ -215,6 +223,8 @@ module MOM_energetic_PBL character*(20), parameter :: ADDITIVE_STRING = "ADDITIVE" !>@} +logical :: report_avg_its = .false. ! Report the average number of ePBL iterations for debugging. + !> A type for conveniently passing around ePBL diagnostics for a column. type, public :: ePBL_column_diags ; private !>@{ Local column copies of energy change diagnostics, all in [R Z3 T-3 ~> W m-2]. @@ -755,6 +765,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs ! manner giving a usable guess. When it does fail, it is due to convection ! within the boundary layer. Likely, a new method e.g. surface_disconnect, ! can improve this. + real :: dMLD_min ! The change in diagnosed mixed layer depth when the guess is min_MLD [Z ~> m] + real :: dMLD_max ! The change in diagnosed mixed layer depth when the guess is max_MLD [Z ~> m] logical :: FIRST_OBL ! Flag for computing "found" Mixing layer depth logical :: OBL_converged ! Flag for convergence of MLD integer :: OBL_it ! Iteration counter @@ -829,6 +841,8 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs max_MLD = 0.0 ; do k=1,nz ; max_MLD = max_MLD + h(k)*GV%H_to_Z ; enddo !min_MLD will initialize as 0. min_MLD = 0.0 + ! Set values of the wrong signs to indicate that these changes are not based on valid estimates + ! dMLD_min = -1.0*US%m_to_Z ; dMLD_max = 1.0*US%m_to_Z ! If no first guess is provided for MLD, try the middle of the water column if (MLD_guess <= min_MLD) MLD_guess = 0.5 * (min_MLD + max_MLD) @@ -1408,17 +1422,37 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs !New method uses ML_DEPTH as computed in ePBL routine MLD_found = MLD_output - if (MLD_found - CS%MLD_tol > MLD_guess) then - min_MLD = MLD_guess - elseif (abs(MLD_guess - MLD_found) < CS%MLD_tol) then + if (MLD_found - MLD_guess > CS%MLD_tol) then + min_MLD = MLD_guess ; dMLD_min = MLD_found - MLD_guess + elseif (abs(MLD_found - MLD_guess) < CS%MLD_tol) then OBL_converged = .true. ! Break convergence loop - else - max_MLD = MLD_guess ! We know this guess was too deep + else ! We know this guess was too deep + max_MLD = MLD_guess ; dMLD_max = MLD_found - MLD_guess ! < -CS%MLD_tol endif - ! For next pass, guess average of minimum and maximum values. - !### We should try using the false position method instead of simple bisection. - MLD_guess = 0.5*(min_MLD + max_MLD) + if (.not.OBL_converged) then ; if (CS%MLD_bisection) then + ! For the next pass, guess the average of the minimum and maximum values. + MLD_guess = 0.5*(min_MLD + max_MLD) + else ! Try using the false position method or the returned value instead of simple bisection. + ! Taking the occasional step with MLD_output empirically step helps to converge faster. + if ((dMLD_min > 0.0) .and. (dMLD_max < 0.0) .and. (OBL_it > 2) .and. (mod(OBL_it-1,4)>0)) then + ! Both bounds have valid change estimates and are probably in the range of possible outputs. + MLD_Guess = (dMLD_min*max_MLD - dMLD_max*min_MLD) / (dMLD_min - dMLD_max) + elseif ((MLD_found > min_MLD) .and. (MLD_found < max_MLD)) then + ! The output MLD_found is an interesting guess, as it likely to bracket the true solution + ! along with the previous value of MLD_guess and to be close to the solution. + MLD_guess = MLD_found + else ! Bisect if the other guesses would be out-of-bounds. This does not happen much. + MLD_guess = 0.5*(min_MLD + max_MLD) + endif + endif ; endif + endif + if ((OBL_converged) .or. (OBL_it==CS%Max_MLD_Its)) then + if (report_avg_its) then + CS%sum_its(1) = CS%sum_its(1) + real_to_EFP(real(OBL_it)) + CS%sum_its(2) = CS%sum_its(2) + real_to_EFP(1.0) + endif + exit endif enddo ! Iteration loop for converged boundary layer thickness. if (CS%Use_LT) then @@ -2131,16 +2165,22 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) endif call get_param(param_file, mdl, "MLD_ITERATION_GUESS", CS%MLD_ITERATION_GUESS, & - "A logical that specifies whether or not to use the "//& - "previous timestep MLD as a first guess in the MLD iteration. "//& - "The default is false to facilitate reproducibility.", default=.false.) + "If true, use the previous timestep MLD as a first guess in the MLD iteration, "//& + "otherwise use half the ocean depth as the first guess of the boundary layer "//& + "depth. The default is false to facilitate reproducibility.", & + default=.false., do_not_log=.not.CS%Use_MLD_iteration) call get_param(param_file, mdl, "EPBL_MLD_TOLERANCE", CS%MLD_tol, & "The tolerance for the iteratively determined mixed "//& "layer depth. This is only used with USE_MLD_ITERATION.", & - units="meter", default=1.0, scale=US%m_to_Z) + units="meter", default=1.0, scale=US%m_to_Z, do_not_log=.not.CS%Use_MLD_iteration) + call get_param(param_file, mdl, "EPBL_MLD_BISECTION", CS%MLD_bisection, & + "If true, use bisection with the iterative determination of the self-consistent "//& + "mixed layer depth. Otherwise use the false position after a maximum and minimum "//& + "bound have been evaluated and the returned value or bisection before this.", & + default=.true., do_not_log=.not.CS%Use_MLD_iteration) !### The default should become false. call get_param(param_file, mdl, "EPBL_MLD_MAX_ITS", CS%max_MLD_its, & "The maximum number of iterations that can be used to find a self-consistent "//& - "mixed layer depth. For now, due to the use of bisection, the maximum number "//& + "mixed layer depth. If EPBL_MLD_BISECTION is true, the maximum number "//& "iteractions needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.", & default=20, do_not_log=.not.CS%Use_MLD_iteration) if (.not.CS%Use_MLD_iteration) CS%Max_MLD_Its = 1 @@ -2339,6 +2379,10 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "If true, temperature and salinity are used as state "//& "variables.", default=.true.) + if (report_avg_its) then + CS%sum_its(1) = real_to_EFP(0.0) ; CS%sum_its(2) = real_to_EFP(0.0) + endif + if (max(CS%id_TKE_wind, CS%id_TKE_MKE, CS%id_TKE_conv, & CS%id_TKE_mixing, CS%id_TKE_mech_decay, CS%id_TKE_forcing, & CS%id_TKE_conv_decay) > 0) then @@ -2370,6 +2414,9 @@ subroutine energetic_PBL_end(CS) type(energetic_PBL_CS), pointer :: CS !< Energetic_PBL control structure that !! will be deallocated in this subroutine. + character(len=256) :: mesg + real :: avg_its + if (.not.associated(CS)) return if (allocated(CS%ML_depth)) deallocate(CS%ML_depth) @@ -2387,6 +2434,14 @@ subroutine energetic_PBL_end(CS) if (allocated(CS%Mixing_Length)) deallocate(CS%Mixing_Length) if (allocated(CS%Velocity_Scale)) deallocate(CS%Velocity_Scale) + if (report_avg_its) then + call EFP_sum_across_PEs(CS%sum_its, 2) + + avg_its = EFP_to_real(CS%sum_its(1)) / EFP_to_real(CS%sum_its(2)) + write (mesg,*) "Average ePBL iterations = ", avg_its + call MOM_mesg(mesg) + endif + deallocate(CS) end subroutine energetic_PBL_end From f90a4151d477f9145aa72ab19fff1430f410122b Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Sun, 10 May 2020 07:05:05 -0400 Subject: [PATCH 15/40] Corrected dOxygenization of two comments --- src/parameterizations/vertical/MOM_energetic_PBL.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index 2a66c60c88..e2fda25f6d 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -52,7 +52,7 @@ module MOM_energetic_PBL !! surface boundary layer to constrain the mixing lengths. logical :: MLD_iteration_guess !< False to default to guessing half the !! ocean depth for the first iteration. - logical :: MLD_bisection !! If true, use bisection with the iterative determination of the + logical :: MLD_bisection !< If true, use bisection with the iterative determination of the !! self-consistent mixed layer depth. Otherwise use the false position !! after a maximum and minimum bound have been evaluated and the !! returned value from the previous guess or bisection before this. @@ -223,7 +223,7 @@ module MOM_energetic_PBL character*(20), parameter :: ADDITIVE_STRING = "ADDITIVE" !>@} -logical :: report_avg_its = .false. ! Report the average number of ePBL iterations for debugging. +logical :: report_avg_its = .false. !< Report the average number of ePBL iterations for debugging. !> A type for conveniently passing around ePBL diagnostics for a column. type, public :: ePBL_column_diags ; private From badd68186080135bc086c0f06c80f74357e77286 Mon Sep 17 00:00:00 2001 From: Hemant Khatri Date: Mon, 11 May 2020 13:44:57 -0400 Subject: [PATCH 16/40] Changed the division by dt to multiplication form --- src/core/MOM_barotropic.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 9eb7a9d2e0..c2a9b3c8b1 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -2315,13 +2315,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif if (CS%id_ubtdt > 0) then do j=js,je ; do I=is-1,ie - ubt_dt(I,j) = (ubt_wtd(I,j) - ubt_st(I,j))/dt + ubt_dt(I,j) = (ubt_wtd(I,j) - ubt_st(I,j))*Idt enddo ; enddo call post_data(CS%id_ubtdt, ubt_dt(IsdB:IedB,jsd:jed), CS%diag) endif if (CS%id_vbtdt > 0) then do J=js-1,je ; do i=is,ie - vbt_dt(i,J) = (vbt_wtd(i,J) - vbt_st(i,J))/dt + vbt_dt(i,J) = (vbt_wtd(i,J) - vbt_st(i,J))*Idt enddo ; enddo call post_data(CS%id_vbtdt, vbt_dt(isd:ied,JsdB:JedB), CS%diag) endif From 6ca31848ea77c8868fb179ef2f136a091f7540a6 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 11 May 2020 14:33:49 -0400 Subject: [PATCH 17/40] Replaced array syntax additions in MOM_ice_shelf Replaced array syntax arithmetic with explicit loops in shelf_calc_flux, so that uninitialized values in halo points could not trigger model failures. All answers are bitwise identical. --- src/ice_shelf/MOM_ice_shelf.F90 | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index 891d6b3ea7..8d012377a5 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -607,15 +607,16 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) ISS%tflux_ocn(i,j) = 0.0 endif -! haline_driving(:,:) = sfc_state%sss(i,j) - Sbdry(i,j) +! haline_driving(i,j) = sfc_state%sss(i,j) - Sbdry(i,j) enddo ! i-loop enddo ! j-loop - ! ISS%water_flux = net liquid water into the ocean [R Z T-1 ~> kg m-2 s-1] - fluxes%iceshelf_melt(:,:) = ISS%water_flux(:,:) * CS%flux_factor do j=js,je ; do i=is,ie + ! ISS%water_flux = net liquid water into the ocean [R Z T-1 ~> kg m-2 s-1] + fluxes%iceshelf_melt(i,j) = ISS%water_flux(i,j) * CS%flux_factor + if ((sfc_state%ocean_mass(i,j) > CS%col_mass_melt_threshold) .and. & (ISS%area_shelf_h(i,j) > 0.0) .and. (CS%isthermo)) then @@ -653,11 +654,10 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) ISS%water_flux(i,j) = 0.0 fluxes%iceshelf_melt(i,j) = 0.0 endif ! area_shelf_h - enddo ; enddo ! i- and j-loops - ! mass flux [R Z L2 T-1 ~> kg s-1], part of ISOMIP diags. - mass_flux(:,:) = 0.0 - mass_flux(:,:) = ISS%water_flux(:,:) * ISS%area_shelf_h(:,:) + ! mass flux [R Z L2 T-1 ~> kg s-1], part of ISOMIP diags. + mass_flux(i,j) = ISS%water_flux(i,j) * ISS%area_shelf_h(i,j) + enddo ; enddo ! i- and j-loops if (CS%active_shelf_dynamics .or. CS%override_shelf_movement) then call cpu_clock_begin(id_clock_pass) @@ -690,7 +690,7 @@ subroutine shelf_calc_flux(sfc_state, fluxes, Time, time_step, CS, forces) ! advect the ice shelf, and advance the front. Calving will be in here somewhere as well.. ! when we decide on how to do it call update_ice_shelf(CS%dCS, ISS, G, US, US%s_to_T*time_step, Time, & - sfc_state%ocean_mass(:,:), coupled_GL) + sfc_state%ocean_mass, coupled_GL) endif @@ -1735,7 +1735,9 @@ subroutine update_shelf_mass(G, US, CS, ISS, Time) call time_interp_external(CS%id_read_mass, Time, ISS%mass_shelf) ! This should only be done if time_interp_external did an update. - ISS%mass_shelf(:,:) = US%kg_m3_to_R*US%m_to_Z * ISS%mass_shelf(:,:) ! Rescale after time_interp + do j=js,je ; do i=is,ie + ISS%mass_shelf(i,j) = US%kg_m3_to_R*US%m_to_Z * ISS%mass_shelf(i,j) ! Rescale after time_interp + enddo ; enddo do j=js,je ; do i=is,ie ISS%area_shelf_h(i,j) = 0.0 From 9f7587a57fe3ac3b045a433bdf04fdeedc8b2aaa Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 11 May 2020 14:38:13 -0400 Subject: [PATCH 18/40] Replaced array syntax in MOM_internal_tides Replaced array syntax sums in MOM_internal_tides with explicit loops. Also documented internal variable units in MOM_internal_tides, and noted an incorrect expression in PPM_angular_advect with a comment (with '###') and a suggested correction. All answers are bitwise identical because the bug was noted but not corrected. --- .../lateral/MOM_internal_tides.F90 | 225 +++++++++--------- 1 file changed, 115 insertions(+), 110 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index dda892dc3e..6145fb1dce 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -558,7 +558,7 @@ subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, & ! Output 2-D energy loss (summed over angles) for each freq and mode do m=1,CS%NMode ; do fr=1,CS%Nfreq if (CS%id_itidal_loss_mode(fr,m) > 0 .or. CS%id_allprocesses_loss_mode(fr,m) > 0) then - itidal_loss_mode(:,:) = 0.0 ! wave-drag processes (could do others as well) + itidal_loss_mode(:,:) = 0.0 ! wave-drag processes (could do others as well) allprocesses_loss_mode(:,:) = 0.0 ! all processes summed together do a=1,CS%nAngle ; do j=js,je ; do i=is,ie itidal_loss_mode(i,j) = itidal_loss_mode(i,j) + CS%TKE_itidal_loss(i,j,a,fr,m) @@ -886,12 +886,16 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) !! across angles [R Z3 T-2 ~> J m-2]. ! Local variables real :: flux - real :: u_ang - real :: Angle_size - real :: I_Angle_size - real :: I_dt + real :: u_ang ! Angular propagation speed [Rad T-1 ~> Rad s-1] + real :: Angle_size ! The size of each orientation wedge in radians [Rad] + real :: I_Angle_size ! The inverse of the the orientation wedges [Rad-1] + real :: I_dt ! The inverse of the timestep [T-1 ~> s-1] + real :: aR, aL ! Left and right edge estimates of energy density [R Z3 T-2 rad-1 ~> J m-2 rad-1] + real :: dMx, dMn + real :: Ep, Ec, Em ! Mean angular energy density for three successive wedges in angular + ! orientation [R Z3 T-2 rad-1 ~> J m-2 rad-1] + real :: dA, mA, a6 ! Difference, mean, and curvature of energy density [R Z3 T-2 rad-1 ~> J m-2 rad-1] integer :: a - real :: aR, aL, dMx, dMn, Ep, Ec, Em, dA, mA, a6 I_dt = 1 / dt Angle_size = (8.0*atan(1.0)) / (real(NAngle)) @@ -902,9 +906,12 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) u_ang = CFL_ang(A)*Angle_size*I_dt if (u_ang >= 0.0) then ! Implementation of PPM-H3 - Ep = En2d(a+1)*I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - Ec = En2d(a) *I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - Em = En2d(a-1)*I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) + ! Convert wedge-integrated energy density into angular energy densities for three successive + ! wedges around the source wedge for this flux [R Z3 T-2 rad-1 ~> J m-2 rad-1]. + Ep = En2d(a+1)*I_Angle_size + Ec = En2d(a) *I_Angle_size + Em = En2d(a-1)*I_Angle_size + ! Calculate and bound edge values of energy density. aL = ( 5.*Ec + ( 2.*Em - Ep ) )/6. ! H3 estimate aL = max( min(Ec,Em), aL) ; aL = min( max(Ec,Em), aL) ! Bound aR = ( 5.*Ec + ( 2.*Ep - Em ) )/6. ! H3 estimate @@ -918,17 +925,21 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) aR = 3.*Ec - 2.*aL !? endif a6 = 6.*Ec - 3. * (aR + aL) ! Curvature - ! CALCULATE FLUX RATE (Jm-2/s) + ! Calculate angular flux rate [R Z3 T-3 ~> W m-2] flux = u_ang*( aR + 0.5 * CFL_ang(A) * ( ( aL - aR ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - !flux = u_ang*( aR - 0.5 * CFL_ang(A) * ( ( aR - aL ) - a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - ! CALCULATE AMOUNT FLUXED (Jm-2) + ! The following expression copied from tracer_advect is equivalent. + ! flux = u_ang*( aR - 0.5 * CFL_ang(A) * ( ( aR - aL ) - a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) + ! Calculate amount of energy fluxed between wedges [R Z3 T-2 ~> J m-2] Flux_En(A) = dt * flux !Flux_En(A) = (dt * I_Angle_size) * flux else ! Implementation of PPM-H3 - Ep = En2d(a+2)*I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - Ec = En2d(a+1)*I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) - Em = En2d(a) *I_Angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad) + ! Convert wedge-integrated energy density into angular energy densities for three successive + ! wedges around the source wedge for this flux [R Z3 T-2 rad-1 ~> J m-2 rad-1]. + Ep = En2d(a+2)*I_Angle_size + Ec = En2d(a+1)*I_Angle_size + Em = En2d(a) *I_Angle_size + ! Calculate and bound edge values of energy density. aL = ( 5.*Ec + ( 2.*Em - Ep ) )/6. ! H3 estimate aL = max( min(Ec,Em), aL) ; aL = min( max(Ec,Em), aL) ! Bound aR = ( 5.*Ec + ( 2.*Ep - Em ) )/6. ! H3 estimate @@ -942,10 +953,12 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) aR = 3.*Ec - 2.*aL endif a6 = 6.*Ec - 3. * (aR + aL) ! Curvature - ! CALCULATE FLUX RATE (Jm-2/s) + ! Calculate angular flux rate [R Z3 T-3 ~> W m-2] + !### This expression is wrong, because it was just copied from above. The correct one is below flux = u_ang*( aR + 0.5 * CFL_ang(A) * ( ( aL - aR ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - !flux = u_ang*( aL + 0.5 * CFL_ang(A) * ( ( aR - aL ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - ! CALCULATE AMOUNT FLUXED (Jm-2) + ! This is the correct expression; note that CFL_ang is negative here, so it looks a bit odd. + !flux = u_ang*( aL - 0.5 * CFL_ang(A) * ( ( aR - aL ) + a6 * ( 1. + 2./3. * CFL_ang(A) ) ) ) + ! Calculate amount of energy fluxed between wedges [R Z3 T-2 ~> J m-2] Flux_En(A) = dt * flux !Flux_En(A) = (dt * I_Angle_size) * flux endif @@ -1014,7 +1027,7 @@ subroutine propagate(En, cn, freq, dt, G, US, CS, NAngle) ! FIND AVERAGE GROUP VELOCITY (SPEED) AT CELL CORNERS ! NOTE: THIS HAS NOT BE ADAPTED FOR REFLECTION YET (BDM)!! ! Fix indexing here later - speed(:,:) = 0 + speed(:,:) = 0.0 do J=jsh-1,jeh ; do I=ish-1,ieh f2 = G%CoriolisBu(I,J)**2 speed(I,J) = 0.25*(cn(i,j) + cn(i+1,j) + cn(i+1,j+1) + cn(i,j+1)) * & @@ -1058,21 +1071,21 @@ subroutine propagate(En, cn, freq, dt, G, US, CS, NAngle) ! Apply propagation in x-direction (reflection included) LB%jsh = jsh ; LB%jeh = jeh ; LB%ish = ish ; LB%ieh = ieh - call propagate_x(En(:,:,:), speed_x, Cgx_av(:), dCgx(:), dt, G, US, CS%nAngle, CS, LB) + call propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, US, CS%nAngle, CS, LB) ! Check for energy conservation on computational domain (for debugging) - !call sum_En(G,CS,En(:,:,:),'post-propagate_x') + !call sum_En(G, CS, En, 'post-propagate_x') ! Update halos - call pass_var(En(:,:,:),G%domain) + call pass_var(En, G%domain) ! Apply propagation in y-direction (reflection included) ! LB%jsh = js ; LB%jeh = je ; LB%ish = is ; LB%ieh = ie ! Use if no teleport LB%jsh = jsh ; LB%jeh = jeh ; LB%ish = ish ; LB%ieh = ieh - call propagate_y(En(:,:,:), speed_y, Cgy_av(:), dCgy(:), dt, G, US, CS%nAngle, CS, LB) + call propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, US, CS%nAngle, CS, LB) ! Check for energy conservation on computational domain (for debugging) - !call sum_En(G,CS,En(:,:,:),'post-propagate_y') + !call sum_En(G, CS, En, 'post-propagate_y') endif end subroutine propagate @@ -1084,7 +1097,7 @@ subroutine propagate_corner_spread(En, energized_wedge, NAngle, speed, dt, G, CS type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. real, dimension(G%isd:G%ied,G%jsd:G%jed), & intent(inout) :: En !< The energy density integrated over an angular - !! band [R Z3 T-2 ~> J m-2], intent in/out. + !! band [R Z3 T-2 ~> J m-2]. real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed), & intent(in) :: speed !< The magnitude of the group velocity at the cell !! corner points [L T-1 ~> m s-1]. @@ -1351,7 +1364,7 @@ subroutine propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, US, Nangle, CS, LB) !! discretized wave energy spectrum. real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), & intent(inout) :: En !< The energy density integrated over an angular - !! band [R Z3 T-2 ~> J m-2], intent in/out. + !! band [R Z3 T-2 ~> J m-2]. real, dimension(G%IsdB:G%IedB,G%jsd:G%jed), & intent(in) :: speed_x !< The magnitude of the group velocity at the !! Cu points [L T-1 ~> m s-1]. @@ -1404,18 +1417,18 @@ subroutine propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, US, Nangle, CS, LB) enddo ! a-loop ! Only reflect newly arrived energy; existing energy in incident wedge is not reflected - ! and will eventually propagate out of cell. (Thid code only reflects if En > 0) - call reflect(Fdt_m(:,:,:), Nangle, CS, G, LB) - call teleport(Fdt_m(:,:,:), Nangle, CS, G, LB) - call reflect(Fdt_p(:,:,:), Nangle, CS, G, LB) - call teleport(Fdt_p(:,:,:), Nangle, CS, G, LB) + ! and will eventually propagate out of cell. (This code only reflects if En > 0.) + call reflect(Fdt_m, Nangle, CS, G, LB) + call teleport(Fdt_m, Nangle, CS, G, LB) + call reflect(Fdt_p, Nangle, CS, G, LB) + call teleport(Fdt_p, Nangle, CS, G, LB) - ! Update reflected energy (Jm-2) - do j=jsh,jeh ; do i=ish,ieh + ! Update reflected energy [R Z3 T-2 ~> J m-2] + do a=1,Nangle ; do j=jsh,jeh ; do i=ish,ieh ! if ((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0) & ! for debugging ! call MOM_error(FATAL, "propagate_x: OutFlux>Available") - En(i,j,:) = En(i,j,:) + G%IareaT(i,j)*(Fdt_m(i,j,:) + Fdt_p(i,j,:)) - enddo ; enddo + En(i,j,a) = En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a)) + enddo ; enddo ; enddo end subroutine propagate_x @@ -1426,7 +1439,7 @@ subroutine propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, US, Nangle, CS, LB) !! discretized wave energy spectrum. real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), & intent(inout) :: En !< The energy density integrated over an angular - !! band [R Z3 T-2 ~> J m-2], intent in/out. + !! band [R Z3 T-2 ~> J m-2]. real, dimension(G%isd:G%ied,G%JsdB:G%JedB), & intent(in) :: speed_y !< The magnitude of the group velocity at the !! Cv points [L T-1 ~> m s-1]. @@ -1486,13 +1499,13 @@ subroutine propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, US, Nangle, CS, LB) enddo ! a-loop ! Only reflect newly arrived energy; existing energy in incident wedge is not reflected - ! and will eventually propagate out of cell. (Thid code only reflects if En > 0) - call reflect(Fdt_m(:,:,:), Nangle, CS, G, LB) - call teleport(Fdt_m(:,:,:), Nangle, CS, G, LB) - call reflect(Fdt_p(:,:,:), Nangle, CS, G, LB) - call teleport(Fdt_p(:,:,:), Nangle, CS, G, LB) + ! and will eventually propagate out of cell. (This code only reflects if En > 0.) + call reflect(Fdt_m, Nangle, CS, G, LB) + call teleport(Fdt_m, Nangle, CS, G, LB) + call reflect(Fdt_p, Nangle, CS, G, LB) + call teleport(Fdt_p, Nangle, CS, G, LB) - ! Update reflected energy (Jm-2) + ! Update reflected energy [R Z3 T-2 ~> J m-2] do a=1,Nangle ; do j=jsh,jeh ; do i=ish,ieh ! if ((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0) & ! for debugging ! call MOM_error(FATAL, "propagate_y: OutFlux>Available", .true.) @@ -1521,8 +1534,7 @@ subroutine zonal_flux_En(u, h, hL, hR, uh, dt, G, US, j, ish, ieh, vol_CFL) !! the cell areas when estimating the CFL number. ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]. - real :: curv_3 ! A measure of the thickness curvature over a grid length, - ! with the same units as h_in. + real :: curv_3 ! A measure of the energy density curvature over a grid length [R Z3 T-2 ~> J m-2] integer :: i do I=ish-1,ieh @@ -1566,8 +1578,7 @@ subroutine merid_flux_En(v, h, hL, hR, vh, dt, G, US, J, ish, ieh, vol_CFL) !! the CFL number. ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim]. - real :: curv_3 ! A measure of the thickness curvature over a grid length, - ! with the same units as h_in. + real :: curv_3 ! A measure of the energy density curvature over a grid length [R Z3 T-2 ~> J m-2] integer :: i do i=ish,ieh @@ -1603,18 +1614,18 @@ subroutine reflect(En, NAngle, CS, G, LB) type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds. ! Local variables real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c - ! angle of boudary wrt equator + ! angle of boundary wrt equator [rad] real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl ! fraction of wave energy reflected - ! values should collocate with angle_c + ! values should collocate with angle_c [nondim] logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge ! tags of cells with double reflection - real :: TwoPi ! 2*pi - real :: Angle_size ! size of beam wedge (rad) - real :: angle_wall ! angle of coast/ridge/shelf wrt equator - real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator - real :: angle_r ! angle of reflected ray wrt equator + real :: TwoPi ! 2*pi = 6.2831853... [nondim] + real :: Angle_size ! size of beam wedge [rad] + real :: angle_wall ! angle of coast/ridge/shelf wrt equator [rad] + real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator [rad] + real :: angle_r ! angle of reflected ray wrt equator [rad] real, dimension(1:Nangle) :: En_reflected integer :: i, j, a, a_r, na !integer :: isd, ied, jsd, jed ! start and end local indices on data domain @@ -1623,7 +1634,6 @@ subroutine reflect(En, NAngle, CS, G, LB) ! (values exclude halos) integer :: ish, ieh, jsh, jeh ! start and end local indices on data domain ! leaving out outdated halo points (march in) - integer :: id_g, jd_g ! global (decomp-invar) indices !isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec @@ -1643,59 +1653,54 @@ subroutine reflect(En, NAngle, CS, G, LB) ridge = CS%refl_dbl En_reflected(:) = 0.0 - !do j=jsc-1,jec+1 - do j=jsh,jeh - !do i=isc-1,iec+1 - do i=ish,ieh - ! jd_g = j + G%jdg_offset ; id_g = i + G%idg_offset - ! redistribute energy in angular space if ray will hit boundary - ! i.e., if energy is in a reflecting cell - if (angle_c(i,j) /= CS%nullangle) then - do a=1,NAngle - if (En(i,j,a) > 0.0) then - ! if ray is incident, keep specified boundary angle - if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then - angle_wall = angle_c(i,j) - ! if ray is not incident but in ridge cell, use complementary angle - elseif (ridge(i,j)) then - angle_wall = angle_c(i,j) + 0.5*TwoPi - if (angle_wall > TwoPi) then - angle_wall = angle_wall - TwoPi*floor(abs(angle_wall)/TwoPi) - elseif (angle_wall < 0.0) then - angle_wall = angle_wall + TwoPi*ceiling(abs(angle_wall)/TwoPi) - endif - ! if ray is not incident and not in a ridge cell, keep specified angle - else - angle_wall = angle_c(i,j) - endif - ! do reflection - if (sin(angle_i(a) - angle_wall) >= 0.0) then - angle_r = 2.0*angle_wall - angle_i(a) - if (angle_r > TwoPi) then - angle_r = angle_r - TwoPi*floor(abs(angle_r)/TwoPi) - elseif (angle_r < 0.0) then - angle_r = angle_r + TwoPi*ceiling(abs(angle_r)/TwoPi) - endif - a_r = nint(angle_r/Angle_size) + 1 - do while (a_r > Nangle) ; a_r = a_r - Nangle ; enddo - if (a /= a_r) then - En_reflected(a_r) = part_refl(i,j)*En(i,j,a) - En(i,j,a) = (1.0-part_refl(i,j))*En(i,j,a) - endif - endif + do j=jsh,jeh ; do i=ish,ieh + ! redistribute energy in angular space if ray will hit boundary + ! i.e., if energy is in a reflecting cell + if (angle_c(i,j) /= CS%nullangle) then + do a=1,NAngle ; if (En(i,j,a) > 0.0) then + if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then + ! if ray is incident, keep specified boundary angle + angle_wall = angle_c(i,j) + elseif (ridge(i,j)) then + ! if ray is not incident but in ridge cell, use complementary angle + angle_wall = angle_c(i,j) + 0.5*TwoPi + if (angle_wall > TwoPi) then + angle_wall = angle_wall - TwoPi*floor(abs(angle_wall)/TwoPi) + elseif (angle_wall < 0.0) then + angle_wall = angle_wall + TwoPi*ceiling(abs(angle_wall)/TwoPi) endif - enddo ! a-loop - En(i,j,:) = En(i,j,:) + En_reflected(:) - En_reflected(:) = 0.0 - endif - enddo ! i-loop - enddo ! j-loop + else + ! if ray is not incident and not in a ridge cell, keep specified angle + angle_wall = angle_c(i,j) + endif + + ! do reflection + if (sin(angle_i(a) - angle_wall) >= 0.0) then + angle_r = 2.0*angle_wall - angle_i(a) + if (angle_r > TwoPi) then + angle_r = angle_r - TwoPi*floor(abs(angle_r)/TwoPi) + elseif (angle_r < 0.0) then + angle_r = angle_r + TwoPi*ceiling(abs(angle_r)/TwoPi) + endif + a_r = nint(angle_r/Angle_size) + 1 + do while (a_r > Nangle) ; a_r = a_r - Nangle ; enddo + if (a /= a_r) then + En_reflected(a_r) = part_refl(i,j)*En(i,j,a) + En(i,j,a) = (1.0-part_refl(i,j))*En(i,j,a) + endif + endif + endif ; enddo ! a-loop + do a=1,NAngle + En(i,j,a) = En(i,j,a) + En_reflected(a) + En_reflected(a) = 0.0 + enddo ! a-loop + endif + enddo ; enddo ! i- and j-loops ! Check to make sure no energy gets onto land (only run for debugging) ! do a=1,NAngle ; do j=jsc,jec ; do i=isc,iec ! if (En(i,j,a) > 0.001 .and. G%mask2dT(i,j) == 0) then - ! jd_g = j + G%jdg_offset ; id_g = i + G%idg_offset - ! write (mesg,*) 'En=', En(i,j,a), 'a=', a, 'ig_g=',id_g, 'jg_g=',jd_g + ! write (mesg,*) 'En=', En(i,j,a), 'a=', a, 'ig_g=',i+G%idg_offset, 'jg_g=',j+G%jdg_offset ! call MOM_error(FATAL, "reflect: Energy detected out of bounds: "//trim(mesg), .true.) ! endif ! enddo ; enddo ; enddo @@ -1717,17 +1722,17 @@ subroutine teleport(En, NAngle, CS, G, LB) type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds. ! Local variables real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c - ! angle of boudary wrt equator + ! angle of boundary wrt equator [rad] real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl ! fraction of wave energy reflected - ! values should collocate with angle_c + ! values should collocate with angle_c [nondim] logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: pref_cell ! flag for partial reflection logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge - ! tags of cells with double reflection - real :: TwoPi ! size of beam wedge (rad) - real :: Angle_size ! size of beam wedge (rad) - real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator + ! tags of cells with double reflection + real :: TwoPi ! 2*pi = 6.2831853... [nondim] + real :: Angle_size ! size of beam wedge [rad] + real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator [rad] real, dimension(1:NAngle) :: cos_angle, sin_angle real :: En_tele ! energy to be "teleported" [R Z3 T-2 ~> J m-2] character(len=160) :: mesg ! The text of an error message @@ -2295,8 +2300,8 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) CS%TKE_itidal_loss(:,:,:,:,:) = 0.0 allocate(CS%TKE_Froude_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode)) CS%TKE_Froude_loss(:,:,:,:,:) = 0.0 - allocate(CS%tot_leak_loss(isd:ied,jsd:jed)) ; CS%tot_leak_loss(:,:) = 0.0 - allocate(CS%tot_quad_loss(isd:ied,jsd:jed) ) ; CS%tot_quad_loss(:,:) = 0.0 + allocate(CS%tot_leak_loss(isd:ied,jsd:jed)) ; CS%tot_leak_loss(:,:) = 0.0 + allocate(CS%tot_quad_loss(isd:ied,jsd:jed) ) ; CS%tot_quad_loss(:,:) = 0.0 allocate(CS%tot_itidal_loss(isd:ied,jsd:jed)) ; CS%tot_itidal_loss(:,:) = 0.0 allocate(CS%tot_Froude_loss(isd:ied,jsd:jed)) ; CS%tot_Froude_loss(:,:) = 0.0 From 34196a2eb4cbd684cd02ec9e3b0a18e1e423eef1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 11 May 2020 14:47:21 -0400 Subject: [PATCH 19/40] Added array-syntax notation for full-array copies Added array-syntax notation for full-array copies in offline_diabatic_ale. All answers are bitwise identical. --- src/tracer/MOM_offline_main.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index 65f83ecfea..a0f7b23346 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -722,9 +722,9 @@ subroutine offline_diabatic_ale(fluxes, Time_start, Time_end, CS, h_pre, eatr, e ! Add diurnal cycle for shortwave radiation (only used if run in ocean-only mode) if (CS%diurnal_SW .and. CS%read_sw) then - sw(:,:) = fluxes%sw - sw_vis(:,:) = fluxes%sw_vis_dir - sw_nir(:,:) = fluxes%sw_nir_dir + sw(:,:) = fluxes%sw(:,:) + sw_vis(:,:) = fluxes%sw_vis_dir(:,:) + sw_nir(:,:) = fluxes%sw_nir_dir(:,:) call offline_add_diurnal_SW(fluxes, CS%G, Time_start, Time_end) endif @@ -738,9 +738,9 @@ subroutine offline_diabatic_ale(fluxes, Time_start, Time_end, CS, h_pre, eatr, e CS%G, CS%GV, CS%US, CS%tv, CS%optics, CS%tracer_flow_CSp, CS%debug) if (CS%diurnal_SW .and. CS%read_sw) then - fluxes%sw(:,:) = sw - fluxes%sw_vis_dir(:,:) = sw_vis - fluxes%sw_nir_dir(:,:) = sw_nir + fluxes%sw(:,:) = sw(:,:) + fluxes%sw_vis_dir(:,:) = sw_vis(:,:) + fluxes%sw_nir_dir(:,:) = sw_nir(:,:) endif if (CS%debug) then From e68a66489d692ac21a746a9684dd716edcf84e31 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 11 May 2020 14:50:14 -0400 Subject: [PATCH 20/40] Removed duplicate dimension declarations Removed duplicated dimension declarations for two ppoly variables in bulk_average in MOM_lateral_boundary_diffusion. When I first saw these declarations, they were confusing to me, as I was unsure at first whether they are actually declaring 4-d arrays or 2-d arrays (it is the latter). Also removed unneeded full array syntax in 2 subroutine calls, and merged a pair of do-loop statements with common loop contents. All answers are bitwise identical. --- src/tracer/MOM_lateral_boundary_diffusion.F90 | 38 +++++++++---------- 1 file changed, 18 insertions(+), 20 deletions(-) diff --git a/src/tracer/MOM_lateral_boundary_diffusion.F90 b/src/tracer/MOM_lateral_boundary_diffusion.F90 index 8b9be533d5..f244931376 100644 --- a/src/tracer/MOM_lateral_boundary_diffusion.F90 +++ b/src/tracer/MOM_lateral_boundary_diffusion.F90 @@ -171,16 +171,14 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) ! for diagnostics if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0) then - tendency(:,:,:) = 0.0 + tendency(:,:,:) = 0.0 endif - do j = G%jsc-1, G%jec+1 - ! Interpolate state to interface - do i = G%isc-1, G%iec+1 - call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), tracer%t(i,j,:), ppoly0_coefs(i,j,:,:), & - ppoly0_E(i,j,:,:), ppoly_S, remap_method, GV%H_subroundoff, GV%H_subroundoff) - enddo - enddo + ! Interpolate state to interface + do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 + call build_reconstructions_1d( CS%remap_CS, G%ke, h(i,j,:), tracer%t(i,j,:), ppoly0_coefs(i,j,:,:), & + ppoly0_E(i,j,:,:), ppoly_S, remap_method, GV%H_subroundoff, GV%H_subroundoff) + enddo ; enddo ! Diffusive fluxes in the i-direction uFlx(:,:,:) = 0. vFlx(:,:,:) = 0. @@ -253,41 +251,41 @@ subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS) if (tracer%id_lbd_dfy>0) call post_data(tracer%id_lbd_dfy, vFlx*Idt, CS%diag) if (tracer%id_lbd_dfx_2d>0) then uwork_2d(:,:) = 0. - do k=1,GV%ke; do j=G%jsc,G%jec; do I=G%isc-1,G%iec + do k=1,GV%ke ; do j=G%jsc,G%jec ; do I=G%isc-1,G%iec uwork_2d(I,j) = uwork_2d(I,j) + (uFlx(I,j,k) * Idt) - enddo; enddo; enddo + enddo ; enddo ; enddo call post_data(tracer%id_lbd_dfx_2d, uwork_2d, CS%diag) endif if (tracer%id_lbd_dfy_2d>0) then vwork_2d(:,:) = 0. - do k=1,GV%ke; do J=G%jsc-1,G%jec; do i=G%isc,G%iec + do k=1,GV%ke ; do J=G%jsc-1,G%jec ; do i=G%isc,G%iec vwork_2d(i,J) = vwork_2d(i,J) + (vFlx(i,J,k) * Idt) - enddo; enddo; enddo + enddo ; enddo ; enddo call post_data(tracer%id_lbd_dfy_2d, vwork_2d, CS%diag) endif ! post tendency of tracer content if (tracer%id_lbdxy_cont > 0) then - call post_data(tracer%id_lbdxy_cont, tendency(:,:,:), CS%diag) + call post_data(tracer%id_lbdxy_cont, tendency, CS%diag) endif ! post depth summed tendency for tracer content if (tracer%id_lbdxy_cont_2d > 0) then tendency_2d(:,:) = 0. - do j = G%jsc,G%jec ; do i = G%isc,G%iec - do k = 1, GV%ke + do j=G%jsc,G%jec ; do i=G%isc,G%iec + do k=1,GV%ke tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k) enddo enddo ; enddo - call post_data(tracer%id_lbdxy_cont_2d, tendency_2d(:,:), CS%diag) + call post_data(tracer%id_lbdxy_cont_2d, tendency_2d, CS%diag) endif ! post tendency of tracer concentration; this step must be ! done after posting tracer content tendency, since we alter ! the tendency array and its units. if (tracer%id_lbdxy_conc > 0) then - do k = 1, GV%ke ; do j = G%jsc,G%jec ; do i = G%isc,G%iec + do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + GV%H_subroundoff ) enddo ; enddo ; enddo call post_data(tracer%id_lbdxy_conc, tendency, CS%diag) @@ -306,9 +304,9 @@ real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coe real, dimension(nk) :: h !< Layer thicknesses [H ~> m or kg m-2] real :: hBLT !< Depth of the boundary layer [H ~> m or kg m-2] real, dimension(nk) :: phi !< Scalar quantity - real, dimension(nk,2) :: ppoly0_E(:,:) !< Edge value of polynomial - real, dimension(nk,deg+1) :: ppoly0_coefs(:,:) !< Coefficients of polynomial - integer :: method !< Remapping scheme to use + real, dimension(nk,2) :: ppoly0_E !< Edge value of polynomial + real, dimension(nk,deg+1) :: ppoly0_coefs !< Coefficients of polynomial + integer :: method !< Remapping scheme to use integer :: k_top !< Index of the first layer within the boundary real :: zeta_top !< Fraction of the layer encompassed by the bottom boundary layer From 21c33b10a5742a274e8a6399c7cc22488c8874f6 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 11 May 2020 14:54:53 -0400 Subject: [PATCH 21/40] Merge scaling factors when reading Nikurashin input Combined scaling factors when reading internal tide TKE input from a file for use with the Nikurashin mixing scheme, while also eliminating an array-syntax multiplication. All answers are bitwise identical. --- src/parameterizations/vertical/MOM_tidal_mixing.F90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 27b316e144..e7d5bcc476 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -509,8 +509,7 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, diag, CS) filename) call safe_alloc_ptr(CS%TKE_Niku,is,ie,js,je) ; CS%TKE_Niku(:,:) = 0.0 call MOM_read_data(filename, 'TKE_input', CS%TKE_Niku, G%domain, timelevel=1, & ! ??? timelevel -aja - scale=US%W_m2_to_RZ3_T3) - CS%TKE_Niku(:,:) = Niku_scale * CS%TKE_Niku(:,:) + scale=Niku_scale*US%W_m2_to_RZ3_T3) call get_param(param_file, mdl, "GAMMA_NIKURASHIN",CS%Gamma_lee, & "The fraction of the lee wave energy that is dissipated "//& @@ -781,7 +780,7 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kv) ! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable - do k = 1,G%ke+1 + do k=1,G%ke+1 N2_int_i(k) = US%s_to_T**2 * N2_int(i,k) enddo @@ -876,7 +875,7 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, US, CS, N2_int, Kd_lay, Kv) CVmix_tidal_params_user = CS%CVMix_tidal_params) ! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable - do k = 1,G%ke+1 + do k=1,G%ke+1 N2_int_i(k) = US%s_to_T**2 * N2_int(i,k) enddo From a6445be6f40fb7b2ba657428ab37949a08497491 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 11 May 2020 15:51:00 -0400 Subject: [PATCH 22/40] +Fixed pointer indexing in update_ALE_sponge_field Added grid type arguments to calls to update_ALE_sponge_field so that the internal array pointers set by this routine will use the same indexing conventions as the rest of the MOM6 code. Also added comments describing some arguments and other variables and got rid of some unneeded line continuations in MOM.F90. All answers are bitwise identical, but there are two new arguments to update_ALE_sponge_field. --- src/core/MOM.F90 | 19 ++-- .../vertical/MOM_ALE_sponge.F90 | 106 +++++++++--------- 2 files changed, 61 insertions(+), 64 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index aff4860a21..25ed1b1f6e 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1013,7 +1013,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & Time_local + real_to_time(US%T_to_s*(bbl_time_int-dt)), CS%diag) ! Calculate the BBL properties and store them inside visc (u,h). call cpu_clock_begin(id_clock_BBL_visc) - call set_viscous_BBL(CS%u(:,:,:), CS%v(:,:,:), CS%h, CS%tv, CS%visc, G, GV, US, & + call set_viscous_BBL(CS%u, CS%v, CS%h, CS%tv, CS%visc, G, GV, US, & CS%set_visc_CSp, symmetrize=.true.) call cpu_clock_end(id_clock_BBL_visc) if (showCallTree) call callTree_wayPoint("done with set_viscous_BBL (step_MOM)") @@ -2204,7 +2204,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & allocate(CS%tv%frazil(isd:ied,jsd:jed)) ; CS%tv%frazil(:,:) = 0.0 endif if (bound_salinity) then - allocate(CS%tv%salt_deficit(isd:ied,jsd:jed)) ; CS%tv%salt_deficit(:,:)=0.0 + allocate(CS%tv%salt_deficit(isd:ied,jsd:jed)) ; CS%tv%salt_deficit(:,:) = 0.0 endif if (bulkmixedlayer .or. use_temperature) then @@ -2369,20 +2369,17 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (associated(sponge_in_CSp)) then ! TODO: Implementation and testing of non-ALE spong rotation - call MOM_error(FATAL, "Index rotation of non-ALE sponge is not yet " & - // "implemented.") + call MOM_error(FATAL, "Index rotation of non-ALE sponge is not yet implemented.") endif if (associated(ALE_sponge_in_CSp)) then - call rotate_ALE_sponge(ALE_sponge_in_CSp, G_in, CS%ALE_sponge_CSp, G, & - turns, param_file) - call update_ALE_sponge_field(CS%ALE_sponge_CSp, T_in, CS%T) - call update_ALE_sponge_field(CS%ALE_sponge_CSp, S_in, CS%S) + call rotate_ALE_sponge(ALE_sponge_in_CSp, G_in, CS%ALE_sponge_CSp, G, turns, param_file) + call update_ALE_sponge_field(CS%ALE_sponge_CSp, T_in, G, GV, CS%T) + call update_ALE_sponge_field(CS%ALE_sponge_CSp, S_in, G, GV, CS%S) endif if (associated(OBC_in)) & - call rotate_OBC_init(OBC_in, G, GV, US, param_file, CS%tv, restart_CSp, & - CS%OBC) + call rotate_OBC_init(OBC_in, G, GV, US, param_file, CS%tv, restart_CSp, CS%OBC) deallocate(u_in) deallocate(v_in) @@ -2427,7 +2424,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & G => CS%G if (CS%debug .or. CS%G%symmetric) then call clone_MOM_domain(CS%G%Domain, CS%G%Domain_aux, symmetric=.false.) - else ; CS%G%Domain_aux => CS%G%Domain ;endif + else ; CS%G%Domain_aux => CS%G%Domain ; endif G%ke = GV%ke endif diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index b791535ed1..fe1ccab53d 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -171,7 +171,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ integer :: i, j, k, col, total_sponge_cols, total_sponge_cols_u, total_sponge_cols_v character(len=10) :: remapScheme if (associated(CS)) then - call MOM_error(WARNING, "initialize_sponge called with an associated "// & + call MOM_error(WARNING, "initialize_ALE_sponge_fixed called with an associated "// & "control structure.") return endif @@ -260,14 +260,14 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ if (CS%sponge_uv) then - allocate(data_hu(G%isdB:G%iedB,G%jsd:G%jed,nz_data)); data_hu(:,:,:)=0.0 - allocate(data_hv(G%isd:G%ied,G%jsdB:G%jedB,nz_data)); data_hv(:,:,:)=0.0 - allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)); Iresttime_u(:,:)=0.0 - allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)); Iresttime_v(:,:)=0.0 + allocate(data_hu(G%isdB:G%iedB,G%jsd:G%jed,nz_data)) ; data_hu(:,:,:) = 0.0 + allocate(data_hv(G%isd:G%ied,G%jsdB:G%jedB,nz_data)) ; data_hv(:,:,:) = 0.0 + allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)) ; Iresttime_u(:,:) = 0.0 + allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)) ; Iresttime_v(:,:) = 0.0 ! u points CS%num_col_u = 0 ; !CS%fldno_u = 0 - do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB + do j=CS%jsc,CS%jec ; do I=CS%iscB,CS%iecB data_hu(I,j,:) = 0.5 * (data_h(i,j,:) + data_h(i+1,j,:)) Iresttime_u(I,j) = 0.5 * (Iresttime(i,j) + Iresttime(i+1,j)) if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) & @@ -276,9 +276,9 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ if (CS%num_col_u > 0) then - allocate(CS%Iresttime_col_u(CS%num_col_u)) ; CS%Iresttime_col_u = 0.0 - allocate(CS%col_i_u(CS%num_col_u)) ; CS%col_i_u = 0 - allocate(CS%col_j_u(CS%num_col_u)) ; CS%col_j_u = 0 + allocate(CS%Iresttime_col_u(CS%num_col_u)) ; CS%Iresttime_col_u(:) = 0.0 + allocate(CS%col_i_u(CS%num_col_u)) ; CS%col_i_u(:) = 0 + allocate(CS%col_j_u(CS%num_col_u)) ; CS%col_j_u(:) = 0 ! pass indices, restoring time to the CS structure col = 1 @@ -286,7 +286,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ if ((Iresttime_u(I,j)>0.0) .and. (G%mask2dCu(I,j)>0)) then CS%col_i_u(col) = i ; CS%col_j_u(col) = j CS%Iresttime_col_u(col) = Iresttime_u(i,j) - col = col +1 + col = col + 1 endif enddo ; enddo @@ -323,7 +323,7 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, param_file, CS, data_h, nz_ if ((Iresttime_v(i,J)>0.0) .and. (G%mask2dCv(i,J)>0)) then CS%col_i_v(col) = i ; CS%col_j_v(col) = j CS%Iresttime_col_v(col) = Iresttime_v(i,j) - col = col +1 + col = col + 1 endif enddo ; enddo @@ -415,7 +415,7 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) character(len=10) :: remapScheme if (associated(CS)) then - call MOM_error(WARNING, "initialize_sponge called with an associated "// & + call MOM_error(WARNING, "initialize_ALE_sponge_varying called with an associated "// & "control structure.") return endif @@ -486,8 +486,8 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, param_file, CS) call log_param(param_file, mdl, "!Total sponge columns at h points", total_sponge_cols, & "The total number of columns where sponges are applied at h points.") if (CS%sponge_uv) then - allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)); Iresttime_u(:,:)=0.0 - allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)); Iresttime_v(:,:)=0.0 + allocate(Iresttime_u(G%isdB:G%iedB,G%jsd:G%jed)) ; Iresttime_u(:,:) = 0.0 + allocate(Iresttime_v(G%isd:G%ied,G%jsdB:G%jedB)) ; Iresttime_v(:,:) = 0.0 ! u points CS%num_col_u = 0 ; !CS%fldno_u = 0 do j=CS%jsc,CS%jec; do I=CS%iscB,CS%iecB @@ -578,7 +578,7 @@ subroutine set_up_ALE_sponge_field_fixed(sp_val, G, f_ptr, CS) if (CS%fldno > MAX_FIELDS_) then write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease & &the number of fields to be damped in the call to & - &initialize_sponge." )') CS%fldno + &initialize_ALE_sponge." )') CS%fldno call MOM_error(FATAL,"set_up_ALE_sponge_field: "//mesg) endif @@ -605,8 +605,8 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, type(time_type), intent(in) :: Time !< The current model time type(ocean_grid_type), intent(in) :: G !< Grid structure (in). type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & target, intent(in) :: f_ptr !< Pointer to the field to be damped (in). type(ALE_sponge_CS), pointer :: CS !< Sponge control structure (in/out). @@ -634,7 +634,7 @@ subroutine set_up_ALE_sponge_field_varying(filename, fieldname, Time, G, GV, US, if (CS%fldno > MAX_FIELDS_) then write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease & &the number of fields to be damped in the call to & - &initialize_sponge." )') CS%fldno + &initialize_ALE_sponge." )') CS%fldno call MOM_error(FATAL,"set_up_ALE_sponge_field: "//mesg) endif ! get a unique time interp id for this field. If sponge data is ongrid, then setup @@ -788,11 +788,11 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure (in). type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in) real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]. type(ALE_sponge_CS), pointer :: CS !< A pointer to the control structure for this module - !! that is set by a previous call to initialize_sponge (in). + !! that is set by a previous call to initialize_ALE_sponge (in). type(time_type), optional, intent(in) :: Time !< The current model date real :: damp ! The timestep times the local damping coefficient [nondim]. @@ -833,8 +833,8 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) nz_data = CS%Ref_val(m)%nz_data allocate(sp_val(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data)) - sp_val(:,:,:)=0.0 - mask_z(:,:,:)=0.0 + sp_val(:,:,:) = 0.0 + mask_z(:,:,:) = 0.0 call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, 1.0, G, sp_val, mask_z, z_in, & z_edges_in, missing_value, .true., .false., .false., & spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, & @@ -1003,17 +1003,20 @@ end subroutine apply_ALE_sponge !> Rotate the ALE sponge fields from the input to the model index map. subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) - type(ALE_sponge_CS), intent(in) :: sponge_in - type(ocean_grid_type), intent(in) :: G_in - type(ALE_sponge_CS), pointer :: sponge - type(ocean_grid_type), intent(in) :: G - integer, intent(in) :: turns - type(param_file_type), intent(in) :: param_file + type(ALE_sponge_CS), intent(in) :: sponge_in !< The control structure for this module with the + !! original grid rotation + type(ocean_grid_type), intent(in) :: G_in !< The ocean's grid structure with the original rotation. + type(ALE_sponge_CS), pointer :: sponge !< A pointer to the control that will be set up with + !! the new grid rotation + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure with the new rotation. + integer, intent(in) :: turns !< The number of 90-degree turns between grids + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. ! First part: Index construction ! 1. Reconstruct Iresttime(:,:) from sponge_in ! 2. rotate Iresttime(:,:) - ! 3. Call initialize_sponge using new grid and rotated Iresttime(:,:) + ! 3. Call initialize_ALE_sponge using new grid and rotated Iresttime(:,:) ! All the index adjustment should follow from the Iresttime rotation real, dimension(:,:), allocatable :: Iresttime_in, Iresttime @@ -1040,15 +1043,13 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) endif ! Re-populate the 2D Iresttime and data_h arrays on the original grid - do c = 1, sponge_in%num_col + do c=1,sponge_in%num_col c_i = sponge_in%col_i(c) c_j = sponge_in%col_j(c) Iresttime_in(c_i, c_j) = sponge_in%Iresttime_col(c) - if (fixed_sponge) then - do k = 1, nz_data - data_h(c_i, c_j, k) = sponge_in%Ref_h%p(k,c) - enddo - endif + if (fixed_sponge) then ; do k=1,nz_data + data_h(c_i, c_j, k) = sponge_in%Ref_h%p(k,c) + enddo ; endif enddo call rotate_array(Iresttime_in, turns, Iresttime) @@ -1076,19 +1077,13 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) allocate(sp_val(G%isd:G%ied, G%jsd:G%jed, nz_data)) endif - do n = 1, sponge_in%fldno + do n=1,sponge_in%fldno ! Assume that tracers are pointers and are remapped in other functions(?) sp_ptr => sponge_in%var(n)%p sp_val_in(:,:,:) = 0.0 - do c = 1, sponge_in%num_col - c_i = sponge_in%col_i(c) - c_j = sponge_in%col_j(c) - if (fixed_sponge) then - do k = 1, nz_data - sp_val_in(c_i, c_j, k) = sponge_in%Ref_val(n)%p(k,c) - enddo - endif - enddo + if (fixed_sponge) then ; do c=1,sponge_in%num_col ; do k=1,nz_data + sp_val_in(sponge_in%col_i(c), sponge_in%col_j(c), k) = sponge_in%Ref_val(n)%p(k,c) + enddo ; enddo ; endif call rotate_array(sp_val_in, turns, sp_val) if (fixed_sponge) then @@ -1144,17 +1139,22 @@ end subroutine rotate_ALE_sponge ! TODO: This function solely exists to replace field pointers in the sponge ! after rotation. This function is part of a temporary solution until ! something more robust is developed. -subroutine update_ALE_sponge_field(sponge, p_old, p_new) - type(ALE_sponge_CS), pointer :: sponge - real, dimension(:,:,:), target, intent(in) :: p_old - real, dimension(:,:,:), target, intent(in) :: p_new +subroutine update_ALE_sponge_field(sponge, p_old, G, GV, p_new) + type(ALE_sponge_CS), pointer :: sponge !< A pointer to the control structure for this module + !! that is set by a previous call to initialize_ALE_sponge. + real, dimension(:,:,:), & + target, intent(in) :: p_old !< The previous array of target values + type(ocean_grid_type), intent(in) :: G !< The updated ocean grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + target, intent(in) :: p_new !< The new array of target values integer :: n - do n = 1, sponge%fldno - if (associated(sponge%var(n)%p, p_old)) & - sponge%var(n)%p => p_new + do n=1,sponge%fldno + if (associated(sponge%var(n)%p, p_old)) sponge%var(n)%p => p_new enddo + end subroutine update_ALE_sponge_field @@ -1163,7 +1163,7 @@ end subroutine update_ALE_sponge_field !> This subroutine deallocates any memory associated with the ALE_sponge module. subroutine ALE_sponge_end(CS) type(ALE_sponge_CS), pointer :: CS !< A pointer to the control structure that is - !! set by a previous call to initialize_sponge. + !! set by a previous call to initialize_ALE_sponge. integer :: m From b2453e4fa78b7769cad724d10edab94837ebe8aa Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 11 May 2020 15:55:30 -0400 Subject: [PATCH 23/40] Added array-syntax notation for a full-array copy Added array-syntax notation for a full-array copy in ISOMIP_Tracer.F90. All answers are bitwise identical. --- src/tracer/ISOMIP_tracer.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index c9bf98f7ff..5503287c50 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -287,7 +287,7 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G if (.not.associated(CS)) return - melt(:,:) = fluxes%iceshelf_melt + melt(:,:) = fluxes%iceshelf_melt(:,:) ! max. melt mmax = MAXVAL(melt(is:ie,js:je)) From aeac463d7bb086edf28a4c50e5a73c113be184d5 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 11 May 2020 17:52:16 -0400 Subject: [PATCH 24/40] Initialized dMLD_min and dMLD_max in ePBL_column Initialized dMLD_min and dMLD_max in ePBL_column, and corrected a comment in response to helpful reviews from Brandon Reichl and Andrew Shao. Because these two arrays are not used until after the 3rd iteration, this may not matter to the solution, although it should help make the code clearer and avoids unused variables. All answers in the MOM6-examples test cases are bitwise identical. --- src/parameterizations/vertical/MOM_energetic_PBL.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index e2fda25f6d..25e1f80ff0 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -839,10 +839,10 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs !/The following lines are for the iteration over MLD ! max_MLD will initialized as ocean bottom depth max_MLD = 0.0 ; do k=1,nz ; max_MLD = max_MLD + h(k)*GV%H_to_Z ; enddo - !min_MLD will initialize as 0. + ! min_MLD will be initialized to 0. min_MLD = 0.0 ! Set values of the wrong signs to indicate that these changes are not based on valid estimates - ! dMLD_min = -1.0*US%m_to_Z ; dMLD_max = 1.0*US%m_to_Z + dMLD_min = -1.0*US%m_to_Z ; dMLD_max = 1.0*US%m_to_Z ! If no first guess is provided for MLD, try the middle of the water column if (MLD_guess <= min_MLD) MLD_guess = 0.5 * (min_MLD + max_MLD) @@ -1434,7 +1434,7 @@ subroutine ePBL_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, abs ! For the next pass, guess the average of the minimum and maximum values. MLD_guess = 0.5*(min_MLD + max_MLD) else ! Try using the false position method or the returned value instead of simple bisection. - ! Taking the occasional step with MLD_output empirically step helps to converge faster. + ! Taking the occasional step with MLD_output empirically helps to converge faster. if ((dMLD_min > 0.0) .and. (dMLD_max < 0.0) .and. (OBL_it > 2) .and. (mod(OBL_it-1,4)>0)) then ! Both bounds have valid change estimates and are probably in the range of possible outputs. MLD_Guess = (dMLD_min*max_MLD - dMLD_max*min_MLD) / (dMLD_min - dMLD_max) From 4b33dd03250c1e9e8d1c126c82c36e89b4f9e3d9 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 13 May 2020 18:46:40 -0400 Subject: [PATCH 25/40] (*)Corrected PPM_angular_advect Corrected the negative CFL branch of PPM_angular_advect in MOM_internal_tides. Simultaneously there was some revision to match other equivalent PPM advection schemes in the MOM6 code and to replace some divisions by a multiplication by a reciprocal. The previous version was sufficiently wrong that it could not ever have been used in any scientifically meaningful solutions, including anything in MOM6-examples. Accordingly, the PPM_angular_advect code was changed without a flag to retain the previous answers. All answers in the MOM6-examples test cases are bitwise identical, and output files are unchanged. --- .../lateral/MOM_internal_tides.F90 | 49 +++++++++---------- 1 file changed, 23 insertions(+), 26 deletions(-) diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 6145fb1dce..a0f1631d6d 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -894,7 +894,8 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) real :: dMx, dMn real :: Ep, Ec, Em ! Mean angular energy density for three successive wedges in angular ! orientation [R Z3 T-2 rad-1 ~> J m-2 rad-1] - real :: dA, mA, a6 ! Difference, mean, and curvature of energy density [R Z3 T-2 rad-1 ~> J m-2 rad-1] + real :: dA, curv_3 ! Difference and curvature of energy density [R Z3 T-2 rad-1 ~> J m-2 rad-1] + real, parameter :: oneSixth = 1.0/6.0 ! One sixth [nondim] integer :: a I_dt = 1 / dt @@ -912,23 +913,21 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) Ec = En2d(a) *I_Angle_size Em = En2d(a-1)*I_Angle_size ! Calculate and bound edge values of energy density. - aL = ( 5.*Ec + ( 2.*Em - Ep ) )/6. ! H3 estimate + aL = ( 5.*Ec + ( 2.*Em - Ep ) ) * oneSixth ! H3 estimate aL = max( min(Ec,Em), aL) ; aL = min( max(Ec,Em), aL) ! Bound - aR = ( 5.*Ec + ( 2.*Ep - Em ) )/6. ! H3 estimate + aR = ( 5.*Ec + ( 2.*Ep - Em ) ) * oneSixth ! H3 estimate aR = max( min(Ec,Ep), aR) ; aR = min( max(Ec,Ep), aR) ! Bound - dA = aR - aL ; mA = 0.5*( aR + aL ) + dA = aR - aL if ((Ep-Ec)*(Ec-Em) <= 0.) then - aL = Ec ; aR = Ec ! PCM for local extremum - elseif ( dA*(Ec-mA) > (dA*dA)/6. ) then - aL = 3.*Ec - 2.*aR !? - elseif ( dA*(Ec-mA) < - (dA*dA)/6. ) then - aR = 3.*Ec - 2.*aL !? + aL = Ec ; aR = Ec ! use PCM for local extremum + elseif ( 3.0*dA*(2.*Ec - (aR + aL)) > (dA*dA) ) then + aL = 3.*Ec - 2.*aR ! Flatten the profile to move the extremum to the left edge + elseif ( 3.0*dA*(2.*Ec - (aR + aL)) < - (dA*dA) ) then + aR = 3.*Ec - 2.*aL ! Flatten the profile to move the extremum to the right edge endif - a6 = 6.*Ec - 3. * (aR + aL) ! Curvature + curv_3 = (aR + aL) - 2.0*Ec ! Curvature ! Calculate angular flux rate [R Z3 T-3 ~> W m-2] - flux = u_ang*( aR + 0.5 * CFL_ang(A) * ( ( aL - aR ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - ! The following expression copied from tracer_advect is equivalent. - ! flux = u_ang*( aR - 0.5 * CFL_ang(A) * ( ( aR - aL ) - a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) + flux = u_ang*( aR + CFL_ang(A) * ( 0.5*(aL - aR) + curv_3 * (CFL_ang(A) - 1.5) ) ) ! Calculate amount of energy fluxed between wedges [R Z3 T-2 ~> J m-2] Flux_En(A) = dt * flux !Flux_En(A) = (dt * I_Angle_size) * flux @@ -940,24 +939,22 @@ subroutine PPM_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang) Ec = En2d(a+1)*I_Angle_size Em = En2d(a) *I_Angle_size ! Calculate and bound edge values of energy density. - aL = ( 5.*Ec + ( 2.*Em - Ep ) )/6. ! H3 estimate + aL = ( 5.*Ec + ( 2.*Em - Ep ) ) * oneSixth ! H3 estimate aL = max( min(Ec,Em), aL) ; aL = min( max(Ec,Em), aL) ! Bound - aR = ( 5.*Ec + ( 2.*Ep - Em ) )/6. ! H3 estimate + aR = ( 5.*Ec + ( 2.*Ep - Em ) ) * oneSixth ! H3 estimate aR = max( min(Ec,Ep), aR) ; aR = min( max(Ec,Ep), aR) ! Bound - dA = aR - aL ; mA = 0.5*( aR + aL ) + dA = aR - aL if ((Ep-Ec)*(Ec-Em) <= 0.) then - aL = Ec ; aR = Ec ! PCM for local extremum - elseif ( dA*(Ec-mA) > (dA*dA)/6. ) then - aL = 3.*Ec - 2.*aR - elseif ( dA*(Ec-mA) < - (dA*dA)/6. ) then - aR = 3.*Ec - 2.*aL + aL = Ec ; aR = Ec ! use PCM for local extremum + elseif ( 3.0*dA*(2.*Ec - (aR + aL)) > (dA*dA) ) then + aL = 3.*Ec - 2.*aR ! Flatten the profile to move the extremum to the left edge + elseif ( 3.0*dA*(2.*Ec - (aR + aL)) < - (dA*dA) ) then + aR = 3.*Ec - 2.*aL ! Flatten the profile to move the extremum to the right edge endif - a6 = 6.*Ec - 3. * (aR + aL) ! Curvature + curv_3 = (aR + aL) - 2.0*Ec ! Curvature ! Calculate angular flux rate [R Z3 T-3 ~> W m-2] - !### This expression is wrong, because it was just copied from above. The correct one is below - flux = u_ang*( aR + 0.5 * CFL_ang(A) * ( ( aL - aR ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) ) - ! This is the correct expression; note that CFL_ang is negative here, so it looks a bit odd. - !flux = u_ang*( aL - 0.5 * CFL_ang(A) * ( ( aR - aL ) + a6 * ( 1. + 2./3. * CFL_ang(A) ) ) ) + ! Note that CFL_ang is negative here, so it looks odd compared with equivalent expressions. + flux = u_ang*( aL - CFL_ang(A) * ( 0.5*(aR - aL) + curv_3 * (-CFL_ang(A) - 1.5) ) ) ! Calculate amount of energy fluxed between wedges [R Z3 T-2 ~> J m-2] Flux_En(A) = dt * flux !Flux_En(A) = (dt * I_Angle_size) * flux From 435fa9b3a3bd4389a521f4dddfa9c6244ba92763 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 14 May 2020 15:52:45 -0400 Subject: [PATCH 26/40] Makefile: reduced logging, -k support The Makefile has been modified to reduce the amount of output during testing. Output is generally omitted on a successful test. When a test fails, we only display a small portion of the total output. We also now run all tests, even if they fail, in order to give a complete profile of the test failures. Regression testing rules have been integrated into the general rules. Finally, Travis config has been modified to further reduce output, and to run all of tests (make -k). --- .testing/Makefile | 212 +++++++++++++++++++++++++++++----------------- .travis.yml | 5 +- 2 files changed, 137 insertions(+), 80 deletions(-) diff --git a/.testing/Makefile b/.testing/Makefile index 99672268c3..f0ecb09a0a 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -119,28 +119,28 @@ build/%/MOM6: build/%/Makefile $(FMS)/lib/libfms.a build/%/Makefile: build/%/path_names cp $(MKMF_TEMPLATE) $(@D) cd $(@D) && $(MKMF) \ - -t $(notdir $(MKMF_TEMPLATE)) \ - -o '-I ../../$(DEPS)/fms/build' \ - -p MOM6 \ - -l '../../$(DEPS)/fms/lib/libfms.a' \ - -c $(MKMF_CPP) \ - path_names + -t $(notdir $(MKMF_TEMPLATE)) \ + -o '-I ../../$(DEPS)/fms/build' \ + -p MOM6 \ + -l '../../$(DEPS)/fms/lib/libfms.a' \ + -c $(MKMF_CPP) \ + path_names # NOTE: These path_names rules could be merged build/target/path_names: $(LIST_PATHS) $(TARGET_CODEBASE) $(TARGET_SOURCE) mkdir -p $(@D) cd $(@D) && $(LIST_PATHS) -l \ - ../../$(TARGET_CODEBASE)/src \ - ../../$(TARGET_CODEBASE)/config_src/solo_driver \ - ../../$(TARGET_CODEBASE)/$(GRID_SRC) + ../../$(TARGET_CODEBASE)/src \ + ../../$(TARGET_CODEBASE)/config_src/solo_driver \ + ../../$(TARGET_CODEBASE)/$(GRID_SRC) build/%/path_names: $(LIST_PATHS) $(MOM_SOURCE) mkdir -p $(@D) cd $(@D) && $(LIST_PATHS) -l \ - ../../../src \ - ../../../config_src/solo_driver \ - ../../../$(GRID_SRC) + ../../../src \ + ../../../config_src/solo_driver \ + ../../../$(GRID_SRC) # Target repository for regression tests $(TARGET_CODEBASE): @@ -158,10 +158,10 @@ $(FMS)/lib/libfms.a: $(FMS)/build/Makefile $(FMS)/build/Makefile: $(FMS)/build/path_names cp $(MKMF_TEMPLATE) $(@D) cd $(@D) && $(MKMF) \ - -t $(notdir $(MKMF_TEMPLATE)) \ - -p ../lib/libfms.a \ - -c $(MKMF_CPP) \ - path_names + -t $(notdir $(MKMF_TEMPLATE)) \ + -p ../lib/libfms.a \ + -c $(MKMF_CPP) \ + path_names $(FMS)/build/path_names: $(LIST_PATHS) $(FMS)/src $(FMS_SOURCE) mkdir -p $(@D) @@ -202,18 +202,38 @@ test.repros: $(foreach c,$(CONFIGS),$(c).repro $(c).repro.diag) test.openmps: $(foreach c,$(CONFIGS),$(c).openmp $(c).openmp.diag) test.nans: $(foreach c,$(CONFIGS),$(c).nan $(c).nan.diag) test.dims: $(foreach c,$(CONFIGS),$(foreach d,$(DIMS),$(c).dim.$(d) $(c).dim.$(d).diag)) - test.regressions: $(foreach c,$(CONFIGS),$(c).regression $(c).regression.diag) - ! ls -1 results/*/*.reg -define CMP_RULE -.PRECIOUS: $(foreach b,$(2),results/%/ocean.stats.$(b)) -%.$(1): $(foreach b,$(2),results/%/ocean.stats.$(b)) - cmp $$^ || diff $$^ +# Color highlights for test results +RED=\033[0;31m +GREEN=\033[0;32m +RESET=\033[0m + +DONE=${GREEN}DONE${RESET} +PASS=${GREEN}PASS${RESET} +FAIL=${RED}FAIL${RESET} -.PRECIOUS: $(foreach b,$(2),results/%/chksum_diag.$(b)) -%.$(1).diag: $(foreach b,$(2),results/%/chksum_diag.$(b)) - cmp $$^ || diff $$^ +# Comparison rules +# $(1): Test type (grid, layout, &c.) +# $(2): Comparison targets (symmetric asymmetric, symmetric layout, &c.) +define CMP_RULE +.PRECIOUS: $(foreach b,$(2),work/%/$(b)/ocean.stats) +%.$(1): $(foreach b,$(2),work/%/$(b)/ocean.stats) + @cmp $$^ || !( \ + mkdir -p results/$$*; \ + (diff $$^ | tee results/$$*/ocean.stats.$(1).diff | head) ; \ + echo -e "${FAIL}: Solutions $$*.$(1) have changed." \ + ) + @echo -e "${PASS}: Solutions $$*.$(1) agree." + +.PRECIOUS: $(foreach b,$(2),work/%/$(b)/chksum_diag) +%.$(1).diag: $(foreach b,$(2),work/%/$(b)/chksum_diag) + @cmp $$^ || !( \ + mkdir -p results/$$*; \ + (diff $$^ | tee results/$$*/chksum_diag.$(1).diff | head) ; \ + echo -e "${FAIL}: Diagnostics $$*.$(1).diag have changed." \ + ) + @echo -e "${PASS}: Diagnostics $$*.$(1).diag agree." endef $(eval $(call CMP_RULE,grid,symmetric asymmetric)) @@ -223,29 +243,31 @@ $(eval $(call CMP_RULE,repro,symmetric repro)) $(eval $(call CMP_RULE,openmp,symmetric openmp)) $(eval $(call CMP_RULE,nan,symmetric nan)) $(foreach d,$(DIMS),$(eval $(call CMP_RULE,dim.$(d),symmetric dim.$(d)))) +$(eval $(call CMP_RULE,regression,symmetric target)) # Custom comparison rules -.PRECIOUS: $(foreach b,symmetric restart target,results/%/ocean.stats.$(b)) - # Restart tests only compare the final stat record -%.restart: $(foreach b,symmetric restart,results/%/ocean.stats.$(b)) - cmp $(foreach f,$^,<(tr -s ' ' < $(f) | cut -d ' ' -f3- | tail -n 1)) \ - || diff $^ +.PRECIOUS: $(foreach b,symmetric restart target,work/%/$(b)/ocean.stats) +%.restart: $(foreach b,symmetric restart,work/%/$(b)/ocean.stats) + #cmp $(foreach f,$^,<(tr -s ' ' < $(f) | cut -d ' ' -f3- | tail -n 1)) \ + # || diff $^ + @cmp $(foreach f,$^,<(tr -s ' ' < $(f) | cut -d ' ' -f3- | tail -n 1)) \ + || !( \ + mkdir -p results/$*; \ + (diff $$^ | tee results/$*/chksum_diag.restart.diff | head) ; \ + echo -e "${FAIL}: Diagnostics $*.restart.diag have changed." \ + ) + @echo -e "${PASS}: Diagnostics $*.restart.diag agree." # TODO: chksum_diag parsing of restart files -# All regression tests must be completed when considering answer changes -%.regression: $(foreach b,symmetric target,results/%/ocean.stats.$(b)) - cmp $^ || (diff $^ > $<.reg || true) - -%.regression.diag: $(foreach b,symmetric target,results/%/chksum_diag.$(b)) - cmp $^ || (diff $^ > $<.reg || true) #--- # Test run output files # Generalized MPI environment variable support +# XXX: Using `-env` in the MPICH test can erroneously producing an `nv` file. # $(1): Environment variables ifeq ($(shell $(MPIRUN) -x tmp=1 true 2> /dev/null ; echo $$?), 0) MPIRUN_CMD=$(MPIRUN) $(if $(1),-x $(1),) @@ -255,7 +277,8 @@ else MPIRUN_CMD=$(1) $(MPIRUN) endif -# Rule to build results//{ocean.stats,chksum_diag}. + +# Rule to build work//{ocean.stats,chksum_diag}. # $(1): Test configuration name # $(2): Executable type # $(3): Enable coverage flag @@ -263,25 +286,28 @@ endif # $(5): Environment variables # $(6): Number of MPI ranks define STAT_RULE -results/%/ocean.stats.$(1): build/$(2)/MOM6 +work/%/$(1)/ocean.stats work/%/$(1)/chksum_diag: build/$(2)/MOM6 + @echo "Running test $$*.$(1)..." if [ $(3) ]; then find build/$(2) -name *.gcda -exec rm -f '{}' \; ; fi - mkdir -p work/$$*/$(1) - cp -rL $$*/* work/$$*/$(1) - cd work/$$*/$(1) && if [ -f Makefile ]; then make; fi - mkdir -p work/$$*/$(1)/RESTART - echo -e "$(4)" > work/$$*/$(1)/MOM_override - cd work/$$*/$(1) && $$(call MPIRUN_CMD,$(5)) -n $(6) ../../../$$< 2> debug.out > std.out \ - || ! sed 's/^/$$*.$(1): /' std.out debug.out \ - && sed 's/^/$$*.$(1): /' std.out mkdir -p $$(@D) - cp work/$$*/$(1)/ocean.stats $$@ + cp -rL $$*/* $$(@D) + cd $$(@D) && if [ -f Makefile ]; then make; fi + mkdir -p $$(@D)/RESTART + echo -e "$(4)" > $$(@D)/MOM_override + cd $$(@D) \ + && $$(call MPIRUN_CMD,$(5)) -n $(6) ../../../$$< 2> std.err > std.out \ + || !( \ + mkdir -p ../../../results/$$*/ ; \ + cat std.out | tee ../../../results/$$*/std.$(1).out | tail ; \ + cat std.err | tee ../../../results/$$*/std.$(1).err | tail ; \ + rm ocean.stats chksum_diag ; \ + echo -e "${FAIL}: $$*.$(1) failed at runtime." \ + ) + @echo -e "${DONE}: $$*.$(1); no runtime errors." if [ $(3) ]; then cd .. && bash <(curl -s https://codecov.io/bash) -n $$@; fi - -results/%/chksum_diag.$(1): results/%/ocean.stats.$(1) - mkdir -p $$(@D) - cp work/$$*/$(1)/chksum_diag $$@ endef + # Define $(,) as comma escape character , := , @@ -300,50 +326,80 @@ $(eval $(call STAT_RULE,dim.z,symmetric,,Z_RESCALE_POWER=11,,1)) $(eval $(call STAT_RULE,dim.q,symmetric,,Q_RESCALE_POWER=11,,1)) $(eval $(call STAT_RULE,dim.r,symmetric,,R_RESCALE_POWER=11,,1)) + # Restart tests require significant preprocessing, and are handled separately. -results/%/ocean.stats.restart: build/symmetric/MOM6 - rm -rf work/$*/restart - mkdir -p work/$*/restart - cp -rL $*/* work/$*/restart +work/%/restart/ocean.stats: build/symmetric/MOM6 + rm -rf $(@D) + mkdir -p $(@D) + cp -rL $*/* $(@D) cd work/$*/restart && if [ -f Makefile ]; then make; fi - mkdir -p work/$*/restart/RESTART + mkdir -p $(@D)/RESTART # Generate the half-period input namelist - # TODO: Assumes runtime set by DAYMAX, will fail if set by input.nml - cd work/$*/restart \ - && daymax=$$(grep DAYMAX MOM_input | cut -d '!' -f 1 | cut -d '=' -f 2 | xargs) \ - && timeunit=$$(grep TIMEUNIT MOM_input | cut -d '!' -f 1 | cut -d '=' -f 2 | xargs) \ - && if [ -z "$${timeunit}" ]; then timeunit="8.64e4"; fi \ - && printf -v timeunit_int "%.f" "$${timeunit}" \ - && halfperiod=$$(printf "%.f" $$(bc <<< "scale=10; 0.5 * $${daymax} * $${timeunit_int}")) \ - && printf "\n&ocean_solo_nml\n seconds = $${halfperiod}\n/\n" >> input.nml + # TODO: Assumes that runtime set by DAYMAX, will fail if set by input.nml + cd $(@D) \ + && daymax=$$(grep DAYMAX MOM_input | cut -d '!' -f 1 | cut -d '=' -f 2 | xargs) \ + && timeunit=$$(grep TIMEUNIT MOM_input | cut -d '!' -f 1 | cut -d '=' -f 2 | xargs) \ + && if [ -z "$${timeunit}" ]; then timeunit="8.64e4"; fi \ + && printf -v timeunit_int "%.f" "$${timeunit}" \ + && halfperiod=$$(printf "%.f" $$(bc <<< "scale=10; 0.5 * $${daymax} * $${timeunit_int}")) \ + && printf "\n&ocean_solo_nml\n seconds = $${halfperiod}\n/\n" >> input.nml # Run the first half-period - cd work/$*/restart && $(MPIRUN) -n 1 ../../../$< 2> debug1.out > std1.out \ - || ! sed 's/^/$*.restart1: /' std1.out debug1.out \ - && sed 's/^/$*.restart1: /' std1.out + cd $(@D) && $(MPIRUN) -n 1 ../../../$< 2> std1.err > std1.out \ + || !( \ + cat std1.out | tee ../../../results/$*/std.restart1.out | tail ; \ + cat std1.err | tee ../../../results/$*/std.restart1.err | tail ; \ + echo -e "${FAIL}: $*.restart failed at runtime." \ + ) # Setup the next inputs - cd work/$*/restart && rm -rf INPUT && mv RESTART INPUT - mkdir work/$*/restart/RESTART - cd work/$*/restart && sed -i -e "s/input_filename *= *'n'/input_filename = 'r'/g" input.nml + cd $(@D) && rm -rf INPUT && mv RESTART INPUT + mkdir $(@D)/RESTART + cd $(@D) && sed -i -e "s/input_filename *= *'n'/input_filename = 'r'/g" input.nml # Run the second half-period - cd work/$*/restart && $(MPIRUN) -n 1 ../../../$< 2> debug2.out > std2.out \ - || ! sed 's/^/$*.restart2: /' std2.out debug2.out \ - && sed 's/^/$*.restart2: /' std2.out - # Archive the results and cleanup - mkdir -p $(@D) - cp work/$*/restart/ocean.stats $@ + cd $(@D) && $(MPIRUN) -n 1 ../../../$< 2> std2.err > std2.out \ + || !( \ + cat std2.out | tee ../../../results/$*/std.restart2.out | tail ; \ + cat std2.err | tee ../../../results/$*/std.restart2.err | tail ; \ + echo -e "${FAIL}: $*.restart failed at runtime." \ + ) # TODO: Restart checksum diagnostics +#--- +# Not a true rule; only call this after `make test` to summarize test results. +.PHONY: test.summary +test.summary: + @if ls results/*/* &> /dev/null; then \ + if ls results/*/std.*.err &> /dev/null; then \ + echo "The following tests failed to complete:" ; \ + ls results/*/std.*.out \ + | awk '{split($$0,a,"/"); split(a[3],t,"."); print " ",a[2],":",t[2]}' ; \ + fi; \ + if ls results/*/ocean.stats.*.diff &> /dev/null; then \ + echo "The following tests report solution regressions:" ; \ + ls results/*/ocean.stats.*.diff \ + | awk '{split($$0,a,"/"); split(a[3],t,"."); print " ",a[2],":",t[3]}' ; \ + fi; \ + if ls results/*/chksum_diag.*.diff &> /dev/null; then \ + echo "The following tests report diagnostic regressions:" ; \ + ls results/*/chksum_diag.*.diff \ + | awk '{split($$0,a,"/"); split(a[3],t,"."); print " ",a[2],":",t[2]}' ; \ + fi; \ + false ; \ + else \ + echo -e "${PASS}: All tests passed!"; \ + fi + + #---- +# NOTE: These tests assert that we are in the .testing directory. + .PHONY: clean clean: clean.stats - @# Assert that we are in .testing for recursive delete @[ $$(basename $$(pwd)) = .testing ] rm -rf build .PHONY: clean.stats clean.stats: - @# Assert that we are in .testing for recursive delete @[ $$(basename $$(pwd)) = .testing ] rm -rf work results diff --git a/.travis.yml b/.travis.yml index ac7cab1b14..34bbe0dcce 100644 --- a/.travis.yml +++ b/.travis.yml @@ -40,7 +40,8 @@ jobs: - make all - echo -en 'travis_fold:end:script.1\\r' - echo 'Running tests...' && echo -en 'travis_fold:start:script.2\\r' - - make test + - make -k -s test + - make test.summary - echo -en 'travis_fold:end:script.2\\r' # NOTE: Code coverage upload is here to reduce load imbalance @@ -58,5 +59,5 @@ jobs: - make build.regressions - echo -en 'travis_fold:end:script.1\\r' - echo 'Running tests...' && echo -en 'travis_fold:start:script.2\\r' - - make test.regressions + - make -k -s test.regressions - echo -en 'travis_fold:end:script.2\\r' From 663799dd98de9c74f8055a706df908227790fa72 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 14 May 2020 16:45:17 -0400 Subject: [PATCH 27/40] Travis: remove test folding Test logging is now much shorter, so folding is less important. --- .travis.yml | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/.travis.yml b/.travis.yml index 34bbe0dcce..6b0b4c2a5e 100644 --- a/.travis.yml +++ b/.travis.yml @@ -39,10 +39,8 @@ jobs: - echo 'Build executables...' && echo -en 'travis_fold:start:script.1\\r' - make all - echo -en 'travis_fold:end:script.1\\r' - - echo 'Running tests...' && echo -en 'travis_fold:start:script.2\\r' - make -k -s test - make test.summary - - echo -en 'travis_fold:end:script.2\\r' # NOTE: Code coverage upload is here to reduce load imbalance - if: type = pull_request @@ -58,6 +56,5 @@ jobs: - echo 'Build executables...' && echo -en 'travis_fold:start:script.1\\r' - make build.regressions - echo -en 'travis_fold:end:script.1\\r' - - echo 'Running tests...' && echo -en 'travis_fold:start:script.2\\r' - make -k -s test.regressions - - echo -en 'travis_fold:end:script.2\\r' + - make test.summary From 3567d8397170cef7d005811a7a95a51efc6b0eda Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 14 May 2020 17:34:10 -0400 Subject: [PATCH 28/40] Makefile: codecov reporting CodeCov reporting log is now saved to results/ rather than piped to stdout, further reducing test logging output. --- .testing/Makefile | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/.testing/Makefile b/.testing/Makefile index f0ecb09a0a..bcfea91e40 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -304,7 +304,12 @@ work/%/$(1)/ocean.stats work/%/$(1)/chksum_diag: build/$(2)/MOM6 echo -e "${FAIL}: $$*.$(1) failed at runtime." \ ) @echo -e "${DONE}: $$*.$(1); no runtime errors." - if [ $(3) ]; then cd .. && bash <(curl -s https://codecov.io/bash) -n $$@; fi + if [ $(3) ]; then \ + mkdir -p results/$$* ; \ + bash <(curl -s https://codecov.io/bash) -n $$@ \ + > results/$$*/codecov.$(1).out \ + 2> results/$$*/codecov.$(1).err ; \ + fi endef From c233e02fe5f6720ddd29724f77e3be27561dd203 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 14 May 2020 20:40:12 -0400 Subject: [PATCH 29/40] Makefile: copy codecov output to work, not results. We were saving codecov output to results, but this was breaking the `test.summary` test, which only returns true if `results` is empty. This change should resolve this issue. --- .testing/Makefile | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.testing/Makefile b/.testing/Makefile index bcfea91e40..66a116a32a 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -307,8 +307,8 @@ work/%/$(1)/ocean.stats work/%/$(1)/chksum_diag: build/$(2)/MOM6 if [ $(3) ]; then \ mkdir -p results/$$* ; \ bash <(curl -s https://codecov.io/bash) -n $$@ \ - > results/$$*/codecov.$(1).out \ - 2> results/$$*/codecov.$(1).err ; \ + > work/$$*/codecov.$(1).out \ + 2> work/$$*/codecov.$(1).err ; \ fi endef From 42a9eaffa6c0c11182f25495176da83e938d97cc Mon Sep 17 00:00:00 2001 From: Hemant Khatri Date: Sat, 16 May 2020 22:19:36 -0400 Subject: [PATCH 30/40] Changed the array intilization to standard halos instead of wide halos for variables used in computing barotropic velocity tendency terms --- src/core/MOM_barotropic.F90 | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index c2a9b3c8b1..ee90dc1d67 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -477,10 +477,14 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! sums less than one due to viscous losses. Nondimensional. real, dimension(SZIB_(G),SZJ_(G)) :: & av_rem_u, & ! The weighted average of visc_rem_u, nondimensional. - tmp_u ! A temporary array at u points. + tmp_u, & ! A temporary array at u points. + ubt_st, & ! The zonal barotropic velocity at the start of timestep [L T-1 ~> m s-1]. + ubt_dt ! The zonal barotropic velocity tendency [L T-2 ~> m s-2]. real, dimension(SZI_(G),SZJB_(G)) :: & av_rem_v, & ! The weighted average of visc_rem_v, nondimensional. - tmp_v ! A temporary array at v points. + tmp_v, & ! A temporary array at v points. + vbt_st, & ! The meridional barotropic velocity at the start of timestep [L T-1 ~> m s-1]. + vbt_dt ! The meridional barotropic velocity tendency [L T-2 ~> m s-2]. real, dimension(SZI_(G),SZJ_(G)) :: & e_anom ! The anomaly in the sea surface height or column mass ! averaged between the beginning and end of the time step, @@ -491,8 +495,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! or [T-1 H-1 ~> s-1 m-1 or m2 s-1 kg-1] real, dimension(SZIBW_(CS),SZJW_(CS)) :: & ubt, & ! The zonal barotropic velocity [L T-1 ~> m s-1]. - ubt_st, & ! The zonal barotropic velocity at the start of timestep [L T-1 ~> m s-1]. - ubt_dt, & ! The zonal barotropic velocity tendency [L T-2 ~> m s-2]. + ! ubt_st, & ! The zonal barotropic velocity at the start of timestep [L T-1 ~> m s-1]. + ! ubt_dt, & ! The zonal barotropic velocity tendency [L T-2 ~> m s-2]. bt_rem_u, & ! The fraction of the barotropic zonal velocity that remains ! after a time step, the remainder being lost to bottom drag. ! bt_rem_u is a nondimensional number between 0 and 1. @@ -526,8 +530,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! spacing [H L ~> m2 or kg m-1]. real, dimension(SZIW_(CS),SZJBW_(CS)) :: & vbt, & ! The meridional barotropic velocity [L T-1 ~> m s-1]. - vbt_st, & ! The meridional barotropic velocity at the start of timestep [L T-1 ~> m s-1]. - vbt_dt, & ! The meridional barotropic velocity tendency [L T-2 ~> m s-2]. +! vbt_st, & ! The meridional barotropic velocity at the start of timestep [L T-1 ~> m s-1]. +! vbt_dt, & ! The meridional barotropic velocity tendency [L T-2 ~> m s-2]. bt_rem_v, & ! The fraction of the barotropic meridional velocity that ! remains after a time step, the rest being lost to bottom ! drag. bt_rem_v is a nondimensional number between 0 and 1. @@ -1553,11 +1557,15 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif if (CS%id_ubtdt > 0) then - ubt_st(:,:) = ubt(:,:) + do j=js-1,je+1 ; do I=is-1,ie + ubt_st(I,j) = ubt(I,j) + enddo ; enddo endif if (CS%id_vbtdt > 0) then - vbt_st(:,:) = vbt(:,:) + do J=js-1,je ; do i=is-1,ie+1 + vbt_st(i,J) = vbt(i,J) + enddo ; enddo endif if (query_averaging_enabled(CS%diag)) then From e261b421f9d3c435923219b477633dff21157107 Mon Sep 17 00:00:00 2001 From: Hemant Khatri Date: Sat, 16 May 2020 22:25:23 -0400 Subject: [PATCH 31/40] Removed unused variables --- src/core/MOM_barotropic.F90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index ee90dc1d67..7bb03dd2cf 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -495,8 +495,6 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! or [T-1 H-1 ~> s-1 m-1 or m2 s-1 kg-1] real, dimension(SZIBW_(CS),SZJW_(CS)) :: & ubt, & ! The zonal barotropic velocity [L T-1 ~> m s-1]. - ! ubt_st, & ! The zonal barotropic velocity at the start of timestep [L T-1 ~> m s-1]. - ! ubt_dt, & ! The zonal barotropic velocity tendency [L T-2 ~> m s-2]. bt_rem_u, & ! The fraction of the barotropic zonal velocity that remains ! after a time step, the remainder being lost to bottom drag. ! bt_rem_u is a nondimensional number between 0 and 1. @@ -530,8 +528,6 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ! spacing [H L ~> m2 or kg m-1]. real, dimension(SZIW_(CS),SZJBW_(CS)) :: & vbt, & ! The meridional barotropic velocity [L T-1 ~> m s-1]. -! vbt_st, & ! The meridional barotropic velocity at the start of timestep [L T-1 ~> m s-1]. -! vbt_dt, & ! The meridional barotropic velocity tendency [L T-2 ~> m s-2]. bt_rem_v, & ! The fraction of the barotropic meridional velocity that ! remains after a time step, the rest being lost to bottom ! drag. bt_rem_v is a nondimensional number between 0 and 1. From 6ef10ffc8fb4425b48125a457bce22c4e0af52d5 Mon Sep 17 00:00:00 2001 From: Hemant Khatri Date: Mon, 18 May 2020 09:54:54 -0400 Subject: [PATCH 32/40] Removed white space --- src/core/MOM_barotropic.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 7bb03dd2cf..f63ab42669 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -1557,7 +1557,6 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, ubt_st(I,j) = ubt(I,j) enddo ; enddo endif - if (CS%id_vbtdt > 0) then do J=js-1,je ; do i=is-1,ie+1 vbt_st(i,J) = vbt(i,J) From 2befa56093270b5eb318d6089812904b394f63bf Mon Sep 17 00:00:00 2001 From: Hemant Khatri Date: Mon, 18 May 2020 15:01:27 -0400 Subject: [PATCH 33/40] Trailing white spaces removed --- src/core/MOM_barotropic.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index f63ab42669..e6f11987a2 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -1553,8 +1553,8 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, endif if (CS%id_ubtdt > 0) then - do j=js-1,je+1 ; do I=is-1,ie - ubt_st(I,j) = ubt(I,j) + do j=js-1,je+1 ; do I=is-1,ie + ubt_st(I,j) = ubt(I,j) enddo ; enddo endif if (CS%id_vbtdt > 0) then From c907001567da3376db30a25365d8974cebb679a7 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 19 May 2020 17:15:25 -0400 Subject: [PATCH 34/40] Rotate: fixed-time ALE sponge bugfixes The sp_val and sp_val_in fields were only allocated for fixed-time sponges (fixed_sponge=.true.) but were being referenced regardless of this flag in a few cases. On x86 architectures this did not create any errors, but was detected on an Arm64 CPU. This patch correctly moves all code associated with these variables into a flag-enabled block. --- .../vertical/MOM_ALE_sponge.F90 | 32 +++++++++++-------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index fe1ccab53d..70914a69ad 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -1047,9 +1047,11 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) c_i = sponge_in%col_i(c) c_j = sponge_in%col_j(c) Iresttime_in(c_i, c_j) = sponge_in%Iresttime_col(c) - if (fixed_sponge) then ; do k=1,nz_data - data_h(c_i, c_j, k) = sponge_in%Ref_h%p(k,c) - enddo ; endif + if (fixed_sponge) then + do k=1,nz_data + data_h(c_i, c_j, k) = sponge_in%Ref_h%p(k,c) + enddo + endif enddo call rotate_array(Iresttime_in, turns, Iresttime) @@ -1080,15 +1082,22 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) do n=1,sponge_in%fldno ! Assume that tracers are pointers and are remapped in other functions(?) sp_ptr => sponge_in%var(n)%p - sp_val_in(:,:,:) = 0.0 - if (fixed_sponge) then ; do c=1,sponge_in%num_col ; do k=1,nz_data - sp_val_in(sponge_in%col_i(c), sponge_in%col_j(c), k) = sponge_in%Ref_val(n)%p(k,c) - enddo ; enddo ; endif - - call rotate_array(sp_val_in, turns, sp_val) if (fixed_sponge) then + sp_val_in(:,:,:) = 0.0 + do c=1,sponge_in%num_col + c_i = sponge_in%col_i(c) + c_j = sponge_in%col_j(c) + do k=1,nz_data + sp_val_in(c_i, c_j, k) = sponge_in%Ref_val(n)%p(k,c) + enddo + enddo + + call rotate_array(sp_val_in, turns, sp_val) + ! NOTE: This points sp_val with the unrotated field. See note below. call set_up_ALE_sponge_field(sp_val, G, sp_ptr, sponge) + + deallocate(sp_val_in) else ! We don't want to repeat FMS init in set_up_ALE_sponge_field_varying() ! (time_interp_external_init, init_external_field, etc), so we manually @@ -1118,11 +1127,6 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) endif enddo - if (fixed_sponge) then - deallocate(sp_val_in) - deallocate(sp_val) - endif - ! TODO: var_u and var_v sponge dampling is not yet supported. if (associated(sponge_in%var_u%p) .or. associated(sponge_in%var_v%p)) & call MOM_error(FATAL, "Rotation of ALE sponge velocities is not yet " & From 5c0dc48a04d1e96f5eb0c82797f59f929789ac95 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 19 May 2020 17:20:40 -0400 Subject: [PATCH 35/40] Rotate: ALE sponge whitespace fix --- src/parameterizations/vertical/MOM_ALE_sponge.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 70914a69ad..438236fc2a 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -1048,7 +1048,7 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) c_j = sponge_in%col_j(c) Iresttime_in(c_i, c_j) = sponge_in%Iresttime_col(c) if (fixed_sponge) then - do k=1,nz_data + do k = 1, nz_data data_h(c_i, c_j, k) = sponge_in%Ref_h%p(k,c) enddo endif @@ -1084,10 +1084,10 @@ subroutine rotate_ALE_sponge(sponge_in, G_in, sponge, G, turns, param_file) sp_ptr => sponge_in%var(n)%p if (fixed_sponge) then sp_val_in(:,:,:) = 0.0 - do c=1,sponge_in%num_col + do c = 1, sponge_in%num_col c_i = sponge_in%col_i(c) c_j = sponge_in%col_j(c) - do k=1,nz_data + do k = 1, nz_data sp_val_in(c_i, c_j, k) = sponge_in%Ref_val(n)%p(k,c) enddo enddo From fd4a47d7b70ed31a6a2486d271758a163a123d5b Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 20 May 2020 11:01:32 -0400 Subject: [PATCH 36/40] Travis: Arm64 tests; git depth restored; Makefile There are four minor changes in this patch: - Arm64 testing has been enabled on Travis. As a non-x86 architecture, this will provide more robust testing of the code. - Travis MPI library switch to MPICH There were issues with using OpenMPI on Arm64, particularly with MPI_Init. Switching to MPICH appeared to resolve those issues, and MPICH has a much simpler implementation, so we've switched all jobs to it. - default git depth has been restored Our regression testing method no longer requires a complete history for comparison. - Makefile fix: dimension now reported in summary The test.summary output now includes the dimension type in the list of failed tests. --- .testing/Makefile | 6 +++--- .travis.yml | 23 +++++++++++++++++------ 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/.testing/Makefile b/.testing/Makefile index 66a116a32a..d38189667b 100644 --- a/.testing/Makefile +++ b/.testing/Makefile @@ -378,17 +378,17 @@ test.summary: if ls results/*/std.*.err &> /dev/null; then \ echo "The following tests failed to complete:" ; \ ls results/*/std.*.out \ - | awk '{split($$0,a,"/"); split(a[3],t,"."); print " ",a[2],":",t[2]}' ; \ + | awk '{split($$0,a,"/"); split(a[3],t,"."); v=t[2]; if(length(t)>3) v=v"."t[3]; print a[2],":",v}'; \ fi; \ if ls results/*/ocean.stats.*.diff &> /dev/null; then \ echo "The following tests report solution regressions:" ; \ ls results/*/ocean.stats.*.diff \ - | awk '{split($$0,a,"/"); split(a[3],t,"."); print " ",a[2],":",t[3]}' ; \ + | awk '{split($$0,a,"/"); split(a[3],t,"."); v=t[3]; if(length(t)>4) v=v"."t[4]; print a[2],":",v}'; \ fi; \ if ls results/*/chksum_diag.*.diff &> /dev/null; then \ echo "The following tests report diagnostic regressions:" ; \ ls results/*/chksum_diag.*.diff \ - | awk '{split($$0,a,"/"); split(a[3],t,"."); print " ",a[2],":",t[2]}' ; \ + | awk '{split($$0,a,"/"); split(a[3],t,"."); v=t[2]; if(length(t)>3) v=v"."t[3]; print a[2],":",v}'; \ fi; \ false ; \ else \ diff --git a/.travis.yml b/.travis.yml index 6b0b4c2a5e..4ceab0a438 100644 --- a/.travis.yml +++ b/.travis.yml @@ -5,19 +5,16 @@ language: c dist: bionic -# --depth flag is breaking our merge, try disabling it -# NOTE: We may be able to go back to depth=50 in production -git: - depth: false - addons: apt: sources: - ubuntu-toolchain-r-test packages: - - tcsh pkg-config netcdf-bin libnetcdf-dev libnetcdff-dev openmpi-bin libopenmpi-dev gfortran + - tcsh pkg-config netcdf-bin libnetcdf-dev libnetcdff-dev gfortran + - mpich libmpich-dev - doxygen graphviz flex bison cmake - python-numpy python-netcdf4 + - bc jobs: include: @@ -58,3 +55,17 @@ jobs: - echo -en 'travis_fold:end:script.1\\r' - make -k -s test.regressions - make test.summary + + - arch: arm64 + env: + - JOB="Configuration testing" + - DO_REGRESSION_TESTS=false + - DO_REPRO_TESTS=false + - MKMF_TEMPLATE=linux-ubuntu-xenial-gnu.mk + script: + - cd .testing + - echo 'Build executables...' && echo -en 'travis_fold:start:script.1\\r' + - make all + - echo -en 'travis_fold:end:script.1\\r' + - make -k -s test + - make test.summary From 2baaea7cd580aa9539d4eb4f0dbf5f2ccbae7821 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 20 May 2020 11:21:45 -0400 Subject: [PATCH 37/40] Travis: Renaming jobs to reflect architecture --- .travis.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.travis.yml b/.travis.yml index 4ceab0a438..6bf509ce8c 100644 --- a/.travis.yml +++ b/.travis.yml @@ -28,7 +28,7 @@ jobs: - test ! -s doxy_errors - env: - - JOB="Configuration testing" + - JOB="x86 Configuration testing" - DO_REGRESSION_TESTS=false - MKMF_TEMPLATE=linux-ubuntu-xenial-gnu.mk script: @@ -42,7 +42,7 @@ jobs: # NOTE: Code coverage upload is here to reduce load imbalance - if: type = pull_request env: - - JOB="Regression testing" + - JOB="x86 Regression testing" - DO_REGRESSION_TESTS=true - REPORT_COVERAGE=true - MKMF_TEMPLATE=linux-ubuntu-xenial-gnu.mk @@ -58,7 +58,7 @@ jobs: - arch: arm64 env: - - JOB="Configuration testing" + - JOB="ARM64 Configuration testing" - DO_REGRESSION_TESTS=false - DO_REPRO_TESTS=false - MKMF_TEMPLATE=linux-ubuntu-xenial-gnu.mk From 7369c39dd13a87e2f6077e5a40b2cb415a111174 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 20 May 2020 20:11:54 -0400 Subject: [PATCH 38/40] +Add just_read_params arg to entrain_diffusive_init Added the new optional argument just_read_params to entrain_diffusive_init, which if present and true causes this routine to read the parameters used by the entrain_diffusive module without logging them. Also set this new parameter based on whether ALE remapping is being used. This prevents the diagnostics and parameters from this module from being offered or logged when they are not actually available, and it will prevent the parameter CORRECT_DENSITY from being logged in ALE-configured runs. All answers are bitwise identical, but there are changes in some MOM_parameter_doc and available_diags files. --- .../vertical/MOM_diabatic_driver.F90 | 14 ++---- .../vertical/MOM_entrain_diffusive.F90 | 47 ++++++++++--------- 2 files changed, 30 insertions(+), 31 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index d753afc97b..0bd3138670 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -2273,11 +2273,8 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e endif - ! This block sets ea, eb from Kd or Kd_int. - ! Otherwise, call entrainment_diffusive() which sets ea and eb - ! based on KD and target densities (ie. does remapping as well). - ! When not using ALE, calculate layer entrainments/detrainments from - ! diffusivities and differences between layer and target densities + ! Calculate layer entrainments and detrainments from diffusivities and differences between + ! layer and target densities (i.e. do remapping as well as diffusion). call cpu_clock_begin(id_clock_entrain) ! Calculate appropriately limited diapycnal mass fluxes to account ! for diapycnal diffusion and advection. Sets: ea, eb. Changes: kb @@ -3688,7 +3685,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di ! False. CS%use_CVMix_conv = CVMix_conv_init(Time, G, GV, US, param_file, diag, CS%CVMix_conv_csp) - call entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS%entrain_diffusive_CSp) + call entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS%entrain_diffusive_CSp, & + just_read_params=CS%useALEalgorithm) ! initialize the geothermal heating module if (CS%use_geothermal) & @@ -3763,12 +3761,10 @@ subroutine diabatic_driver_end(CS) call entrain_diffusive_end(CS%entrain_diffusive_CSp) call set_diffusivity_end(CS%set_diff_CSp) - if (CS%useKPP) then + if (CS%useKPP) then deallocate( CS%KPP_buoy_flux ) deallocate( CS%KPP_temp_flux ) deallocate( CS%KPP_salt_flux ) - endif - if (CS%useKPP) then deallocate( CS%KPP_NLTheat ) deallocate( CS%KPP_NLTscalar ) call KPP_end(CS%KPP_CSp) diff --git a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 index 4e30756f7b..1be3421534 100644 --- a/src/parameterizations/vertical/MOM_entrain_diffusive.F90 +++ b/src/parameterizations/vertical/MOM_entrain_diffusive.F90 @@ -2081,7 +2081,7 @@ end subroutine find_maxF_kb !> This subroutine initializes the parameters and memory associated with the !! entrain_diffusive module. -subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS) +subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS, just_read_params) type(time_type), intent(in) :: Time !< The current model time. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -2092,18 +2092,15 @@ subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS) !! output. type(entrain_diffusive_CS), pointer :: CS !< A pointer that is set to point to the control !! structure. -! for this module -! Arguments: Time - The current model time. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) diag - A structure that is used to regulate diagnostic output. -! (in/out) CS - A pointer that is set to point to the control structure -! for this module + logical, optional, intent(in) :: just_read_params !< If present and true, this call will + !! only read parameters logging them or registering + !! any diagnostics + + ! Local variables real :: decay_length, dt, Kd -! This include declares and sets the variable "version". -#include "version_variable.h" + logical :: just_read ! If true, just read parameters but do nothing else. + ! This include declares and sets the variable "version". +# include "version_variable.h" character(len=40) :: mdl = "MOM_entrain_diffusive" ! This module's name. if (associated(CS)) then @@ -2113,37 +2110,43 @@ subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS) endif allocate(CS) + just_read = .false. ; if (present(just_read_params)) just_read = just_read_params + CS%diag => diag CS%bulkmixedlayer = (GV%nkml > 0) ! Set default, read and log parameters - call log_version(param_file, mdl, version, "") + if (.not.just_read) call log_version(param_file, mdl, version, "") call get_param(param_file, mdl, "CORRECT_DENSITY", CS%correct_density, & "If true, and USE_EOS is true, the layer densities are "//& "restored toward their target values by the diapycnal "//& "mixing, as described in Hallberg (MWR, 2000).", & - default=.true.) + default=.true., do_not_log=just_read) call get_param(param_file, mdl, "MAX_ENT_IT", CS%max_ent_it, & "The maximum number of iterations that may be used to "//& - "calculate the interior diapycnal entrainment.", default=5) + "calculate the interior diapycnal entrainment.", default=5, do_not_log=just_read) ! In this module, KD is only used to set the default for TOLERANCE_ENT. [m2 s-1] call get_param(param_file, mdl, "KD", Kd, fail_if_missing=.true.) call get_param(param_file, mdl, "DT", dt, & "The (baroclinic) dynamics time step.", units = "s", & - fail_if_missing=.true.) + fail_if_missing=.true., do_not_log=just_read) ! CS%Tolerance_Ent = MAX(100.0*GV%Angstrom_H,1.0e-4*sqrt(dt*Kd)) ! call get_param(param_file, mdl, "TOLERANCE_ENT", CS%Tolerance_Ent, & "The tolerance with which to solve for entrainment values.", & - units="m", default=MAX(100.0*GV%Angstrom_m,1.0e-4*sqrt(dt*Kd)), scale=GV%m_to_H) + units="m", default=MAX(100.0*GV%Angstrom_m,1.0e-4*sqrt(dt*Kd)), scale=GV%m_to_H, do_not_log=just_read) CS%Rho_sig_off = 1000.0*US%kg_m3_to_R - CS%id_Kd = register_diag_field('ocean_model', 'Kd_effective', diag%axesTL, Time, & - 'Diapycnal diffusivity as applied', 'm2 s-1', conversion=US%Z2_T_to_m2_s) - CS%id_diff_work = register_diag_field('ocean_model', 'diff_work', diag%axesTi, Time, & - 'Work actually done by diapycnal diffusion across each interface', & - 'W m-2', conversion=US%RZ3_T3_to_W_m2) + if (.not.just_read) then + CS%id_Kd = register_diag_field('ocean_model', 'Kd_effective', diag%axesTL, Time, & + 'Diapycnal diffusivity as applied', 'm2 s-1', conversion=US%Z2_T_to_m2_s) + CS%id_diff_work = register_diag_field('ocean_model', 'diff_work', diag%axesTi, Time, & + 'Work actually done by diapycnal diffusion across each interface', & + 'W m-2', conversion=US%RZ3_T3_to_W_m2) + endif + + if (just_read) deallocate(CS) end subroutine entrain_diffusive_init From 64651d5f3ae290928a4c67bac68b0ff119a19402 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 21 May 2020 15:32:09 -0400 Subject: [PATCH 39/40] CodeCov: Set number of reports to 8 This setting will prevent CodeCov from updating the repository until all 8 test suites are completed. The number is fixed to 8 in the file, so we now need to keep these in sync until this can be sorted out some other way. --- .codecov.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.codecov.yml b/.codecov.yml index 576633bf6a..05fe474ab3 100644 --- a/.codecov.yml +++ b/.codecov.yml @@ -6,3 +6,6 @@ coverage: patch: default: threshold: 100% +comment: + # This must be set to the number of test cases (TCs) + after_n_builds: 8 From eb28c73a57763b0a7a4c32704b168b085c547693 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Fri, 22 May 2020 15:35:44 +0000 Subject: [PATCH 40/40] Switch to intel/18 on gaea - The compiler versions used in the gaea pipelines are encoded in the MRS package of scripts. This commit switches to the master branch which will now pull the latest version of MRS. --- .gitlab-ci.yml | 2 -- 1 file changed, 2 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 39b63c8f85..5a05694fef 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -32,13 +32,11 @@ setup: - git clone --recursive http://gitlab.gfdl.noaa.gov/ogrp/Gaea-stats-MOM6-examples.git tests && cd tests # Install / update testing scripts - git clone https://github.com/adcroft/MRS.git MRS - - (cd MRS ; git checkout xanadu-fms) # Update MOM6-examples and submodules - (cd MOM6-examples && git checkout . && git checkout dev/gfdl && git pull && git submodule init && git submodule update) - (cd MOM6-examples/src/MOM6 && git submodule update) - test -d MOM6-examples/src/LM3 || make -f MRS/Makefile.clone clone_gfdl -s - make -f MRS/Makefile.clone MOM6-examples/.datasets -s - #- (cd MOM6-examples/src/mkmf && git pull https://github.com/adcroft/mkmf.git add_coverage_mode) - env > gitlab_session.log # Cache everything under tests to unpack for each subsequent stage - cd ../ ; time tar zcf $CACHE_DIR/tests_$CI_PIPELINE_ID.tgz tests