Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

(+*) Fix bugs in tracer index in tracer reservoirs #480

Merged
merged 12 commits into from
Oct 12, 2023
72 changes: 51 additions & 21 deletions src/core/MOM_open_boundary.F90
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@ module MOM_open_boundary
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : thermo_var_ptrs
use MOM_verticalGrid, only : verticalGrid_type
use MOM_string_functions, only : lowercase

implicit none ; private

Expand Down Expand Up @@ -118,6 +119,8 @@ module MOM_open_boundary
real :: scale !< A scaling factor for converting the units of input
!! data, like [S ppt-1 ~> 1] for salinity.
logical :: is_initialized !< reservoir values have been set when True
integer :: ntr_index = -1 !< index of segment tracer in the global tracer registry
integer :: fd_index = -1 !< index of segment tracer in the input fields
end type OBC_segment_tracer_type

!> Registry type for tracers on segments
Expand Down Expand Up @@ -4647,8 +4650,8 @@ end subroutine segment_tracer_registry_init

!> Register a tracer array that is active on an OBC segment, potentially also specifying how the
!! tracer inflow values are specified.
subroutine register_segment_tracer(tr_ptr, param_file, GV, segment, &
OBC_scalar, OBC_array, scale)
subroutine register_segment_tracer(tr_ptr, ntr_index, param_file, GV, segment, &
OBC_scalar, OBC_array, scale, fd_index)
type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
type(tracer_type), target :: tr_ptr !< A target that can be used to set a pointer to the
!! stored value of tr. This target must be
Expand All @@ -4657,6 +4660,7 @@ subroutine register_segment_tracer(tr_ptr, param_file, GV, segment, &
!! but it also means that any updates to this
!! structure in the calling module will be
!! available subsequently to the tracer registry.
integer, intent(in) :: ntr_index !< index of segment tracer in the global tracer registry
type(param_file_type), intent(in) :: param_file !< file to parse for model parameter values
type(OBC_segment_type), intent(inout) :: segment !< current segment data structure
real, optional, intent(in) :: OBC_scalar !< If present, use scalar value for segment tracer
Expand All @@ -4667,6 +4671,7 @@ subroutine register_segment_tracer(tr_ptr, param_file, GV, segment, &
real, optional, intent(in) :: scale !< A scaling factor that should be used with any
!! data that is read in, to convert it to the internal
!! units of this tracer.
integer, optional, intent(in) :: fd_index !< index of segment tracer in the input field

! Local variables
real :: rescale ! A multiplicative correction to the scaling factor.
Expand All @@ -4690,6 +4695,8 @@ subroutine register_segment_tracer(tr_ptr, param_file, GV, segment, &

segment%tr_Reg%Tr(ntseg)%Tr => tr_ptr
segment%tr_Reg%Tr(ntseg)%name = tr_ptr%name
segment%tr_Reg%Tr(ntseg)%ntr_index = ntr_index
if (present(fd_index)) segment%tr_Reg%Tr(ntseg)%fd_index = fd_index

segment%tr_Reg%Tr(ntseg)%scale = 1.0
if (present(scale)) then
Expand Down Expand Up @@ -4752,7 +4759,7 @@ subroutine register_temp_salt_segments(GV, US, OBC, tr_Reg, param_file)
type(param_file_type), intent(in) :: param_file !< file to parse for model parameter values

! Local variables
integer :: n
integer :: n, ntr_id
character(len=32) :: name
type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list
type(tracer_type), pointer :: tr_ptr => NULL()
Expand All @@ -4767,12 +4774,12 @@ subroutine register_temp_salt_segments(GV, US, OBC, tr_Reg, param_file)
call MOM_error(FATAL,"register_temp_salt_segments: tracer array was previously allocated")

name = 'temp'
call tracer_name_lookup(tr_Reg, tr_ptr, name)
call register_segment_tracer(tr_ptr, param_file, GV, segment, &
call tracer_name_lookup(tr_Reg, ntr_id, tr_ptr, name)
call register_segment_tracer(tr_ptr, ntr_id, param_file, GV, segment, &
OBC_array=segment%temp_segment_data_exists, scale=US%degC_to_C)
name = 'salt'
call tracer_name_lookup(tr_Reg, tr_ptr, name)
call register_segment_tracer(tr_ptr, param_file, GV, segment, &
call tracer_name_lookup(tr_Reg, ntr_id, tr_ptr, name)
call register_segment_tracer(tr_ptr, ntr_id, param_file, GV, segment, &
OBC_array=segment%salt_segment_data_exists, scale=US%ppt_to_S)
enddo

Expand Down Expand Up @@ -4825,8 +4832,8 @@ subroutine register_obgc_segments(GV, OBC, tr_Reg, param_file, tr_name)
type(param_file_type), intent(in) :: param_file !< file to parse for model parameter values
character(len=*), intent(in) :: tr_name!< Tracer name
! Local variables
integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, nz, nf
integer :: i, j, k, n
integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, nz, nf, ntr_id, fd_id
integer :: i, j, k, n, m
type(OBC_segment_type), pointer :: segment => NULL() ! pointer to segment type list
type(tracer_type), pointer :: tr_ptr => NULL()

Expand All @@ -4835,8 +4842,13 @@ subroutine register_obgc_segments(GV, OBC, tr_Reg, param_file, tr_name)
do n=1, OBC%number_of_segments
segment=>OBC%segment(n)
if (.not. segment%on_pe) cycle
call tracer_name_lookup(tr_Reg, tr_ptr, tr_name)
call register_segment_tracer(tr_ptr, param_file, GV, segment, OBC_array=.True.)
call tracer_name_lookup(tr_Reg, ntr_id, tr_ptr, tr_name)
! get the obgc field index
fd_id = -1
do m=1,segment%num_fields
if (lowercase(segment%field(m)%name) == lowercase(tr_name)) fd_id = m
enddo
call register_segment_tracer(tr_ptr, ntr_id, param_file, GV, segment, OBC_array=.True., fd_index=fd_id)
enddo

end subroutine register_obgc_segments
Expand Down Expand Up @@ -5336,8 +5348,9 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
real :: fac1 ! The denominator of the expression for tracer updates [nondim]
real :: I_scale ! The inverse of the scaling factor for the tracers.
! For salinity the units would be [ppt S-1 ~> 1]
integer :: i, j, k, m, n, ntr, nz
integer :: i, j, k, m, n, ntr, nz, ntr_id, fd_id
integer :: ishift, idir, jshift, jdir
real :: resrv_lfac_out, resrv_lfac_in
real :: b_in, b_out ! The 0 and 1 switch for tracer reservoirs
! 1 if the length scale of reservoir is zero [nondim]
real :: a_in, a_out ! The 0 and 1(-1) switch for reservoir source weights
Expand Down Expand Up @@ -5365,7 +5378,16 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
! Can keep this or take it out, either way
if (G%mask2dT(I+ishift,j) == 0.0) cycle
! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep
do m=1,ntr
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
if(fd_id == -1) then
resrv_lfac_out = 1.0
resrv_lfac_in = 1.0
else
resrv_lfac_out = segment%field(fd_id)%resrv_lfac_out
resrv_lfac_in = segment%field(fd_id)%resrv_lfac_in
endif
I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale
if (allocated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz
! Calculate weights. Both a and u_L are nodim. Adding them together has no meaning.
Expand All @@ -5374,14 +5396,14 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
! When InvLscale_in is 0 and inflow, only nudged data is applied to reservoirs
a_out = b_out * max(0.0, sign(1.0, idir*uhr(I,j,k)))
a_in = b_in * min(0.0, sign(1.0, idir*uhr(I,j,k)))
u_L_out = max(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_out*segment%field(m)%resrv_lfac_out / &
u_L_out = max(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_out*resrv_lfac_out / &
((h(i+ishift,j,k) + GV%H_subroundoff)*G%dyCu(I,j)))
u_L_in = min(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_in*segment%field(m)%resrv_lfac_in / &
u_L_in = min(0.0, (idir*uhr(I,j,k))*segment%Tr_InvLscale_in*resrv_lfac_in / &
((h(i+ishift,j,k) + GV%H_subroundoff)*G%dyCu(I,j)))
fac1 = (1.0 - (a_out - a_in)) + ((u_L_out + a_out) - (u_L_in + a_in))
segment%tr_Reg%Tr(m)%tres(I,j,k) = (1.0/fac1) * &
((1.0-a_out+a_in)*segment%tr_Reg%Tr(m)%tres(I,j,k)+ &
((u_L_out+a_out)*Reg%Tr(m)%t(I+ishift,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)
enddo ; endif
Expand All @@ -5400,20 +5422,28 @@ subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
! Can keep this or take it out, either way
if (G%mask2dT(i,j+jshift) == 0.0) cycle
! Update the reservoir tracer concentration implicitly using a Backward-Euler timestep
do m=1,ntr
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
if(fd_id == -1) then
resrv_lfac_out = 1.0
resrv_lfac_in = 1.0
else
resrv_lfac_out = segment%field(fd_id)%resrv_lfac_out
resrv_lfac_in = segment%field(fd_id)%resrv_lfac_in
endif
I_scale = 1.0 ; if (segment%tr_Reg%Tr(m)%scale /= 0.0) I_scale = 1.0 / segment%tr_Reg%Tr(m)%scale
if (allocated(segment%tr_Reg%Tr(m)%tres)) then ; do k=1,nz
a_out = b_out * max(0.0, sign(1.0, jdir*vhr(i,J,k)))
a_in = b_in * min(0.0, sign(1.0, jdir*vhr(i,J,k)))
v_L_out = max(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_out*segment%field(m)%resrv_lfac_out / &
v_L_out = max(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_out*resrv_lfac_out / &
((h(i,j+jshift,k) + GV%H_subroundoff)*G%dxCv(i,J)))
v_L_in = min(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_in*segment%field(m)%resrv_lfac_in / &
v_L_in = min(0.0, (jdir*vhr(i,J,k))*segment%Tr_InvLscale_in*resrv_lfac_in / &
((h(i,j+jshift,k) + GV%H_subroundoff)*G%dxCv(i,J)))
fac1 = 1.0 + (v_L_out-v_L_in)
fac1 = (1.0 - (a_out - a_in)) + ((v_L_out + a_out) - (v_L_in + a_in))
segment%tr_Reg%Tr(m)%tres(i,J,k) = (1.0/fac1) * &
((1.0-a_out+a_in)*segment%tr_Reg%Tr(m)%tres(i,J,k) + &
((v_L_out+a_out)*Reg%Tr(m)%t(i,J+jshift,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)
enddo ; endif
Expand Down
54 changes: 30 additions & 24 deletions src/tracer/MOM_tracer_advect.F90
Original file line number Diff line number Diff line change
Expand Up @@ -381,7 +381,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
real :: a6 ! Curvature of the reconstruction tracer values [conc]
logical :: do_i(SZIB_(G),SZJ_(G)) ! If true, work on given points.
logical :: usePLMslope
integer :: i, j, m, n, i_up, stencil
integer :: i, j, m, n, i_up, stencil, ntr_id
type(OBC_segment_type), pointer :: segment=>NULL()
logical, dimension(SZJ_(G),SZK_(GV)) :: domore_u_initial

Expand Down Expand Up @@ -442,18 +442,19 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
if (segment%is_E_or_W) then
if (j>=segment%HI%jsd .and. j<=segment%HI%jed) then
I = segment%HI%IsdB
do m = 1,ntr ! replace tracers with OBC values
do m = 1,segment%tr_Reg%ntseg ! replace tracers with OBC values
ntr_id = segment%tr_reg%Tr(m)%ntr_index
if (allocated(segment%tr_Reg%Tr(m)%tres)) then
if (segment%direction == OBC_DIRECTION_W) then
T_tmp(i,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
T_tmp(i,ntr_id) = segment%tr_Reg%Tr(m)%tres(i,j,k)
else
T_tmp(i+1,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
T_tmp(i+1,ntr_id) = segment%tr_Reg%Tr(m)%tres(i,j,k)
endif
else
if (segment%direction == OBC_DIRECTION_W) then
T_tmp(i,m) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
T_tmp(i,ntr_id) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
else
T_tmp(i+1,m) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
T_tmp(i+1,ntr_id) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
endif
endif
enddo
Expand Down Expand Up @@ -586,10 +587,11 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
(uhr(I,j,k) < 0.0) .and. (segment%direction == OBC_DIRECTION_E)) then
uhh(I) = uhr(I,j,k)
! should the reservoir evolve for this case Kate ?? - Nope
do m=1,ntr
do m=1,segment%tr_Reg%ntseg
ntr_id = segment%tr_reg%Tr(m)%ntr_index
if (allocated(segment%tr_Reg%Tr(m)%tres)) then
flux_x(I,j,m) = uhh(I)*segment%tr_Reg%Tr(m)%tres(I,j,k)
else ; flux_x(I,j,m) = uhh(I)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
flux_x(I,j,ntr_id) = uhh(I)*segment%tr_Reg%Tr(m)%tres(I,j,k)
else ; flux_x(I,j,ntr_id) = uhh(I)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
enddo
endif
endif
Expand All @@ -609,10 +611,11 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
if ((uhr(I,j,k) > 0.0) .and. (G%mask2dT(i,j) < 0.5) .or. &
(uhr(I,j,k) < 0.0) .and. (G%mask2dT(i+1,j) < 0.5)) then
uhh(I) = uhr(I,j,k)
do m=1,ntr
do m=1,segment%tr_Reg%ntseg
ntr_id = segment%tr_reg%Tr(m)%ntr_index
if (allocated(segment%tr_Reg%Tr(m)%tres)) then
flux_x(I,j,m) = uhh(I)*segment%tr_Reg%Tr(m)%tres(I,j,k)
else; flux_x(I,j,m) = uhh(I)*segment%tr_Reg%Tr(m)%OBC_inflow_conc; endif
flux_x(I,j,ntr_id) = uhh(I)*segment%tr_Reg%Tr(m)%tres(I,j,k)
else; flux_x(I,j,ntr_id) = uhh(I)*segment%tr_Reg%Tr(m)%OBC_inflow_conc; endif
enddo
endif
endif
Expand Down Expand Up @@ -754,7 +757,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
logical :: do_j_tr(SZJ_(G)) ! If true, calculate the tracer profiles.
logical :: do_i(SZIB_(G), SZJ_(G)) ! If true, work on given points.
logical :: usePLMslope
integer :: i, j, j2, m, n, j_up, stencil
integer :: i, j, j2, m, n, j_up, stencil, ntr_id
type(OBC_segment_type), pointer :: segment=>NULL()
logical :: domore_v_initial(SZJB_(G)) ! Initial state of domore_v

Expand Down Expand Up @@ -823,18 +826,19 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
if (segment%is_N_or_S) then
if (i>=segment%HI%isd .and. i<=segment%HI%ied) then
J = segment%HI%JsdB
do m = 1,ntr ! replace tracers with OBC values
do m = 1,segment%tr_Reg%ntseg ! replace tracers with OBC values
ntr_id = segment%tr_reg%Tr(m)%ntr_index
if (allocated(segment%tr_Reg%Tr(m)%tres)) then
if (segment%direction == OBC_DIRECTION_S) then
T_tmp(i,m,j) = segment%tr_Reg%Tr(m)%tres(i,j,k)
T_tmp(i,ntr_id,j) = segment%tr_Reg%Tr(m)%tres(i,j,k)
else
T_tmp(i,m,j+1) = segment%tr_Reg%Tr(m)%tres(i,j,k)
T_tmp(i,ntr_id,j+1) = segment%tr_Reg%Tr(m)%tres(i,j,k)
endif
else
if (segment%direction == OBC_DIRECTION_S) then
T_tmp(i,m,j) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
T_tmp(i,ntr_id,j) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
else
T_tmp(i,m,j+1) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
T_tmp(i,ntr_id,j+1) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
endif
endif
enddo
Expand Down Expand Up @@ -968,10 +972,11 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
if ((vhr(i,J,k) > 0.0) .and. (segment%direction == OBC_DIRECTION_S) .or. &
(vhr(i,J,k) < 0.0) .and. (segment%direction == OBC_DIRECTION_N)) then
vhh(i,J) = vhr(i,J,k)
do m=1,ntr
do m=1,segment%tr_Reg%ntseg
ntr_id = segment%tr_reg%Tr(m)%ntr_index
if (allocated(segment%tr_Reg%Tr(m)%tres)) then
flux_y(i,m,J) = vhh(i,J)*OBC%segment(n)%tr_Reg%Tr(m)%tres(i,J,k)
else ; flux_y(i,m,J) = vhh(i,J)*OBC%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
flux_y(i,ntr_id,J) = vhh(i,J)*OBC%segment(n)%tr_Reg%Tr(m)%tres(i,J,k)
else ; flux_y(i,ntr_id,J) = vhh(i,J)*OBC%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
enddo
endif
enddo
Expand All @@ -991,10 +996,11 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
if ((vhr(i,J,k) > 0.0) .and. (G%mask2dT(i,j) < 0.5) .or. &
(vhr(i,J,k) < 0.0) .and. (G%mask2dT(i,j+1) < 0.5)) then
vhh(i,J) = vhr(i,J,k)
do m=1,ntr
do m=1,segment%tr_Reg%ntseg
ntr_id = segment%tr_reg%Tr(m)%ntr_index
if (allocated(segment%tr_Reg%Tr(m)%tres)) then
flux_y(i,m,J) = vhh(i,J)*segment%tr_Reg%Tr(m)%tres(i,J,k)
else ; flux_y(i,m,J) = vhh(i,J)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
flux_y(i,ntr_id,J) = vhh(i,J)*segment%tr_Reg%Tr(m)%tres(i,J,k)
else ; flux_y(i,ntr_id,J) = vhh(i,J)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
enddo
endif
enddo
Expand Down
11 changes: 8 additions & 3 deletions src/tracer/MOM_tracer_registry.F90
Original file line number Diff line number Diff line change
Expand Up @@ -848,16 +848,21 @@ end subroutine tracer_Reg_chkinv


!> Find a tracer in the tracer registry by name.
subroutine tracer_name_lookup(Reg, tr_ptr, name)
subroutine tracer_name_lookup(Reg, n, tr_ptr, name)
type(tracer_registry_type), pointer :: Reg !< pointer to tracer registry
type(tracer_type), pointer :: tr_ptr !< target or pointer to the tracer array
character(len=32), intent(in) :: name !< tracer name
integer, intent(out) :: n !< index to tracer registery

integer n
do n=1,Reg%ntr
if (lowercase(Reg%Tr(n)%name) == lowercase(name)) tr_ptr => Reg%Tr(n)
if (lowercase(Reg%Tr(n)%name) == lowercase(name)) then
tr_ptr => Reg%Tr(n)
return
endif
enddo

call MOM_error(FATAL,"MOM cannot find registered tracer: "//name)

end subroutine tracer_name_lookup

!> Initialize the tracer registry.
Expand Down
Loading