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