Skip to content

Commit

Permalink
+Add optional argument du_cor to continuity_PPM
Browse files Browse the repository at this point in the history 
  Added the new optional arguments du_cor and dv_cor to continuity_PPM and
zonal_mass_flux or meridional_mass_flux to return the barotropic velocity
increments that make the summed barotropic transports match uhbt or vhbt.  Also
flipped the sign of the now unused du_aux and dv_aux optional arguments.  All
answers are bitwise identical, but there are new optional arguments to three
publicly visible routines.
  • Loading branch information
@Hallberg-NOAA
Hallberg-NOAA committed Dec 12, 2023
1 parent 480da18 commit 7028b4c
Showing 1 changed file with 42 additions and 18 deletions.
60 changes: 42 additions & 18 deletions src/core/MOM_continuity_PPM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@ module MOM_continuity_PPM
!! based on Lin (1994).
subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhbt, vhbt, &
visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont, &
uhbt_aux, vhbt_aux, ddu, ddv)
uhbt_aux, vhbt_aux, ddu, ddv, du_cor, dv_cor)
type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
Expand Down Expand Up @@ -142,12 +142,18 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhb
!! perhaps at a different time [H L2 T-1 ~> m3 s-1 or kg s-1].
real, dimension(SZIB_(G),SZJ_(G)), &
optional, intent(out) :: ddu !< The difference between the barotropic zonal velocity
!! increment required to give a transport of uhbt_aux
!! and the one that gives uhbt [L T-1 ~> m s-1]
!! increment required to give a transport of uhbt
!! and the one that gives uhbt_aux [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJB_(G)), &
optional, intent(out) :: ddv !< The difference between the barotropic meridional velocity
!! increment required to give a transport of vhbt_aux
!! and the one that gives vhbt [L T-1 ~> m s-1]
!! increment required to give a transport of vhbt
!! and the one that gives vhbt_aux [L T-1 ~> m s-1]
real, dimension(SZIB_(G),SZJ_(G)), &
optional, intent(out) :: du_cor !< The zonal velocity increments from u that give uhbt
!! as the depth-integrated transports [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJB_(G)), &
optional, intent(out) :: dv_cor !< The meridional velocity increments from v that give vhbt
!! as the depth-integrated transports [L T-1 ~> m s-1].

! Local variables
real :: h_W(SZI_(G),SZJ_(G),SZK_(GV)) ! West edge thicknesses in the zonal PPM reconstruction [H ~> m or kg m-2]
Expand All @@ -174,29 +180,29 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhb
LB = set_continuity_loop_bounds(G, CS, i_stencil=.false., j_stencil=.true.)
call zonal_edge_thickness(hin, h_W, h_E, G, GV, US, CS, OBC, LB)
call zonal_mass_flux(u, hin, h_W, h_E, uh, dt, G, GV, US, CS, OBC, pbv%por_face_areaU, &
LB, uhbt, visc_rem_u, u_cor, BT_cont, uhbt_aux, ddu)
LB, uhbt, visc_rem_u, u_cor, BT_cont, uhbt_aux, ddu, du_cor)
call continuity_zonal_convergence(h, uh, dt, G, GV, LB, hin)

! Now advect meridionally, using the updated thicknesses to determine the fluxes.
LB = set_continuity_loop_bounds(G, CS, i_stencil=.false., j_stencil=.false.)
call meridional_edge_thickness(h, h_S, h_N, G, GV, US, CS, OBC, LB)
call meridional_mass_flux(v, h, h_S, h_N, vh, dt, G, GV, US, CS, OBC, pbv%por_face_areaV, &
LB, vhbt, visc_rem_v, v_cor, BT_cont, vhbt_aux, ddv)
LB, vhbt, visc_rem_v, v_cor, BT_cont, vhbt_aux, ddv, dv_cor)
call continuity_merdional_convergence(h, vh, dt, G, GV, LB, hmin=h_min)

