diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index df4f16409c..b632487cdf 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -719,7 +719,7 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n, u, v, OBC, Reg, dt, dzReg ! we have to keep track of the total dzInterface if for some reason ! we're using the old remapping algorithm for u/v real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzInterface, dzIntTotal - real :: h_neglect, h_neglect_edge + real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2] nz = GV%ke diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index e46bfc7c37..5bc52f9fe3 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -1271,11 +1271,11 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (Htot_avg*CS%dx_Cv(i,J) <= 0.0) then CS%IDatv(i,J) = 0.0 elseif (integral_BT_cont) then - CS%IDatv(i,J) = CS%dx_Cv(i,J) / (max(find_dvhbt_dv(vbt(i,J)*dt, BTCL_v(i,J)), & - CS%dx_Cv(i,J)*Htot_avg) ) * GV%Z_to_H + CS%IDatv(i,J) = GV%Z_to_H * CS%dx_Cv(i,J) / (max(find_dvhbt_dv(vbt(i,J)*dt, BTCL_v(i,J)), & + CS%dx_Cv(i,J)*Htot_avg) ) elseif (use_BT_cont) then ! Reconsider the max and whether there should be some scaling. - CS%IDatv(i,J) = CS%dx_Cv(i,J) / (max(find_dvhbt_dv(vbt(i,J), BTCL_v(i,J)), & - CS%dx_Cv(i,J)*Htot_avg) ) * GV%Z_to_H + CS%IDatv(i,J) = GV%Z_to_H * CS%dx_Cv(i,J) / (max(find_dvhbt_dv(vbt(i,J), BTCL_v(i,J)), & + CS%dx_Cv(i,J)*Htot_avg) ) else CS%IDatv(i,J) = GV%Z_to_H / Htot_avg endif diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 47061f9447..9d87d0bc5d 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -310,6 +310,9 @@ module MOM_open_boundary real :: ramp_value !< If ramp is True, where we are on the ramp from !! zero to one [nondim]. type(time_type) :: ramp_start_time !< Time when model was started. + logical :: answers_2018 !< If true, use the order of arithmetic and expressions for remapping + !! that recover the answers from the end of 2018. Otherwise, use more + !! robust and accurate forms of mathematically equivalent expressions. end type ocean_OBC_type !> Control structure for open boundaries that read from files. @@ -358,9 +361,9 @@ subroutine open_boundary_config(G, US, param_file, OBC) character(len=200) :: config1 ! String for OBC_USER_CONFIG real :: Lscale_in, Lscale_out ! parameters controlling tracer values at the boundaries [L ~> m] character(len=128) :: inputdir - logical :: answers_2018, default_2018_answers logical :: check_reconstruction, check_remapping, force_bounds_in_subcell character(len=64) :: remappingScheme + logical :: default_2018_answers ! This include declares and sets the variable "version". # include "version_variable.h" @@ -608,7 +611,7 @@ subroutine open_boundary_config(G, US, param_file, OBC) call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & "This sets the default value for the various _2018_ANSWERS parameters.", & default=.false.) - call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", answers_2018, & + call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", OBC%answers_2018, & "If true, use the order of arithmetic and expressions that recover the "//& "answers from the end of 2018. Otherwise, use updated and more robust "//& "forms of the same expressions.", default=default_2018_answers) @@ -616,7 +619,7 @@ subroutine open_boundary_config(G, US, param_file, OBC) allocate(OBC%remap_CS) call initialize_remapping(OBC%remap_CS, remappingScheme, boundary_extrapolation = .false., & check_reconstruction=check_reconstruction, check_remapping=check_remapping, & - force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=answers_2018) + force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=OBC%answers_2018) endif ! OBC%number_of_segments > 0 @@ -3730,6 +3733,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) real, allocatable :: normal_trans_bt(:,:) ! barotropic transport [H L2 T-1 ~> m3 s-1] integer :: turns ! Number of index quarter turns real :: time_delta ! Time since tidal reference date [T ~> s] + real :: h_neglect, h_neglect_edge ! Small thicknesses [H ~> m or kg m-2] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -3742,6 +3746,14 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) if (OBC%add_tide_constituents) time_delta = US%s_to_T * time_type_to_real(Time - OBC%time_ref) + if (.not. OBC%answers_2018) then + h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff + elseif (GV%Boussinesq) then + h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10 + else + h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10 + endif + do n = 1, OBC%number_of_segments segment => OBC%segment(n) @@ -4008,21 +4020,21 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), & - GV%H_subroundoff, GV%H_subroundoff) + h_neglect, h_neglect_edge) elseif (G%mask2dCu(I,j)>0.) then h_stack(:) = h(i+ishift,j,:) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), & - GV%H_subroundoff, GV%H_subroundoff) + h_neglect, h_neglect_edge) elseif (G%mask2dCu(I,j+1)>0.) then h_stack(:) = h(i+ishift,j+1,:) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src,segment%field(m)%dz_src(I,j,:), & segment%field(m)%buffer_src(I,J,:), & GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), & - GV%H_subroundoff, GV%H_subroundoff) + h_neglect, h_neglect_edge) endif enddo else @@ -4038,7 +4050,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(I,j,:), & segment%field(m)%buffer_src(I,j,:), & GV%ke, h(i+ishift,j,:), segment%field(m)%buffer_dst(I,j,:), & - GV%H_subroundoff, GV%H_subroundoff) + h_neglect, h_neglect_edge) endif enddo endif @@ -4058,21 +4070,21 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), & - GV%H_subroundoff, GV%H_subroundoff) + h_neglect, h_neglect_edge) elseif (G%mask2dCv(i,J)>0.) then h_stack(:) = h(i,j+jshift,:) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), & - GV%H_subroundoff, GV%H_subroundoff) + h_neglect, h_neglect_edge) elseif (G%mask2dCv(i+1,J)>0.) then h_stack(:) = h(i+1,j+jshift,:) call remapping_core_h(OBC%remap_CS, & segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), & segment%field(m)%buffer_src(I,J,:), & GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), & - GV%H_subroundoff, GV%H_subroundoff) + h_neglect, h_neglect_edge) endif enddo else @@ -4088,7 +4100,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(i,J,:), & segment%field(m)%buffer_src(i,J,:), & GV%ke, h(i,j+jshift,:), segment%field(m)%buffer_dst(i,J,:), & - GV%H_subroundoff, GV%H_subroundoff) + h_neglect, h_neglect_edge) endif enddo endif diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index a4c59be85c..140cf3a4d6 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -134,6 +134,9 @@ module MOM_oda_driver_mod type(INC_CS) :: INC_CS !< A Structure containing integer file handles for bias adjustment integer :: id_inc_t !< A diagnostic handle for the temperature climatological adjustment integer :: id_inc_s !< A diagnostic handle for the salinity climatological adjustment + logical :: answers_2018 !< If true, use the order of arithmetic and expressions for remapping + !! that recover the answers from the end of 2018. Otherwise, use more + !! robust and accurate forms of mathematically equivalent expressions. end type ODA_CS @@ -396,6 +399,7 @@ subroutine set_prior_tracer(Time, G, GV, h, tv, CS) integer :: isg, ieg, jsg, jeg, idg_offset, jdg_offset integer :: id logical :: used, symmetric + real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2] ! return if not time for analysis if (Time < CS%Time) return @@ -407,6 +411,14 @@ subroutine set_prior_tracer(Time, G, GV, h, tv, CS) call set_PElist(CS%filter_pelist) !call MOM_mesg('Setting prior') + if (.not. CS%answers_2018) then + h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff + elseif (GV%Boussinesq) then + h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10 + else + h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10 + endif + ! computational domain for the analysis grid isc=CS%Grid%isc;iec=CS%Grid%iec;jsc=CS%Grid%jsc;jec=CS%Grid%jec ! array extents for the ensemble member @@ -415,9 +427,9 @@ subroutine set_prior_tracer(Time, G, GV, h, tv, CS) ! remap temperature and salinity from the ensemble member to the analysis grid do j=G%jsc,G%jec ; do i=G%isc,G%iec call remapping_core_h(CS%remapCS, GV%ke, h(i,j,:), tv%T(i,j,:), & - CS%nk, CS%h(i,j,:), T(i,j,:), GV%H_subroundoff, GV%H_subroundoff) + CS%nk, CS%h(i,j,:), T(i,j,:), h_neglect, h_neglect_edge) call remapping_core_h(CS%remapCS, GV%ke, h(i,j,:), tv%S(i,j,:), & - CS%nk, CS%h(i,j,:), S(i,j,:), GV%H_subroundoff, GV%H_subroundoff) + CS%nk, CS%h(i,j,:), S(i,j,:), h_neglect, h_neglect_edge) enddo ; enddo ! cast ensemble members to the analysis domain do m=1,CS%ensemble_size @@ -652,6 +664,7 @@ subroutine apply_oda_tracer_increments(dt, Time_end, G, GV, tv, h, CS) real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: T !< The updated temperature [degC] real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: S !< The updated salinity [g kg-1] real :: missing_value + real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2] if (.not. associated(CS)) return if (CS%assim_method == NO_ASSIM .and. (.not. CS%do_bias_adjustment)) return @@ -668,12 +681,20 @@ subroutine apply_oda_tracer_increments(dt, Time_end, G, GV, tv, h, CS) S = S + CS%tv_bc%S endif + if (.not. CS%answers_2018) then + h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff + elseif (GV%Boussinesq) then + h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10 + else + h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10 + endif + isc=G%isc; iec=G%iec; jsc=G%jsc; jec=G%jec do j=jsc,jec; do i=isc,iec call remapping_core_h(CS%remapCS, CS%nk, CS%h(i,j,:), T(i,j,:), & - G%ke, h(i,j,:), T_inc(i,j,:), GV%H_subroundoff, GV%H_subroundoff) + G%ke, h(i,j,:), T_inc(i,j,:), h_neglect, h_neglect_edge) call remapping_core_h(CS%remapCS, CS%nk, CS%h(i,j,:), S(i,j,:), & - G%ke, h(i,j,:), S_inc(i,j,:), GV%H_subroundoff, GV%H_subroundoff) + G%ke, h(i,j,:), S_inc(i,j,:), h_neglect, h_neglect_edge) enddo; enddo