Skip to content

Commit

Permalink
More cleaning up for issue NOAA-EMC#54
Browse files Browse the repository at this point in the history
  • Loading branch information
@kshedstrom
kshedstrom committed Apr 7, 2022
1 parent 4b218c0 commit 3c7bd19
Show file tree
Hide file tree
Showing 4 changed files with 53 additions and 20 deletions.
2 changes: 1 addition & 1 deletion src/ALE/MOM_ALE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
8 changes: 4 additions & 4 deletions src/core/MOM_barotropic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
34 changes: 23 additions & 11 deletions src/core/MOM_open_boundary.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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"

Expand Down Expand Up @@ -608,15 +611,15 @@ 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)

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

Expand Down Expand Up @@ -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
Expand All @@ -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)

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
29 changes: 25 additions & 4 deletions src/ocean_data_assim/MOM_oda_driver.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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


Expand Down

0 comments on commit 3c7bd19

Please sign in to comment.