Skip to content

Commit

Permalink
(+*) Fix bugs in tracer index in tracer reservoirs (#480)
Browse files Browse the repository at this point in the history 
* Fix nudged OBCs for tracers

* Fix index in tracer reservoirs

* fix index in the function update_segment_tracer_reservoirs

* Fix index bugs in OBC tracer reservoirs

* Fix tracer index bugs at open boundaries
  • Loading branch information
@WenhaoChen89
WenhaoChen89 authored Oct 12, 2023
1 parent bd4c87c commit 95d6e93
Show file tree
Hide file tree
Showing 6 changed files with 104 additions and 63 deletions.
72 changes: 51 additions & 21 deletions src/core/MOM_open_boundary.F90
View file
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
View file
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
View file
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

0 comments on commit 95d6e93

Please sign in to comment.