From 4d9ed4fd1c2bf4dbb93022c8ec4ef63f157ae56c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 15 Oct 2021 15:42:00 -0400 Subject: [PATCH] +Make 37 optional arguments in src/core mandatory Made 37 optional arguments that are always present in calls into non-optional arguments to routines in the src/core directory. Many of these are pointer arguments related to things like open boundary conditions, so if these types are not to be used, they are simply not allocated. In several cases, this required the order of the arguments to be shifted around, but the various types of the arguments should prevent the model from compiling if the calls (e.g., in user-modified code that is not in the github repository) are not changed equivalently. Also eliminated 3 internal arguments in MOM_barotropic.F90 that are always hard-coded to the same values (the maximize argument to BT_cont_to_face_areas()) or are never used (the guess arguments to uhbt_to_ubt() and vhbt_to_vbt()), along with the code that they would exercise. All answers and output are bitwise identical. --- src/core/MOM.F90 | 11 +- src/core/MOM_PressureForce.F90 | 5 +- src/core/MOM_PressureForce_FV.F90 | 16 +- src/core/MOM_PressureForce_Montgomery.F90 | 22 +-- src/core/MOM_barotropic.F90 | 214 ++++++++-------------- src/core/MOM_continuity.F90 | 7 +- src/core/MOM_continuity_PPM.F90 | 88 +++++---- src/core/MOM_dynamics_split_RK2.F90 | 41 ++--- src/core/MOM_dynamics_unsplit.F90 | 12 +- src/core/MOM_dynamics_unsplit_RK2.F90 | 10 +- src/core/MOM_open_boundary.F90 | 2 +- src/core/MOM_transcribe_grid.F90 | 4 +- 12 files changed, 176 insertions(+), 256 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 73a7cd58a7..e8c770d247 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -3056,8 +3056,8 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS) type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: ssh !< time mean surface height [m] - real, dimension(:,:), optional, pointer :: p_atm !< Ocean surface pressure [R L2 T-2 ~> Pa] - logical, optional, intent(in) :: use_EOS !< If true, calculate the density for + real, dimension(:,:), pointer :: p_atm !< Ocean surface pressure [R L2 T-2 ~> Pa] + logical, intent(in) :: use_EOS !< If true, calculate the density for !! the SSH correction using the equation of state. real :: Rho_conv(SZI_(G)) ! The density used to convert surface pressure to @@ -3069,9 +3069,8 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec EOSdom(:) = EOS_domain(G%HI) - if (present(p_atm)) then ; if (associated(p_atm)) then - calc_rho = associated(tv%eqn_of_state) - if (present(use_EOS) .and. calc_rho) calc_rho = use_EOS + if (associated(p_atm)) then + calc_rho = use_EOS .and. associated(tv%eqn_of_state) ! Correct the output sea surface height for the contribution from the ice pressure. do j=js,je if (calc_rho) then @@ -3087,7 +3086,7 @@ subroutine adjust_ssh_for_p_atm(tv, G, GV, US, ssh, p_atm, use_EOS) enddo endif enddo - endif ; endif + endif end subroutine adjust_ssh_for_p_atm diff --git a/src/core/MOM_PressureForce.F90 b/src/core/MOM_PressureForce.F90 index 2316bb9725..dbc01dcc27 100644 --- a/src/core/MOM_PressureForce.F90 +++ b/src/core/MOM_PressureForce.F90 @@ -50,8 +50,7 @@ subroutine PressureForce(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, e intent(out) :: PFv !< Meridional pressure force acceleration [L T-2 ~> m s-2] type(PressureForce_CS), pointer :: CS !< Pressure force control structure type(ALE_CS), pointer :: ALE_CSp !< ALE control structure - real, dimension(:,:), & - optional, pointer :: p_atm !< The pressure at the ice-ocean or + real, dimension(:,:), pointer :: p_atm !< The pressure at the ice-ocean or !! atmosphere-ocean interface [R L2 T-2 ~> Pa]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & optional, intent(out) :: pbce !< The baroclinic pressure anomaly in each layer @@ -89,7 +88,7 @@ subroutine PressureForce_init(Time, G, GV, US, param_file, diag, CS, tides_CSp) type(param_file_type), intent(in) :: param_file !< Parameter file handles type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure type(PressureForce_CS), pointer :: CS !< Pressure force control structure - type(tidal_forcing_CS), optional, pointer :: tides_CSp !< Tide control structure + type(tidal_forcing_CS), pointer :: tides_CSp !< Tide control structure #include "version_variable.h" character(len=40) :: mdl = "MOM_PressureForce" ! This module's name. diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90 index 23e58272ed..1963d3f2c5 100644 --- a/src/core/MOM_PressureForce_FV.F90 +++ b/src/core/MOM_PressureForce_FV.F90 @@ -86,7 +86,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: PFv !< Meridional acceleration [L T-2 ~> m s-2] type(PressureForce_FV_CS), pointer :: CS !< Finite volume PGF control structure type(ALE_CS), pointer :: ALE_CSp !< ALE control structure - real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean + real, dimension(:,:), pointer :: p_atm !< The pressure at the ice-ocean !! or atmosphere-ocean interface [R L2 T-2 ~> Pa]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(out) :: pbce !< The baroclinic pressure !! anomaly in each layer due to eta anomalies @@ -167,8 +167,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_ "MOM_PressureForce_FV_nonBouss: The Stanley parameterization is not yet"//& "implemented in non-Boussinesq mode.") - use_p_atm = .false. - if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif + use_p_atm = associated(p_atm) use_EOS = associated(tv%eqn_of_state) use_ALE = .false. if (associated(ALE_CSp)) use_ALE = CS%reconstruct .and. use_EOS @@ -425,7 +424,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: PFv !< Meridional acceleration [L T-2 ~> m s-2] type(PressureForce_FV_CS), pointer :: CS !< Finite volume PGF control structure type(ALE_CS), pointer :: ALE_CSp !< ALE control structure - real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean + real, dimension(:,:), pointer :: p_atm !< The pressure at the ice-ocean !! or atmosphere-ocean interface [R L2 T-2 ~> Pa]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(out) :: pbce !< The baroclinic pressure !! anomaly in each layer due to eta anomalies @@ -501,8 +500,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm if (.not.associated(CS)) call MOM_error(FATAL, & "MOM_PressureForce_FV_Bouss: Module must be initialized before it is used.") - use_p_atm = .false. - if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif + use_p_atm = associated(p_atm) use_EOS = associated(tv%eqn_of_state) do i=Isq,Ieq+1 ; p0(i) = 0.0 ; enddo use_ALE = .false. @@ -808,7 +806,7 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS type(param_file_type), intent(in) :: param_file !< Parameter file handles type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure type(PressureForce_FV_CS), pointer :: CS !< Finite volume PGF control structure - type(tidal_forcing_CS), optional, pointer :: tides_CSp !< Tides control structure + type(tidal_forcing_CS), pointer :: tides_CSp !< Tides control structure ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl ! This module's name. @@ -821,9 +819,7 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS else ; allocate(CS) ; endif CS%diag => diag ; CS%Time => Time - if (present(tides_CSp)) then - if (associated(tides_CSp)) CS%tides_CSp => tides_CSp - endif + if (associated(tides_CSp)) CS%tides_CSp => tides_CSp mdl = "MOM_PressureForce_FV" call log_version(param_file, mdl, version, "") diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90 index 05e68aef12..e832f72158 100644 --- a/src/core/MOM_PressureForce_Montgomery.F90 +++ b/src/core/MOM_PressureForce_Montgomery.F90 @@ -71,7 +71,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: PFv !< Meridional acceleration due to pressure gradients !! (equal to -dM/dy) [L T-2 ~> m s-2]. type(PressureForce_Mont_CS), pointer :: CS !< Control structure for Montgomery potential PGF - real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean or + real, dimension(:,:), pointer :: p_atm !< The pressure at the ice-ocean or !! atmosphere-ocean [R L2 T-2 ~> Pa]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & optional, intent(out) :: pbce !< The baroclinic pressure anomaly in @@ -133,9 +133,8 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) - use_p_atm = .false. - if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif - is_split = .false. ; if (present(pbce)) is_split = .true. + use_p_atm = associated(p_atm) + is_split = present(pbce) use_EOS = associated(tv%eqn_of_state) if (.not.associated(CS)) call MOM_error(FATAL, & @@ -368,7 +367,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: PFv !< Meridional acceleration due to pressure gradients !! (equal to -dM/dy) [L T-2 ~> m s2]. type(PressureForce_Mont_CS), pointer :: CS !< Control structure for Montgomery potential PGF - real, dimension(:,:), optional, pointer :: p_atm !< The pressure at the ice-ocean or + real, dimension(:,:), pointer :: p_atm !< The pressure at the ice-ocean or !! atmosphere-ocean [R L2 T-2 ~> Pa]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), optional, intent(out) :: pbce !< The baroclinic pressure anomaly in !! each layer due to free surface height anomalies @@ -421,9 +420,8 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce, Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB EOSdom(1) = Isq - (G%isd-1) ; EOSdom(2) = G%iec+1 - (G%isd-1) - use_p_atm = .false. - if (present(p_atm)) then ; if (associated(p_atm)) use_p_atm = .true. ; endif - is_split = .false. ; if (present(pbce)) is_split = .true. + use_p_atm = associated(p_atm) + is_split = present(pbce) use_EOS = associated(tv%eqn_of_state) if (.not.associated(CS)) call MOM_error(FATAL, & @@ -826,8 +824,8 @@ subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, tides_ type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< Parameter file handles type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure - type(PressureForce_Mont_CS), pointer :: CS !< Montgomery PGF control structure - type(tidal_forcing_CS), optional, pointer :: tides_CSp !< Tides control structure + type(PressureForce_Mont_CS), pointer :: CS !< Montgomery PGF control structure + type(tidal_forcing_CS), pointer :: tides_CSp !< Tides control structure ! Local variables logical :: use_temperature, use_EOS @@ -842,9 +840,7 @@ subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, tides_ else ; allocate(CS) ; endif CS%diag => diag ; CS%Time => Time - if (present(tides_CSp)) then - if (associated(tides_CSp)) CS%tides_CSp => tides_CSp - endif + if (associated(tides_CSp)) CS%tides_CSp => tides_CSp mdl = "MOM_PressureForce_Mont" call log_version(param_file, mdl, version, "") diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index f912ff3275..131b7f705d 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -409,8 +409,8 @@ module MOM_barotropic subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, & eta_PF_in, U_Cor, V_Cor, accel_layer_u, accel_layer_v, & eta_out, uhbtav, vhbtav, G, GV, US, CS, & - visc_rem_u, visc_rem_v, etaav, ADp, OBC, BT_cont, eta_PF_start, & - taux_bot, tauy_bot, uh0, vh0, u_uh0, v_vh0) + visc_rem_u, visc_rem_v, ADp, OBC, BT_cont, eta_PF_start, & + taux_bot, tauy_bot, uh0, vh0, u_uh0, v_vh0, etaav) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -460,28 +460,28 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, !! viscosity is applied, in the zonal direction. Nondimensional !! between 0 (at the bottom) and 1 (far above the bottom). real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: visc_rem_v !< Ditto for meridional direction [nondim]. - real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: etaav !< The free surface height or column mass - !! averaged over the barotropic integration [H ~> m or kg m-2]. - type(accel_diag_ptrs), optional, pointer :: ADp !< Acceleration diagnostic pointers - type(ocean_OBC_type), optional, pointer :: OBC !< The open boundary condition structure. - type(BT_cont_type), optional, pointer :: BT_cont !< A structure with elements that describe + type(accel_diag_ptrs), pointer :: ADp !< Acceleration diagnostic pointers + type(ocean_OBC_type), pointer :: OBC !< The open boundary condition structure. + type(BT_cont_type), pointer :: BT_cont !< A structure with elements that describe !! the effective open face areas as a function of barotropic !! flow. - real, dimension(:,:), optional, pointer :: eta_PF_start !< The eta field consistent with the pressure + real, dimension(:,:), pointer :: eta_PF_start !< The eta field consistent with the pressure !! gradient at the start of the barotropic stepping !! [H ~> m or kg m-2]. - real, dimension(:,:), optional, pointer :: taux_bot !< The zonal bottom frictional stress from + real, dimension(:,:), pointer :: taux_bot !< The zonal bottom frictional stress from !! ocean to the seafloor [R L Z T-2 ~> Pa]. - real, dimension(:,:), optional, pointer :: tauy_bot !< The meridional bottom frictional stress + real, dimension(:,:), pointer :: tauy_bot !< The meridional bottom frictional stress !! from ocean to the seafloor [R L Z T-2 ~> Pa]. - real, dimension(:,:,:), optional, pointer :: uh0 !< The zonal layer transports at reference + real, dimension(:,:,:), pointer :: uh0 !< The zonal layer transports at reference !! velocities [H L2 T-1 ~> m3 s-1 or kg s-1]. - real, dimension(:,:,:), optional, pointer :: u_uh0 !< The velocities used to calculate + real, dimension(:,:,:), pointer :: u_uh0 !< The velocities used to calculate !! uh0 [L T-1 ~> m s-1] - real, dimension(:,:,:), optional, pointer :: vh0 !< The zonal layer transports at reference + real, dimension(:,:,:), pointer :: vh0 !< The zonal layer transports at reference !! velocities [H L2 T-1 ~> m3 s-1 or kg s-1]. - real, dimension(:,:,:), optional, pointer :: v_vh0 !< The velocities used to calculate + real, dimension(:,:,:), pointer :: v_vh0 !< The velocities used to calculate !! vh0 [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: etaav !< The free surface height or column mass + !! averaged over the barotropic integration [H ~> m or kg m-2]. ! Local variables real :: ubt_Cor(SZIB_(G),SZJ_(G)) ! The barotropic velocities that had been @@ -709,12 +709,10 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, Idt = 1.0 / dt accel_underflow = CS%vel_underflow * Idt - use_BT_cont = .false. - if (present(BT_cont)) use_BT_cont = (associated(BT_cont)) + use_BT_cont = associated(BT_cont) integral_BT_cont = use_BT_cont .and. CS%integral_BT_cont - interp_eta_PF = .false. - if (present(eta_PF_start)) interp_eta_PF = (associated(eta_PF_start)) + interp_eta_PF = associated(eta_PF_start) project_velocity = CS%BT_project_velocity @@ -728,11 +726,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, find_PF = (do_ave .and. ((CS%id_PFu_bt > 0) .or. (CS%id_PFv_bt > 0))) find_Cor = (do_ave .and. ((CS%id_Coru_bt > 0) .or. (CS%id_Corv_bt > 0))) - add_uh0 = .false. - if (present(uh0)) add_uh0 = associated(uh0) - if (add_uh0 .and. .not.(present(vh0) .and. present(u_uh0) .and. & - present(v_vh0))) call MOM_error(FATAL, & - "btstep: vh0, u_uh0, and v_vh0 must be present if uh0 is used.") + add_uh0 = associated(uh0) if (add_uh0 .and. .not.(associated(vh0) .and. associated(u_uh0) .and. & associated(v_vh0))) call MOM_error(FATAL, & "btstep: vh0, u_uh0, and v_vh0 must be associated if uh0 is used.") @@ -745,7 +739,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, apply_OBCs = .false. ; CS%BT_OBC%apply_u_OBCs = .false. ; CS%BT_OBC%apply_v_OBCs = .false. apply_OBC_open = .false. apply_OBC_flather = .false. - if (present(OBC)) then ; if (associated(OBC)) then + if (associated(OBC)) then CS%BT_OBC%apply_u_OBCs = OBC%open_u_BCs_exist_globally .or. OBC%specified_u_BCs_exist_globally CS%BT_OBC%apply_v_OBCs = OBC%open_v_BCs_exist_globally .or. OBC%specified_v_BCs_exist_globally apply_OBC_flather = open_boundary_query(OBC, apply_Flather_OBC=.true.) @@ -756,7 +750,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (apply_OBC_flather .and. .not.GV%Boussinesq) call MOM_error(FATAL, & "btstep: Flather open boundary conditions have not yet been "// & "implemented for a non-Boussinesq model.") - endif ; endif + endif num_cycles = 1 if (CS%use_wide_halos) & @@ -1074,9 +1068,9 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, call set_local_BT_cont_types(BT_cont, BTCL_u, BTCL_v, G, US, MS, CS%BT_Domain, 1+ievf-ie) else if (CS%Nonlinear_continuity) then - call find_face_areas(Datu, Datv, G, GV, US, CS, MS, eta, 1) + call find_face_areas(Datu, Datv, G, GV, US, CS, MS, 1, eta) else - call find_face_areas(Datu, Datv, G, GV, US, CS, MS, halo=1) + call find_face_areas(Datu, Datv, G, GV, US, CS, MS, 1) endif endif @@ -1130,10 +1124,10 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (integral_BT_cont) then call adjust_local_BT_cont_types(ubt, uhbt, vbt, vhbt, BTCL_u, BTCL_v, & - G, US, MS, halo=1+ievf-ie, dt_baroclinic=dt) + G, US, MS, 1+ievf-ie, dt_baroclinic=dt) else call adjust_local_BT_cont_types(ubt, uhbt, vbt, vhbt, BTCL_u, BTCL_v, & - G, US, MS, halo=1+ievf-ie) + G, US, MS, 1+ievf-ie) endif endif if (integral_BT_cont) then @@ -1279,17 +1273,15 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, else BT_force_v(i,J) = 0.0 endif ; enddo ; enddo - if (present(taux_bot) .and. present(tauy_bot)) then if (associated(taux_bot) .and. associated(tauy_bot)) then - !$OMP parallel do default(shared) - do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then - BT_force_u(I,j) = BT_force_u(I,j) - taux_bot(I,j) * mass_to_Z * CS%IDatu(I,j) - endif ; enddo ; enddo - !$OMP parallel do default(shared) - do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then - BT_force_v(i,J) = BT_force_v(i,J) - tauy_bot(i,J) * mass_to_Z * CS%IDatv(i,J) - endif ; enddo ; enddo - endif + !$OMP parallel do default(shared) + do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0.0) then + BT_force_u(I,j) = BT_force_u(I,j) - taux_bot(I,j) * mass_to_Z * CS%IDatu(I,j) + endif ; enddo ; enddo + !$OMP parallel do default(shared) + do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.0) then + BT_force_v(i,J) = BT_force_v(i,J) - tauy_bot(i,J) * mass_to_Z * CS%IDatv(i,J) + endif ; enddo ; enddo endif ! bc_accel_u & bc_accel_v are only available on the potentially @@ -1563,7 +1555,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, associated(forces%rigidity_ice_v)) H_min_dyn = GV%Z_to_H * CS%Dmin_dyn_psurf if (ice_is_rigid .and. use_BT_cont) & - call BT_cont_to_face_areas(BT_cont, Datu, Datv, G, US, MS, 0, .true.) + call BT_cont_to_face_areas(BT_cont, Datu, Datv, G, US, MS, halo=0) if (ice_is_rigid) then !$OMP parallel do default(shared) private(Idt_max2,H_eff_dx2,dyn_coef_max,ice_strength) do j=js,je ; do i=is,ie @@ -1776,7 +1768,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if ((.not.use_BT_cont) .and. CS%Nonlinear_continuity .and. & (CS%Nonlin_cont_update_period > 0)) then if ((n>1) .and. (mod(n-1,CS%Nonlin_cont_update_period) == 0)) & - call find_face_areas(Datu, Datv, G, GV, US, CS, MS, eta, 1+iev-ie) + call find_face_areas(Datu, Datv, G, GV, US, CS, MS, 1+iev-ie, eta) endif if (integral_BT_cont) then @@ -2681,33 +2673,33 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (CS%id_frhatv1 > 0) CS%frhatv1(:,:,:) = CS%frhatv(:,:,:) endif - if ((present(ADp)) .and. (associated(ADp%diag_hfrac_u))) then + if (associated(ADp%diag_hfrac_u)) then do k=1,nz ; do j=js,je ; do I=is-1,ie ADp%diag_hfrac_u(I,j,k) = CS%frhatu(I,j,k) enddo ; enddo ; enddo endif - if ((present(ADp)) .and. (associated(ADp%diag_hfrac_v))) then + if (associated(ADp%diag_hfrac_v)) then do k=1,nz ; do J=js-1,je ; do i=is,ie ADp%diag_hfrac_v(i,J,k) = CS%frhatv(i,J,k) enddo ; enddo ; enddo endif - if ((present(ADp)) .and. (present(BT_cont)) .and. (associated(ADp%diag_hu))) then + if (use_BT_cont .and. associated(ADp%diag_hu)) then do k=1,nz ; do j=js,je ; do I=is-1,ie ADp%diag_hu(I,j,k) = BT_cont%h_u(I,j,k) enddo ; enddo ; enddo endif - if ((present(ADp)) .and. (present(BT_cont)) .and. (associated(ADp%diag_hv))) then + if (use_BT_cont .and. associated(ADp%diag_hv)) then do k=1,nz ; do J=js-1,je ; do i=is,ie ADp%diag_hv(i,J,k) = BT_cont%h_v(i,J,k) enddo ; enddo ; enddo endif - if (present(ADp) .and. (associated(ADp%visc_rem_u))) then + if (associated(ADp%visc_rem_u)) then do k=1,nz ; do j=js,je ; do I=is-1,ie ADp%visc_rem_u(I,j,k) = visc_rem_u(I,j,k) enddo ; enddo ; enddo endif - if (present(ADp) .and. (associated(ADp%visc_rem_u))) then + if (associated(ADp%visc_rem_u)) then do k=1,nz ; do J=js-1,je ; do i=is,ie ADp%visc_rem_v(i,J,k) = visc_rem_v(i,J,k) enddo ; enddo ; enddo @@ -2790,11 +2782,11 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add) if (present(BT_cont)) use_BT_cont = (associated(BT_cont)) if (use_BT_cont) then - call BT_cont_to_face_areas(BT_cont, Datu, Datv, G, US, MS, 0, .true.) + call BT_cont_to_face_areas(BT_cont, Datu, Datv, G, US, MS, halo=0) elseif (CS%Nonlinear_continuity .and. present(eta)) then - call find_face_areas(Datu, Datv, G, GV, US, CS, MS, eta=eta, halo=0) + call find_face_areas(Datu, Datv, G, GV, US, CS, MS, 0, eta=eta) else - call find_face_areas(Datu, Datv, G, GV, US, CS, MS, halo=0, add_max=add_SSH) + call find_face_areas(Datu, Datv, G, GV, US, CS, MS, 0, add_max=add_SSH) endif det_de = 0.0 @@ -3593,7 +3585,7 @@ end function find_duhbt_du !> This function inverts the transport function to determine the barotopic !! velocity that is consistent with a given transport, or if INTEGRAL_BT_CONT=True !! this finds the time-integrated velocity that is consistent with a time-integrated transport. -function uhbt_to_ubt(uhbt, BTC, guess) result(ubt) +function uhbt_to_ubt(uhbt, BTC) result(ubt) real, intent(in) :: uhbt !< The barotropic zonal transport that should be inverted for, !! [H L2 T-1 ~> m3 s-1 or kg s-1] or the time-integrated !! transport [H L2 ~> m3 or kg]. @@ -3601,8 +3593,6 @@ function uhbt_to_ubt(uhbt, BTC, guess) result(ubt) !! barotropic transports to be calculated consistently with the !! layers' continuity equations. The dimensions of some !! of the elements in this type vary depending on INTEGRAL_BT_CONT. - real, optional, intent(in) :: guess !< A guess at what ubt will be [L T-1 ~> m s-1] or [L ~> m]. - !! The result is not allowed to be dramatically larger than guess. real :: ubt !< The result - The velocity that gives uhbt transport [L T-1 ~> m s-1] !! or the time-integrated velocity [L ~> m]. @@ -3676,18 +3666,6 @@ function uhbt_to_ubt(uhbt, BTC, guess) result(ubt) ubt = BTC%uBT_WW + (uhbt - BTC%uh_WW) / BTC%FA_u_WW endif - if (present(guess)) then - dvel = abs(ubt) - vs1*abs(guess) - if (dvel > 0.0) then ! Limit the velocity - if (dvel < 40.0 * (abs(guess)*(vs2-vs1)) ) then - vsr = vs2 - (vs2-vs1) * exp(-dvel / (abs(guess)*(vs2-vs1))) - else ! The exp is less than 4e-18 anyway in this case, so neglect it. - vsr = vs2 - endif - ubt = SIGN(vsr * guess, ubt) - endif - endif - end function uhbt_to_ubt !> The function find_vhbt determines the meridional transport for a given velocity, or with @@ -3742,7 +3720,7 @@ end function find_dvhbt_dv !> This function inverts the transport function to determine the barotopic !! velocity that is consistent with a given transport, or if INTEGRAL_BT_CONT=True !! this finds the time-integrated velocity that is consistent with a time-integrated transport. -function vhbt_to_vbt(vhbt, BTC, guess) result(vbt) +function vhbt_to_vbt(vhbt, BTC) result(vbt) real, intent(in) :: vhbt !< The barotropic meridional transport that should be !! inverted for [H L2 T-1 ~> m3 s-1 or kg s-1] or the !! time-integrated transport [H L2 ~> m3 or kg]. @@ -3750,8 +3728,6 @@ function vhbt_to_vbt(vhbt, BTC, guess) result(vbt) !! barotropic transports to be calculated consistently !! with the layers' continuity equations. The dimensions of some !! of the elements in this type vary depending on INTEGRAL_BT_CONT. - real, optional, intent(in) :: guess !< A guess at what vbt will be [L T-1 ~> m s-1] or [L ~> m]. - !! The result is not allowed to be dramatically larger than guess. real :: vbt !< The result - The velocity that gives vhbt transport [L T-1 ~> m s-1] !! or the time-integrated velocity [L ~> m]. @@ -3825,40 +3801,25 @@ function vhbt_to_vbt(vhbt, BTC, guess) result(vbt) vbt = BTC%vBT_SS + (vhbt - BTC%vh_SS) / BTC%FA_v_SS endif - if (present(guess)) then - dvel = abs(vbt) - vs1*abs(guess) - if (dvel > 0.0) then ! Limit the velocity - if (dvel < 40.0 * (abs(guess)*(vs2-vs1)) ) then - vsr = vs2 - (vs2-vs1) * exp(-dvel / (abs(guess)*(vs2-vs1))) - else ! The exp is less than 4e-18 anyway in this case, so neglect it. - vsr = vs2 - endif - vbt = SIGN(guess * vsr, vbt) - endif - endif - end function vhbt_to_vbt !> This subroutine sets up reordered versions of the BT_cont type in the !! local_BT_cont types, which have wide halos properly filled in. subroutine set_local_BT_cont_types(BT_cont, BTCL_u, BTCL_v, G, US, MS, BT_Domain, halo, dt_baroclinic) - type(BT_cont_type), intent(inout) :: BT_cont !< The BT_cont_type input to the - !! barotropic solver. - type(memory_size_type), intent(in) :: MS !< A type that describes the - !! memory sizes of the argument - !! arrays. - type(local_BT_cont_u_type), dimension(SZIBW_(MS),SZJW_(MS)), intent(out) :: BTCL_u !< A structure with the u - !! information from BT_cont. - type(local_BT_cont_v_type), dimension(SZIW_(MS),SZJBW_(MS)), intent(out) :: BTCL_v !< A structure with the v - !! information from BT_cont. - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(MOM_domain_type), intent(inout) :: BT_Domain !< The domain to use for updating - !! the halos of wide arrays. - integer, optional, intent(in) :: halo !< The extra halo size to use here. - real, optional, intent(in) :: dt_baroclinic !< The baroclinic time step - !! [T ~> s], which is provided if - !! INTEGRAL_BT_CONTINUITY is true. + type(BT_cont_type), intent(inout) :: BT_cont !< The BT_cont_type input to the barotropic solver + type(memory_size_type), intent(in) :: MS !< A type that describes the memory sizes of + !! the argument arrays + type(local_BT_cont_u_type), dimension(SZIBW_(MS),SZJW_(MS)), & + intent(out) :: BTCL_u !< A structure with the u information from BT_cont + type(local_BT_cont_v_type), dimension(SZIW_(MS),SZJBW_(MS)), & + intent(out) :: BTCL_v !< A structure with the v information from BT_cont + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(MOM_domain_type), intent(inout) :: BT_Domain !< The domain to use for updating the halos + !! of wide arrays + integer, intent(in) :: halo !< The extra halo size to use here + real, optional, intent(in) :: dt_baroclinic !< The baroclinic time step [T ~> s], which + !! is provided if INTEGRAL_BT_CONTINUITY is true. ! Local variables real, dimension(SZIBW_(MS),SZJW_(MS)) :: & @@ -3874,7 +3835,7 @@ subroutine set_local_BT_cont_types(BT_cont, BTCL_u, BTCL_v, G, US, MS, BT_Domain integer :: i, j, is, ie, js, je, hs is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - hs = 1 ; if (present(halo)) hs = max(halo,0) + hs = max(halo,0) dt = 1.0 ; if (present(dt_baroclinic)) dt = dt_baroclinic ! Copy the BT_cont arrays into symmetric, potentially wide haloed arrays. @@ -3999,7 +3960,7 @@ subroutine adjust_local_BT_cont_types(ubt, uhbt, vbt, vhbt, BTCL_u, BTCL_v, & intent(out) :: BTCL_v !< A structure with the v information from BT_cont. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - integer, optional, intent(in) :: halo !< The extra halo size to use here. + integer, intent(in) :: halo !< The extra halo size to use here. real, optional, intent(in) :: dt_baroclinic !< The baroclinic time step [T ~> s], which is !! provided if INTEGRAL_BT_CONTINUITY is true. @@ -4013,7 +3974,7 @@ subroutine adjust_local_BT_cont_types(ubt, uhbt, vbt, vhbt, BTCL_u, BTCL_v, & integer :: i, j, is, ie, js, je, hs is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - hs = 1 ; if (present(halo)) hs = max(halo,0) + hs = max(halo,0) dt = 1.0 ; if (present(dt_baroclinic)) dt = dt_baroclinic !$OMP parallel do default(shared) @@ -4079,9 +4040,9 @@ subroutine adjust_local_BT_cont_types(ubt, uhbt, vbt, vhbt, BTCL_u, BTCL_v, & end subroutine adjust_local_BT_cont_types -!> This subroutine uses the BTCL types to find typical or maximum face +!> This subroutine uses the BT_cont_type to find the maximum face !! areas, which can then be used for finding wave speeds, etc. -subroutine BT_cont_to_face_areas(BT_cont, Datu, Datv, G, US, MS, halo, maximize) +subroutine BT_cont_to_face_areas(BT_cont, Datu, Datv, G, US, MS, halo) type(BT_cont_type), intent(inout) :: BT_cont !< The BT_cont_type input to the !! barotropic solver. type(memory_size_type), intent(in) :: MS !< A type that describes the memory @@ -4091,35 +4052,22 @@ subroutine BT_cont_to_face_areas(BT_cont, Datu, Datv, G, US, MS, halo, maximize) real, dimension(MS%isdw:MS%iedw,MS%jsdw-1:MS%jedw), & intent(out) :: Datv !< The effective meridional face area [H L ~> m2 or kg m-1]. type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type integer, optional, intent(in) :: halo !< The extra halo size to use here. - logical, optional, intent(in) :: maximize !< If present and true, find the - !! maximum face area for any velocity. ! Local variables - logical :: find_max integer :: i, j, is, ie, js, je, hs is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec hs = 1 ; if (present(halo)) hs = max(halo,0) - find_max = .false. ; if (present(maximize)) find_max = maximize - if (find_max) then - do j=js-hs,je+hs ; do I=is-1-hs,ie+hs - Datu(I,j) = max(BT_cont%FA_u_EE(I,j), BT_cont%FA_u_E0(I,j), & - BT_cont%FA_u_W0(I,j), BT_cont%FA_u_WW(I,j)) - enddo ; enddo - do J=js-1-hs,je+hs ; do i=is-hs,ie+hs - Datv(i,J) = max(BT_cont%FA_v_NN(i,J), BT_cont%FA_v_N0(i,J), & - BT_cont%FA_v_S0(i,J), BT_cont%FA_v_SS(i,J)) - enddo ; enddo - else - do j=js-hs,je+hs ; do I=is-1-hs,ie+hs - Datu(I,j) = 0.5 * (BT_cont%FA_u_E0(I,j) + BT_cont%FA_u_W0(I,j)) - enddo ; enddo - do J=js-1-hs,je+hs ; do i=is-hs,ie+hs - Datv(i,J) = 0.5 * (BT_cont%FA_v_N0(i,J) + BT_cont%FA_v_S0(i,J)) - enddo ; enddo - endif + do j=js-hs,je+hs ; do I=is-1-hs,ie+hs + Datu(I,j) = max(BT_cont%FA_u_EE(I,j), BT_cont%FA_u_E0(I,j), & + BT_cont%FA_u_W0(I,j), BT_cont%FA_u_WW(I,j)) + enddo ; enddo + do J=js-1-hs,je+hs ; do i=is-hs,ie+hs + Datv(i,J) = max(BT_cont%FA_v_NN(i,J), BT_cont%FA_v_N0(i,J), & + BT_cont%FA_v_S0(i,J), BT_cont%FA_v_SS(i,J)) + enddo ; enddo end subroutine BT_cont_to_face_areas @@ -4133,7 +4081,7 @@ end subroutine swap !> This subroutine determines the open face areas of cells for calculating !! the barotropic transport. -subroutine find_face_areas(Datu, Datv, G, GV, US, CS, MS, eta, halo, add_max) +subroutine find_face_areas(Datu, Datv, G, GV, US, CS, MS, halo, eta, add_max) type(memory_size_type), intent(in) :: MS !< A type that describes the memory sizes of the argument arrays. real, dimension(MS%isdw-1:MS%iedw,MS%jsdw:MS%jedw), & intent(out) :: Datu !< The open zonal face area [H L ~> m2 or kg m-1]. @@ -4144,10 +4092,10 @@ subroutine find_face_areas(Datu, Datv, G, GV, US, CS, MS, eta, halo, add_max) type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(barotropic_CS), pointer :: CS !< The control structure returned by a previous !! call to barotropic_init. + integer, intent(in) :: halo !< The halo size to use, default = 1. real, dimension(MS%isdw:MS%iedw,MS%jsdw:MS%jedw), & optional, intent(in) :: eta !< The barotropic free surface height anomaly !! or column mass anomaly [H ~> m or kg m-2]. - integer, optional, intent(in) :: halo !< The halo size to use, default = 1. real, optional, intent(in) :: add_max !< A value to add to the maximum depth (used !! to overestimate the external wave speed) [Z ~> m]. @@ -4155,7 +4103,7 @@ subroutine find_face_areas(Datu, Datv, G, GV, US, CS, MS, eta, halo, add_max) real :: H1, H2 ! Temporary total thicknesses [H ~> m or kg m-2]. integer :: i, j, is, ie, js, je, hs is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - hs = 1 ; if (present(halo)) hs = max(halo,0) + hs = max(halo,0) !$OMP parallel default(none) shared(is,ie,js,je,hs,eta,GV,G,CS,Datu,Datv,add_max) & !$OMP private(H1,H2) @@ -4308,12 +4256,10 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure. logical, intent(out) :: calc_dtbt !< If true, the barotropic time step must !! be recalculated before stepping. - type(BT_cont_type), optional, & - pointer :: BT_cont !< A structure with elements that describe the + type(BT_cont_type), pointer :: BT_cont !< A structure with elements that describe the !! effective open face areas as a function of !! barotropic flow. - type(tidal_forcing_CS), optional, & - pointer :: tides_CSp !< A pointer to the control structure of the + type(tidal_forcing_CS), pointer :: tides_CSp !< A pointer to the control structure of the !! tide module. ! This include declares and sets the variable "version". @@ -4370,9 +4316,7 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, CS%module_is_initialized = .true. CS%diag => diag ; CS%Time => Time - if (present(tides_CSp)) then - if (associated(tides_CSp)) CS%tides_CSp => tides_CSp - endif + if (associated(tides_CSp)) CS%tides_CSp => tides_CSp ! Read all relevant parameters and write them to the model log. call get_param(param_file, mdl, "SPLIT", CS%split, default=.true., do_not_log=.true.) @@ -4986,7 +4930,7 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS, enddo ; enddo endif - call find_face_areas(Datu, Datv, G, GV, US, CS, MS, halo=1) + call find_face_areas(Datu, Datv, G, GV, US, CS, MS, 1) if ((CS%bound_BT_corr) .and. .not.(use_BT_Cont_type .and. CS%BT_cont_bounds)) then ! This is not used in most test cases. Were it ever to become more widely used, consider ! replacing maxvel with min(G%dxT(i,j),G%dyT(i,j)) * (CS%maxCFL_BT_cont*Idt) . diff --git a/src/core/MOM_continuity.F90 b/src/core/MOM_continuity.F90 index 480568c545..655055b03d 100644 --- a/src/core/MOM_continuity.F90 +++ b/src/core/MOM_continuity.F90 @@ -39,7 +39,7 @@ module MOM_continuity !> Time steps the layer thicknesses, using a monotonically limited, directionally split PPM scheme, !! based on Lin (1994). -subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, OBC, & +subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vhbt, & visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont) type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure. type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. @@ -60,14 +60,13 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, OBC, real, intent(in) :: dt !< Time increment [T ~> s]. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(continuity_CS), pointer :: CS !< Control structure for mom_continuity. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. 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]. 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]. - type(ocean_OBC_type), & - optional, pointer :: OBC !< Open boundaries control structure. 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 @@ -96,7 +95,7 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, OBC, " one must be present in call to continuity.") if (CS%continuity_scheme == PPM_SCHEME) then - call continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS%PPM_CSp, uhbt, vhbt, OBC, & + call continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS%PPM_CSp, OBC, uhbt, vhbt, & visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont=BT_cont) else call MOM_error(FATAL, "continuity: Unrecognized value of continuity_scheme") diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index d8b6cddaaa..d30e1af0f2 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -73,11 +73,10 @@ module MOM_continuity_PPM !> Time steps the layer thicknesses, using a monotonically limit, directionally split PPM scheme, !! based on Lin (1994). -subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, OBC, & +subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, 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(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. - type(continuity_PPM_CS), pointer :: CS !< Module's control 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)), & @@ -91,15 +90,15 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, O real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & intent(out) :: vh !< Meridional volume flux, 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(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(continuity_PPM_CS), pointer :: CS !< Module's control structure. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control 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]. 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]. - type(ocean_OBC_type), & - optional, pointer :: OBC !< Open boundaries control structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & optional, intent(in) :: visc_rem_u !< The fraction of zonal momentum originally @@ -149,7 +148,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, O ! 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, uhbt, OBC, visc_rem_u, u_cor, BT_cont) + call zonal_mass_flux(u, hin, uh, dt, G, GV, US, CS, LB, OBC, uhbt, visc_rem_u, u_cor, BT_cont) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) @@ -164,7 +163,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, O ! Now advect meridionally, using the updated thicknesses to determine ! the fluxes. - call meridional_mass_flux(v, h, vh, dt, G, GV, US, CS, LB, vhbt, OBC, visc_rem_v, v_cor, BT_cont) + call meridional_mass_flux(v, h, vh, dt, G, GV, US, CS, LB, OBC, vhbt, visc_rem_v, v_cor, BT_cont) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) @@ -180,7 +179,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, O 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, vhbt, OBC, visc_rem_v, v_cor, BT_cont) + call meridional_mass_flux(v, hin, vh, dt, G, GV, US, CS, LB, OBC, vhbt, visc_rem_v, v_cor, BT_cont) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) @@ -192,7 +191,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, O ! 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, uhbt, OBC, visc_rem_u, u_cor, BT_cont) + call zonal_mass_flux(u, h, uh, dt, G, GV, US, CS, LB, OBC, uhbt, visc_rem_u, u_cor, BT_cont) call cpu_clock_begin(id_clock_update) !$OMP parallel do default(shared) @@ -208,7 +207,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, O 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, uhbt, OBC, & +subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, OBC, uhbt, & visc_rem_u, u_cor, BT_cont) type(ocean_grid_type), intent(inout) :: G !< Ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Ocean's vertical grid structure. @@ -223,17 +222,16 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(continuity_PPM_CS), pointer :: CS !< This module's control structure. type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure. - type(ocean_OBC_type), & - optional, pointer :: OBC !< Open boundaries control structure. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control 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]. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & optional, intent(in) :: visc_rem_u !< The fraction of zonal momentum originally in a layer 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. !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). - 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]. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & optional, intent(out) :: u_cor !< The zonal velocitiess (u with a barotropic correction) @@ -272,7 +270,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & local_specified_BC = .false. ; set_BT_cont = .false. ; local_Flather_OBC = .false. local_open_BC = .false. if (present(BT_cont)) set_BT_cont = (associated(BT_cont)) - if (present(OBC)) then ; if (associated(OBC)) then + if (associated(OBC)) then ; if (OBC%OBC_pe) then local_specified_BC = OBC%specified_u_BCs_exist_globally local_Flather_OBC = OBC%Flather_u_BCs_exist_globally local_open_BC = OBC%open_u_BCs_exist_globally @@ -293,7 +291,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & enddo ; enddo else call PPM_reconstruction_x(h_in(:,:,k), h_L(:,:,k), h_R(:,:,k), G, LB, & - 2.0*GV%Angstrom_H, CS%monotonic, simple_2nd=CS%simple_2nd, OBC=OBC) + 2.0*GV%Angstrom_H, CS%monotonic, CS%simple_2nd, OBC) endif do I=ish-1,ieh ; visc_rem(I,k) = 1.0 ; enddo enddo @@ -484,7 +482,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & if (OBC%segment(n)%open .and. OBC%segment(n)%is_E_or_W) then I = OBC%segment(n)%HI%IsdB if (OBC%segment(n)%direction == OBC_DIRECTION_E) then - do J = OBC%segment(n)%HI%Jsd, OBC%segment(n)%HI%Jed + do j = OBC%segment(n)%HI%Jsd, OBC%segment(n)%HI%Jed FA_u = 0.0 do k=1,nz ; FA_u = FA_u + h_in(i,j,k)*G%dy_Cu(I,j) ; enddo BT_cont%FA_u_W0(I,j) = FA_u ; BT_cont%FA_u_E0(I,j) = FA_u @@ -492,7 +490,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & BT_cont%uBT_WW(I,j) = 0.0 ; BT_cont%uBT_EE(I,j) = 0.0 enddo else - do J = OBC%segment(n)%HI%Jsd, OBC%segment(n)%HI%Jed + do j = OBC%segment(n)%HI%Jsd, OBC%segment(n)%HI%Jed FA_u = 0.0 do k=1,nz ; FA_u = FA_u + h_in(i+1,j,k)*G%dy_Cu(I,j) ; enddo BT_cont%FA_u_W0(I,j) = FA_u ; BT_cont%FA_u_E0(I,j) = FA_u @@ -508,10 +506,10 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & if (set_BT_cont) then ; if (allocated(BT_cont%h_u)) then if (present(u_cor)) then call zonal_face_thickness(u_cor, h_in, h_L, h_R, BT_cont%h_u, dt, G, GV, US, LB, & - CS%vol_CFL, CS%marginal_faces, visc_rem_u, OBC) + CS%vol_CFL, CS%marginal_faces, OBC, visc_rem_u) else call zonal_face_thickness(u, h_in, h_L, h_R, BT_cont%h_u, dt, G, GV, US, LB, & - CS%vol_CFL, CS%marginal_faces, visc_rem_u, OBC) + CS%vol_CFL, CS%marginal_faces, OBC, visc_rem_u) endif endif ; endif @@ -601,7 +599,7 @@ end subroutine zonal_flux_layer !> Sets the effective interface thickness at each zonal velocity point. subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL, & - marginal, visc_rem_u, OBC) + marginal, OBC, visc_rem_u) type(ocean_grid_type), intent(inout) :: 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]. @@ -619,13 +617,13 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL, !! of face areas to the cell areas when estimating the CFL number. logical, intent(in) :: marginal !< If true, report the !! marginal face thicknesses; otherwise report transport-averaged thicknesses. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & optional, intent(in) :: visc_rem_u !< Both the fraction of the momentum originally in a layer 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. !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). - type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim] @@ -672,9 +670,7 @@ subroutine zonal_face_thickness(u, h, h_L, h_R, h_u, dt, G, GV, US, LB, vol_CFL, endif local_open_BC = .false. - if (present(OBC)) then ; if (associated(OBC)) then - local_open_BC = OBC%open_u_BCs_exist_globally - endif ; endif + if (associated(OBC)) local_open_BC = OBC%open_u_BCs_exist_globally if (local_open_BC) then do n = 1, OBC%number_of_segments if (OBC%segment(n)%open .and. OBC%segment(n)%is_E_or_W) then @@ -1022,12 +1018,12 @@ 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, vhbt, OBC, & +subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, OBC, vhbt, & visc_rem_v, v_cor, BT_cont) type(ocean_grid_type), intent(inout) :: 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 + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_in !< Layer thickness used to !! calculate fluxes [H ~> m or kg m-2] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: vh !< Volume flux through meridional !! faces = v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1] @@ -1035,16 +1031,16 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, & type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(continuity_PPM_CS), pointer :: CS !< This module's control structure.G type(loop_bounds_type), intent(in) :: LB !< Loop bounds structure. - type(ocean_OBC_type), optional, pointer :: OBC !< Open boundary condition type + 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)), 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)), & optional, intent(in) :: visc_rem_v !< Both the fraction of the momentum !! originally in a layer 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. Nondimensional between !! 0 (at the bottom) and 1 (far above the bottom). - 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)), & optional, intent(out) :: v_cor !< The meridional velocitiess (v with a barotropic correction) @@ -1084,11 +1080,11 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, & local_specified_BC = .false. ; set_BT_cont = .false. ; local_Flather_OBC = .false. local_open_BC = .false. if (present(BT_cont)) set_BT_cont = (associated(BT_cont)) - if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then + if (associated(OBC)) then ; if (OBC%OBC_pe) then local_specified_BC = OBC%specified_v_BCs_exist_globally local_Flather_OBC = OBC%Flather_v_BCs_exist_globally local_open_BC = OBC%open_v_BCs_exist_globally - endif ; endif ; endif + endif ; endif ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = GV%ke CFL_dt = CS%CFL_limit_adjust / dt @@ -1105,7 +1101,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, & enddo ; enddo else call PPM_reconstruction_y(h_in(:,:,k), h_L(:,:,k), h_R(:,:,k), G, LB, & - 2.0*GV%Angstrom_H, CS%monotonic, simple_2nd=CS%simple_2nd, OBC=OBC) + 2.0*GV%Angstrom_H, CS%monotonic, CS%simple_2nd, OBC) endif do i=ish,ieh ; visc_rem(i,k) = 1.0 ; enddo enddo @@ -1315,10 +1311,10 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, & if (set_BT_cont) then ; if (allocated(BT_cont%h_v)) then if (present(v_cor)) then call merid_face_thickness(v_cor, h_in, h_L, h_R, BT_cont%h_v, dt, G, GV, US, LB, & - CS%vol_CFL, CS%marginal_faces, visc_rem_v, OBC) + CS%vol_CFL, CS%marginal_faces, OBC, visc_rem_v) else call merid_face_thickness(v, h_in, h_L, h_R, BT_cont%h_v, dt, G, GV, US, LB, & - CS%vol_CFL, CS%marginal_faces, visc_rem_v, OBC) + CS%vol_CFL, CS%marginal_faces, OBC, visc_rem_v) endif endif ; endif @@ -1412,7 +1408,7 @@ end subroutine merid_flux_layer !> Sets the effective interface thickness at each meridional velocity point. subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL, & - marginal, visc_rem_v, OBC) + marginal, OBC, visc_rem_v) type(ocean_grid_type), intent(inout) :: 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]. @@ -1431,12 +1427,12 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL, !! of face areas to the cell areas when estimating the CFL number. logical, intent(in) :: marginal !< If true, report the marginal !! face thicknesses; otherwise report transport-averaged thicknesses. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), optional, intent(in) :: visc_rem_v !< Both the fraction !! of the momentum originally in a layer 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. !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom). - type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. ! Local variables real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim] @@ -1485,9 +1481,7 @@ subroutine merid_face_thickness(v, h, h_L, h_R, h_v, dt, G, GV, US, LB, vol_CFL, endif local_open_BC = .false. - if (present(OBC)) then ; if (associated(OBC)) then - local_open_BC = OBC%open_v_BCs_exist_globally - endif ; endif + if (associated(OBC)) local_open_BC = OBC%open_v_BCs_exist_globally if (local_open_BC) then do n = 1, OBC%number_of_segments if (OBC%segment(n)%open .and. OBC%segment(n)%is_N_or_S) then @@ -1848,7 +1842,7 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ logical, intent(in) :: simple_2nd !< If true, use the !! arithmetic mean thicknesses as the default edge values !! for a simple 2nd order scheme. - type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. ! Local variables with useful mnemonic names. real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes. @@ -1861,9 +1855,9 @@ subroutine PPM_reconstruction_x(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ type(OBC_segment_type), pointer :: segment => NULL() local_open_BC = .false. - if (present(OBC)) then ; if (associated(OBC)) then + if (associated(OBC)) then local_open_BC = OBC%open_u_BCs_exist_globally - endif ; endif + endif isl = LB%ish-1 ; iel = LB%ieh+1 ; jsl = LB%jsh ; jel = LB%jeh @@ -1983,7 +1977,7 @@ subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ logical, intent(in) :: simple_2nd !< If true, use the !! arithmetic mean thicknesses as the default edge values !! for a simple 2nd order scheme. - type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. + type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure. ! Local variables with useful mnemonic names. real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes. @@ -1996,9 +1990,9 @@ subroutine PPM_reconstruction_y(h_in, h_L, h_R, G, LB, h_min, monotonic, simple_ type(OBC_segment_type), pointer :: segment => NULL() local_open_BC = .false. - if (present(OBC)) then ; if (associated(OBC)) then + if (associated(OBC)) then local_open_BC = OBC%open_v_BCs_exist_globally - endif ; endif + endif isl = LB%ish ; iel = LB%ieh ; jsl = LB%jsh-1 ; jel = LB%jeh+1 diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 0532aeac53..51f12329a5 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -563,8 +563,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! u_accel_bt = layer accelerations due to barotropic solver if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then call cpu_clock_begin(id_clock_continuity) - call continuity(u, v, h, hp, uh_in, vh_in, dt, G, GV, US, CS%continuity_CSp, & - OBC=CS%OBC, visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont) + call continuity(u, v, h, hp, uh_in, vh_in, dt, G, GV, US, CS%continuity_CSp, CS%OBC, & + visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) if (BT_cont_BT_thick) then call btcalc(h, G, GV, CS%barotropic_CSp, CS%BT_cont%h_u, CS%BT_cont%h_v, & @@ -581,12 +581,10 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (calc_dtbt) call set_dtbt(G, GV, US, CS%barotropic_CSp, eta, CS%pbce) if (showCallTree) call callTree_enter("btstep(), MOM_barotropic.F90") ! This is the predictor step call to btstep. - call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, & - u_av, v_av, CS%u_accel_bt, CS%v_accel_bt, eta_pred, CS%uhbt, CS%vhbt, & - G, GV, US, CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, ADp=CS%ADp, & - OBC=CS%OBC, BT_cont=CS%BT_cont, eta_PF_start=eta_PF_start, & - taux_bot=taux_bot, tauy_bot=tauy_bot, & - uh0=uh_ptr, vh0=vh_ptr, u_uh0=u_ptr, v_vh0=v_ptr) + call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, & + CS%u_accel_bt, CS%v_accel_bt, eta_pred, CS%uhbt, CS%vhbt, G, GV, US, & + CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, CS%ADp, CS%OBC, CS%BT_cont, & + eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr) if (showCallTree) call callTree_leave("btstep()") call cpu_clock_end(id_clock_btstep) @@ -650,8 +648,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! uh = u_av * h ! hp = h + dt * div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(up, vp, h, hp, uh, vh, dt, G, GV, US, CS%continuity_CSp, & - CS%uhbt, CS%vhbt, CS%OBC, CS%visc_rem_u, CS%visc_rem_v, & + call continuity(up, vp, h, hp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, & + CS%uhbt, CS%vhbt, CS%visc_rem_u, CS%visc_rem_v, & u_av, v_av, BT_cont=CS%BT_cont) call cpu_clock_end(id_clock_continuity) if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)") @@ -782,13 +780,10 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (showCallTree) call callTree_enter("btstep(), MOM_barotropic.F90") ! This is the corrector step call to btstep. - call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, & - CS%eta_PF, u_av, v_av, CS%u_accel_bt, CS%v_accel_bt, & - eta_pred, CS%uhbt, CS%vhbt, G, GV, US, CS%barotropic_CSp, & - CS%visc_rem_u, CS%visc_rem_v, etaav=eta_av, ADp=CS%ADp, & - OBC=CS%OBC, BT_cont = CS%BT_cont, eta_PF_start=eta_PF_start, & - taux_bot=taux_bot, tauy_bot=tauy_bot, & - uh0=uh_ptr, vh0=vh_ptr, u_uh0=u_ptr, v_vh0=v_ptr) + call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, & + CS%u_accel_bt, CS%v_accel_bt, eta_pred, CS%uhbt, CS%vhbt, G, GV, US, & + CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, CS%ADp, CS%OBC, CS%BT_cont, & + eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr, etaav=eta_av) do j=js,je ; do i=is,ie ; eta(i,j) = eta_pred(i,j) ; enddo ; enddo call cpu_clock_end(id_clock_btstep) @@ -856,8 +851,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! h = h + dt * div . uh ! u_av and v_av adjusted so their mass transports match uhbt and vhbt. call cpu_clock_begin(id_clock_continuity) - call continuity(u, v, h, h, uh, vh, dt, G, GV, US, CS%continuity_CSp, & - CS%uhbt, CS%vhbt, CS%OBC, CS%visc_rem_u, CS%visc_rem_v, u_av, v_av) + call continuity(u, v, h, h, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, & + CS%uhbt, CS%vhbt, CS%visc_rem_u, CS%visc_rem_v, u_av, v_av) call cpu_clock_end(id_clock_continuity) call do_group_pass(CS%pass_h, G%Domain, clock=id_clock_pass) ! Whenever thickness changes let the diag manager know, target grids @@ -1272,8 +1267,6 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param !! diagnostic pointers type(VarMix_CS), pointer :: VarMix !< points to spatially variable viscosities type(MEKE_type), pointer :: MEKE !< points to mesoscale eddy kinetic energy fields -! type(Barotropic_CS), pointer :: Barotropic_CSp !< Pointer to the control structure for -! !! the barotropic module type(thickness_diffuse_CS), pointer :: thickness_diffuse_CSp !< Pointer to the control structure !! used for the isopycnal height diffusive transport. type(ocean_OBC_type), pointer :: OBC !< points to OBC related fields @@ -1286,7 +1279,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param !! the number of times the velocity is !! truncated (this should be 0). logical, intent(out) :: calc_dtbt !< If true, recalculate the barotropic time step - integer, optional, intent(out) :: cont_stencil !< The stencil for thickness + integer, intent(out) :: cont_stencil !< The stencil for thickness !! from the continuity solver. ! local variables @@ -1402,7 +1395,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param grain=CLOCK_ROUTINE) call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp) - if (present(cont_stencil)) cont_stencil = continuity_stencil(CS%continuity_CSp) + cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv_CSp) if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp) call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & @@ -1482,7 +1475,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param if (.not. query_initialized(uh,"uh",restart_CS) .or. & .not. query_initialized(vh,"vh",restart_CS)) then do k=1,nz ; do j=jsd,jed ; do i=isd,ied ; h_tmp(i,j,k) = h(i,j,k) ; enddo ; enddo ; enddo - call continuity(u, v, h, h_tmp, uh, vh, dt, G, GV, US, CS%continuity_CSp, OBC=CS%OBC) + call continuity(u, v, h, h_tmp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC) call pass_var(h_tmp, G%Domain, clock=id_clock_pass_init) do k=1,nz ; do j=jsd,jed ; do i=isd,ied CS%h_av(i,j,k) = 0.5*(h(i,j,k) + h_tmp(i,j,k)) diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 6f33a00768..48d767e1a8 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -265,7 +265,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! uh = u*h ! hp = h + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(u, v, h, hp, uh, vh, dt*0.5, G, GV, US, CS%continuity_CSp, OBC=CS%OBC) + call continuity(u, v, h, hp, uh, vh, dt*0.5, G, GV, US, CS%continuity_CSp, CS%OBC) call cpu_clock_end(id_clock_continuity) call pass_var(hp, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -355,7 +355,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! uh = up * hp ! h_av = hp + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(up, vp, hp, h_av, uh, vh, (0.5*dt), G, GV, US, CS%continuity_CSp, OBC=CS%OBC) + call continuity(up, vp, hp, h_av, uh, vh, (0.5*dt), G, GV, US, CS%continuity_CSp, CS%OBC) call cpu_clock_end(id_clock_continuity) call pass_var(h_av, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -415,7 +415,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & ! uh = upp * hp ! h = hp + dt/2 div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(upp, vpp, hp, h, uh, vh, (dt*0.5), G, GV, US, CS%continuity_CSp, OBC=CS%OBC) + call continuity(upp, vpp, hp, h, uh, vh, (dt*0.5), G, GV, US, CS%continuity_CSp, CS%OBC) call cpu_clock_end(id_clock_continuity) call pass_var(h, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -572,7 +572,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS type(MOM_dyn_unsplit_CS), pointer :: CS !< The control structure set up !! by initialize_dyn_unsplit. type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control - !!structure. + !! structure. type(accel_diag_ptrs), target, intent(inout) :: Accel_diag !< A set of pointers to the various !! accelerations in the momentum equations, which can be used !! for later derived diagnostics, like energy budgets. @@ -601,7 +601,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS integer, target, intent(inout) :: ntrunc !< A target for the variable that !! records the number of times the velocity !! is truncated (this should be 0). - integer, optional, intent(out) :: cont_stencil !< The stencil for thickness + integer, intent(out) :: cont_stencil !< The stencil for thickness !! from the continuity solver. ! This subroutine initializes all of the variables that are used by this @@ -653,7 +653,7 @@ subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, US, param_file, diag, CS Accel_diag%CAu => CS%CAu ; Accel_diag%CAv => CS%CAv call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp) - if (present(cont_stencil)) cont_stencil = continuity_stencil(CS%continuity_CSp) + cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv_CSp) if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp) call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index 18a192cb39..e6fec7f61e 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -281,7 +281,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, call cpu_clock_begin(id_clock_continuity) ! This is a duplicate calculation of the last continuity from the previous step ! and could/should be optimized out. -AJA - call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, GV, US, CS%continuity_CSp, OBC=CS%OBC) + call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, G, GV, US, CS%continuity_CSp, CS%OBC) call cpu_clock_end(id_clock_continuity) call pass_var(hp, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -350,7 +350,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! uh = up[n-1/2] * h[n-1/2] ! h_av = h + dt div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(up, vp, h_in, hp, uh, vh, dt, G, GV, US, CS%continuity_CSp, OBC=CS%OBC) + call continuity(up, vp, h_in, hp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC) call cpu_clock_end(id_clock_continuity) call pass_var(hp, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -405,7 +405,7 @@ subroutine step_MOM_dyn_unsplit_RK2(u_in, v_in, h_in, tv, visc, Time_local, dt, ! uh = up[n] * h[n] (up[n] might be extrapolated to damp GWs) ! h[n+1] = h[n] + dt div . uh call cpu_clock_begin(id_clock_continuity) - call continuity(up, vp, h_in, h_in, uh, vh, dt, G, GV, US, CS%continuity_CSp, OBC=CS%OBC) + call continuity(up, vp, h_in, h_in, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC) call cpu_clock_end(id_clock_continuity) call pass_var(h_in, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) @@ -547,7 +547,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag integer, target, intent(inout) :: ntrunc !< A target for the variable !! that records the number of times the !! velocity is truncated (this should be 0). - integer, optional, intent(out) :: cont_stencil !< The stencil for + integer, intent(out) :: cont_stencil !< The stencil for !! thickness from the continuity solver. ! This subroutine initializes all of the variables that are used by this @@ -615,7 +615,7 @@ subroutine initialize_dyn_unsplit_RK2(u, v, h, Time, G, GV, US, param_file, diag Accel_diag%CAu => CS%CAu ; Accel_diag%CAv => CS%CAv call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp) - if (present(cont_stencil)) cont_stencil = continuity_stencil(CS%continuity_CSp) + cont_stencil = continuity_stencil(CS%continuity_CSp) call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv_CSp) if (use_tides) call tidal_forcing_init(Time, G, param_file, CS%tides_CSp) call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, & diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index f0b1158b22..1601d6dd56 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -160,7 +160,7 @@ module MOM_open_boundary integer :: zphase_index !< Save where zphase is in segment%field. real :: Velocity_nudging_timescale_in !< Nudging timescale on inflow [T ~> s]. real :: Velocity_nudging_timescale_out !< Nudging timescale on outflow [T ~> s]. - logical :: on_pe !< true if segment is located in the computational domain + logical :: on_pe !< true if any portion of the segment is located in this PE's data domain logical :: temp_segment_data_exists !< true if temperature data arrays are present logical :: salt_segment_data_exists !< true if salinity data arrays are present real, pointer, dimension(:,:) :: Cg=>NULL() !< The external gravity wave speed [L T-1 ~> m s-1] diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 index a9626a805c..f176d6671c 100644 --- a/src/core/MOM_transcribe_grid.F90 +++ b/src/core/MOM_transcribe_grid.F90 @@ -170,7 +170,7 @@ end subroutine copy_dyngrid_to_MOM_grid subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) type(ocean_grid_type), intent(in) :: oG !< Ocean grid type type(dyn_horgrid_type), intent(inout) :: dG !< Common horizontal grid type - type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type integer :: isd, ied, jsd, jed ! Common data domains. integer :: IsdB, IedB, JsdB, JedB ! Common data domains. @@ -305,7 +305,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US) call pass_vector(dG%Dopen_u, dG%Dopen_v, dG%Domain, To_All+Scalar_Pair, CGRID_NE) endif - call set_derived_dyn_horgrid(dG, US) + call set_derived_dyn_horgrid(dG, US) end subroutine copy_MOM_grid_to_dyngrid