Skip to content

Commit

Permalink
*Add new update_segment_tracer_reservoirs routine
Browse files Browse the repository at this point in the history 
- The OBC tracer reservoirs were being updated in MOM_tracer_advect -
  twice each! Update them separately after tracer advection.
- The OBC tracer lengthscale was being cubed to get the volume. Change
  that to a lengthscale times the face area where the advection is
  happening.
- Changes answers if the tracer lengthscales were not set to zero.
  • Loading branch information
@kshedstrom
kshedstrom committed Aug 29, 2019
1 parent 9547fce commit 5c97200
Show file tree
Hide file tree
Showing 3 changed files with 96 additions and 90 deletions.
3 changes: 3 additions & 0 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,6 +86,7 @@ module MOM
use MOM_open_boundary, only : ocean_OBC_type, OBC_registry_type
use MOM_open_boundary, only : register_temp_salt_segments
use MOM_open_boundary, only : open_boundary_register_restarts
use MOM_open_boundary, only : update_segment_tracer_reservoirs
use MOM_set_visc, only : set_viscous_BBL, set_viscous_ML, set_visc_init
use MOM_set_visc, only : set_visc_register_restarts, set_visc_CS
use MOM_sponge, only : init_sponge_diags, sponge_CS
Expand Down Expand Up @@ -1089,6 +1090,8 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local)
call tracer_hordiff(h, CS%t_dyn_rel_adv, CS%MEKE, CS%VarMix, G, GV, US, &
CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv)
if (showCallTree) call callTree_waypoint("finished tracer advection/diffusion (step_MOM)")
call update_segment_tracer_reservoirs(G, GV, CS%uhtr, CS%vhtr, h, CS%OBC, &
CS%t_dyn_rel_adv, CS%tracer_Reg)
call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo)

call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
Expand Down
102 changes: 92 additions & 10 deletions src/core/MOM_open_boundary.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@ module MOM_open_boundary
public register_temp_salt_segments
public fill_temp_salt_segments
public open_boundary_register_restarts
public update_segment_tracer_reservoirs

integer, parameter, public :: OBC_NONE = 0 !< Indicates the use of no open boundary
integer, parameter, public :: OBC_SIMPLE = 1 !< Indicates the use of a simple inflow open boundary
Expand Down Expand Up @@ -181,11 +182,11 @@ module MOM_open_boundary
!! can occur [T-1 ~> s-1].
type(segment_tracer_registry_type), pointer :: tr_Reg=> NULL()!< A pointer to the tracer registry for the segment.
type(hor_index_type) :: HI !< Horizontal index ranges
real :: Tr_InvLscale3_out !< An effective inverse length scale cubed [m-3]
real :: Tr_InvLscale3_in !< for restoring the tracer concentration in a
!! ficticious reservior towards interior values
!! when flow is exiting the domain, or towards
!! an externally imposed value when flow is entering
real :: Tr_InvLscale_out !< An effective inverse length scale [m-1]
real :: Tr_InvLscale_in !< for restoring the tracer concentration in a
!! ficticious reservior towards interior values
!! when flow is exiting the domain, or towards
!! an externally imposed value when flow is entering
end type OBC_segment_type

!> Open-boundary data
Expand Down Expand Up @@ -494,10 +495,10 @@ subroutine open_boundary_config(G, US, param_file, OBC)
! tracer-specific in the future for example, in cases where certain tracers are poorly constrained
! by data while others are well constrained - MJH.
do l = 1, OBC%number_of_segments
OBC%segment(l)%Tr_InvLscale3_in=0.0
if (Lscale_in>0.) OBC%segment(l)%Tr_InvLscale3_in = 1.0/(Lscale_in*Lscale_in*Lscale_in)
OBC%segment(l)%Tr_InvLscale3_out=0.0
if (Lscale_out>0.) OBC%segment(l)%Tr_InvLscale3_out = 1.0/(Lscale_out*Lscale_out*Lscale_out)
OBC%segment(l)%Tr_InvLscale_in=0.0
if (Lscale_in>0.) OBC%segment(l)%Tr_InvLscale_in = 1.0/Lscale_in
OBC%segment(l)%Tr_InvLscale_out=0.0
if (Lscale_out>0.) OBC%segment(l)%Tr_InvLscale_out = 1.0/Lscale_out
enddo

