From e368ea72edc25e336399e7f6589bcdadfdca1f19 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 7 Nov 2023 08:48:04 -0500 Subject: [PATCH] +(*)Add continuity_fluxes & continuity_adjust_vel Added the new public interfaces continuity_fluxes and continuity_adjust_vel to make use of the continuity code without actually updating the layer thicknesses, and made the existing routines zonal_mass_flux and meridional_mass_flux in the continuity_PPM module public. Continuity_fluxes is overloaded to provide either the layer thickness fluxes or the vertically summed barotropic fluxes. As a part of this change, the LB arguments to meridional_mass_flux and zonal_mass_flux were made optional and moved toward the end of the list of arguments. The intent of the ocean_grid_type argument to 11 routines in the continuity_PPM module was changed from inout to in. Missing factors of por_face_area[UV] were also added to merid_face_thickness and zonal_face_thickness when the visc_rem arguments are not present or where OBCs are in use. This could change answers, but it seems very unlikely that any impacted cases are in use yet. Answers are bitwise identical all known cases, but there are 4 new public interfaces, and some bugs were fixed in cases that are not likely in use yet. --- src/core/MOM_continuity.F90 | 192 ++++++++++++++++++++++++++++++-- src/core/MOM_continuity_PPM.F90 | 89 +++++++++------ 2 files changed, 239 insertions(+), 42 deletions(-) diff --git a/src/core/MOM_continuity.F90 b/src/core/MOM_continuity.F90 index 76e1bbc623..7a833ccb5c 100644 --- a/src/core/MOM_continuity.F90 +++ b/src/core/MOM_continuity.F90 @@ -3,8 +3,8 @@ module MOM_continuity ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_continuity_PPM, only : continuity_PPM, continuity_PPM_init -use MOM_continuity_PPM, only : continuity_PPM_stencil +use MOM_continuity_PPM, only : continuity_PPM, zonal_mass_flux, meridional_mass_flux +use MOM_continuity_PPM, only : continuity_PPM_stencil, continuity_PPM_init use MOM_continuity_PPM, only : continuity_PPM_CS use MOM_diag_mediator, only : time_type, diag_ctrl use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe @@ -20,7 +20,7 @@ module MOM_continuity #include -public continuity, continuity_init, continuity_stencil +public continuity, continuity_fluxes, continuity_adjust_vel, continuity_init, continuity_stencil !> Control structure for mom_continuity type, public :: continuity_CS ; private @@ -35,10 +35,16 @@ module MOM_continuity integer, parameter :: PPM_SCHEME = 1 !< Enumerated constant to select PPM character(len=20), parameter :: PPM_STRING = "PPM" !< String to select PPM +!> Finds the thickness fluxes from the continuity solver or their vertical sum without +!! actually updating the layer thicknesses. +interface continuity_fluxes + module procedure continuity_3d_fluxes, continuity_2d_fluxes +end interface continuity_fluxes + contains -!> Time steps the layer thicknesses, using a monotonically limited, directionally split PPM scheme, -!! based on Lin (1994). +!> Time steps the layer thicknesses, using a monotonically or positive-definite limited, directionally +!! split PPM scheme, based on Lin (1994). subroutine continuity(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) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure. @@ -63,21 +69,23 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhbt, v type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics real, dimension(SZIB_(G),SZJ_(G)), & - optional, intent(in) :: uhbt !< The vertically summed volume - !! flux through zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1]. + optional, intent(in) :: uhbt !< The vertically summed volume flux through + !! zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1]. real, dimension(SZI_(G),SZJB_(G)), & - optional, intent(in) :: vhbt !< The vertically summed volume - !! flux through meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1]. + optional, intent(in) :: vhbt !< The vertically summed volume flux through + !! meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1]. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & optional, intent(in) :: visc_rem_u !< Both the fraction of !! zonal momentum that remains after a time-step of viscosity, and the fraction of a time-step's !! worth of a barotropic acceleration that a layer experiences after viscosity is applied [nondim]. - !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). + !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). When this column is + !! under an ice shelf, this can also go to 0 at the top due to the no-slip boundary condition there. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & optional, intent(in) :: visc_rem_v !< Both the fraction of !! meridional momentum that remains after a time-step of viscosity, and the fraction of a time-step's !! worth of a barotropic acceleration that a layer experiences after viscosity is applied [nondim]. - !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). + !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). When this column is + !! under an ice shelf, this can also go to 0 at the top due to the no-slip boundary condition there. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & optional, intent(out) :: u_cor !< The zonal velocities that !! give uhbt as the depth-integrated transport [L T-1 ~> m s-1]. @@ -104,6 +112,168 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhbt, v end subroutine continuity +!> Finds the thickness fluxes from the continuity solver without actually updating the +!! layer thicknesses. Because the fluxes in the two directions are calculated based on the +!! input thicknesses, which are not updated between the direcitons, the fluxes returned here +!! are not the same as those that would be returned by a call to continuity. +subroutine continuity_3d_fluxes(u, v, h, uh, vh, dt, G, GV, US, CS, OBC, pbv) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(out) :: uh !< Thickness fluxes through zonal faces, + !! u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(out) :: vh !< Thickness fluxes through meridional faces, + !! v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1]. + real, intent(in) :: dt !< Time increment [T ~> s]. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(continuity_CS), intent(in) :: CS !< Control structure for mom_continuity. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. + type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics + + if (CS%continuity_scheme == PPM_SCHEME) then + call zonal_mass_flux(u, h, uh, dt, G, GV, US, CS%PPM, OBC, pbv%por_face_areaU) + + call meridional_mass_flux(v, h, vh, dt, G, GV, US, CS%PPM, OBC, pbv%por_face_areaV) + else + call MOM_error(FATAL, "continuity: Unrecognized value of continuity_scheme") + endif + +end subroutine continuity_3d_fluxes + +!> Find the vertical sum of the thickness fluxes from the continuity solver without actually +!! updating the layer thicknesses. Because the fluxes in the two directions are calculated +!! based on the input thicknesses, which are not updated between the directions, the fluxes +!! returned here are not the same as those that would be returned by a call to continuity. +subroutine continuity_2d_fluxes(u, v, h, uhbt, vhbt, dt, G, GV, US, CS, OBC, pbv) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. + real, dimension(SZIB_(G),SZJ_(G)), & + intent(out) :: uhbt !< Vertically summed thickness flux through + !! zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1]. + real, dimension(SZI_(G),SZJB_(G)), & + intent(out) :: vhbt !< Vertically summed thickness flux through + !! meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1]. + real, intent(in) :: dt !< Time increment [T ~> s]. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(continuity_CS), intent(in) :: CS !< Control structure for mom_continuity. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. + type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics + + ! Local variables + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uh ! Thickness fluxes through zonal faces, + ! u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vh ! Thickness fluxes through meridional faces, + ! v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1]. + integer :: i, j, k + + uh(:,:,:) = 0.0 + vh(:,:,:) = 0.0 + + if (CS%continuity_scheme == PPM_SCHEME) then + call zonal_mass_flux(u, h, uh, dt, G, GV, US, CS%PPM, OBC, pbv%por_face_areaU) + + call meridional_mass_flux(v, h, vh, dt, G, GV, US, CS%PPM, OBC, pbv%por_face_areaV) + else + call MOM_error(FATAL, "continuity: Unrecognized value of continuity_scheme") + endif + + uhbt(:,:) = 0.0 + vhbt(:,:) = 0.0 + + do k=1,GV%ke ; do j=G%jsc,G%jec ; do I=G%isc-1,G%iec + uhbt(I,j) = uhbt(I,j) + uh(I,j,k) + enddo ; enddo ; enddo + + do k=1,GV%ke ; do J=G%jsc-1,G%jec ; do i=G%isc,G%iec + vhbt(I,j) = vhbt(I,j) + vh(I,j,k) + enddo ; enddo ; enddo + +end subroutine continuity_2d_fluxes + + +!> Correct the velocities to give the specified depth-integrated transports by applying a +!! barotropic acceleration (subject to viscous drag) to the velocities. +subroutine continuity_adjust_vel(u, v, h, dt, G, GV, US, CS, OBC, pbv, uhbt, vhbt, visc_rem_u, visc_rem_v) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: u !< Zonal velocity, which will be adjusted to + !! give uhbt as the depth-integrated + !! transport [L T-1 ~> m s-1]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(inout) :: v !< Meridional velocity, which will be adjusted + !! to give vhbt as the depth-integrated + !! transport [L T-1 ~> m s-1]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. + real, intent(in) :: dt !< Time increment [T ~> s]. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(continuity_CS), intent(in) :: CS !< Control structure for mom_continuity. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. + type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics + real, dimension(SZIB_(G),SZJ_(G)), & + intent(in) :: uhbt !< The vertically summed thickness flux through + !! zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1]. + real, dimension(SZI_(G),SZJB_(G)), & + intent(in) :: vhbt !< The vertically summed thickness flux through + !! meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1]. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + optional, intent(in) :: visc_rem_u !< Both the fraction of the zonal momentum + !! that remains after a time-step of viscosity, and + !! the fraction of a time-step's worth of a barotropic + !! acceleration that a layer experiences after viscosity + !! is applied [nondim]. This goes between 0 (at the + !! bottom) and 1 (far above the bottom). When this + !! column is under an ice shelf, this also goes to 0 + !! at the top due to the no-slip boundary condition there. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + optional, intent(in) :: visc_rem_v !< Both the fraction of the meridional momentum + !! that remains after a time-step of viscosity, and + !! the fraction of a time-step's worth of a barotropic + !! acceleration that a layer experiences after viscosity + !! is applied [nondim]. This goes between 0 (at the + !! bottom) and 1 (far above the bottom). When this + !! column is under an ice shelf, this also goes to 0 + !! at the top due to the no-slip boundary condition there. + + ! Local variables + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_in !< Input zonal velocity [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_in !< Input meridional velocity [L T-1 ~> m s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uh !< Volume flux through zonal faces = + !! u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vh !< Volume flux through meridional faces = + !! v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1]. + + ! It might not be necessary to separate the input velocity array from the adjusted velocities, + ! but it seems safer to do so, even if it might be less efficient. + u_in(:,:,:) = u(:,:,:) + v_in(:,:,:) = v(:,:,:) + + if (CS%continuity_scheme == PPM_SCHEME) then + call zonal_mass_flux(u_in, h, uh, dt, G, GV, US, CS%PPM, OBC, pbv%por_face_areaU, & + uhbt=uhbt, visc_rem_u=visc_rem_u, u_cor=u) + + call meridional_mass_flux(v_in, h, vh, dt, G, GV, US, CS%PPM, OBC, pbv%por_face_areaV, & + vhbt=vhbt, visc_rem_v=visc_rem_v, v_cor=v) + else + call MOM_error(FATAL, "continuity: Unrecognized value of continuity_scheme") + endif + +end subroutine continuity_adjust_vel + !> Initializes continuity_cs subroutine continuity_init(Time, G, GV, US, param_file, diag, CS) type(time_type), target, intent(in) :: Time !< Current model time. diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 73c6503242..7f96e675cd 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -19,6 +19,7 @@ module MOM_continuity_PPM #include public continuity_PPM, continuity_PPM_init, continuity_PPM_stencil +public zonal_mass_flux, meridional_mass_flux !>@{ CPU time clock IDs integer :: id_clock_update, id_clock_correct @@ -72,7 +73,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) - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. + 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)), & intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. @@ -147,8 +148,8 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhb ! First, advect zonally. LB%ish = G%isc ; LB%ieh = G%iec LB%jsh = G%jsc-stencil ; LB%jeh = G%jec+stencil - call zonal_mass_flux(u, hin, uh, dt, G, GV, US, CS, LB, OBC, & - pbv%por_face_areaU, uhbt, visc_rem_u, u_cor, BT_cont) + call zonal_mass_flux(u, hin, uh, dt, G, GV, US, CS, OBC, pbv%por_face_areaU, & + LB, uhbt, visc_rem_u, u_cor, BT_cont) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) @@ -163,8 +164,8 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhb ! Now advect meridionally, using the updated thicknesses to determine ! the fluxes. - call meridional_mass_flux(v, h, vh, dt, G, GV, US, CS, LB, OBC, & - pbv%por_face_areaV, vhbt, visc_rem_v, v_cor, BT_cont) + call meridional_mass_flux(v, h, vh, dt, G, GV, US, CS, OBC, pbv%por_face_areaV, & + LB, vhbt, visc_rem_v, v_cor, BT_cont) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) @@ -180,8 +181,8 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhb LB%ish = G%isc-stencil ; LB%ieh = G%iec+stencil LB%jsh = G%jsc ; LB%jeh = G%jec - call meridional_mass_flux(v, hin, vh, dt, G, GV, US, CS, LB, OBC, & - pbv%por_face_areaV, vhbt, visc_rem_v, v_cor, BT_cont) + call meridional_mass_flux(v, hin, vh, dt, G, GV, US, CS, OBC, pbv%por_face_areaV, & + LB, vhbt, visc_rem_v, v_cor, BT_cont) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) @@ -193,8 +194,8 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhb ! Now advect zonally, using the updated thicknesses to determine ! the fluxes. LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec - call zonal_mass_flux(u, h, uh, dt, G, GV, US, CS, LB, OBC, & - pbv%por_face_areaU, uhbt, visc_rem_u, u_cor, BT_cont) + call zonal_mass_flux(u, h, uh, dt, G, GV, US, CS, OBC, pbv%por_face_areaU, & + LB, uhbt, visc_rem_u, u_cor, BT_cont) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) @@ -210,9 +211,9 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhb end subroutine continuity_PPM !> Calculates the mass or volume fluxes through the zonal faces, and other related quantities. -subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, por_face_areaU, uhbt, & - visc_rem_u, u_cor, BT_cont) - type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. +subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, OBC, por_face_areaU, & + LB_in, uhbt, visc_rem_u, u_cor, BT_cont) + 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)), & intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. @@ -224,10 +225,11 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, por_face_are real, intent(in) :: dt !< Time increment [T ~> s]. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(continuity_PPM_CS), intent(in) :: CS !< This module's control structure. - type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure. type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. real, dimension(SZIB_(G), SZJ_(G), SZK_(G)), & intent(in) :: por_face_areaU !< fractional open area of U-faces [nondim] + type(loop_bounds_type), & + optional, intent(in) :: LB_in !< Loop bounds structure. real, dimension(SZIB_(G),SZJ_(G)), & optional, intent(in) :: uhbt !< The summed volume flux through zonal faces !! [H L2 T-1 ~> m3 s-1 or kg s-1]. @@ -265,6 +267,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, por_face_are real :: I_dt ! 1.0 / dt [T-1 ~> s-1]. real :: du_lim ! The velocity change that give a relative CFL of 1 [L T-1 ~> m s-1]. real :: dx_E, dx_W ! Effective x-grid spacings to the east and west [L ~> m]. + type(loop_bounds_type) :: LB integer :: i, j, k, ish, ieh, jsh, jeh, n, nz integer :: l_seg logical :: local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC @@ -279,6 +282,12 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, por_face_are local_Flather_OBC = OBC%Flather_u_BCs_exist_globally local_open_BC = OBC%open_u_BCs_exist_globally endif ; endif + + if (present(LB_in)) then + LB = LB_in + else + LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec + endif ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke CFL_dt = CS%CFL_limit_adjust / dt @@ -523,7 +532,7 @@ end subroutine zonal_mass_flux !> Evaluates the zonal mass or volume fluxes in a layer. subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, US, j, & ish, ieh, do_I, vol_CFL, por_face_areaU, OBC) - type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. + type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure. real, dimension(SZIB_(G)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. real, dimension(SZIB_(G)), intent(in) :: visc_rem !< Both the fraction of the !! momentum originally in a layer that remains after a time-step @@ -606,7 +615,7 @@ end subroutine zonal_flux_layer !! back these thicknesses to account for viscosity and fractional open areas. subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL, & marginal, OBC, por_face_areaU, visc_rem_u) - type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. + 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)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness used to @@ -617,7 +626,7 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL, !! reconstruction [H ~> m or kg m-2]. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h_u !< Effective thickness at zonal faces, !! scaled down to account for the effects of - !! viscoity and the fractional open area + !! viscosity and the fractional open area !! [H ~> m or kg m-2]. real, intent(in) :: dt !< Time increment [T ~> s]. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -680,6 +689,11 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL, do k=1,nz ; do j=jsh,jeh ; do I=ish-1,ieh h_u(I,j,k) = h_u(I,j,k) * (visc_rem_u(I,j,k) * por_face_areaU(I,j,k)) enddo ; enddo ; enddo + else + !$OMP parallel do default(shared) + do k=1,nz ; do j=jsh,jeh ; do I=ish-1,ieh + h_u(I,j,k) = h_u(I,j,k) * por_face_areaU(I,j,k) + enddo ; enddo ; enddo endif local_open_BC = .false. @@ -695,7 +709,7 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL, enddo enddo ; else ; do k=1,nz do j = OBC%segment(n)%HI%jsd, OBC%segment(n)%HI%jed - h_u(I,j,k) = h(i,j,k) + h_u(I,j,k) = h(i,j,k) * por_face_areaU(I,j,k) enddo enddo ; endif else @@ -705,7 +719,7 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL, enddo enddo ; else ; do k=1,nz do j = OBC%segment(n)%HI%jsd, OBC%segment(n)%HI%jed - h_u(I,j,k) = h(i+1,j,k) + h_u(I,j,k) = h(i+1,j,k) * por_face_areaU(I,j,k) enddo enddo ; endif endif @@ -721,7 +735,7 @@ subroutine zonal_flux_adjust(u, h_in, h_L, h_R, uhbt, uh_tot_0, duhdu_tot_0, & du, du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, & j, ish, ieh, do_I_in, por_face_areaU, uh_3d, OBC) - type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. + 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)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to @@ -873,7 +887,7 @@ end subroutine zonal_flux_adjust subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, & du_max_CFL, du_min_CFL, dt, G, GV, US, CS, visc_rem, & visc_rem_max, j, ish, ieh, do_I, por_face_areaU) - type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. + 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)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to @@ -1036,9 +1050,9 @@ subroutine set_zonal_BT_cont(u, h_in, h_L, h_R, BT_cont, uh_tot_0, duhdu_tot_0, 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, vh, dt, G, GV, US, CS, LB, OBC, por_face_areaV, vhbt, & - visc_rem_v, v_cor, BT_cont) - type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. +subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, OBC, por_face_areaV, & + LB_in, vhbt, visc_rem_v, v_cor, BT_cont) + 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] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to @@ -1048,11 +1062,12 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, por_fac real, intent(in) :: dt !< Time increment [T ~> s]. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(continuity_PPM_CS), intent(in) :: CS !< This module's control structure.G - type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure. type(ocean_OBC_type), pointer :: OBC !< Open boundary condition type !! specifies whether, where, and what open boundary conditions are used. real, dimension(SZI_(G),SZJB_(G),SZK_(G)), & intent(in) :: por_face_areaV !< fractional open area of V-faces [nondim] + type(loop_bounds_type), & + optional, intent(in) :: LB_in !< Loop bounds structure. real, dimension(SZI_(G),SZJB_(G)), optional, intent(in) :: vhbt !< The summed volume flux through !< meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1]. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & @@ -1090,6 +1105,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, por_fac real :: I_dt ! 1.0 / dt [T-1 ~> s-1]. real :: dv_lim ! The velocity change that give a relative CFL of 1 [L T-1 ~> m s-1]. real :: dy_N, dy_S ! Effective y-grid spacings to the north and south [L ~> m]. + type(loop_bounds_type) :: LB integer :: i, j, k, ish, ieh, jsh, jeh, n, nz integer :: l_seg logical :: local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC @@ -1104,6 +1120,12 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, por_fac local_Flather_OBC = OBC%Flather_v_BCs_exist_globally local_open_BC = OBC%open_v_BCs_exist_globally endif ; endif + + if (present(LB_in)) then + LB = LB_in + else + LB%ish = G%isc ; LB%ieh = G%iec ; LB%jsh = G%jsc ; LB%jeh = G%jec + endif ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke CFL_dt = CS%CFL_limit_adjust / dt @@ -1343,7 +1365,7 @@ end subroutine meridional_mass_flux !> Evaluates the meridional mass or volume fluxes in a layer. subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, & ish, ieh, do_I, vol_CFL, por_face_areaV, OBC) - type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. + type(ocean_grid_type), intent(in) :: G !< Ocean's grid structure. real, dimension(SZI_(G)), intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]. real, dimension(SZI_(G)), intent(in) :: visc_rem !< Both the fraction of the !! momentum originally in a layer that remains after a time-step @@ -1432,7 +1454,7 @@ end subroutine merid_flux_layer !! back these thicknesses to account for viscosity and fractional open areas. subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL, & marginal, OBC, por_face_areaV, visc_rem_v) - type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. + 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]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness used to calculate fluxes, @@ -1443,7 +1465,7 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL, !! [H ~> m or kg m-2]. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: h_v !< Effective thickness at meridional faces, !! scaled down to account for the effects of - !! viscoity and the fractional open area + !! viscosity and the fractional open area !! [H ~> m or kg m-2]. real, intent(in) :: dt !< Time increment [T ~> s]. type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure. @@ -1508,6 +1530,11 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL, do k=1,nz ; do J=jsh-1,jeh ; do i=ish,ieh h_v(i,J,k) = h_v(i,J,k) * (visc_rem_v(i,J,k) * por_face_areaV(i,J,k)) enddo ; enddo ; enddo + else + !$OMP parallel do default(shared) + do k=1,nz ; do J=jsh-1,jeh ; do i=ish,ieh + h_v(i,J,k) = h_v(i,J,k) * por_face_areaV(i,J,k) + enddo ; enddo ; enddo endif local_open_BC = .false. @@ -1523,7 +1550,7 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL, enddo enddo ; else ; do k=1,nz do i = OBC%segment(n)%HI%isd, OBC%segment(n)%HI%ied - h_v(i,J,k) = h(i,j,k) + h_v(i,J,k) = h(i,j,k) * por_face_areaV(i,J,k) enddo enddo ; endif else @@ -1533,7 +1560,7 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL, enddo enddo ; else ; do k=1,nz do i = OBC%segment(n)%HI%isd, OBC%segment(n)%HI%ied - h_v(i,J,k) = h(i,j+1,k) + h_v(i,J,k) = h(i,j+1,k) * por_face_areaV(i,J,k) enddo enddo ; endif endif @@ -1547,7 +1574,7 @@ end subroutine merid_face_thickness subroutine meridional_flux_adjust(v, h_in, h_L, h_R, vhbt, vh_tot_0, dvhdv_tot_0, & dv, dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, & j, ish, ieh, do_I_in, por_face_areaV, vh_3d, OBC) - type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. + 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]. @@ -1699,7 +1726,7 @@ end subroutine meridional_flux_adjust subroutine set_merid_BT_cont(v, h_in, h_L, h_R, BT_cont, vh_tot_0, dvhdv_tot_0, & dv_max_CFL, dv_min_CFL, dt, G, GV, US, CS, visc_rem, & visc_rem_max, j, ish, ieh, do_I, por_face_areaV) - type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. + 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]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to calculate fluxes,