else ! .not. x_first
! First advect meridionally, with loop bounds that accomodate the subsequent zonal advection.
LB = set_continuity_loop_bounds(G, CS, i_stencil=.true., j_stencil=.false.)
call meridional_edge_thickness(hin, h_S, h_N, G, GV, US, CS, OBC, LB)
call meridional_mass_flux(v, hin, h_S, h_N, vh, dt, G, GV, US, CS, OBC, pbv%por_face_areaV, &
LB, vhbt, visc_rem_v, v_cor, BT_cont, vhbt_aux, ddv)
LB, vhbt, visc_rem_v, v_cor, BT_cont, vhbt_aux, ddv, dv_cor)
call continuity_merdional_convergence(h, vh, dt, G, GV, LB, hin)

! Now advect zonally, using the updated thicknesses to determine the fluxes.
LB = set_continuity_loop_bounds(G, CS, i_stencil=.false., j_stencil=.false.)
call zonal_edge_thickness(h, h_W, h_E, G, GV, US, CS, OBC, LB)
call zonal_mass_flux(u, h, h_W, h_E, uh, dt, G, GV, US, CS, OBC, pbv%por_face_areaU, &
LB, uhbt, visc_rem_u, u_cor, BT_cont, uhbt_aux, ddu)
LB, uhbt, visc_rem_u, u_cor, BT_cont, uhbt_aux, ddu, du_cor)
call continuity_zonal_convergence(h, uh, dt, G, GV, LB, hmin=h_min)
endif

Expand Down Expand Up @@ -526,7 +532,7 @@ end subroutine meridional_edge_thickness