endif ! OBC%number_of_segments > 0
Expand Down Expand Up @@ -4041,7 +4042,7 @@ subroutine flood_fill2(G, color, cin, cout, cland)
end subroutine flood_fill2

!> Register OBC segment data for restarts
subroutine open_boundary_register_restarts(HI, GV, OBC_CS,restart_CSp)
subroutine open_boundary_register_restarts(HI, GV, OBC_CS, restart_CSp)
type(hor_index_type), intent(in) :: HI !< Horizontal indices
type(verticalGrid_type), pointer :: GV !< Container for vertical grid information
type(ocean_OBC_type), pointer :: OBC_CS !< OBC data structure, data intent(inout)
Expand Down Expand Up @@ -4080,6 +4081,87 @@ subroutine open_boundary_register_restarts(HI, GV, OBC_CS,restart_CSp)

end subroutine open_boundary_register_restarts

!> Update the OBC tracer reservoirs after the tracers have been updated.
subroutine update_segment_tracer_reservoirs(G, GV, uhr, vhr, h, OBC, dt, Reg)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uhr !< accumulated volume/mass flux through
!! the zonal face [H L2 ~> m3 or kg]
real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vhr !< accumulated volume/mass flux through
!! the meridional face [H L2 ~> m3 or kg]
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< layer thickness after advection
!! [H ~> m or kg m-2]
type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
real, intent(in) :: dt !< time increment [s]
type(tracer_registry_type), pointer :: Reg !< pointer to tracer registry
! Local variables
integer :: i, j, k, m, n, ntr, nz
integer :: ishift, idir, jshift, jdir
type(OBC_segment_type), pointer :: segment=>NULL()
real :: u_L_in, u_L_out
real :: v_L_in, v_L_out
real :: fac1

nz = GV%ke
ntr = Reg%ntr
if (associated(OBC)) then ; if (OBC%OBC_pe) then
do n=1,OBC%number_of_segments
segment=>OBC%segment(n)
if (.not. associated(segment%tr_Reg)) cycle
if (segment%is_E_or_W) then
do j=segment%HI%jsd,segment%HI%jed
I = segment%HI%IsdB

ishift=0 ! ishift+I corresponds to the nearest interior tracer cell index
idir=1 ! idir switches the sign of the flow so that positive is into the reservoir
if (segment%direction == OBC_DIRECTION_W) then
ishift=1
idir=-1
endif
! update the reservoir tracer concentration implicitly
! using Backward-Euler timestep
do m=1,ntr
if (associated(segment%tr_Reg%Tr(m)%tres)) then
do k=1,nz
u_L_in=max((idir*uhr(I,j,k))*segment%Tr_InvLscale_in/(h(i+ishift,j,k)*G%dyCu(I,j)),0.)
u_L_out=min((idir*uhr(I,j,k))*segment%Tr_InvLscale_out/(h(i+ishift,j,k)*G%dyCu(I,j)),0.)
fac1=1.0+dt*(u_L_in-u_L_out)
segment%tr_Reg%Tr(m)%tres(I,j,k)= (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(I,j,k) + &
dt*(u_L_in*Reg%Tr(m)%t(I+ishift,j,k) - &
u_L_out*segment%tr_Reg%Tr(m)%t(I,j,k)))
enddo
endif
enddo
enddo
else
do i=segment%HI%isd,segment%HI%ied
J = segment%HI%JsdB
jshift=0 ! jshift+J corresponds to the nearest interior tracer cell index
jdir=1 ! jdir switches the sign of the flow so that positive is into the reservoir
if (segment%direction == OBC_DIRECTION_S) then
jshift=1
jdir=-1
endif
! update the reservoir tracer concentration implicitly
! using Backward-Euler timestep
do m=1,ntr
if (associated(segment%tr_Reg%Tr(m)%tres)) then
do k=1,nz
v_L_in=max((jdir*vhr(i,J,k))*segment%Tr_InvLscale_in/(h(i,j+jshift,k)*G%dxCv(i,J)),0.)
v_L_out=min((jdir*vhr(i,J,k))*segment%Tr_InvLscale_out/(h(i,j+jshift,k)*G%dxCv(i,J)),0.)
fac1=1.0+dt*(v_L_in-v_L_out)
segment%tr_Reg%Tr(m)%tres(i,J,k)= (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,J,k) + &
dt*(v_L_in*Reg%Tr(m)%t(i,J+jshift,k) - &
v_L_out*segment%tr_Reg%Tr(m)%t(i,J,k)))
enddo
endif
enddo
enddo
endif
enddo
endif; endif
end subroutine update_segment_tracer_reservoirs

