diff --git a/config_src/infra/FMS1/MOM_constants.F90 b/config_src/infra/FMS1/MOM_constants.F90 index 2db177e08c..a632267a7f 100644 --- a/config_src/infra/FMS1/MOM_constants.F90 +++ b/config_src/infra/FMS1/MOM_constants.F90 @@ -3,12 +3,16 @@ module MOM_constants ! This file is part of MOM6. See LICENSE.md for the license. -use constants_mod, only : HLV, HLF +use constants_mod, only : FMS_HLV => HLV +use constants_mod, only : FMS_HLF => HLF implicit none ; private -!> The constant offset for converting temperatures in Kelvin to Celsius real, public, parameter :: CELSIUS_KELVIN_OFFSET = 273.15 -public :: HLV, HLF + !< The constant offset for converting temperatures in Kelvin to Celsius [K] +real, public, parameter :: HLV = real(FMS_HLV, kind=kind(1.0)) + !< Latent heat of vaporization [J kg-1] +real, public, parameter :: HLF = real(FMS_HLF, kind=kind(1.0)) + !< Latent heat of fusion [J kg-1] end module MOM_constants diff --git a/config_src/infra/FMS2/MOM_constants.F90 b/config_src/infra/FMS2/MOM_constants.F90 index 2db177e08c..a632267a7f 100644 --- a/config_src/infra/FMS2/MOM_constants.F90 +++ b/config_src/infra/FMS2/MOM_constants.F90 @@ -3,12 +3,16 @@ module MOM_constants ! This file is part of MOM6. See LICENSE.md for the license. -use constants_mod, only : HLV, HLF +use constants_mod, only : FMS_HLV => HLV +use constants_mod, only : FMS_HLF => HLF implicit none ; private -!> The constant offset for converting temperatures in Kelvin to Celsius real, public, parameter :: CELSIUS_KELVIN_OFFSET = 273.15 -public :: HLV, HLF + !< The constant offset for converting temperatures in Kelvin to Celsius [K] +real, public, parameter :: HLV = real(FMS_HLV, kind=kind(1.0)) + !< Latent heat of vaporization [J kg-1] +real, public, parameter :: HLF = real(FMS_HLF, kind=kind(1.0)) + !< Latent heat of fusion [J kg-1] end module MOM_constants diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 7ee90d746a..059e515051 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2883,6 +2883,26 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, & ! reservoirs are used. call open_boundary_register_restarts(HI, GV, US, CS%OBC, CS%tracer_Reg, & param_file, restart_CSp, use_temperature) + if (turns /= 0) then + if (CS%OBC%radiation_BCs_exist_globally) then + OBC_in%rx_normal => CS%OBC%rx_normal + OBC_in%ry_normal => CS%OBC%ry_normal + endif + if (CS%OBC%oblique_BCs_exist_globally) then + OBC_in%rx_oblique_u => CS%OBC%rx_oblique_u + OBC_in%ry_oblique_u => CS%OBC%ry_oblique_u + OBC_in%rx_oblique_v => CS%OBC%rx_oblique_v + OBC_in%ry_oblique_v => CS%OBC%ry_oblique_v + OBC_in%cff_normal_u => CS%OBC%cff_normal_u + OBC_in%cff_normal_v => CS%OBC%cff_normal_v + endif + if (any(CS%OBC%tracer_x_reservoirs_used)) then + OBC_in%tres_x => CS%OBC%tres_x + endif + if (any(CS%OBC%tracer_y_reservoirs_used)) then + OBC_in%tres_y => CS%OBC%tres_y + endif + endif endif if (present(waves_CSp)) then diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index af2beca1fb..aa502bdd4b 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -742,6 +742,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB integer :: ioff, joff integer :: l_seg + real :: factor(SZI_(G),SZJ_(G)) ! If non-zero, work on given points. if (.not.CS%module_is_initialized) call MOM_error(FATAL, & "btstep: Module MOM_barotropic must be initialized before it is used.") @@ -2457,17 +2458,69 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, haloshift=iev-ie, unscale=US%L_to_m**2*GV%H_to_m) endif + do j=jsv,jev + do i=isv,iev + factor(i,j) = CS%IareaT(i,j) + enddo + enddo + + ! Update factor so that nothing changes outside of the OBC (problem for interior OBCs only) + if (associated(OBC)) then ; if (OBC%OBC_pe) then +! do j=jsv,jev +! if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then +! do i=isv,iev-1 ; if (OBC%segnum_u(I,j) /= OBC_NONE) then +! if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then +! factor(i+1,j) = 0.0 +! elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then +! factor(i,j) = 0.0 +! endif +! endif ; enddo +! endif + do i=isv,iev + if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then + do j=jsv,jev + if (OBC%segnum_u(I-1,j) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_u(I-1,j))%direction == OBC_DIRECTION_E) then + factor(i,j) = 0.0 + endif + endif + if (OBC%segnum_u(I,j) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + factor(i,j) = 0.0 + endif + endif + enddo + endif + enddo + do j=jsv,jev + if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then + do i=isv,iev + if (OBC%segnum_v(i,J-1) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then + factor(i,j) = 0.0 + endif + endif + if (OBC%segnum_v(i,J) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + factor(i,j) = 0.0 + endif + endif + enddo + endif + enddo + endif ; endif + if (integral_BT_cont) then !$OMP do do j=jsv,jev ; do i=isv,iev - eta(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + CS%IareaT(i,j) * & + eta(i,j) = (eta_IC(i,j) + n*eta_src(i,j)) + factor(i,j) * & ((uhbt_int(I-1,j) - uhbt_int(I,j)) + (vhbt_int(i,J-1) - vhbt_int(i,J))) eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n) enddo ; enddo else !$OMP do do j=jsv,jev ; do i=isv,iev - eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * CS%IareaT(i,j)) * & + eta(i,j) = (eta(i,j) + eta_src(i,j)) + (dtbt * factor(i,j)) * & ((uhbt(I-1,j) - uhbt(I,j)) + (vhbt(i,J-1) - vhbt(i,J))) eta_wtd(i,j) = eta_wtd(i,j) + eta(i,j) * wt_eta(n) enddo ; enddo @@ -5270,7 +5323,7 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS) character(len=40) :: mdl = "MOM_barotropic" ! This module's name. integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB real :: am2, ak1 !< Bandwidth parameters of the M2 and K1 streaming filters [nondim] - real :: om2, ok1 !< Target frequencies of the M2 and K1 streaming filters [T-1 ~> s-1] + real :: om2, ok1 !< Target frequencies of the M2 and K1 streaming filters [rad T-1 ~> rad s-1] isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed IsdB = HI%IsdB ; IedB = HI%IedB ; JsdB = HI%JsdB ; JedB = HI%JedB @@ -5301,13 +5354,13 @@ subroutine register_barotropic_restarts(HI, GV, US, param_file, CS, restart_CS) "Frequency of the M2 tidal constituent. "//& "This is only used if TIDES and TIDE_M2"// & " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and M2"// & - " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=tidal_frequency("M2"), & + " is in OBC_TIDE_CONSTITUENTS.", units="rad s-1", default=tidal_frequency("M2"), & scale=US%T_to_s, do_not_log=.true.) call get_param(param_file, mdl, "TIDE_K1_FREQ", ok1, & "Frequency of the K1 tidal constituent. "//& "This is only used if TIDES and TIDE_K1"// & " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and K1"// & - " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=tidal_frequency("K1"), & + " is in OBC_TIDE_CONSTITUENTS.", units="rad s-1", default=tidal_frequency("K1"), & scale=US%T_to_s, do_not_log=.true.) ALLOC_(CS%ubtav(IsdB:IedB,jsd:jed)) ; CS%ubtav(:,:) = 0.0 diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index f89c8953ab..b15d28194e 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -320,7 +320,8 @@ module MOM_open_boundary logical :: add_tide_constituents = .false. !< If true, add tidal constituents to the boundary elevation !! and velocity. Will be set to true if n_tide_constituents > 0. character(len=2), allocatable, dimension(:) :: tide_names !< Names of tidal constituents to add to the boundary data. - real, allocatable, dimension(:) :: tide_frequencies !< Angular frequencies of chosen tidal constituents [T-1 ~> s-1]. + real, allocatable, dimension(:) :: tide_frequencies !< Angular frequencies of chosen tidal + !! constituents [rad T-1 ~> rad s-1]. real, allocatable, dimension(:) :: tide_eq_phases !< Equilibrium phases of chosen tidal constituents [rad]. real, allocatable, dimension(:) :: tide_fn !< Amplitude modulation of boundary tides by nodal cycle [nondim]. real, allocatable, dimension(:) :: tide_un !< Phase modulation of boundary tides by nodal cycle [rad]. @@ -352,24 +353,24 @@ module MOM_open_boundary type(remapping_CS), pointer :: remap_h_CS=> NULL() !< ALE remapping control structure for !! thickness-based fields on segments type(OBC_registry_type), pointer :: OBC_Reg => NULL() !< Registry type for boundaries - real, allocatable :: rx_normal(:,:,:) !< Array storage for normal phase speed for EW radiation OBCs in units of + real, pointer :: rx_normal(:,:,:) !< Array storage for normal phase speed for EW radiation OBCs in units of !! grid points per timestep [nondim] - real, allocatable :: ry_normal(:,:,:) !< Array storage for normal phase speed for NS radiation OBCs in units of + real, pointer :: ry_normal(:,:,:) !< Array storage for normal phase speed for NS radiation OBCs in units of !! grid points per timestep [nondim] - real, allocatable :: rx_oblique_u(:,:,:) !< X-direction oblique boundary condition radiation speeds squared + real, pointer :: rx_oblique_u(:,:,:) !< X-direction oblique boundary condition radiation speeds squared !! at u points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: ry_oblique_u(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared + real, pointer :: ry_oblique_u(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared !! at u points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: rx_oblique_v(:,:,:) !< X-direction oblique boundary condition radiation speeds squared + real, pointer :: rx_oblique_v(:,:,:) !< X-direction oblique boundary condition radiation speeds squared !! at v points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: ry_oblique_v(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared + real, pointer :: ry_oblique_v(:,:,:) !< Y-direction oblique boundary condition radiation speeds squared !! at v points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: cff_normal_u(:,:,:) !< Denominator for normalizing EW oblique boundary condition radiation + real, pointer :: cff_normal_u(:,:,:) !< Denominator for normalizing EW oblique boundary condition radiation !! rates at u points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: cff_normal_v(:,:,:) !< Denominator for normalizing NS oblique boundary condition radiation + real, pointer :: cff_normal_v(:,:,:) !< Denominator for normalizing NS oblique boundary condition radiation !! rates at v points for restarts [L2 T-2 ~> m2 s-2] - real, allocatable :: tres_x(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] - real, allocatable :: tres_y(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] + real, pointer :: tres_x(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] + real, pointer :: tres_y(:,:,:,:) !< Array storage of tracer reservoirs for restarts, in unscaled units [conc] logical :: debug !< If true, write verbose checksums for debugging purposes. real :: silly_h !< A silly value of thickness outside of the domain that can be used to test !! the independence of the OBCs to this external data [Z ~> m]. @@ -705,10 +706,13 @@ subroutine open_boundary_config(G, US, param_file, OBC) "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", OBC%om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) + call get_param(param_file, mdl, "OBC_REMAPPING_USE_OM4_SUBCELLS", OBC%om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for neutral diffusion. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=OBC%om4_remap_via_sub_cells) endif ! OBC%number_of_segments > 0 @@ -1231,7 +1235,7 @@ subroutine initialize_obc_tides(OBC, US, param_file) "This is only used if TIDES and TIDE_"//trim(OBC%tide_names(c))// & " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and "//trim(OBC%tide_names(c))//& " is in OBC_TIDE_CONSTITUENTS.", & - units="s-1", default=tidal_frequency(trim(OBC%tide_names(c))), scale=US%T_to_s) + units="rad s-1", default=tidal_frequency(trim(OBC%tide_names(c))), scale=US%T_to_s) ! Find equilibrium phase if needed if (OBC%add_eq_phase) then @@ -1476,7 +1480,7 @@ subroutine setup_u_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_y) "Timescales in days for nudging along a segment, "//& "for inflow, then outflow. Setting both to zero should "//& "behave like SIMPLE obcs for the baroclinic velocities.", & - fail_if_missing=.true., default=0., units="days", scale=86400.0*US%s_to_T) + fail_if_missing=.true., units="days", scale=86400.0*US%s_to_T) OBC%segment(l_seg)%Velocity_nudging_timescale_in = tnudge(1) OBC%segment(l_seg)%Velocity_nudging_timescale_out = tnudge(2) deallocate(tnudge) @@ -1617,7 +1621,7 @@ subroutine setup_v_point_obc(OBC, G, US, segment_str, l_seg, PF, reentrant_x) "Timescales in days for nudging along a segment, "//& "for inflow, then outflow. Setting both to zero should "//& "behave like SIMPLE obcs for the baroclinic velocities.", & - fail_if_missing=.true., default=0., units="days", scale=86400.0*US%s_to_T) + fail_if_missing=.true., units="days", scale=86400.0*US%s_to_T) OBC%segment(l_seg)%Velocity_nudging_timescale_in = tnudge(1) OBC%segment(l_seg)%Velocity_nudging_timescale_out = tnudge(2) deallocate(tnudge) @@ -1948,15 +1952,15 @@ subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CS) call create_group_pass(OBC%pass_oblique, OBC%cff_normal_u, OBC%cff_normal_v, G%Domain, To_All+Scalar_Pair) call do_group_pass(OBC%pass_oblique, G%Domain) endif - if (allocated(OBC%tres_x) .and. allocated(OBC%tres_y)) then + if (associated(OBC%tres_x) .and. associated(OBC%tres_y)) then do m=1,OBC%ntr call pass_vector(OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%Domain, To_All+Scalar_Pair) enddo - elseif (allocated(OBC%tres_x)) then + elseif (associated(OBC%tres_x)) then do m=1,OBC%ntr call pass_var(OBC%tres_x(:,:,:,m), G%Domain, position=EAST_FACE) enddo - elseif (allocated(OBC%tres_y)) then + elseif (associated(OBC%tres_y)) then do m=1,OBC%ntr call pass_var(OBC%tres_y(:,:,:,m), G%Domain, position=NORTH_FACE) enddo @@ -2001,16 +2005,16 @@ subroutine open_boundary_dealloc(OBC) if (allocated(OBC%segment)) deallocate(OBC%segment) if (allocated(OBC%segnum_u)) deallocate(OBC%segnum_u) if (allocated(OBC%segnum_v)) deallocate(OBC%segnum_v) - if (allocated(OBC%rx_normal)) deallocate(OBC%rx_normal) - if (allocated(OBC%ry_normal)) deallocate(OBC%ry_normal) - if (allocated(OBC%rx_oblique_u)) deallocate(OBC%rx_oblique_u) - if (allocated(OBC%ry_oblique_u)) deallocate(OBC%ry_oblique_u) - if (allocated(OBC%rx_oblique_v)) deallocate(OBC%rx_oblique_v) - if (allocated(OBC%ry_oblique_v)) deallocate(OBC%ry_oblique_v) - if (allocated(OBC%cff_normal_u)) deallocate(OBC%cff_normal_u) - if (allocated(OBC%cff_normal_v)) deallocate(OBC%cff_normal_v) - if (allocated(OBC%tres_x)) deallocate(OBC%tres_x) - if (allocated(OBC%tres_y)) deallocate(OBC%tres_y) + if (associated(OBC%rx_normal)) nullify(OBC%rx_normal) + if (associated(OBC%ry_normal)) nullify(OBC%ry_normal) + if (associated(OBC%rx_oblique_u)) nullify(OBC%rx_oblique_u) + if (associated(OBC%ry_oblique_u)) nullify(OBC%ry_oblique_u) + if (associated(OBC%rx_oblique_v)) nullify(OBC%rx_oblique_v) + if (associated(OBC%ry_oblique_v)) nullify(OBC%ry_oblique_v) + if (associated(OBC%cff_normal_u)) nullify(OBC%cff_normal_u) + if (associated(OBC%cff_normal_v)) nullify(OBC%cff_normal_v) + if (associated(OBC%tres_x)) nullify(OBC%tres_x) + if (associated(OBC%tres_y)) nullify(OBC%tres_y) if (associated(OBC%remap_z_CS)) deallocate(OBC%remap_z_CS) if (associated(OBC%remap_h_CS)) deallocate(OBC%remap_h_CS) deallocate(OBC) @@ -3369,7 +3373,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US, haloshift=0, symmetric=sym, unscale=1.0/US%L_T_to_m_s**2) endif if (OBC%ntr == 0) return - if (.not. allocated (OBC%tres_x) .or. .not. allocated (OBC%tres_y)) return + if (.not. associated (OBC%tres_x) .or. .not. associated (OBC%tres_y)) return do m=1,OBC%ntr write(var_num,'(I3.3)') m call uvchksum("radiation_OBCs: OBC%tres_[xy]_"//var_num, OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%HI, & @@ -5053,7 +5057,9 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) integer :: i, j integer :: l_seg logical :: fatal_error = .False. - real :: min_depth ! The minimum depth for ocean points [Z ~> m] + real :: min_depth ! The minimum depth for ocean points [Z ~> m] + real :: mask_depth ! The masking depth for ocean points [Z ~> m] + real :: Dmask ! The depth for masking in the same units as G%bathyT [Z ~> m]. integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2 character(len=256) :: mesg ! Message for error messages. real, allocatable, dimension(:,:) :: color, color2 ! For sorting inside from outside, @@ -5063,6 +5069,12 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, & units="m", default=0.0, scale=US%m_to_Z, do_not_log=.true.) + call get_param(param_file, mdl, "MASKING_DEPTH", mask_depth, & + units="m", default=-9999.0, scale=US%m_to_Z, do_not_log=.true.) + + Dmask = mask_depth + if (mask_depth == -9999.0*US%m_to_Z) Dmask = min_depth + ! The reference depth on a dyn_horgrid is 0, otherwise would need: min_depth = min_depth - G%Z_ref allocate(color(G%isd:G%ied, G%jsd:G%jed), source=0.0) @@ -5153,7 +5165,7 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) &"the masking of the outside grid points.")') i, j call MOM_error(WARNING,"MOM mask_outside_OBCs: "//mesg, all_print=.true.) endif - if (color(i,j) == cout) G%bathyT(i,j) = min_depth + if (color(i,j) == cout) G%bathyT(i,j) = Dmask enddo ; enddo if (fatal_error) call MOM_error(FATAL, & "MOM_open_boundary: inconsistent OBC segments.") @@ -5455,8 +5467,8 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) if (G%mask2dT(I+ishift,j) == 0.0) cycle ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep do m=1,segment%tr_Reg%ntseg - ntr_id = segment%tr_reg%Tr(m)%ntr_index - fd_id = segment%tr_reg%Tr(m)%fd_index + ntr_id = segment%tr_Reg%Tr(m)%ntr_index + fd_id = segment%tr_Reg%Tr(m)%fd_index if (fd_id == -1) then resrv_lfac_out = 1.0 resrv_lfac_in = 1.0 @@ -5481,7 +5493,7 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) ((1.0-a_out+a_in)*segment%tr_Reg%Tr(m)%tres(I,j,k)+ & ((u_L_out+a_out)*Reg%Tr(ntr_id)%t(I+ishift,j,k) - & (u_L_in+a_in)*segment%tr_Reg%Tr(m)%t(I,j,k))) - if (allocated(OBC%tres_x)) OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k) + if (associated(OBC%tres_x)) OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k) enddo ; endif enddo enddo @@ -5499,8 +5511,8 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) if (G%mask2dT(i,j+jshift) == 0.0) cycle ! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep do m=1,segment%tr_Reg%ntseg - ntr_id = segment%tr_reg%Tr(m)%ntr_index - fd_id = segment%tr_reg%Tr(m)%fd_index + ntr_id = segment%tr_Reg%Tr(m)%ntr_index + fd_id = segment%tr_Reg%Tr(m)%fd_index if (fd_id == -1) then resrv_lfac_out = 1.0 resrv_lfac_in = 1.0 @@ -5521,7 +5533,7 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg) ((1.0-a_out+a_in)*segment%tr_Reg%Tr(m)%tres(i,J,k) + & ((v_L_out+a_out)*Reg%Tr(ntr_id)%t(i,J+jshift,k) - & (v_L_in+a_in)*segment%tr_Reg%Tr(m)%t(i,J,k))) - if (allocated(OBC%tres_y)) OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) + if (associated(OBC%tres_y)) OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) enddo ; endif enddo enddo @@ -5597,7 +5609,7 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) ! Update tracer concentrations segment%tr_Reg%Tr(m)%tres(I,j,:) = tr_column(:) - if (allocated(OBC%tres_x)) then ; do k=1,nz + if (associated(OBC%tres_x)) then ; do k=1,nz OBC%tres_x(I,j,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(I,j,k) enddo ; endif @@ -5664,7 +5676,7 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell) ! Update tracer concentrations segment%tr_Reg%Tr(m)%tres(i,J,:) = tr_column(:) - if (allocated(OBC%tres_y)) then ; do k=1,nz + if (associated(OBC%tres_y)) then ; do k=1,nz OBC%tres_y(i,J,k,m) = I_scale * segment%tr_Reg%Tr(m)%tres(i,J,k) enddo ; endif @@ -6047,13 +6059,13 @@ subroutine rotate_OBC_init(OBC_in, G, GV, US, param_file, tv, restart_CS, OBC) "If true, Temperature and salinity are used as state "//& "variables.", default=.true., do_not_log=.true.) + if (use_temperature) & + call fill_temp_salt_segments(G, GV, US, OBC, tv) + do l = 1, OBC%number_of_segments call rotate_OBC_segment_data(OBC_in%segment(l), OBC%segment(l), G%HI%turns) enddo - if (use_temperature) & - call fill_temp_salt_segments(G, GV, US, OBC, tv) - call open_boundary_init(G, GV, US, param_file, OBC, restart_CS) end subroutine rotate_OBC_init @@ -6076,6 +6088,14 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%field(n)%handle = segment_in%field(n)%handle segment%field(n)%dz_handle = segment_in%field(n)%dz_handle + if (allocated(segment_in%field(n)%buffer_dst)) then + call allocate_rotated_array(segment_in%field(n)%buffer_dst, & + lbound(segment_in%field(n)%buffer_dst), turns, & + segment%field(n)%buffer_dst) + call rotate_array(segment_in%field(n)%buffer_dst, turns, & + segment%field(n)%buffer_dst) + endif + if (modulo(turns, 2) /= 0) then select case (segment_in%field(n)%name) case ('U') @@ -6086,6 +6106,7 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%field(n)%name = 'Vphase' case ('V') segment%field(n)%name = 'U' + segment%field(n)%buffer_dst(:,:,:) = -segment%field(n)%buffer_dst(:,:,:) case ('Vamp') segment%field(n)%name = 'Uamp' case ('Vphase') @@ -6122,6 +6143,93 @@ subroutine rotate_OBC_segment_data(segment_in, segment, turns) segment%field(n)%value = segment_in%field(n)%value enddo + if (allocated(segment_in%SSH)) & + call rotate_array(segment_in%SSH, turns, segment%SSH) + if (allocated(segment_in%cg)) & + call rotate_array(segment_in%cg, turns, segment%cg) + if (allocated(segment_in%htot)) & + call rotate_array(segment_in%htot, turns, segment%htot) + if (allocated(segment_in%dztot)) & + call rotate_array(segment_in%dztot, turns, segment%dztot) + if (allocated(segment_in%h)) & + call rotate_array(segment_in%h, turns, segment%h) + if (allocated(segment_in%normal_vel)) & + call rotate_array(segment_in%normal_vel, turns, segment%normal_vel) + if (allocated(segment_in%normal_trans)) & + call rotate_array(segment_in%normal_trans, turns, segment%normal_trans) + if (allocated(segment_in%normal_vel_bt)) & + call rotate_array(segment_in%normal_vel_bt, turns, segment%normal_vel_bt) + if (allocated(segment_in%tangential_vel)) & + call rotate_array(segment_in%tangential_vel, turns, segment%tangential_vel) + if (allocated(segment_in%tangential_grad)) & + call rotate_array(segment_in%tangential_grad, turns, segment%tangential_grad) + if (allocated(segment_in%grad_normal)) & + call rotate_array(segment_in%grad_normal, turns, segment%grad_normal) + if (allocated(segment_in%grad_tan)) & + call rotate_array(segment_in%grad_tan, turns, segment%grad_tan) + if (allocated(segment_in%grad_gradient)) & + call rotate_array(segment_in%grad_gradient, turns, segment%grad_gradient) + if (allocated(segment_in%rx_norm_rad)) & + call rotate_array(segment_in%rx_norm_rad, turns, segment%ry_norm_rad) + if (allocated(segment_in%ry_norm_rad)) & + call rotate_array(segment_in%ry_norm_rad, turns, segment%rx_norm_rad) + if (allocated(segment_in%rx_norm_obl)) & + call rotate_array(segment_in%rx_norm_obl, turns, segment%ry_norm_obl) + if (allocated(segment_in%ry_norm_obl)) & + call rotate_array(segment_in%ry_norm_obl, turns, segment%rx_norm_obl) + if (allocated(segment_in%cff_normal)) & + call rotate_array(segment_in%cff_normal, turns, segment%cff_normal) + if (allocated(segment_in%nudged_normal_vel)) & + call rotate_array(segment_in%nudged_normal_vel, turns, segment%nudged_normal_vel) + if (allocated(segment_in%nudged_tangential_vel)) & + call rotate_array(segment_in%nudged_tangential_vel, turns, segment%nudged_tangential_vel) + if (allocated(segment_in%nudged_tangential_grad)) & + call rotate_array(segment_in%nudged_tangential_grad, turns, segment%nudged_tangential_grad) + if (associated(segment_in%tr_Reg)) then + do n = 1, segment_in%tr_Reg%ntseg + call rotate_array(segment_in%tr_Reg%tr(n)%tres, turns, segment%tr_Reg%tr(n)%tres) + ! Testing this to see if it works for contant tres values. Probably wrong otherwise. + segment%tr_Reg%Tr(n)%is_initialized=.true. + enddo + endif + + ! Only doing this here for 1 turns, get flow directed in the opposite direction. + ! Currents that are now E_W were N_S and the turns should change their sign. + if (turns == 1 .and. segment%is_E_or_W) then + if (allocated(segment_in%normal_vel)) & + segment%normal_vel(:,:,:) = - segment%normal_vel(:,:,:) + if (allocated(segment_in%normal_trans)) & + segment%normal_trans(:,:,:) = - segment%normal_trans(:,:,:) + if (allocated(segment_in%normal_vel_bt)) & + segment%normal_vel_bt(:,:) = - segment%normal_vel_bt(:,:) + if (allocated(segment_in%tangential_vel)) & + segment%tangential_vel(:,:,:) = - segment%tangential_vel(:,:,:) + if (allocated(segment_in%tangential_grad)) & + segment%tangential_grad(:,:,:) = - segment%tangential_grad(:,:,:) + if (allocated(segment_in%grad_normal)) & + segment%grad_normal(:,:,:) = - segment%grad_normal(:,:,:) + if (allocated(segment_in%grad_tan)) & + segment%grad_tan(:,:,:) = - segment%grad_tan(:,:,:) + if (allocated(segment_in%grad_gradient)) & + segment%grad_gradient(:,:,:) = - segment%grad_gradient(:,:,:) + if (allocated(segment%rx_norm_rad)) & + segment%rx_norm_rad(:,:,:) = - segment%rx_norm_rad(:,:,:) + if (allocated(segment%ry_norm_rad)) & + segment%ry_norm_rad(:,:,:) = - segment%ry_norm_rad(:,:,:) + if (allocated(segment%rx_norm_obl)) & + segment%rx_norm_obl(:,:,:) = - segment%rx_norm_obl(:,:,:) + if (allocated(segment%ry_norm_obl)) & + segment%ry_norm_obl(:,:,:) = - segment%ry_norm_obl(:,:,:) + if (allocated(segment_in%cff_normal)) & + segment%cff_normal(:,:,:) = - segment%cff_normal(:,:,:) + if (allocated(segment_in%nudged_normal_vel)) & + segment%nudged_normal_vel(:,:,:) = - segment%nudged_normal_vel(:,:,:) + if (allocated(segment_in%nudged_tangential_vel)) & + segment%nudged_tangential_vel(:,:,:) = - segment%nudged_tangential_vel(:,:,:) + if (allocated(segment_in%nudged_tangential_grad)) & + segment%nudged_tangential_grad(:,:,:) = - segment%nudged_tangential_grad(:,:,:) + endif + segment%temp_segment_data_exists = segment_in%temp_segment_data_exists segment%salt_segment_data_exists = segment_in%salt_segment_data_exists end subroutine rotate_OBC_segment_data diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 677fdfe6dc..abce27909b 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1616,10 +1616,13 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, & "If true, use a more robust estimate of the first mode wave speed as the "//& "starting point for iterations.", default=.true.) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) + call get_param(param_file, mdl, "INTWAVE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//& "See REMAPPING_USE_OM4_SUBCELLS for details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & "This sets the default value for the various _ANSWER_DATE parameters.", & default=99991231) diff --git a/src/diagnostics/MOM_harmonic_analysis.F90 b/src/diagnostics/MOM_harmonic_analysis.F90 index 1e3b9895cb..0013a8ab83 100644 --- a/src/diagnostics/MOM_harmonic_analysis.F90 +++ b/src/diagnostics/MOM_harmonic_analysis.F90 @@ -173,6 +173,7 @@ end subroutine HA_register !> This subroutine accumulates the temporal basis functions in FtF. !! The tidal constituents are those used in MOM_tidal_forcing, plus the mean (of zero frequency). +!! Only the main diagonal and entries below it are calculated, which are needed for Cholesky decomposition. subroutine HA_accum_FtF(Time, CS) type(time_type), intent(in) :: Time !< The current model time type(harmonic_analysis_CS), intent(inout) :: CS !< Control structure of the MOM_harmonic_analysis module @@ -191,27 +192,31 @@ subroutine HA_accum_FtF(Time, CS) nc = CS%nc now = CS%US%s_to_T * time_type_to_real(Time - CS%time_ref) - ! Accumulate FtF - CS%FtF(1,1) = CS%FtF(1,1) + 1.0 !< For the zero frequency + !< First entry, corresponding to the zero frequency constituent (mean) + CS%FtF(1,1) = CS%FtF(1,1) + 1.0 + do c=1,nc icos = 2*c isin = 2*c+1 cosomegat = cos(CS%freq(c) * now + CS%phase0(c)) sinomegat = sin(CS%freq(c) * now + CS%phase0(c)) + + ! First column, corresponding to the zero frequency constituent (mean) CS%FtF(icos,1) = CS%FtF(icos,1) + cosomegat CS%FtF(isin,1) = CS%FtF(isin,1) + sinomegat - CS%FtF(1,icos) = CS%FtF(icos,1) - CS%FtF(1,isin) = CS%FtF(isin,1) - do cc=c,nc + + do cc=1,c iccos = 2*cc issin = 2*cc+1 ccosomegat = cos(CS%freq(cc) * now + CS%phase0(cc)) ssinomegat = sin(CS%freq(cc) * now + CS%phase0(cc)) + + ! Interior of the matrix, corresponding to the products of cosine and sine terms CS%FtF(icos,iccos) = CS%FtF(icos,iccos) + cosomegat * ccosomegat CS%FtF(icos,issin) = CS%FtF(icos,issin) + cosomegat * ssinomegat CS%FtF(isin,iccos) = CS%FtF(isin,iccos) + sinomegat * ccosomegat CS%FtF(isin,issin) = CS%FtF(isin,issin) + sinomegat * ssinomegat - enddo ! cc=c,nc + enddo ! cc=1,c enddo ! c=1,nc end subroutine HA_accum_FtF @@ -276,14 +281,18 @@ subroutine HA_accum_FtSSH(key, data, Time, G, CS) is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je - ! Accumulate FtF and FtSSH + !< First entry, corresponding to the zero frequency constituent (mean) + do j=js,je ; do i=is,ie + ha1%FtSSH(i,j,1) = ha1%FtSSH(i,j,1) + (data(i,j) - ha1%ref(i,j)) + enddo ; enddo + + !< The remaining entries do c=1,nc icos = 2*c isin = 2*c+1 cosomegat = cos(CS%freq(c) * now + CS%phase0(c)) sinomegat = sin(CS%freq(c) * now + CS%phase0(c)) do j=js,je ; do i=is,ie - ha1%FtSSH(i,j,1) = ha1%FtSSH(i,j,1) + (data(i,j) - ha1%ref(i,j)) ha1%FtSSH(i,j,icos) = ha1%FtSSH(i,j,icos) + (data(i,j) - ha1%ref(i,j)) * cosomegat ha1%FtSSH(i,j,isin) = ha1%FtSSH(i,j,isin) + (data(i,j) - ha1%ref(i,j)) * sinomegat enddo ; enddo @@ -315,7 +324,7 @@ subroutine HA_write(ha1, Time, G, CS) ! Local variables real, dimension(:,:,:), allocatable :: FtSSHw !< An array containing the harmonic constants [A] integer :: year, month, day, hour, minute, second - integer :: nc, k, is, ie, js, je + integer :: nc, i, j, k, is, ie, js, je character(len=255) :: filename !< Output file name type(MOM_infra_file) :: cdf !< The file handle for output harmonic constants @@ -348,6 +357,11 @@ subroutine HA_write(ha1, Time, G, CS) call create_MOM_file(cdf, trim(filename), cdf_vars, & 2*nc+1, cdf_fields, SINGLE_FILE, 86400.0, G=G) + ! Add the initial field back to the mean state + do j=js,je ; do i=is,ie + FtSSHw(i,j,1) = FtSSHw(i,j,1) + ha1%ref(i,j) + enddo ; enddo + ! Write data call MOM_write_field(cdf, cdf_fields(1), G%domain, FtSSHw(:,:,1), 0.0) do k=1,nc @@ -362,75 +376,68 @@ subroutine HA_write(ha1, Time, G, CS) end subroutine HA_write -!> This subroutine computes the harmonic constants (stored in FtSSHw) using the dot products of the temporal +!> This subroutine computes the harmonic constants (stored in x) using the dot products of the temporal !! basis functions accumulated in FtF, and the dot products of the SSH (or other fields) with the temporal basis !! functions accumulated in FtSSH. The system is solved by Cholesky decomposition, !! -!! FtF * FtSSHw = FtSSH, => FtFw * (FtFw' * FtSSHw) = FtSSH, +!! FtF * x = FtSSH, => L * (L' * x) = FtSSH, => L * y = FtSSH, !! -!! where FtFw is a lower triangular matrix, and the prime denotes matrix transpose. +!! where L is the lower triangular matrix, y = L' * x, and x is the solution vector. !! -subroutine HA_solver(ha1, nc, FtF, FtSSHw) +subroutine HA_solver(ha1, nc, FtF, x) type(HA_type), pointer, intent(in) :: ha1 !< Control structure for the current field integer, intent(in) :: nc !< Number of harmonic constituents real, dimension(:,:), intent(in) :: FtF !< Accumulator of (F' * F) for all fields [nondim] - real, dimension(:,:,:), allocatable, intent(out) :: FtSSHw !< Work array for Cholesky decomposition [A] + real, dimension(ha1%is:ha1%ie,ha1%js:ha1%je,2*nc+1), & + intent(out) :: x !< Solution vector of harmonic constants [A] ! Local variables - real :: tmp0 !< Temporary variable for Cholesky decomposition [nondim] - real, dimension(:), allocatable :: tmp1 !< Temporary variable for Cholesky decomposition [nondim] - real, dimension(:,:), allocatable :: tmp2 !< Temporary variable for Cholesky decomposition [A] - real, dimension(:,:), allocatable :: FtFw !< Lower triangular matrix for Cholesky decomposition [nondim] - integer :: k, m, n, is, ie, js, je - - is = ha1%is ; ie = ha1%ie ; js = ha1%js ; je = ha1%je - - allocate(tmp1(1:2*nc+1), source=0.0) - allocate(tmp2(is:ie,js:je), source=0.0) - allocate(FtFw(1:2*nc+1,1:2*nc+1), source=0.0) - allocate(FtSSHw(is:ie,js:je,2*nc+1), source=0.0) - - ! Construct FtFw - FtFw(:,:) = 0.0 + real :: tmp0 !< Temporary variable for Cholesky decomposition [nondim] + real, dimension(2*nc+1,2*nc+1) :: L !< Lower triangular matrix of Cholesky decomposition [nondim] + real, dimension(2*nc+1) :: tmp1 !< Inverse of the diagonal entries of L [nondim] + real, dimension(ha1%is:ha1%ie,ha1%js:ha1%je) :: tmp2 !< 2D temporary array involving FtSSH [A] + real, dimension(ha1%is:ha1%ie,ha1%js:ha1%je,2*nc+1) :: y !< 3D temporary array, i.e., L' * x [A] + integer :: k, m, n + + ! Cholesky decomposition do m=1,2*nc+1 + + ! First, calculate the diagonal entries tmp0 = 0.0 - do k=1,m-1 - tmp0 = tmp0 + FtFw(m,k) * FtFw(m,k) + do k=1,m-1 ! This loop operates along the m-th row + tmp0 = tmp0 + L(m,k) * L(m,k) enddo - FtFw(m,m) = sqrt(FtF(m,m) - tmp0) - tmp1(m) = 1 / FtFw(m,m) - do k=m+1,2*nc+1 + L(m,m) = sqrt(FtF(m,m) - tmp0) ! This is the m-th diagonal entry + + ! Now calculate the off-diagonal entries + tmp1(m) = 1 / L(m,m) + do k=m+1,2*nc+1 ! This loop operates along the column below the m-th diagonal entry tmp0 = 0.0 do n=1,m-1 - tmp0 = tmp0 + FtFw(k,n) * FtFw(m,n) + tmp0 = tmp0 + L(k,n) * L(m,n) enddo - FtFw(k,m) = (FtF(k,m) - tmp0) * tmp1(m) + L(k,m) = (FtF(k,m) - tmp0) * tmp1(m) ! This is the k-th off-diagonal entry below the m-th diagonal entry enddo enddo - ! Solve for (FtFw' * FtSSHw) - FtSSHw(:,:,:) = ha1%FtSSH(:,:,:) + ! Solve for y from L * y = FtSSH do k=1,2*nc+1 tmp2(:,:) = 0.0 do m=1,k-1 - tmp2(:,:) = tmp2(:,:) + FtFw(k,m) * FtSSHw(:,:,m) + tmp2(:,:) = tmp2(:,:) + L(k,m) * y(:,:,m) enddo - FtSSHw(:,:,k) = (FtSSHw(:,:,k) - tmp2(:,:)) * tmp1(k) + y(:,:,k) = (ha1%FtSSH(:,:,k) - tmp2(:,:)) * tmp1(k) enddo - ! Solve for FtSSHw + ! Solve for x from L' * x = y do k=2*nc+1,1,-1 tmp2(:,:) = 0.0 do m=k+1,2*nc+1 - tmp2(:,:) = tmp2(:,:) + FtSSHw(:,:,m) * FtFw(m,k) + tmp2(:,:) = tmp2(:,:) + L(m,k) * x(:,:,m) enddo - FtSSHw(:,:,k) = (FtSSHw(:,:,k) - tmp2(:,:)) * tmp1(k) + x(:,:,k) = (y(:,:,k) - tmp2(:,:)) * tmp1(k) enddo - deallocate(tmp1) - deallocate(tmp2) - deallocate(FtFw) - end subroutine HA_solver !> \namespace harmonic_analysis @@ -441,7 +448,7 @@ end subroutine HA_solver !! step, and x is a 2*nc-by-1 vector containing the constant coefficients of the sine/cosine for each constituent !! (i.e., the harmonic constants). At each grid point, the harmonic constants are computed using least squares, !! -!! (F' * F) * x = F' * SSH_in, +!! (F' * F) * x = F' * SSH_in, => FtF * x = FtSSH, !! !! where the prime denotes matrix transpose, and SSH_in is the sea surface height (or other fields) determined by !! the model. The dot products (F' * F) and (F' * SSH_in) are computed by accumulating the sums as the model is diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 398241b98c..f5ff19630b 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -21,9 +21,9 @@ module MOM_sum_output use MOM_io, only : attribute_info, set_attribute_info, delete_attribute_info use MOM_io, only : APPEND_FILE, SINGLE_FILE, WRITEONLY_FILE use MOM_spatial_means, only : array_global_min_max -use MOM_time_manager, only : time_type, get_time, get_date, set_time, operator(>) +use MOM_time_manager, only : time_type, get_time, get_date, set_time use MOM_time_manager, only : operator(+), operator(-), operator(*), operator(/) -use MOM_time_manager, only : operator(/=), operator(<=), operator(>=), operator(<) +use MOM_time_manager, only : operator(/=), operator(<=), operator(>=), operator(<), operator(>) use MOM_time_manager, only : get_calendar_type, time_type_to_real, NO_CALENDAR use MOM_tracer_flow_control, only : tracer_flow_control_CS, call_tracer_stocks use MOM_unit_scaling, only : unit_scale_type @@ -489,7 +489,7 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci CS%write_energy_time = CS%Start_time + CS%energysavedays * & (1 + (day - CS%Start_time) / CS%energysavedays) endif - elseif (day + (dt_force/2) <= CS%write_energy_time) then + elseif (day + (dt_force/2) < CS%write_energy_time) then return ! Do not write this step else ! Determine the next write time before proceeding if (CS%energysave_geometric) then diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index c28e2e5896..3bb73e4c57 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -3280,10 +3280,12 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & "This sets the default value for the various _ANSWER_DATE parameters.", & default=99991231) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(param_file, mdl, "DIAG_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for diagnostics. "//& "See REMAPPING_USE_OM4_SUBCELLS for details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, & "The vintage of the expressions and order of arithmetic to use for remapping. "//& "Values below 20190101 result in the use of older, less accurate expressions "//& diff --git a/src/framework/MOM_document.F90 b/src/framework/MOM_document.F90 index f32573815f..eceb87d7d4 100644 --- a/src/framework/MOM_document.F90 +++ b/src/framework/MOM_document.F90 @@ -303,14 +303,16 @@ subroutine doc_param_real(doc, varname, desc, units, val, default, debuggingPara end subroutine doc_param_real !> This subroutine handles parameter documentation for arrays of reals. -subroutine doc_param_real_array(doc, varname, desc, units, vals, default, debuggingParam, like_default) +subroutine doc_param_real_array(doc, varname, desc, units, vals, default, defaults, & + debuggingParam, like_default) type(doc_type), pointer :: doc !< A pointer to a structure that controls where the !! documentation occurs and its formatting character(len=*), intent(in) :: varname !< The name of the parameter being documented character(len=*), intent(in) :: desc !< A description of the parameter being documented character(len=*), intent(in) :: units !< The units of the parameter being documented real, intent(in) :: vals(:) !< The array of values to record - real, optional, intent(in) :: default !< The default value of this parameter + real, optional, intent(in) :: default !< A uniform default value of this parameter + real, optional, intent(in) :: defaults(:) !< The element-wise default values of this parameter logical, optional, intent(in) :: debuggingParam !< If present and true, this is a debugging parameter. logical, optional, intent(in) :: like_default !< If present and true, log this parameter as though !! it has the default value, even if there is no default. @@ -334,6 +336,11 @@ subroutine doc_param_real_array(doc, varname, desc, units, vals, default, debugg do i=1,size(vals) ; if (vals(i) /= default) equalsDefault = .false. ; enddo mesg = trim(mesg)//" default = "//trim(real_string(default)) endif + if (present(defaults)) then + equalsDefault = .true. + do i=1,size(vals) ; if (vals(i) /= defaults(i)) equalsDefault = .false. ; enddo + mesg = trim(mesg)//" default = "//trim(real_array_string(defaults)) + endif if (present(like_default)) then ; if (like_default) equalsDefault = .true. ; endif if (mesgHasBeenDocumented(doc, varName, mesg)) return ! Avoid duplicates diff --git a/src/framework/MOM_file_parser.F90 b/src/framework/MOM_file_parser.F90 index 22d3789ea5..fc496ac1b5 100644 --- a/src/framework/MOM_file_parser.F90 +++ b/src/framework/MOM_file_parser.F90 @@ -1464,7 +1464,7 @@ end subroutine log_param_real !> Log the name and values of an array of real model parameter in documentation files. subroutine log_param_real_array(CS, modulename, varname, value, desc, & - units, default, debuggingParam, like_default, unscale) + units, default, defaults, debuggingParam, like_default, unscale) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters character(len=*), intent(in) :: modulename !< The name of the calling module @@ -1473,7 +1473,8 @@ subroutine log_param_real_array(CS, modulename, varname, value, desc, & character(len=*), optional, intent(in) :: desc !< A description of this variable; if not !! present, this parameter is not written to a doc file character(len=*), intent(in) :: units !< The units of this parameter - real, optional, intent(in) :: default !< The default value of the parameter + real, optional, intent(in) :: default !< A uniform default value of the parameter + real, optional, intent(in) :: defaults(:) !< The element-wise defaults of the parameter logical, optional, intent(in) :: debuggingParam !< If present and true, this parameter is !! logged in the debugging parameter file logical, optional, intent(in) :: like_default !< If present and true, log this parameter as @@ -1498,7 +1499,7 @@ subroutine log_param_real_array(CS, modulename, varname, value, desc, & write(myunits(1:240),'(A)') trim(units) if (present(desc)) & - call doc_param(CS%doc, varname, desc, myunits, log_val, default, & + call doc_param(CS%doc, varname, desc, myunits, log_val, default, defaults, & debuggingParam=debuggingParam, like_default=like_default) end subroutine log_param_real_array @@ -1835,7 +1836,7 @@ end subroutine get_param_real !> This subroutine reads the values of an array of real model parameters from a parameter file !! and logs them in documentation files. subroutine get_param_real_array(CS, modulename, varname, value, desc, units, & - default, fail_if_missing, do_not_read, do_not_log, debuggingParam, & + default, defaults, fail_if_missing, do_not_read, do_not_log, debuggingParam, & scale, unscaled) type(param_file_type), intent(in) :: CS !< The control structure for the file_parser module, !! it is also a structure to parse for run-time parameters @@ -1846,7 +1847,8 @@ subroutine get_param_real_array(CS, modulename, varname, value, desc, units, & character(len=*), optional, intent(in) :: desc !< A description of this variable; if not !! present, this parameter is not written to a doc file character(len=*), intent(in) :: units !< The units of this parameter - real, optional, intent(in) :: default !< The default value of the parameter + real, optional, intent(in) :: default !< A uniform default value of the parameter + real, optional, intent(in) :: defaults(:) !< The element-wise defaults of the parameter logical, optional, intent(in) :: fail_if_missing !< If present and true, a fatal error occurs !! if this variable is not found in the parameter file logical, optional, intent(in) :: do_not_read !< If present and true, do not read a @@ -1865,14 +1867,22 @@ subroutine get_param_real_array(CS, modulename, varname, value, desc, units, & do_read = .true. ; if (present(do_not_read)) do_read = .not.do_not_read do_log = .true. ; if (present(do_not_log)) do_log = .not.do_not_log + if (present(defaults)) then + if (present(default)) call MOM_error(FATAL, & + "get_param_real_array: Only one of default and defaults can be specified at a time.") + if (size(defaults) /= size(value)) call MOM_error(FATAL, & + "get_param_real_array: The size of defaults nad value are not the same.") + endif + if (do_read) then if (present(default)) value(:) = default + if (present(defaults)) value(:) = defaults(:) call read_param_real_array(CS, varname, value, fail_if_missing) endif if (do_log) then call log_param_real_array(CS, modulename, varname, value, desc, & - units, default, debuggingParam) + units, default, defaults, debuggingParam) endif if (present(unscaled)) unscaled(:) = value(:) diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index edf08da3aa..eabd376512 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -1327,10 +1327,9 @@ subroutine compute_global_grid_integrals(G, US) G%areaT_global = reproducing_sum(tmpForSumming) if (G%areaT_global == 0.0) & - call MOM_error(FATAL, "compute_global_grid_integrals: "//& - "zero ocean area (check topography?)") + call MOM_error(FATAL, "compute_global_grid_integrals: zero ocean area (check topography?)") - G%IareaT_global = 1.0 / (G%areaT_global) + G%IareaT_global = 1.0 / G%areaT_global end subroutine compute_global_grid_integrals ! ----------------------------------------------------------------------------- diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 769d60d51d..41b407d6a1 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -278,6 +278,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & " \t uniform - uniform thickness layers evenly distributed \n"//& " \t\t between the surface and MAXIMUM_DEPTH. \n"//& " \t list - read a list of positive interface depths. \n"//& + " \t param - use thicknesses from parameter THICKNESS_INIT_VALUES. \n"//& " \t DOME - use a slope and channel configuration for the \n"//& " \t\t DOME sill-overflow test case. \n"//& " \t ISOMIP - use a configuration for the \n"//& @@ -318,6 +319,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & just_read=just_read) case ("list"); call initialize_thickness_list(dz, depth_tot, G, GV, US, PF, & just_read=just_read) + case ("param"); call initialize_thickness_param(dz, depth_tot, G, GV, US, PF, & + just_read=just_read) case ("DOME"); call DOME_initialize_thickness(dz, depth_tot, G, GV, PF, & just_read=just_read) case ("ISOMIP"); call ISOMIP_initialize_thickness(dz, depth_tot, G, GV, US, PF, tv, & @@ -1011,6 +1014,68 @@ subroutine initialize_thickness_list(h, depth_tot, G, GV, US, param_file, just_r call callTree_leave(trim(mdl)//'()') end subroutine initialize_thickness_list +!> Initializes thickness based on a run-time parameter with nominal thickness +!! for each layer +subroutine initialize_thickness_param(h, depth_tot, G, GV, US, param_file, just_read) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: h !< The thickness that is being initialized [Z ~> m] + real, dimension(SZI_(G),SZJ_(G)), & + intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + logical, intent(in) :: just_read !< If true, this call will only read + !! parameters without changing h. + ! Local variables + character(len=40) :: mdl = "initialize_thickness_param" ! This subroutine's name. + real :: e0(SZK_(GV)+1) ! The resting interface heights [Z ~> m], usually + ! negative because it is positive upward. + real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface, + ! positive upward [Z ~> m]. + real :: dz(SZK_(GV)) ! The nominal initial layer thickness [Z ~> m], usually + real :: h0_def(SZK_(GV)) ! Uniform default values for dz [Z ~> m], usually + integer :: i, j, k, is, ie, js, je, nz + + call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90") + if (G%max_depth<=0.) call MOM_error(FATAL, "initialize_thickness_param: "// & + "MAXIMUM_DEPTH has a nonsensical value! Was it set?") + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + h0_def(:) = ( G%max_depth / real(nz) ) * US%Z_to_m + call get_param(param_file, mdl, "THICKNESS_INIT_VALUES", dz, & + "A list of nominal thickness for each layer to initialize with", & + units="m", scale=US%m_to_Z, defaults=h0_def, do_not_log=just_read) + if (just_read) return ! This subroutine has no run-time parameters. + + e0(nz+1) = -G%max_depth + do k=nz, 1, -1 + e0(K) = e0(K+1) + dz(k) + enddo + + do j=js,je ; do i=is,ie + ! This sets the initial thickness (in m) of the layers. The + ! thicknesses are set to insure that: 1. each layer is at least an + ! Angstrom thick, and 2. the interfaces are where they should be + ! based on the resting depths and interface height perturbations, + ! as long at this doesn't interfere with 1. + eta1D(nz+1) = -depth_tot(i,j) + do k=nz,1,-1 + eta1D(K) = e0(K) + if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then + eta1D(K) = eta1D(K+1) + GV%Angstrom_Z + h(i,j,k) = GV%Angstrom_Z + else + h(i,j,k) = eta1D(K) - eta1D(K+1) + endif + enddo + enddo ; enddo + + call callTree_leave(trim(mdl)//'()') +end subroutine initialize_thickness_param + !> Search density space for location of layers (not implemented!) subroutine initialize_thickness_search call MOM_error(FATAL," MOM_state_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED") @@ -2339,26 +2404,6 @@ subroutine set_velocity_depth_max(G) enddo ; enddo end subroutine set_velocity_depth_max -!> Subroutine to pre-compute global integrals of grid quantities for -!! later use in reporting diagnostics -subroutine compute_global_grid_integrals(G, US) - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - ! Local variables - real, dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming ! Masked and unscaled areas for sums [m2] - real :: area_scale ! A conversion factor to prepare for reproducing sums [m2 L-2 ~> 1] - integer :: i,j - - area_scale = US%L_to_m**2 - tmpForSumming(:,:) = 0. - G%areaT_global = 0.0 ; G%IareaT_global = 0.0 - do j=G%jsc,G%jec ; do i=G%isc,G%iec - tmpForSumming(i,j) = area_scale*G%areaT(i,j) * G%mask2dT(i,j) - enddo ; enddo - G%areaT_global = reproducing_sum(tmpForSumming) - G%IareaT_global = 1. / (G%areaT_global) -end subroutine compute_global_grid_integrals - !> This subroutine sets the 4 bottom depths at velocity points to be the !! minimum of the adjacent depths. subroutine set_velocity_depth_min(G) @@ -2577,10 +2622,12 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=just_read.or.(.not.GV%Boussinesq)) + call get_param(PF, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(PF, mdl, "Z_INIT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for initialization. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701) endif call get_param(PF, mdl, "HOR_REGRID_ANSWER_DATE", hor_regrid_answer_date, & diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index 6e3da385ce..615ce07f0d 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -138,10 +138,12 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, US, PF, src_file, src_var_ "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=.not.GV%Boussinesq) + call get_param(PF, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(PF, mdl, "Z_INIT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for initialization. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701) endif call get_param(PF, mdl, "HOR_REGRID_ANSWER_DATE", hor_regrid_answer_date, & diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index c7101ac6b7..899dcbbbf0 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -3610,10 +3610,12 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS) "mode speeds are not calculated but are simply reported as 0. This must be "//& "non-negative for the wave_speeds routine to be used.", & units="m s-1", default=0.01, scale=US%m_s_to_L_T) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(param_file, mdl, "INTWAVE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//& "See REMAPPING_USE_OM4_SUBCELLS for details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) call get_param(param_file, mdl, "UNIFORM_TEST_CG", CS%uniform_test_cg, & "If positive, a uniform group velocity of internal tide for test case", & default=-1., units="m s-1", scale=US%m_s_to_L_T) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 8f388dc263..ac6efe2268 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -1834,10 +1834,12 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, & "If true, use a more robust estimate of the first mode wave speed as the "//& "starting point for iterations.", default=.true.) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(param_file, mdl, "EBT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//& "See REMAPPING_USE_OM4_SUBCELLS for details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) call wave_speed_init(CS%wave_speed, GV, use_ebt_mode=CS%Resoln_use_ebt, & mono_N2_depth=N2_filter_depth, remap_answer_date=remap_answer_date, & better_speed_est=better_speed_est, min_speed=wave_speed_min, & diff --git a/src/parameterizations/lateral/MOM_streaming_filter.F90 b/src/parameterizations/lateral/MOM_streaming_filter.F90 index a91f6661f2..f7a4f736c7 100644 --- a/src/parameterizations/lateral/MOM_streaming_filter.F90 +++ b/src/parameterizations/lateral/MOM_streaming_filter.F90 @@ -15,7 +15,7 @@ module MOM_streaming_filter !> The control structure for storing the filter infomation of a particular field type, public :: Filter_CS ; private real :: a, & !< Parameter that determines the bandwidth [nondim] - om, & !< Target frequency of the filter [T-1 ~> s-1] + om, & !< Target frequency of the filter [rad T-1 ~> rad s-1] old_time = -1.0 !< The time of the previous accumulating step [T ~> s] real, allocatable, dimension(:,:) :: s1, & !< Dummy variable [A] u1 !< Filtered data [A] @@ -29,7 +29,7 @@ module MOM_streaming_filter !> This subroutine registers each of the fields to be filtered. subroutine Filt_register(a, om, grid, HI, CS) real, intent(in) :: a !< Parameter that determines the bandwidth [nondim] - real, intent(in) :: om !< Target frequency of the filter [T-1 ~> s-1] + real, intent(in) :: om !< Target frequency of the filter [rad T-1 ~> rad s-1] character(len=*), intent(in) :: grid !< Horizontal grid location: h, u, or v type(hor_index_type), intent(in) :: HI !< Horizontal index type structure type(Filter_CS), intent(out) :: CS !< Control structure for the current field diff --git a/src/parameterizations/lateral/MOM_tidal_forcing.F90 b/src/parameterizations/lateral/MOM_tidal_forcing.F90 index 177becf84f..43885cccc3 100644 --- a/src/parameterizations/lateral/MOM_tidal_forcing.F90 +++ b/src/parameterizations/lateral/MOM_tidal_forcing.F90 @@ -50,7 +50,7 @@ module MOM_tidal_forcing !! and bottom geopotential anomalies [nondim]. integer :: nc !< The number of tidal constituents in use. real, dimension(MAX_CONSTITUENTS) :: & - freq, & !< The frequency of a tidal constituent [T-1 ~> s-1]. + freq, & !< The frequency of a tidal constituent [rad T-1 ~> rad s-1]. phase0, & !< The phase of a tidal constituent at time 0 [rad]. amp, & !< The amplitude of a tidal constituent at time 0 [Z ~> m]. love_no !< The Love number of a tidal constituent at time 0 [nondim]. @@ -151,7 +151,7 @@ end function eq_phase !! Values used here are from previous versions of MOM. function tidal_frequency(constit) character (len=2), intent(in) :: constit !> Constituent to look up - real :: tidal_frequency !> Angular frequency [s-1] + real :: tidal_frequency !> Angular frequency [rad s-1] select case (constit) case ("M2") @@ -246,7 +246,7 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS, HA_CS) phase, & ! The phase of some tidal constituent [radians]. lat_rad, lon_rad ! Latitudes and longitudes of h-points [radians]. real :: deg_to_rad ! A conversion factor from degrees to radians [radian degree-1] - real, dimension(MAX_CONSTITUENTS) :: freq_def ! Default frequency for each tidal constituent [s-1] + real, dimension(MAX_CONSTITUENTS) :: freq_def ! Default frequency for each tidal constituent [rad s-1] real, dimension(MAX_CONSTITUENTS) :: phase0_def ! Default reference phase for each tidal constituent [rad] real, dimension(MAX_CONSTITUENTS) :: amp_def ! Default amplitude for each tidal constituent [m] real, dimension(MAX_CONSTITUENTS) :: love_def ! Default love number for each constituent [nondim] @@ -480,7 +480,8 @@ subroutine tidal_forcing_init(Time, G, US, param_file, CS, HA_CS) "Frequency of the "//trim(CS%const_name(c))//" tidal constituent. "//& "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// & " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and "//trim(CS%const_name(c))// & - " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=freq_def(c), scale=US%T_to_s) + " is in OBC_TIDE_CONSTITUENTS.", units="rad s-1", default=freq_def(c), & + scale=US%T_to_s) call get_param(param_file, mdl, "TIDE_"//trim(CS%const_name(c))//"_AMP", CS%amp(c), & "Amplitude of the "//trim(CS%const_name(c))//" tidal constituent. "//& "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// & diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 0dfead633c..c00c72ea3c 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -236,10 +236,12 @@ subroutine initialize_ALE_sponge_fixed(Iresttime, G, GV, param_file, CS, data_h, "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=.not.GV%Boussinesq) if (.not.GV%Boussinesq) CS%remap_answer_date = max(CS%remap_answer_date, 20230701) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(param_file, mdl, "SPONGE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for ALE sponge. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) call get_param(param_file, mdl, "HOR_REGRID_ANSWER_DATE", CS%hor_regrid_answer_date, & "The vintage of the order of arithmetic for horizontal regridding. "//& @@ -520,10 +522,12 @@ subroutine initialize_ALE_sponge_varying(Iresttime, G, GV, US, param_file, CS, I "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(param_file, mdl, "SPONGE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for ALE sponge. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) call get_param(param_file, mdl, "HOR_REGRID_ANSWER_DATE", CS%hor_regrid_answer_date, & "The vintage of the order of arithmetic for horizontal regridding. "//& "Dates before 20190101 give the same answers as the code did in late 2018, "//& diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index c5297a3cf0..e984d5831c 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -3256,13 +3256,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di 'Mixed layer depth (delta rho = 0.125)', 'm', conversion=US%Z_to_m) call get_param(param_file, mdl, "MLD_EN_VALS", CS%MLD_En_vals, & "The energy values used to compute MLDs. If not set (or all set to 0.), the "//& - "default will overwrite to 25., 2500., 250000.", & - units='J/m2', default=0., scale=US%W_m2_to_RZ3_T3*US%s_to_T) - if ((CS%MLD_En_vals(1)==0.).and.(CS%MLD_En_vals(2)==0.).and.(CS%MLD_En_vals(3)==0.)) then - CS%MLD_En_vals = (/ 25.*US%W_m2_to_RZ3_T3*US%s_to_T, & - 2500.*US%W_m2_to_RZ3_T3*US%s_to_T, & - 250000.*US%W_m2_to_RZ3_T3*US%s_to_T /) - endif + "default will overwrite to 25., 2500., 250000.", units='J/m2', & + defaults=(/25., 2500., 250000./), scale=US%W_m2_to_RZ3_T3*US%s_to_T) write(EN1,'(F10.2)') CS%MLD_En_vals(1)*US%RZ3_T3_to_W_m2*US%T_to_s write(EN2,'(F10.2)') CS%MLD_En_vals(2)*US%RZ3_T3_to_W_m2*US%T_to_s write(EN3,'(F10.2)') CS%MLD_En_vals(3)*US%RZ3_T3_to_W_m2*US%T_to_s diff --git a/src/tracer/MOM_hor_bnd_diffusion.F90 b/src/tracer/MOM_hor_bnd_diffusion.F90 index 6d8fe881d1..e2718590fb 100644 --- a/src/tracer/MOM_hor_bnd_diffusion.F90 +++ b/src/tracer/MOM_hor_bnd_diffusion.F90 @@ -143,10 +143,13 @@ logical function hor_bnd_diffusion_init(Time, G, GV, US, param_file, diag, diaba "for vertical remapping for all variables. "//& "It can be one of the following schemes: "//& trim(remappingSchemesDoc), default=remappingDefaultScheme) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) + call get_param(param_file, mdl, "HBD_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for horizontal boundary diffusion. "//& "See REMAPPING_USE_OM4_SUBCELLS for details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) ! GMM, TODO: add HBD params to control optional arguments in initialize_remapping. call initialize_remapping( CS%remap_CS, string, boundary_extrapolation=boundary_extrap, & diff --git a/src/tracer/MOM_neutral_diffusion.F90 b/src/tracer/MOM_neutral_diffusion.F90 index 1aaf7409d2..feb5dde247 100644 --- a/src/tracer/MOM_neutral_diffusion.F90 +++ b/src/tracer/MOM_neutral_diffusion.F90 @@ -238,10 +238,12 @@ logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, "that were in use at the end of 2018. Higher values result in the use of more "//& "robust and accurate forms of mathematically equivalent expressions.", & default=default_answer_date, do_not_log=.not.GV%Boussinesq) + call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & + do_not_log=.true., default=.true.) call get_param(param_file, mdl, "NDIFF_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, & "If true, use the OM4 remapping-via-subcells algorithm for neutral diffusion. "//& "See REMAPPING_USE_OM4_SUBCELLS for more details. "//& - "We recommend setting this option to false.", default=.true.) + "We recommend setting this option to false.", default=om4_remap_via_sub_cells) if (.not.GV%Boussinesq) CS%remap_answer_date = max(CS%remap_answer_date, 20230701) call initialize_remapping( CS%remap_CS, string, boundary_extrapolation=boundary_extrap, & om4_remap_via_sub_cells=om4_remap_via_sub_cells, & diff --git a/src/tracer/MOM_tracer_advect.F90 b/src/tracer/MOM_tracer_advect.F90 index 0e129b2d03..e11ee6e5a5 100644 --- a/src/tracer/MOM_tracer_advect.F90 +++ b/src/tracer/MOM_tracer_advect.F90 @@ -648,6 +648,19 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, & endif enddo + ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only) + if (associated(OBC)) then ; if (OBC%OBC_pe) then + if (OBC%specified_u_BCs_exist_globally .or. OBC%open_u_BCs_exist_globally) then + do i=is,ie-1 ; if (OBC%segnum_u(I,j) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then + do_i(i+1,j) = .false. + elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + do_i(i,j) = .false. + endif + endif ; enddo + endif + endif ; endif + ! update tracer concentration from i-flux and save some diagnostics do m=1,ntr @@ -1039,6 +1052,24 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, & else ; do_i(i,j) = .false. ; endif enddo + ! Update do_i so that nothing changes outside of the OBC (problem for interior OBCs only) + if (associated(OBC)) then ; if (OBC%OBC_pe) then + if (OBC%specified_v_BCs_exist_globally .or. OBC%open_v_BCs_exist_globally) then + do i=is,ie + if (OBC%segnum_v(i,J-1) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J-1))%direction == OBC_DIRECTION_N) then + do_i(i,j) = .false. + endif + endif + if (OBC%segnum_v(i,J) /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + do_i(i,j) = .false. + endif + endif + enddo + endif + endif ; endif + ! update tracer and save some diagnostics do m=1,ntr do i=is,ie ; if (do_i(i,j)) then diff --git a/src/user/user_change_diffusivity.F90 b/src/user/user_change_diffusivity.F90 index 9a56c12b9c..1a1881a42b 100644 --- a/src/user/user_change_diffusivity.F90 +++ b/src/user/user_change_diffusivity.F90 @@ -230,14 +230,15 @@ subroutine user_change_diff_init(Time, G, GV, US, param_file, diag, CS) "applied. The four values specify the latitudes at "//& "which the extra diffusivity starts to increase from 0, "//& "hits its full value, starts to decrease again, and is "//& - "back to 0.", units="degrees_N", default=-1.0e9) + "back to 0.", units="degrees_N", defaults=(/-1.0e9,-1.0e9,-1.0e9,-1.0e9/)) call get_param(param_file, mdl, "USER_KD_ADD_RHO_RANGE", CS%rho_range(:), & "Four successive values that define a range of potential "//& "densities over which the user-given extra diffusivity "//& "is applied. The four values specify the density at "//& "which the extra diffusivity starts to increase from 0, "//& "hits its full value, starts to decrease again, and is "//& - "back to 0.", units="kg m-3", default=-1.0e9, scale=US%kg_m3_to_R) + "back to 0.", units="kg m-3", defaults=(/-1.0e9,-1.0e9,-1.0e9,-1.0e9/),& + scale=US%kg_m3_to_R) call get_param(param_file, mdl, "USER_KD_ADD_USE_ABS_LAT", CS%use_abs_lat, & "If true, use the absolute value of latitude when "//& "checking whether a point fits into range of latitudes.", &