!> Calculates the mass or volume fluxes through the zonal faces, and other related quantities.
subroutine zonal_mass_flux(u, h_in, h_W, h_E, uh, dt, G, GV, US, CS, OBC, por_face_areaU, &
LB_in, uhbt, visc_rem_u, u_cor, BT_cont, uhbt_aux, ddu)
LB_in, uhbt, visc_rem_u, u_cor, BT_cont, uhbt_aux, ddu, du_cor)
type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
Expand Down Expand Up @@ -568,8 +574,11 @@ subroutine zonal_mass_flux(u, h_in, h_W, h_E, uh, dt, G, GV, US, CS, OBC, por_fa
!! perhaps at a different time [H L2 T-1 ~> m3 s-1 or kg s-1].
real, dimension(SZIB_(G),SZJ_(G)), &
optional, intent(out) :: ddu !< The difference between the barotropic zonal velocity
!! increment required to give a transport of uhbt_aux
!! and the one that gives uhbt [L T-1 ~> m s-1]
!! increment required to give a transport of uhbt
!! and the one that gives uhbt_aux [L T-1 ~> m s-1]
real, dimension(SZIB_(G),SZJ_(G)), &
optional, intent(out) :: du_cor !< The zonal velocity increments from u that give uhbt
!! as the depth-integrated transports [L T-1 ~> m s-1].

! Local variables
real, dimension(SZIB_(G),SZK_(GV)) :: duhdu ! Partial derivative of uh with u [H L ~> m2 or kg m-1].
Expand Down Expand Up @@ -615,6 +624,7 @@ subroutine zonal_mass_flux(u, h_in, h_W, h_E, uh, dt, G, GV, US, CS, OBC, por_fa
if (present(ddu)) ddu(:,:) = 0.0
if (present(ddu) .and. .not.(present(uhbt) .and. present(uhbt_aux))) &
call MOM_error(FATAL, "zonal_mass_flux: uhbt and uhbt_aux must be present to calculate ddu.")
if (present(du_cor)) du_cor(:,:) = 0.0

if (present(LB_in)) then
LB = LB_in
Expand Down Expand Up @@ -762,7 +772,7 @@ subroutine zonal_mass_flux(u, h_in, h_W, h_E, uh, dt, G, GV, US, CS, OBC, por_fa
du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, &
j, ish, ieh, do_I, por_face_areaU, OBC=OBC)
! Specified OBC points give ddu = 0 and do not need special handling.
do I=ish-1,ieh ; ddu(I,j) = du_aux(I) - du(I) ; enddo
do I=ish-1,ieh ; ddu(I,j) = du(I) - du_aux(I) ; enddo
endif

if (present(u_cor)) then ; do k=1,nz
Expand All @@ -772,6 +782,10 @@ subroutine zonal_mass_flux(u, h_in, h_W, h_E, uh, dt, G, GV, US, CS, OBC, por_fa
endif ; enddo ; endif
enddo ; endif ! u-corrected

if (present(du_cor)) then
do I=ish-1,ieh ; du_cor(I,j) = du(I) ; enddo
endif

endif

if (set_BT_cont) then
Expand Down Expand Up @@ -1432,7 +1446,7 @@ end subroutine set_zonal_BT_cont

!> Calculates the mass or volume fluxes through the meridional faces, and other related quantities.
subroutine meridional_mass_flux(v, h_in, h_S, h_N, vh, dt, G, GV, US, CS, OBC, por_face_areaV, &
LB_in, vhbt, visc_rem_v, v_cor, BT_cont, vhbt_aux, ddv)
LB_in, vhbt, visc_rem_v, v_cor, BT_cont, vhbt_aux, ddv, dv_cor)
type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure.
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]
Expand Down Expand Up @@ -1473,8 +1487,12 @@ subroutine meridional_mass_flux(v, h_in, h_S, h_N, vh, dt, G, GV, US, CS, OBC, p
real, dimension(SZI_(G),SZJB_(G)), &
optional, intent(out) :: ddv !< The difference between the barotropic
!! meridional velocity increment required to
!! give a transport of vhbt_aux and the one
!! that gives vhbt [L T-1 ~> m s-1]
!! give a transport of vhbt and the one
!! that gives vhbt_aux [L T-1 ~> m s-1]
real, dimension(SZI_(G),SZJB_(G)), &
optional, intent(out) :: dv_cor !< The meridional velocity increments from v
!! that give vhbt as the depth-integrated
!! transports [L T-1 ~> m s-1].

! Local variables
real, dimension(SZI_(G),SZK_(GV)) :: &
Expand Down Expand Up @@ -1521,6 +1539,7 @@ subroutine meridional_mass_flux(v, h_in, h_S, h_N, vh, dt, G, GV, US, CS, OBC, p
if (present(ddv)) ddv(:,:) = 0.0
if (present(ddv) .and. .not.(present(vhbt) .and. present(vhbt_aux))) &
call MOM_error(FATAL, "meridional_mass_flux: vhbt and vhbt_aux must be present to calculate ddv.")
if (present(dv_cor)) dv_cor(:,:) = 0.0

if (present(LB_in)) then
LB = LB_in
Expand Down Expand Up @@ -1665,7 +1684,7 @@ subroutine meridional_mass_flux(v, h_in, h_S, h_N, vh, dt, G, GV, US, CS, OBC, p
dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, &
j, ish, ieh, do_I, por_face_areaV, OBC=OBC)
! Specified OBC points give ddv = 0 and do not need special handling.
do i=ish,ieh ; ddv(i,J) = dv_aux(i) - dv(i) ; enddo
do i=ish,ieh ; ddv(i,J) = dv(i) - dv_aux(i) ; enddo
endif

if (present(v_cor)) then ; do k=1,nz
Expand All @@ -1674,6 +1693,11 @@ subroutine meridional_mass_flux(v, h_in, h_S, h_N, vh, dt, G, GV, US, CS, OBC, p
v_cor(i,J,k) = OBC%segment(OBC%segnum_v(i,J))%normal_vel(i,J,k)
endif ; enddo ; endif
enddo ; endif ! v-corrected

if (present(dv_cor)) then
do i=ish,ieh ; dv_cor(i,J) = dv(i) ; enddo
endif

endif

if (set_BT_cont) then
Expand Down Expand Up @@ -2079,7 +2103,7 @@ subroutine meridional_flux_adjust(v, h_in, h_S, h_N, vhbt, vh_tot_0, dvhdv_tot_0
dv_min, & ! Lower limit on dv correction based on CFL limits and previous iterations [L T-1 ~> m s-1]
dv_max ! Upper limit on dv correction based on CFL limits and previous iterations [L T-1 ~> m s-1]
real :: dv_prev ! The previous value of dv [L T-1 ~> m s-1].
real :: ddv ! The change in dv from the previous iteration [L T-1 ~> m s-1].
real :: ddv ! The change in dv from the previous iteration [L T-1 ~> m s-1].
real :: tol_eta ! The tolerance for the current iteration [H ~> m or kg m-2].
real :: tol_vel ! The tolerance for velocity in the current iteration [L T-1 ~> m s-1].
integer :: i, k, nz, itt, max_itts = 20
Expand Down

0 comments on commit 7028b4c

Please sign in to comment.