!> Adjust interface heights to fit the bathymetry and diagnose layer thickness.
!!
!! If the bottom most interface is below the topography then the bottom-most
Expand Down
81 changes: 1 addition & 80 deletions src/tracer/MOM_tracer_advect.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -380,7 +380,6 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
real :: fac1,u_L_in,u_L_out ! terms used for time-stepping OBC reservoirs
type(OBC_segment_type), pointer :: segment=>NULL()
integer :: ishift, idir
real :: dt ! the inverse of Idt, needed for time-stepping of tracer reservoirs
logical :: usePLMslope

Expand Down Expand Up @@ -439,27 +438,6 @@ 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

ishift=0 ! ishift+I corresponds to the nearest interior tracer cell index
idir=1 ! idir switches the sign of the flow so that positive is into the reservoir
if (segment%direction == OBC_DIRECTION_W) then
ishift=1
idir=-1
endif
! update the reservoir tracer concentration implicitly
! using Backward-Euler timestep
do m=1,ntr
if (associated(segment%tr_Reg%Tr(m)%tres)) then
uhh(I)=uhr(I,j,k)
u_L_in=max(idir*uhh(I)*segment%Tr_InvLscale3_in,0.)
u_L_out=min(idir*uhh(I)*segment%Tr_InvLscale3_out,0.)
fac1=1.0+dt*(u_L_in-u_L_out)
segment%tr_Reg%Tr(m)%tres(I,j,k)= (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(I,j,k) + &
dt*(u_L_in*Tr(m)%t(I+ishift,j,k) - &
u_L_out*segment%tr_Reg%Tr(m)%t(I,j,k)))
endif
enddo

do m = 1,ntr ! replace tracers with OBC values
if (associated(segment%tr_Reg%Tr(m)%tres)) then
if (segment%direction == OBC_DIRECTION_W) then
Expand Down Expand Up @@ -624,7 +602,7 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
endif
endif
enddo
endif
endif

if (OBC%open_u_BCs_exist_globally) then
do n=1,OBC%number_of_segments
Expand All @@ -633,25 +611,6 @@ subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
if (segment%is_E_or_W .and. (j >= segment%HI%jsd .and. j<= segment%HI%jed)) then
if (segment%specified) cycle
if (.not. associated(segment%tr_Reg)) cycle
ishift=0 ! ishift+I corresponds to the nearest interior tracer cell index
idir=1 ! idir switches the sign of the flow so that positive is into the reservoir
if (segment%direction == OBC_DIRECTION_W) then
ishift = 1
idir = -1
endif
! update the reservoir tracer concentration implicitly
! using Backward-Euler timestep
do m=1,ntr
if (associated(segment%tr_Reg%Tr(m)%tres)) then
uhh(I) = uhr(I,j,k)
u_L_in = max(idir*uhh(I)*segment%Tr_InvLscale3_in,0.)
u_L_out = min(idir*uhh(I)*segment%Tr_InvLscale3_out,0.)
fac1 = 1.0+dt*(u_L_in-u_L_out)
segment%tr_Reg%Tr(m)%tres(I,j,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(I,j,k) + &
dt*(u_L_in*Tr(m)%t(I+ishift,j,k) - &
u_L_out*segment%tr_Reg%Tr(m)%t(I,j,k)))
endif
enddo

