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.

Sign up for GitHub

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

*Add new update_segment_tracer_reservoirs routine #987

Merged
merged 3 commits into from
Aug 30, 2019
Merged

*Add new update_segment_tracer_reservoirs routine #987

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5c97200
*Add new update_segment_tracer_reservoirs routine
kshedstrom Aug 29, 2019
2294f80
Fixed the OBC checking for data and dying if missing.
kshedstrom Aug 30, 2019
66a6478
Merge branch 'dev/gfdl' into dev/esmg
adcroft Aug 30, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions src/core/MOM.F90
View file
Open in desktop
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
110 changes: 98 additions & 12 deletions src/core/MOM_open_boundary.F90
View file
Open in desktop
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 @@ -611,7 +612,13 @@ subroutine initialize_segment_data(G, OBC, PF)
! needs documentation !! Yet, unsafe for now, causes grief for
! MOM_parameter_docs in circle_obcs on two processes.
! call get_param(PF, mdl, segnam, segstr, 'xyz')
! Clear out any old values
segstr = ''
call get_param(PF, mdl, segnam, segstr)
if (segstr == '') then
write(mesg,'("No OBC_SEGMENT_XXX_DATA string for OBC segment ",I3)') n
call MOM_error(FATAL, mesg)
endif

call parse_segment_data_str(trim(segstr), fields=fields, num_fields=num_fields)
if (num_fields == 0) then
Expand Down Expand Up @@ -771,8 +778,7 @@ subroutine initialize_segment_data(G, OBC, PF)
segment%t_values_needed .or. segment%s_values_needed .or. &
segment%z_values_needed .or. segment%g_values_needed) then
write(mesg,'("Values needed for OBC segment ",I3)') n
! call MOM_error(FATAL, mesg)
call MOM_error(WARNING, mesg)
call MOM_error(FATAL, mesg)
endif
enddo

Expand Down Expand Up @@ -4041,7 +4047,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 +4086,86 @@ 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
Open in desktop
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