From 76a87b41d3626b27aef80e6c8082abd016862463 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 10 Sep 2018 19:32:27 -0400 Subject: [PATCH 1/4] +Added an optional halo argument to make_frazil Added an optional argument to make_frazil and adjust_salt that will cause them to work on an extended region beyond the computational domain. All answers are bitwise identical, but there are new optional arguments to publicly visible subroutines. --- src/parameterizations/vertical/MOM_diabatic_aux.F90 | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index f662eda365..5fa4125fcb 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -77,7 +77,7 @@ module MOM_diabatic_aux !! This subroutine warms any water that is colder than the (currently !! surface) freezing point up to the freezing point and accumulates !! the required heat (in J m-2) in tv%frazil. -subroutine make_frazil(h, tv, G, GV, CS, p_surf) +subroutine make_frazil(h, tv, G, GV, CS, p_surf, halo) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & @@ -88,6 +88,7 @@ subroutine make_frazil(h, tv, G, GV, CS, p_surf) !! call to diabatic_aux_init. real, dimension(SZI_(G),SZJ_(G)), & optional, intent(in) :: p_surf !< The pressure at the ocean surface, in Pa. + integer, optional, intent(in) :: halo !< Halo width over which to calculate frazil ! Frazil formation keeps the temperature above the freezing point. ! This subroutine warms any water that is colder than the (currently @@ -113,7 +114,11 @@ subroutine make_frazil(h, tv, G, GV, CS, p_surf) logical :: T_fr_set ! True if the freezing point has been calculated for a ! row of points. integer :: i, j, k, is, ie, js, je, nz + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + if (present(halo)) then + is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo + endif call cpu_clock_begin(id_clock_frazil) @@ -309,7 +314,7 @@ end subroutine differential_diffuse_T_S !> This subroutine keeps salinity from falling below a small but positive threshold. !! This usually occurs when the ice model attempts to extract more salt then !! is actually available to it from the ocean. -subroutine adjust_salt(h, tv, G, GV, CS) +subroutine adjust_salt(h, tv, G, GV, CS, halo) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & @@ -318,6 +323,7 @@ subroutine adjust_salt(h, tv, G, GV, CS) !! available thermodynamic fields. type(diabatic_aux_CS), intent(in) :: CS !< The control structure returned by a previous !! call to diabatic_aux_init. + integer, optional, intent(in) :: halo !< Halo width over which to work ! local variables real :: salt_add_col(SZI_(G),SZJ_(G)) !< The accumulated salt requirement @@ -325,6 +331,9 @@ subroutine adjust_salt(h, tv, G, GV, CS) real :: mc !< A layer's mass kg m-2 . integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + if (present(halo)) then + is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo + endif ! call cpu_clock_begin(id_clock_adjust_salt) From fd6fb2e8bd8bf9409aef126f515372075c167524 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 10 Sep 2018 19:34:06 -0400 Subject: [PATCH 2/4] +Added an optional halo argument to geothermal Added an optional argument to geothermal that will cause it to work on an extended region beyond the computational domain, along with a pass_var call on the static geothermal heat flux to enable this to work properly. All answers are bitwise identical, but there is a new optional argument to a publicly visible subroutines. --- src/parameterizations/vertical/MOM_geothermal.F90 | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_geothermal.F90 b/src/parameterizations/vertical/MOM_geothermal.F90 index 360c3a791d..b1fc1fd177 100644 --- a/src/parameterizations/vertical/MOM_geothermal.F90 +++ b/src/parameterizations/vertical/MOM_geothermal.F90 @@ -5,6 +5,7 @@ module MOM_geothermal use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : register_static_field, time_type, diag_ctrl +use MOM_domains, only : pass_var use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_io, only : MOM_read_data, slasher @@ -44,7 +45,7 @@ module MOM_geothermal !! the partial derivative of the coordinate density with temperature is positive !! or very small, the layers are simply heated in place. Any heat that can not !! be applied to the ocean is returned (WHERE)? -subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS) +subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS, halo) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid !! structure. @@ -69,6 +70,7 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS) type(geothermal_CS), pointer :: CS !< The control structure returned by !! a previous call to !! geothermal_init. + integer, optional, intent(in) :: halo !< Halo width over which to work ! Local variables real, dimension(SZI_(G)) :: & heat_rem, & ! remaining heat (H * degC) @@ -105,6 +107,9 @@ subroutine geothermal(h, tv, dt, ea, eb, G, GV, CS) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke + if (present(halo)) then + is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo + endif if (.not. associated(CS)) call MOM_error(FATAL, "MOM_geothermal: "//& "Module must be initialized before it is used.") @@ -377,6 +382,7 @@ subroutine geothermal_init(Time, G, param_file, diag, CS) CS%geo_heat(i,j) = G%mask2dT(i,j) * scale enddo ; enddo endif + call pass_var(CS%geo_heat, G%domain) ! post the static geothermal heating field id = register_static_field('ocean_model', 'geo_heat', diag%axesT1, & From 1cc8d7202223876e80ceb264a907d80faad299ad Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 10 Sep 2018 19:34:28 -0400 Subject: [PATCH 3/4] +Added optional halo_TS arg to set_diffusivity_init Added an optional argument to set_diffusivity_init to return the size of the temperature and salinity halos that are expected to be valid to work with the options in set_diffusivity. Also corrected the checksum call for Kv_shear_Bu. All answers are bitwise identical, but there is a new optional argument to a publicly visible subroutine. --- .../vertical/MOM_set_diffusivity.F90 | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 8d3206303c..7ef7f972ec 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -354,7 +354,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (associated(visc%Kv_shear)) visc%Kv_shear(:,:,:) = 0.0 ! needed for other parameterizations if (CS%debug) then call hchksum(visc%Kd_shear, "after calc_KS_vert visc%Kd_shear",G%HI) - call Bchksum(visc%Kv_shear, "after calc_KS_vert visc%Kv_shear_Bu",G%HI) + call Bchksum(visc%Kv_shear_Bu, "after calc_KS_vert visc%Kv_shear_Bu",G%HI) call Bchksum(visc%TKE_turb, "after calc_KS_vert visc%TKE_turb",G%HI) endif else @@ -1892,7 +1892,7 @@ subroutine set_density_ratios(h, tv, kb, G, GV, CS, j, ds_dsp1, rho_0) end subroutine set_density_ratios subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp, int_tide_CSp, & - tm_CSp) + tm_CSp, halo_TS) type(time_type), intent(in) :: Time !< The current model time type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -1907,6 +1907,8 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp !! structure (BDM) type(tidal_mixing_cs), pointer :: tm_csp !< pointer to tidal mixing control !! structure + integer, optional, intent(out) :: halo_TS !< The halo size of tracer points that must be + !! valid for the calculations in set_diffusivity. ! local variables real :: decay_length @@ -2228,6 +2230,11 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp if (CS%use_CVMix_ddiff) & id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion via CVMix)', grain=CLOCK_MODULE) + if (present(halo_TS)) then + halo_TS = 0 + if (CS%Vertex_Shear) halo_TS = 1 + endif + end subroutine set_diffusivity_init !> Clear pointers and dealocate memory From 9ede9b61a32119251bc2920388af117a881b222c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Mon, 10 Sep 2018 19:35:05 -0400 Subject: [PATCH 4/4] +(*)Corrected halo data when VERTEX_SHEAR=True Added code to work on extra points or do appropriate halo updates for those calls that modify temperatures, salinities and thicknesses before the call to set_diffusivity in both diabatic and legacy_diabatic. All answers are bitwise identical in the existing MOM6 test cases, but this corrects a problem with answers that do not reproduce across PE layouts when VERTEX_SHEAR=True. --- .../vertical/MOM_diabatic_driver.F90 | 36 +++++++++++-------- 1 file changed, 22 insertions(+), 14 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index e3806fd684..5985c6f054 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -154,7 +154,8 @@ module MOM_diabatic_driver !! fluxes are applied, in m. real :: evap_CFL_limit = 0.8 !< The largest fraction of a layer that can be !! evaporated in one time-step (non-dim). - + integer :: halo_TS_diff = 0 !< The temperature, salinity and thickness halo size that + !! must be valid for the diffusivity calculations. logical :: useKPP = .false. !< use CVMix/KPP diffusivities and non-local transport logical :: salt_reject_below_ML !< If true, add salt below mixed layer (layer mode only) logical :: KPPisPassive !< If true, KPP is in passive mode, not changing answers. @@ -379,7 +380,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & integer :: z_ids(7) ! id numbers of diagnostics to be interpolated to depth integer :: dir_flag ! An integer encoding the directions in which to do halo updates. logical :: showCallTree ! If true, show the call tree - integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m + integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m, halo integer :: ig, jg ! global indices for testing testing itide point source (BDM) logical :: avg_enabled ! for testing internal tides (BDM) @@ -447,9 +448,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif if (associated(fluxes%p_surf_full)) then - call make_frazil(h, tv, G, GV, CS%diabatic_aux_CSp, fluxes%p_surf_full) + call make_frazil(h, tv, G, GV, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff) else - call make_frazil(h, tv, G, GV, CS%diabatic_aux_CSp) + call make_frazil(h, tv, G, GV, CS%diabatic_aux_CSp, halo=CS%halo_TS_diff) endif if (showCallTree) call callTree_waypoint("done with 1st make_frazil (diabatic)") @@ -465,15 +466,16 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%debugConservation) call MOM_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, G) if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then -!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_orig,h,eaml,ebml) - do k=1,nz ; do j=js,je ; do i=is,ie + halo = CS%halo_TS_diff + !$OMP parallel do default(shared) + do k=1,nz ; do j=js-halo,je+halo ; do i=is-halo,ie+halo h_orig(i,j,k) = h(i,j,k) ; eaml(i,j,k) = 0.0 ; ebml(i,j,k) = 0.0 enddo ; enddo ; enddo endif if (CS%use_geothermal) then call cpu_clock_begin(id_clock_geothermal) - call geothermal(h, tv, dt, eaml, ebml, G, GV, CS%geothermal_CSp) + call geothermal(h, tv, dt, eaml, ebml, G, GV, CS%geothermal_CSp, halo=CS%halo_TS_diff) call cpu_clock_end(id_clock_geothermal) if (showCallTree) call callTree_waypoint("geothermal (diabatic)") if (CS%debugConservation) call MOM_state_stats('geothermal', u, v, h, tv%T, tv%S, G) @@ -1258,7 +1260,7 @@ subroutine legacy_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_en integer :: z_ids(7) ! id numbers of diagnostics to be interpolated to depth integer :: dir_flag ! An integer encoding the directions in which to do halo updates. logical :: showCallTree ! If true, show the call tree - integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m + integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m, halo integer :: ig, jg ! global indices for testing testing itide point source (BDM) logical :: avg_enabled ! for testing internal tides (BDM) @@ -1323,9 +1325,9 @@ subroutine legacy_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_en endif if (associated(fluxes%p_surf_full)) then - call make_frazil(h, tv, G, GV, CS%diabatic_aux_CSp, fluxes%p_surf_full) + call make_frazil(h, tv, G, GV, CS%diabatic_aux_CSp, fluxes%p_surf_full, halo=CS%halo_TS_diff) else - call make_frazil(h, tv, G, GV, CS%diabatic_aux_CSp) + call make_frazil(h, tv, G, GV, CS%diabatic_aux_CSp, halo=CS%halo_TS_diff) endif if (showCallTree) call callTree_waypoint("done with 1st make_frazil (diabatic)") @@ -1340,15 +1342,16 @@ subroutine legacy_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_en if (CS%debugConservation) call MOM_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, G) if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then -!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_orig,h,eaml,ebml) - do k=1,nz ; do j=js,je ; do i=is,ie + halo = CS%halo_TS_diff + !$OMP parallel do default(shared) + do k=1,nz ; do j=js-halo,je+halo ; do i=is-halo,ie+halo h_orig(i,j,k) = h(i,j,k) ; eaml(i,j,k) = 0.0 ; ebml(i,j,k) = 0.0 enddo ; enddo ; enddo endif if (CS%use_geothermal) then call cpu_clock_begin(id_clock_geothermal) - call geothermal(h, tv, dt, eaml, ebml, G, GV, CS%geothermal_CSp) + call geothermal(h, tv, dt, eaml, ebml, G, GV, CS%geothermal_CSp, halo=CS%halo_TS_diff) call cpu_clock_end(id_clock_geothermal) if (showCallTree) call callTree_waypoint("geothermal (diabatic)") if (CS%debugConservation) call MOM_state_stats('geothermal', u, v, h, tv%T, tv%S, G) @@ -1478,6 +1481,11 @@ subroutine legacy_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_en ! Sets: Kd, Kd_int, visc%Kd_extra_T, visc%Kd_extra_S ! Also changes: visc%Kd_shear, visc%TKE_turb (not clear that TKE_turb is used as input ???? ! And sets visc%Kv_shear + if ((CS%halo_TS_diff > 0) .and. (CS%ML_mix_first > 0.0)) then + if (associated(tv%T)) call pass_var(tv%T, G%Domain, halo=CS%halo_TS_diff, complete=.false.) + if (associated(tv%T)) call pass_var(tv%S, G%Domain, halo=CS%halo_TS_diff, complete=.false.) + call pass_var(h, G%domain, halo=CS%halo_TS_diff, complete=.true.) + endif call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, CS%optics, visc, dt, G, GV, CS%set_diff_CSp, Kd, Kd_int) call cpu_clock_end(id_clock_set_diffusivity) if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)") @@ -3277,7 +3285,7 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, ! initialize module for setting diffusivities call set_diffusivity_init(Time, G, GV, param_file, diag, CS%set_diff_CSp, diag_to_Z_CSp, & - CS%int_tide_CSp, CS%tidal_mixing_CSp) + CS%int_tide_CSp, CS%tidal_mixing_CSp, CS%halo_TS_diff) ! set up the clocks for this module