! Tracer fluxes are set to prescribed values only for inflows from masked areas.
if ((uhr(I,j,k) > 0.0) .and. (G%mask2dT(i,j) < 0.5) .or. &
Expand Down Expand Up @@ -777,7 +736,6 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
integer :: i, j, j2, m, n, j_up, stencil
real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
real :: fac1,v_L_in,v_L_out ! terms used for time-stepping OBC reservoirs
integer :: jshift, jdir
real :: dt ! The inverse of Idt, needed for segment reservoir time-stepping
type(OBC_segment_type), pointer :: segment=>NULL()
logical :: usePLMslope
Expand Down Expand Up @@ -847,26 +805,6 @@ 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
jshift=0 ! jshift+J corresponds to the nearest interior tracer cell index
jdir=1 ! jdir switches the sign of the flow so that positive is into the reservoir
if (segment%direction == OBC_DIRECTION_S) then
jshift=1
jdir=-1
endif
! update the reservoir tracer concentration implicitly
! using Backward-Euler timestep
do m=1,ntr
if (associated(segment%tr_Reg%Tr(m)%tres)) then
vhh(i,J)=vhr(i,J,k)
v_L_in=max(jdir*vhh(i,J)*segment%Tr_InvLscale3_in,0.)
v_L_out=min(jdir*vhh(i,J)*segment%Tr_InvLscale3_out,0.)
fac1=1.0+dt*(v_L_in-v_L_out)
segment%tr_Reg%Tr(m)%tres(i,J,k)= (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,J,k) + &
dt*(v_L_in*Tr(m)%t(i,J+jshift,k) - &
v_L_out*segment%tr_Reg%Tr(m)%t(i,J,k)))
endif
enddo

do m = 1,ntr ! replace tracers with OBC values
if (associated(segment%tr_Reg%Tr(m)%tres)) then
if (segment%direction == OBC_DIRECTION_S) then
Expand Down Expand Up @@ -1040,24 +978,7 @@ subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
if (segment%specified) cycle
if (.not. associated(segment%tr_Reg)) cycle
if (segment%is_N_or_S .and. (J >= segment%HI%JsdB .and. J<= segment%HI%JedB)) then
jshift = 0 ; jdir = 1
if (segment%direction == OBC_DIRECTION_S) then
jshift = 1 ; jdir = -1
endif
do i=segment%HI%isd,segment%HI%ied
! update the reservoir tracer concentration implicitly
! using Backward-Euler timestep
do m=1,ntr
if (associated(segment%tr_Reg%Tr(m)%tres)) then
vhh(i,J)=vhr(i,J,k)
v_L_in = max(jdir*vhh(i,J)*segment%Tr_InvLscale3_in,0.)
v_L_out = min(jdir*vhh(i,J)*segment%Tr_InvLscale3_out,0.)
fac1 = 1.0 + dt*(v_L_in-v_L_out)
segment%tr_Reg%Tr(m)%tres(i,J,k) = (1.0/fac1)*(segment%tr_Reg%Tr(m)%tres(i,J,k) + &
dt*v_L_in*Tr(m)%t(i,j+jshift,k) - &
dt*v_L_out*segment%tr_Reg%Tr(m)%t(i,j,k))
endif
enddo
! Tracer fluxes are set to prescribed values only for inflows from masked areas.
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
Expand Down

1 comment on commit 5c97200

@Hallberg-NOAA
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe that some of the expressions in update_segment_tracer_reservoirs are dimensionally inconsistent. Specifically, the volume fluxes being passed in are in units of [L2 H ~> m3], not [L2 H T-1 ~> m3 s-1], so there should be no need to multiply any of these expressions by a timestep. @kshedstrom, would you be willing to try modifying the code to avoid the multiplication by the timestep on lines 4128, 4130, 4152, and 4154 of MOM_open_boundaries.F90, and see whether you are happy with the results and agree with my diagnosis?

Please sign in to comment.