From 4c030e613a66e57ef6bdd14504fb5ae877af473c Mon Sep 17 00:00:00 2001 From: Hemant Khatri Date: Wed, 29 Jul 2020 08:44:37 -0400 Subject: [PATCH 1/4] Momentum budget terms multiplied by fractional layer-thicknesses (#1163) Squash merge: * New diagnostics for barotropic momentum budget calculations are created. In these, different budget terms multiplied by fractional layer thicknesses are saved. Thus, these terms can be added over the whole to obtain the barotropic momentum budget. 'btstep' subroutine in MOM_barotropic module is modified to return fractional thicknesses. Pressure force acceleration multiplied by fractional thickness as diagnostics are added in this commit. * More thickness weighted diagnostics * Implemented Andrew Shao's suggestions * More fractional thickness multiplied diagnostics * Some barotropic diagnostics obtained from fractional thickness multiplied momentum budget terms * Removed trailing spaces * All fractional-thickness multiplied diagnostics implemented * define diagnostic variables as pointers * Shorter initialization of 2D arrays * Modifications in vertical friction diagnostics * Diagnostics initialization as pointers to save memory allocation * Vertical friction module * Removed commented lines and trailing spaces * Diagnostic description change * Fractional thickness-weighted diagnostics for acceleration due to relative vorticity and gradient of kinetic energy * Modifications in vertical friction diagnostics (fractional thickness-weighted) * Modifications in diag_hfrac_u and diag_hfrac_v calls * Made allocation of hfrac at u/v points optional. These arrays are only allocated if any of the relevant diagnostics are called. * Fractional-thickness weighted 3D diagnostics are now commented out as we do not the proper grid remapping option. * Initialization of 2D diagnostic arrays changed from pointer type to allocatable array style. * Trailing spaces removed and line lengths reduced to be under 120 characters * Comment modified Reviewed by Robert Hallberg --- src/core/MOM_CoriolisAdv.F90 | 154 ++++++++++++++ src/core/MOM_barotropic.F90 | 15 +- src/core/MOM_dynamics_split_RK2.F90 | 193 +++++++++++++++++- src/core/MOM_variables.F90 | 2 + src/diagnostics/MOM_diagnostics.F90 | 93 +++++++++ .../lateral/MOM_hor_visc.F90 | 83 +++++++- .../vertical/MOM_vert_friction.F90 | 84 +++++++- 7 files changed, 614 insertions(+), 10 deletions(-) diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index f470338c4e..cf274d32a9 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -74,6 +74,10 @@ module MOM_CoriolisAdv !>@{ Diagnostic IDs integer :: id_rv = -1, id_PV = -1, id_gKEu = -1, id_gKEv = -1 integer :: id_rvxu = -1, id_rvxv = -1 + ! integer :: id_hf_gKEu = -1, id_hf_gKEv = -1 + integer :: id_hf_gKEu_2d = -1, id_hf_gKEv_2d = -1 + ! integer :: id_hf_rvxu = -1, id_hf_rvxv = -1 + integer :: id_hf_rvxu_2d = -1, id_hf_rvxv_2d = -1 !>@} end type CoriolisAdv_CS @@ -211,6 +215,16 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) real :: QUHeff,QVHeff ! More temporary variables [H L2 T-1 s-1 ~> m3 s-2 or kg s-2]. integer :: i, j, k, n, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz +! Diagnostics for fractional thickness-weighted terms + real, allocatable, dimension(:,:) :: & + hf_gKEu_2d, hf_gKEv_2d, & ! Depth sum of hf_gKEu, hf_gKEv [L T-2 ~> m s-2]. + hf_rvxu_2d, hf_rvxv_2d ! Depth sum of hf_rvxu, hf_rvxv [L T-2 ~> m s-2]. + !real, allocatable, dimension(:,:,:) :: & + ! hf_gKEu, hf_gKEv, & ! accel. due to KE gradient x fract. thickness [L T-2 ~> m s-2]. + ! hf_rvxu, hf_rvxv ! accel. due to RV x fract. thickness [L T-2 ~> m s-2]. + ! 3D diagnostics hf_gKEu etc. are commented because there is no clarity on proper remapping grid option. + ! The code is retained for degugging purposes in the future. + ! To work, the following fields must be set outside of the usual ! is to ie range before this subroutine is called: ! v(is-1:ie+2,js-1:je+1), u(is-1:ie+1,js-1:je+2), h(is-1:ie+2,js-1:je+2), @@ -828,6 +842,82 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) if (CS%id_gKEv>0) call post_data(CS%id_gKEv, AD%gradKEv, CS%diag) if (CS%id_rvxu > 0) call post_data(CS%id_rvxu, AD%rv_x_u, CS%diag) if (CS%id_rvxv > 0) call post_data(CS%id_rvxv, AD%rv_x_v, CS%diag) + + ! Diagnostics for terms multiplied by fractional thicknesses + + ! 3D diagnostics hf_gKEu etc. are commented because there is no clarity on proper remapping grid option. + ! The code is retained for degugging purposes in the future. + !if (CS%id_hf_gKEu > 0) then + ! allocate(hf_gKEu(G%IsdB:G%IedB,G%jsd:G%jed,G%ke)) + ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq + ! hf_gKEu(I,j,k) = AD%gradKEu(I,j,k) * AD%diag_hfrac_u(I,j,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_gKEu, hf_gKEu, CS%diag) + !endif + + !if (CS%id_hf_gKEv > 0) then + ! allocate(hf_gKEv(G%isd:G%ied,G%JsdB:G%JedB,G%ke)) + ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + ! hf_gKEv(i,J,k) = AD%gradKEv(i,J,k) * AD%diag_hfrac_v(i,J,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_gKEv, hf_gKEv, CS%diag) + !endif + + if (CS%id_hf_gKEu_2d > 0) then + allocate(hf_gKEu_2d(G%IsdB:G%IedB,G%jsd:G%jed)) + hf_gKEu_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + hf_gKEu_2d(I,j) = hf_gKEu_2d(I,j) + AD%gradKEu(I,j,k) * AD%diag_hfrac_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_gKEu_2d, hf_gKEu_2d, CS%diag) + deallocate(hf_gKEu_2d) + endif + + if (CS%id_hf_gKEv_2d > 0) then + allocate(hf_gKEv_2d(G%isd:G%ied,G%JsdB:G%JedB)) + hf_gKEv_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + hf_gKEv_2d(i,J) = hf_gKEv_2d(i,J) + AD%gradKEv(i,J,k) * AD%diag_hfrac_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_gKEv_2d, hf_gKEv_2d, CS%diag) + deallocate(hf_gKEv_2d) + endif + + !if (CS%id_hf_rvxv > 0) then + ! allocate(hf_rvxv(G%IsdB:G%IedB,G%jsd:G%jed,G%ke)) + ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq + ! hf_rvxv(I,j,k) = AD%rv_x_v(I,j,k) * AD%diag_hfrac_u(I,j,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_rvxv, hf_rvxv, CS%diag) + !endif + + !if (CS%id_hf_rvxu > 0) then + ! allocate(hf_rvxu(G%isd:G%ied,G%JsdB:G%JedB,G%ke)) + ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + ! hf_rvxu(i,J,k) = AD%rv_x_u(i,J,k) * AD%diag_hfrac_v(i,J,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_rvxu, hf_rvxu, CS%diag) + !endif + + if (CS%id_hf_rvxv_2d > 0) then + allocate(hf_rvxv_2d(G%IsdB:G%IedB,G%jsd:G%jed)) + hf_rvxv_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + hf_rvxv_2d(I,j) = hf_rvxv_2d(I,j) + AD%rv_x_v(I,j,k) * AD%diag_hfrac_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_rvxv_2d, hf_rvxv_2d, CS%diag) + deallocate(hf_rvxv_2d) + endif + + if (CS%id_hf_rvxu_2d > 0) then + allocate(hf_rvxu_2d(G%isd:G%ied,G%JsdB:G%JedB)) + hf_rvxu_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + hf_rvxu_2d(i,J) = hf_rvxu_2d(i,J) + AD%rv_x_u(i,J,k) * AD%diag_hfrac_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_rvxu_2d, hf_rvxu_2d, CS%diag) + deallocate(hf_rvxu_2d) + endif endif end subroutine CorAdCalc @@ -1087,6 +1177,70 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS) 'Zonal Acceleration from Relative Vorticity', 'm-1 s-2', conversion=US%L_T2_to_m_s2) if (CS%id_rvxv > 0) call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz) + !CS%id_hf_gKEu = register_diag_field('ocean_model', 'hf_gKEu', diag%axesCuL, Time, & + ! 'Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', & + ! 'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if (CS%id_hf_gKEu > 0) then + ! call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz) + ! call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + !endif + + !CS%id_hf_gKEv = register_diag_field('ocean_model', 'hf_gKEv', diag%axesCvL, Time, & + ! 'Fractional Thickness-weighted Meridional Acceleration from Grad. Kinetic Energy', & + ! 'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if (CS%id_hf_gKEv > 0) then + ! call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz) + ! call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + !endif + + CS%id_hf_gKEu_2d = register_diag_field('ocean_model', 'hf_gKEu_2d', diag%axesCuL, Time, & + 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', & + 'm-1 s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_hf_gKEu_2d > 0) then + call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + endif + + CS%id_hf_gKEv_2d = register_diag_field('ocean_model', 'hf_gKEv_2d', diag%axesCvL, Time, & + 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Grad. Kinetic Energy', & + 'm-1 s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_hf_gKEv_2d > 0) then + call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + endif + + !CS%id_hf_rvxu = register_diag_field('ocean_model', 'hf_rvxu', diag%axesCvL, Time, & + ! 'Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', & + ! 'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if (CS%id_hf_rvxu > 0) then + ! call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz) + ! call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + !endif + + !CS%id_hf_rvxv = register_diag_field('ocean_model', 'hf_rvxv', diag%axesCuL, Time, & + ! 'Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', & + ! 'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if (CS%id_hf_rvxv > 0) then + ! call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz) + ! call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + !endif + + CS%id_hf_rvxu_2d = register_diag_field('ocean_model', 'hf_rvxu_2d', diag%axesCvL, Time, & + 'Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', & + 'm-1 s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_hf_rvxu_2d > 0) then + call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + endif + + CS%id_hf_rvxv_2d = register_diag_field('ocean_model', 'hf_rvxv_2d', diag%axesCuL, Time, & + 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', & + 'm-1 s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_hf_rvxv_2d > 0) then + call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + endif + end subroutine CoriolisAdv_init !> Destructor for coriolisadv_cs diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 7f34f18998..db85bb40e2 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -27,6 +27,7 @@ module MOM_barotropic use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : BT_cont_type, alloc_bt_cont_type use MOM_verticalGrid, only : verticalGrid_type +use MOM_variables, only : accel_diag_ptrs implicit none ; private @@ -405,7 +406,7 @@ 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, OBC, BT_cont, eta_PF_start, & + visc_rem_u, visc_rem_v, etaav, ADp, OBC, BT_cont, eta_PF_start, & taux_bot, tauy_bot, uh0, vh0, u_uh0, v_vh0) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -458,6 +459,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, real, dimension(SZI_(G),SZJB_(G),SZK_(G)), 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 !! the effective open face areas as a function of barotropic @@ -2583,6 +2585,17 @@ 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 + 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 + 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 (G%nonblocking_updates) then if (find_etaav) call complete_group_pass(CS%pass_etaav, G%Domain) call complete_group_pass(CS%pass_ubta_uhbta, G%Domain) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 276c3c330f..64a9c18b97 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -160,10 +160,16 @@ module MOM_dynamics_split_RK2 integer :: id_umo_2d = -1, id_vmo_2d = -1 integer :: id_PFu = -1, id_PFv = -1 integer :: id_CAu = -1, id_CAv = -1 + ! integer :: id_hf_PFu = -1, id_hf_PFv = -1 + integer :: id_hf_PFu_2d = -1, id_hf_PFv_2d = -1 + ! integer :: id_hf_CAu = -1, id_hf_CAv = -1 + integer :: id_hf_CAu_2d = -1, id_hf_CAv_2d = -1 ! Split scheme only. integer :: id_uav = -1, id_vav = -1 integer :: id_u_BT_accel = -1, id_v_BT_accel = -1 + ! integer :: id_hf_u_BT_accel = -1, id_hf_v_BT_accel = -1 + integer :: id_hf_u_BT_accel_2d = -1, id_hf_v_BT_accel_2d = -1 !>@} type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the @@ -318,6 +324,19 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s u_av, & ! The zonal velocity time-averaged over a time step [L T-1 ~> m s-1]. v_av, & ! The meridional velocity time-averaged over a time step [L T-1 ~> m s-1]. h_av ! The layer thickness time-averaged over a time step [H ~> m or kg m-2]. + + ! real, allocatable, dimension(:,:,:) :: & + ! hf_PFu, hf_PFv, & ! Pressure force accel. x fract. thickness [L T-2 ~> m s-2]. + ! hf_CAu, hf_CAv, & ! Coriolis force accel. x fract. thickness [L T-2 ~> m s-2]. + ! hf_u_BT_accel, hf_v_BT_accel ! barotropic correction accel. x fract. thickness [L T-2 ~> m s-2]. + ! 3D diagnostics hf_PFu etc. are commented because there is no clarity on proper remapping grid option. + ! The code is retained for degugging purposes in the future. + + real, allocatable, dimension(:,:) :: & + hf_PFu_2d, hf_PFv_2d, & ! Depth integeral of hf_PFu, hf_PFv [L T-2 ~> m s-2]. + hf_CAu_2d, hf_CAv_2d, & ! Depth integeral of hf_CAu, hf_CAv [L T-2 ~> m s-2]. + hf_u_BT_accel_2d, hf_v_BT_accel_2d ! Depth integeral of hf_u_BT_accel, hf_v_BT_accel + real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s]. logical :: dyn_p_surf @@ -532,7 +551,7 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s ! 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, & + 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) @@ -682,7 +701,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s call cpu_clock_begin(id_clock_horvisc) call horizontal_viscosity(u_av, v_av, h_av, CS%diffu, CS%diffv, & MEKE, Varmix, G, GV, US, CS%hor_visc_CSp, & - OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp) + OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, & + ADp=CS%ADp) call cpu_clock_end(id_clock_horvisc) if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)") @@ -733,8 +753,8 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s 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, OBC=CS%OBC, & - BT_cont = CS%BT_cont, eta_PF_start=eta_PF_start, & + 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) do j=js,je ; do i=is,ie ; eta(i,j) = eta_pred(i,j) ; enddo ; enddo @@ -860,6 +880,109 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s if (CS%id_u_BT_accel > 0) call post_data(CS%id_u_BT_accel, CS%u_accel_bt, CS%diag) if (CS%id_v_BT_accel > 0) call post_data(CS%id_v_BT_accel, CS%v_accel_bt, CS%diag) + ! Diagnostics for terms multiplied by fractional thicknesses + + ! 3D diagnostics hf_PFu etc. are commented because there is no clarity on proper remapping grid option. + ! The code is retained for degugging purposes in the future. + !if (CS%id_hf_PFu > 0) then + ! allocate(hf_PFu(G%IsdB:G%IedB,G%jsd:G%jed,G%ke)) + ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq + ! hf_PFu(I,j,k) = CS%PFu(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_PFu, hf_PFu, CS%diag) + !endif + !if (CS%id_hf_PFv > 0) then + ! allocate(hf_PFv(G%isd:G%ied,G%JsdB:G%JedB,G%ke)) + ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + ! hf_PFv(i,J,k) = CS%PFv(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_PFv, hf_PFv, CS%diag) + !endif + if (CS%id_hf_PFu_2d > 0) then + allocate(hf_PFu_2d(G%IsdB:G%IedB,G%jsd:G%jed)) + hf_PFu_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + hf_PFu_2d(I,j) = hf_PFu_2d(I,j) + CS%PFu(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_PFu_2d, hf_PFu_2d, CS%diag) + deallocate(hf_PFu_2d) + endif + if (CS%id_hf_PFv_2d > 0) then + allocate(hf_PFv_2d(G%isd:G%ied,G%JsdB:G%JedB)) + hf_PFv_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + hf_PFv_2d(i,J) = hf_PFv_2d(i,J) + CS%PFv(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_PFv_2d, hf_PFv_2d, CS%diag) + deallocate(hf_PFv_2d) + endif + + !if (CS%id_hf_CAu > 0) then + ! allocate(hf_CAu(G%IsdB:G%IedB,G%jsd:G%jed,G%ke)) + ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq + ! hf_CAu(I,j,k) = CS%CAu(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_CAu, hf_CAu, CS%diag) + !endif + !if (CS%id_hf_CAv > 0) then + ! allocate(hf_CAv(G%isd:G%ied,G%JsdB:G%JedB,G%ke)) + ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + ! hf_CAv(i,J,k) = CS%CAv(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_CAv, hf_CAv, CS%diag) + !endif + if (CS%id_hf_CAu_2d > 0) then + allocate(hf_CAu_2d(G%IsdB:G%IedB,G%jsd:G%jed)) + hf_CAu_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + hf_CAu_2d(I,j) = hf_CAu_2d(I,j) + CS%CAu(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_CAu_2d, hf_CAu_2d, CS%diag) + deallocate(hf_CAu_2d) + endif + if (CS%id_hf_CAv_2d > 0) then + allocate(hf_CAv_2d(G%isd:G%ied,G%JsdB:G%JedB)) + hf_CAv_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + hf_CAv_2d(i,J) = hf_CAv_2d(i,J) + CS%CAv(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_CAv_2d, hf_CAv_2d, CS%diag) + deallocate(hf_CAv_2d) + endif + + !if (CS%id_hf_u_BT_accel > 0) then + ! allocate(hf_u_BT_accel(G%IsdB:G%IedB,G%jsd:G%jed,G%ke)) + ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq + ! hf_u_BT_accel(I,j,k) = CS%u_accel_bt(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_u_BT_accel, hf_u_BT_accel, CS%diag) + !endif + !if (CS%id_hf_v_BT_accel > 0) then + ! allocate(hf_v_BT_accel(G%isd:G%ied,G%JsdB:G%JedB,G%ke)) + ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + ! hf_v_BT_accel(i,J,k) = CS%v_accel_bt(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_v_BT_accel, hf_v_BT_accel, CS%diag) + !endif + if (CS%id_hf_u_BT_accel_2d > 0) then + allocate(hf_u_BT_accel_2d(G%IsdB:G%IedB,G%jsd:G%jed)) + hf_u_BT_accel_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + hf_u_BT_accel_2d(I,j) = hf_u_BT_accel_2d(I,j) + CS%u_accel_bt(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_u_BT_accel_2d, hf_u_BT_accel_2d, CS%diag) + deallocate(hf_u_BT_accel_2d) + endif + if (CS%id_hf_v_BT_accel_2d > 0) then + allocate(hf_v_BT_accel_2d(G%isd:G%ied,G%JsdB:G%JedB)) + hf_v_BT_accel_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + hf_v_BT_accel_2d(i,J) = hf_v_BT_accel_2d(i,J) + CS%v_accel_bt(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_v_BT_accel_2d, hf_v_BT_accel_2d, CS%diag) + deallocate(hf_v_BT_accel_2d) + endif + if (CS%debug) then call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym) call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s) @@ -1110,7 +1233,7 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param 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, & CS%tides_CSp) - call hor_visc_init(Time, G, US, param_file, diag, CS%hor_visc_CSp, MEKE) + call hor_visc_init(Time, G, US, param_file, diag, CS%hor_visc_CSp, MEKE, ADp=CS%ADp) call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, & ntrunc, CS%vertvisc_CSp) if (.not.associated(setVisc_CSp)) call MOM_error(FATAL, & @@ -1232,6 +1355,46 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param CS%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, Time, & 'Meridional Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) + !CS%id_hf_PFu = register_diag_field('ocean_model', 'hf_PFu', diag%axesCuL, Time, & + ! 'Fractional Thickness-weighted Zonal Pressure Force Acceleration', 'm s-2', & + ! v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if(CS%id_hf_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + + !CS%id_hf_PFv = register_diag_field('ocean_model', 'hf_PFv', diag%axesCvL, Time, & + ! 'Fractional Thickness-weighted Meridional Pressure Force Acceleration', 'm s-2', & + ! v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if(CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + + !CS%id_hf_CAu = register_diag_field('ocean_model', 'hf_CAu', diag%axesCuL, Time, & + ! 'Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', & + ! v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if(CS%id_hf_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + + !CS%id_hf_CAv = register_diag_field('ocean_model', 'hf_CAv', diag%axesCvL, Time, & + ! 'Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', 'm s-2', & + ! v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if(CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + + CS%id_hf_PFu_2d = register_diag_field('ocean_model', 'hf_PFu_2d', diag%axesCu1, Time, & + 'Depth-sum Fractional Thickness-weighted Zonal Pressure Force Acceleration', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if(CS%id_hf_PFu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + + CS%id_hf_PFv_2d = register_diag_field('ocean_model', 'hf_PFv_2d', diag%axesCv1, Time, & + 'Depth-sum Fractional Thickness-weighted Meridional Pressure Force Acceleration', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if(CS%id_hf_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + + CS%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, Time, & + 'Depth-sum Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if(CS%id_hf_CAu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + + CS%id_hf_CAv_2d = register_diag_field('ocean_model', 'hf_CAv_2d', diag%axesCv1, Time, & + 'Depth-sum Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if(CS%id_hf_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + CS%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, Time, & 'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=US%L_T_to_m_s) CS%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, Time, & @@ -1242,6 +1405,26 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param CS%id_v_BT_accel = register_diag_field('ocean_model', 'v_BT_accel', diag%axesCvL, Time, & 'Barotropic Anomaly Meridional Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) + !CS%id_hf_u_BT_accel = register_diag_field('ocean_model', 'hf_u_BT_accel', diag%axesCuL, Time, & + ! 'Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', 'm s-2', & + ! v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if(CS%id_hf_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + + !CS%id_hf_v_BT_accel = register_diag_field('ocean_model', 'hf_v_BT_accel', diag%axesCvL, Time, & + ! 'Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', 'm s-2', & + ! v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if(CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + + CS%id_hf_u_BT_accel_2d = register_diag_field('ocean_model', 'hf_u_BT_accel_2d', diag%axesCu1, Time, & + 'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if(CS%id_hf_u_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + + CS%id_hf_v_BT_accel_2d = register_diag_field('ocean_model', 'hf_v_BT_accel_2d', diag%axesCv1, Time, & + 'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if(CS%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE) id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE) id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=CLOCK_MODULE) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 26c2344f44..0b225f0bf7 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -185,6 +185,8 @@ module MOM_variables real, pointer :: rv_x_v(:,:,:) => NULL() !< rv_x_v = rv * v at u [L T-2 ~> m s-2] real, pointer :: rv_x_u(:,:,:) => NULL() !< rv_x_u = rv * u at v [L T-2 ~> m s-2] + real, pointer :: diag_hfrac_u(:,:,:) => NULL() !< Fractional layer thickness at u points + real, pointer :: diag_hfrac_v(:,:,:) => NULL() !< Fractional layer thickness at v points end type accel_diag_ptrs !> Pointers to arrays with transports, which can later be used for derived diagnostics, like energy balances. diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index d51173c16b..3936a788d0 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -71,6 +71,9 @@ module MOM_diagnostics dv_dt => NULL(), & !< net j-acceleration [L T-2 ~> m s-2] dh_dt => NULL(), & !< thickness rate of change [H T-1 ~> m s-1 or kg m-2 s-1] p_ebt => NULL() !< Equivalent barotropic modal structure [nondim] + ! hf_du_dt => NULL(), hf_dv_dt => NULL() !< du_dt, dv_dt x fract. thickness [L T-2 ~> m s-2]. + ! 3D diagnostics hf_du(dv)_dt are commented because there is no clarity on proper remapping grid option. + ! The code is retained for degugging purposes in the future. real, pointer, dimension(:,:,:) :: h_Rlay => NULL() !< Layer thicknesses in potential density !! coordinates [H ~> m or kg m-2] @@ -110,6 +113,8 @@ module MOM_diagnostics integer :: id_u = -1, id_v = -1, id_h = -1 integer :: id_e = -1, id_e_D = -1 integer :: id_du_dt = -1, id_dv_dt = -1 + ! integer :: id_hf_du_dt = -1, id_hf_dv_dt = -1 + integer :: id_hf_du_dt_2d = -1, id_hf_dv_dt_2d = -1 integer :: id_col_ht = -1, id_dh_dt = -1 integer :: id_KE = -1, id_dKEdt = -1 integer :: id_PE_to_KE = -1, id_KE_Coradv = -1 @@ -233,6 +238,9 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & real :: work_2d(SZI_(G),SZJ_(G)) ! A 2-d temporary work array. real :: rho_in_situ(SZI_(G)) ! In situ density [R ~> kg m-3] + real, allocatable, dimension(:,:) :: & + hf_du_dt_2d, hf_dv_dt_2d ! z integeral of hf_du_dt, hf_dv_dt [L T-2 ~> m s-2]. + ! tmp array for surface properties real :: surface_field(SZI_(G),SZJ_(G)) real :: pressure_1d(SZI_(G)) ! Temporary array for pressure when calling EOS [R L2 T-2 ~> Pa] @@ -272,6 +280,44 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & if (CS%id_dh_dt>0) call post_data(CS%id_dh_dt, CS%dh_dt, CS%diag, alt_h = diag_pre_sync%h_state) + !! Diagnostics for terms multiplied by fractional thicknesses + + ! 3D diagnostics hf_du(dv)_dt are commented because there is no clarity on proper remapping grid option. + ! The code is retained for degugging purposes in the future. + !if (CS%id_hf_du_dt > 0) then + ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq + ! CS%hf_du_dt(I,j,k) = CS%du_dt(I,j,k) * ADp%diag_hfrac_u(I,j,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_du_dt, CS%hf_du_dt, CS%diag, alt_h = diag_pre_sync%h_state) + !endif + + !if (CS%id_hf_dv_dt > 0) then + ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + ! CS%hf_dv_dt(i,J,k) = CS%dv_dt(i,J,k) * ADp%diag_hfrac_v(i,J,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_dv_dt, CS%hf_dv_dt, CS%diag, alt_h = diag_pre_sync%h_state) + !endif + + if (CS%id_hf_du_dt_2d > 0) then + allocate(hf_du_dt_2d(G%IsdB:G%IedB,G%jsd:G%jed)) + hf_du_dt_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + hf_du_dt_2d(I,j) = hf_du_dt_2d(I,j) + CS%du_dt(I,j,k) * ADp%diag_hfrac_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_du_dt_2d, hf_du_dt_2d, CS%diag) + deallocate(hf_du_dt_2d) + endif + + if (CS%id_hf_dv_dt_2d > 0) then + allocate(hf_dv_dt_2d(G%isd:G%ied,G%JsdB:G%JedB)) + hf_dv_dt_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + hf_dv_dt_2d(i,J) = hf_dv_dt_2d(i,J) + CS%dv_dt(i,J,k) * ADp%diag_hfrac_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_dv_dt_2d, hf_dv_dt_2d, CS%diag) + deallocate(hf_dv_dt_2d) + endif + call diag_restore_grids(CS%diag) call calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS) @@ -1644,6 +1690,50 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag call register_time_deriv(lbound(MIS%h), MIS%h, CS%dh_dt, CS) endif + !CS%id_hf_du_dt = register_diag_field('ocean_model', 'hf_dudt', diag%axesCuL, Time, & + ! 'Fractional Thickness-weighted Zonal Acceleration', 'm s-2', v_extensive=.true., & + ! conversion=US%L_T2_to_m_s2) + !if (CS%id_hf_du_dt > 0) then + ! call safe_alloc_ptr(CS%hf_du_dt,IsdB,IedB,jsd,jed,nz) + ! if (.not.associated(CS%du_dt)) then + ! call safe_alloc_ptr(CS%du_dt,IsdB,IedB,jsd,jed,nz) + ! call register_time_deriv(lbound(MIS%u), MIS%u, CS%du_dt, CS) + ! endif + ! call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + !endif + + !CS%id_hf_dv_dt = register_diag_field('ocean_model', 'hf_dvdt', diag%axesCvL, Time, & + ! 'Fractional Thickness-weighted Meridional Acceleration', 'm s-2', v_extensive=.true., & + ! conversion=US%L_T2_to_m_s2) + !if (CS%id_hf_dv_dt > 0) then + ! call safe_alloc_ptr(CS%hf_dv_dt,isd,ied,JsdB,JedB,nz) + ! if (.not.associated(CS%dv_dt)) then + ! call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz) + ! call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS) + ! endif + ! call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + !endif + + CS%id_hf_du_dt_2d = register_diag_field('ocean_model', 'hf_dudt_2d', diag%axesCu1, Time, & + 'Depth-sum Fractional Thickness-weighted Zonal Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_hf_du_dt_2d > 0) then + if (.not.associated(CS%du_dt)) then + call safe_alloc_ptr(CS%du_dt,IsdB,IedB,jsd,jed,nz) + call register_time_deriv(lbound(MIS%u), MIS%u, CS%du_dt, CS) + endif + call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + endif + + CS%id_hf_dv_dt_2d = register_diag_field('ocean_model', 'hf_dvdt_2d', diag%axesCv1, Time, & + 'Depth-sum Fractional Thickness-weighted Meridional Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2) + if (CS%id_hf_dv_dt_2d > 0) then + if (.not.associated(CS%dv_dt)) then + call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz) + call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS) + endif + call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + endif + ! layer thickness variables !if (GV%nk_rho_varies > 0) then CS%id_h_Rlay = register_diag_field('ocean_model', 'h_rho', diag%axesTL, Time, & @@ -2178,6 +2268,9 @@ subroutine MOM_diagnostics_end(CS, ADp) if (associated(ADp%du_other)) deallocate(ADp%du_other) if (associated(ADp%dv_other)) deallocate(ADp%dv_other) + if (associated(ADp%diag_hfrac_u)) deallocate(ADp%diag_hfrac_u) + if (associated(ADp%diag_hfrac_v)) deallocate(ADp%diag_hfrac_v) + do m=1,CS%num_time_deriv ; deallocate(CS%prev_val(m)%p) ; enddo deallocate(CS) diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 7c1405308f..4e9897f6eb 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -19,6 +19,7 @@ module MOM_hor_visc use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_NONE use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type +use MOM_variables, only : accel_diag_ptrs implicit none ; private @@ -174,9 +175,16 @@ module MOM_hor_visc type(diag_ctrl), pointer :: diag => NULL() !< structure to regulate diagnostics + ! real, pointer :: hf_diffu(:,:,:) => NULL() ! Zonal hor. visc. accel. x fract. thickness [L T-2 ~> m s-2]. + ! real, pointer :: hf_diffv(:,:,:) => NULL() ! Merdional hor. visc. accel. x fract. thickness [L T-2 ~> m s-2]. + ! 3D diagnostics hf_diffu(diffv) are commented because there is no clarity on proper remapping grid option. + ! The code is retained for degugging purposes in the future. + !>@{ !! Diagnostic id integer :: id_diffu = -1, id_diffv = -1 + ! integer :: id_hf_diffu = -1, id_hf_diffv = -1 + integer :: id_hf_diffu_2d = -1, id_hf_diffv_2d = -1 integer :: id_Ah_h = -1, id_Ah_q = -1 integer :: id_Kh_h = -1, id_Kh_q = -1 integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1 @@ -203,7 +211,7 @@ module MOM_hor_visc !! v[is-2:ie+2,js-2:je+2] !! h[is-1:ie+1,js-1:je+1] subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, & - CS, OBC, BT, TD) + CS, OBC, BT, TD, ADp) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), & @@ -230,6 +238,8 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, !! barotropic velocities. type(thickness_diffuse_CS), optional, pointer :: TD !< Pointer to a structure containing !! thickness diffusivities. + type(accel_diag_ptrs), optional, pointer :: ADp !< Acceleration diagnostic pointers + ! Local variables real, dimension(SZIB_(G),SZJ_(G)) :: & Del2u, & ! The u-compontent of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1] @@ -263,6 +273,9 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, grad_d2vel_mag_h, & ! Magnitude of the Laplacian of the velocity vector, squared [L-2 T-2 ~> m-2 s-2] boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim] + real, allocatable, dimension(:,:) :: hf_diffu_2d ! Depth sum of hf_diffu [L T-2 ~> m s-2] + real, allocatable, dimension(:,:) :: hf_diffv_2d ! Depth sum of hf_diffv [L T-2 ~> m s-2] + real, dimension(SZIB_(G),SZJB_(G)) :: & dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1] dDel2vdx, dDel2udy, & ! Components in the biharmonic equivalent of the shearing strain [L-2 T-1 ~> m-2 s-1] @@ -1307,12 +1320,47 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, call post_data(CS%id_FrictWorkIntz, FrictWorkIntz, CS%diag) endif + ! Diagnostics for terms multiplied by fractional thicknesses + + ! 3D diagnostics hf_diffu(diffv) are commented because there is no clarity on proper remapping grid option. + ! The code is retained for degugging purposes in the future. + !if (present(ADp) .and. (CS%id_hf_diffu > 0)) then + ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq + ! CS%hf_diffu(I,j,k) = diffu(I,j,k) * ADp%diag_hfrac_u(I,j,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_diffu, CS%hf_diffu, CS%diag) + !endif + !if (present(ADp) .and. (CS%id_hf_diffv > 0)) then + ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + ! CS%hf_diffv(i,J,k) = diffv(i,J,k) * ADp%diag_hfrac_v(i,J,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_diffv, CS%hf_diffv, CS%diag) + !endif + if (present(ADp) .and. (CS%id_hf_diffu_2d > 0)) then + allocate(hf_diffu_2d(G%IsdB:G%IedB,G%jsd:G%jed)) + hf_diffu_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + hf_diffu_2d(I,j) = hf_diffu_2d(I,j) + diffu(I,j,k) * ADp%diag_hfrac_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_diffu_2d, hf_diffu_2d, CS%diag) + deallocate(hf_diffu_2d) + endif + if (present(ADp) .and. (CS%id_hf_diffv_2d > 0)) then + allocate(hf_diffv_2d(G%isd:G%ied,G%JsdB:G%JedB)) + hf_diffv_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + hf_diffv_2d(i,J) = hf_diffv_2d(i,J) + diffv(i,J,k) * ADp%diag_hfrac_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_diffv_2d, hf_diffv_2d, CS%diag) + deallocate(hf_diffv_2d) + endif + end subroutine horizontal_viscosity !> Allocates space for and calculates static variables used by horizontal_viscosity(). !! hor_visc_init calculates and stores the values of a number of metric functions that !! are used in horizontal_viscosity(). -subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) +subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE, ADp) type(time_type), intent(in) :: Time !< Current model time. type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -1321,6 +1369,7 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) type(diag_ctrl), target, intent(inout) :: diag !< Structure to regulate diagnostic output. type(hor_visc_CS), pointer :: CS !< Pointer to the control structure for this module type(MEKE_type), pointer :: MEKE !< MEKE data + type(accel_diag_ptrs), optional, pointer :: ADp !< Acceleration diagnostic pointers ! Local variables real, dimension(SZIB_(G),SZJ_(G)) :: u0u, u0v real, dimension(SZI_(G),SZJB_(G)) :: v0u, v0v @@ -2016,6 +2065,36 @@ subroutine hor_visc_init(Time, G, US, param_file, diag, CS, MEKE) CS%id_diffv = register_diag_field('ocean_model', 'diffv', diag%axesCvL, Time, & 'Meridional Acceleration from Horizontal Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2) + !CS%id_hf_diffu = register_diag_field('ocean_model', 'hf_diffu', diag%axesCuL, Time, & + ! 'Fractional Thickness-weighted Zonal Acceleration from Horizontal Viscosity', 'm s-2', & + ! v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if ((CS%id_hf_diffu > 0) .and. (present(ADp))) then + ! call safe_alloc_ptr(CS%hf_diffu,G%IsdB,G%IedB,G%jsd,G%jed,G%ke) + ! call safe_alloc_ptr(ADp%diag_hfrac_u,G%IsdB,G%IedB,G%jsd,G%jed,G%ke) + !endif + + !CS%id_hf_diffv = register_diag_field('ocean_model', 'hf_diffv', diag%axesCvL, Time, & + ! 'Fractional Thickness-weighted Meridional Acceleration from Horizontal Viscosity', 'm s-2', & + ! v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if ((CS%id_hf_diffv > 0) .and. (present(ADp))) then + ! call safe_alloc_ptr(CS%hf_diffv,G%isd,G%ied,G%JsdB,G%JedB,G%ke) + ! call safe_alloc_ptr(ADp%diag_hfrac_v,G%isd,G%ied,G%JsdB,G%JedB,G%ke) + !endif + + CS%id_hf_diffu_2d = register_diag_field('ocean_model', 'hf_diffu_2d', diag%axesCu1, Time, & + 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Horizontal Viscosity', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if ((CS%id_hf_diffu_2d > 0) .and. (present(ADp))) then + call safe_alloc_ptr(ADp%diag_hfrac_u,G%IsdB,G%IedB,G%jsd,G%jed,G%ke) + endif + + CS%id_hf_diffv_2d = register_diag_field('ocean_model', 'hf_diffv_2d', diag%axesCv1, Time, & + 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Horizontal Viscosity', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if ((CS%id_hf_diffv_2d > 0) .and. (present(ADp))) then + call safe_alloc_ptr(ADp%diag_hfrac_v,G%isd,G%ied,G%JsdB,G%JedB,G%ke) + endif + if (CS%biharmonic) then CS%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, Time, & 'Biharmonic Horizontal Viscosity at h Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T, & diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index f03cee72b8..c6a6f37b39 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -124,10 +124,18 @@ module MOM_vert_friction integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1 integer :: id_Ray_u = -1, id_Ray_v = -1, id_taux_bot = -1, id_tauy_bot = -1 integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1 + ! integer :: id_hf_du_dt_visc = -1, id_hf_dv_dt_visc = -1 + integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1 !>@} type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() !< A pointer to the control structure !! for recording accelerations leading to velocity truncations + + ! real, pointer :: hf_du_dt_visc(:,:,:) => NULL() ! Zonal friction accel. x fract. thickness [L T-2 ~> m s-2]. + ! real, pointer :: hf_dv_dt_visc(:,:,:) => NULL() ! Merdional friction accel. x fract. thickness [L T-2 ~> m s-2]. + ! 3D diagnostics hf_du(dv)_dt_visc are commented because there is no clarity on proper remapping grid option. + ! The code is retained for degugging purposes in the future. + end type vertvisc_CS contains @@ -202,11 +210,14 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & real :: surface_stress(SZIB_(G))! The same as stress, unless the wind stress ! stress is applied as a body force [H L T-1 ~> m2 s-1 or kg m-1 s-1]. + real, allocatable, dimension(:,:) :: hf_du_dt_visc_2d ! Depth sum of hf_du_dt_visc [L T-2 ~> m s-2] + real, allocatable, dimension(:,:) :: hf_dv_dt_visc_2d ! Depth sum of hf_dv_dt_visc [L T-2 ~> m s-2] + logical :: do_i(SZIB_(G)) logical :: DoStokesMixing - integer :: i, j, k, is, ie, Isq, Ieq, Jsq, Jeq, nz, n - is = G%isc ; ie = G%iec + integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, n + is = G%isc ; ie = G%iec; js = G%jsc; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = G%ke if (.not.associated(CS)) call MOM_error(FATAL,"MOM_vert_friction(visc): "// & @@ -453,6 +464,41 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & if (present(tauy_bot) .and. (CS%id_tauy_bot > 0)) & call post_data(CS%id_tauy_bot, tauy_bot, CS%diag) + ! Diagnostics for terms multiplied by fractional thicknesses + + ! 3D diagnostics hf_du(dv)_dt_visc are commented because there is no clarity on proper remapping grid option. + ! The code is retained for degugging purposes in the future. + !if (CS%id_hf_du_dt_visc > 0) then + ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq + ! CS%hf_du_dt_visc(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%diag_hfrac_u(I,j,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_du_dt_visc, CS%hf_du_dt_visc, CS%diag) + !endif + !if (CS%id_hf_dv_dt_visc > 0) then + ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + ! CS%hf_dv_dt_visc(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%diag_hfrac_v(i,J,k) + ! enddo ; enddo ; enddo + ! call post_data(CS%id_hf_dv_dt_visc, CS%hf_dv_dt_visc, CS%diag) + !endif + if (CS%id_hf_du_dt_visc_2d > 0) then + allocate(hf_du_dt_visc_2d(G%IsdB:G%IedB,G%jsd:G%jed)) + hf_du_dt_visc_2d(:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + hf_du_dt_visc_2d(I,j) = hf_du_dt_visc_2d(I,j) + ADp%du_dt_visc(I,j,k) * ADp%diag_hfrac_u(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_du_dt_visc_2d, hf_du_dt_visc_2d, CS%diag) + deallocate(hf_du_dt_visc_2d) + endif + if (CS%id_hf_dv_dt_visc_2d > 0) then + allocate(hf_dv_dt_visc_2d(G%isd:G%ied,G%JsdB:G%JedB)) + hf_dv_dt_visc_2d(:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + hf_dv_dt_visc_2d(i,J) = hf_dv_dt_visc_2d(i,J) + ADp%dv_dt_visc(i,J,k) * ADp%diag_hfrac_v(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_hf_dv_dt_visc_2d, hf_dv_dt_visc_2d, CS%diag) + deallocate(hf_dv_dt_visc_2d) + endif + end subroutine vertvisc !> Calculate the fraction of momentum originally in a layer that remains @@ -1760,6 +1806,40 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & Time, 'Meridional Bottom Stress from Ocean to Earth', 'Pa', & conversion=US%RZ_to_kg_m2*US%L_T2_to_m_s2) + !CS%id_hf_du_dt_visc = register_diag_field('ocean_model', 'hf_du_dt_visc', diag%axesCuL, Time, & + ! 'Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity', 'm s-2', & + ! v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if (CS%id_hf_du_dt_visc > 0) then + ! call safe_alloc_ptr(CS%hf_du_dt_visc,IsdB,IedB,jsd,jed,nz) + ! call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) + ! call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + !endif + + !CS%id_hf_dv_dt_visc = register_diag_field('ocean_model', 'hf_dv_dt_visc', diag%axesCvL, Time, & + ! 'Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity', 'm s-2', & + ! v_extensive=.true., conversion=US%L_T2_to_m_s2) + !if (CS%id_hf_dv_dt_visc > 0) then + ! call safe_alloc_ptr(CS%hf_dv_dt_visc,isd,ied,JsdB,JedB,nz) + ! call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) + ! call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + !endif + + CS%id_hf_du_dt_visc_2d = register_diag_field('ocean_model', 'hf_du_dt_visc_2d', diag%axesCu1, Time, & + 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if (CS%id_hf_du_dt_visc_2d > 0) then + call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) + endif + + CS%id_hf_dv_dt_visc_2d = register_diag_field('ocean_model', 'hf_dv_dt_visc_2d', diag%axesCv1, Time, & + 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity', 'm s-2', & + conversion=US%L_T2_to_m_s2) + if (CS%id_hf_dv_dt_visc_2d > 0) then + call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz) + endif + if ((len_trim(CS%u_trunc_file) > 0) .or. (len_trim(CS%v_trunc_file) > 0)) & call PointAccel_init(MIS, Time, G, param_file, diag, dirs, CS%PointAccel_CSp) From 16a0a06566866cdf576320b0909841847a37430f Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Wed, 29 Jul 2020 21:25:37 -0400 Subject: [PATCH 2/4] OBC: H-dimensionality fixes This patch fixes three dimensionality errors in the OBC segments. - We add a missing GV%m_to_H conversion for time-dependent eta segments. - The `adjustSegmentEtaToFitBathymetry` function depends on the `segment%Htot` field when computing the dz_src cell spacings. While the calculation primarily assumes that all quantities scale as Z, the segment%Htot field scales as H, which was causing dimensionality errors. We resolve this by converting Htot from Z to H whenever it is used in the calculation. - Segment barotropic velocitys based on transports were limited to a hard-coded thickness of 1e-12, which was not scaled. We have added H-dimensional scaling to these constants. --- src/core/MOM_open_boundary.F90 | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 31f037c66e..35d559ab52 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -3947,7 +3947,8 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) G%dyCu(I,j) normal_trans_bt(I,j) = normal_trans_bt(I,j) + segment%normal_trans(I,j,k) enddo - segment%normal_vel_bt(I,j) = normal_trans_bt(I,j) / (max(segment%Htot(I,j),1.e-12) * G%dyCu(I,j)) + segment%normal_vel_bt(I,j) = normal_trans_bt(I,j) & + / (max(segment%Htot(I,j), 1.e-12 * GV%m_to_H) * G%dyCu(I,j)) if (associated(segment%nudged_normal_vel)) segment%nudged_normal_vel(I,j,:) = segment%normal_vel(I,j,:) enddo elseif (trim(segment%field(m)%name) == 'V' .and. segment%is_N_or_S) then @@ -3960,7 +3961,8 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) G%dxCv(i,J) normal_trans_bt(i,J) = normal_trans_bt(i,J) + segment%normal_trans(i,J,k) enddo - segment%normal_vel_bt(i,J) = normal_trans_bt(i,J) / (max(segment%Htot(i,J),1.e-12) * G%dxCv(i,J)) + segment%normal_vel_bt(i,J) = normal_trans_bt(i,J) & + / (max(segment%Htot(i,J), 1.e-12 * GV%m_to_H) * G%dxCv(i,J)) if (associated(segment%nudged_normal_vel)) segment%nudged_normal_vel(i,J,:) = segment%normal_vel(i,J,:) enddo elseif (trim(segment%field(m)%name) == 'V' .and. segment%is_E_or_W .and. & @@ -4028,13 +4030,14 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time) if (OBC%ramp) then do j=js_obc2,je_obc do i=is_obc2,ie_obc - segment%eta(i,j) = OBC%ramp_value * segment%field(m)%buffer_dst(i,j,1) + segment%eta(i,j) = GV%m_to_H * OBC%ramp_value & + * segment%field(m)%buffer_dst(i,j,1) enddo enddo else do j=js_obc2,je_obc do i=is_obc2,ie_obc - segment%eta(i,j) = segment%field(m)%buffer_dst(i,j,1) + segment%eta(i,j) = GV%m_to_H * segment%field(m)%buffer_dst(i,j,1) enddo enddo endif @@ -4883,8 +4886,8 @@ subroutine adjustSegmentEtaToFitBathymetry(G, GV, US, segment,fld) ! The normal slope at the boundary is zero by a ! previous call to open_boundary_impose_normal_slope do k=nz+1,1,-1 - if (-eta(i,j,k) > segment%Htot(i,j) + hTolerance) then - eta(i,j,k) = -segment%Htot(i,j) + if (-eta(i,j,k) > segment%Htot(i,j)*GV%H_to_Z + hTolerance) then + eta(i,j,k) = -segment%Htot(i,j)*GV%H_to_Z contractions = contractions + 1 endif enddo @@ -4902,10 +4905,10 @@ subroutine adjustSegmentEtaToFitBathymetry(G, GV, US, segment,fld) ! The whole column is dilated to accommodate deeper topography than ! the bathymetry would indicate. - if (-eta(i,j,nz+1) < segment%Htot(i,j) - hTolerance) then + if (-eta(i,j,nz+1) < (segment%Htot(i,j) * GV%H_to_Z) - hTolerance) then dilations = dilations + 1 ! expand bottom-most cell only - eta(i,j,nz+1) = -segment%Htot(i,j) + eta(i,j,nz+1) = -(segment%Htot(i,j) * GV%H_to_Z) segment%field(fld)%dz_src(i,j,nz)= eta(i,j,nz)-eta(i,j,nz+1) ! if (eta(i,j,1) <= eta(i,j,nz+1)) then ! do k=1,nz ; segment%field(fld)%dz_src(i,j,k) = (eta(i,j,1) + G%bathyT(i,j)) / real(nz) ; enddo From 46714ecfc2a5e1a6d798758c300f10d8ae3899ca Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 30 Jul 2020 10:56:15 -0400 Subject: [PATCH 3/4] OBC: Removal of segment zero The current OBC segment list includes a "segment zero" which corresponds to when the segment maps (segnum_u, segnum_v) point to OBC_NONE, which is set to zero. The purpose of this "segment zero" seems to be for avoiding invalid references, so that OBC%segment(OBC%segnum_[uv](i,j)) always returns a valid result, included when set to OBC_NONE. This results in reading and setting values which are unused or unneeded. This patch replaces these redundant accesses to segment zero with various conditional blocks for avoiding these accesses. --- src/core/MOM_barotropic.F90 | 15 ++- src/core/MOM_continuity_PPM.F90 | 109 +++++++++++++----- src/core/MOM_isopycnal_slopes.F90 | 49 ++++---- src/core/MOM_open_boundary.F90 | 9 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 36 ++++-- .../vertical/MOM_set_diffusivity.F90 | 44 ++++--- 6 files changed, 178 insertions(+), 84 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 7f34f18998..ae43868341 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -687,6 +687,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB integer :: ioff, joff + integer :: l_seg if (.not.associated(CS)) call MOM_error(FATAL, & "btstep: Module MOM_barotropic must be initialized before it is used.") @@ -2324,9 +2325,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (CS%BT_OBC%apply_u_OBCs) then ! copy back the value for u-points on the boundary. !GOMP parallel do default(shared) do j=js,je ; do I=is-1,ie - if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then + l_seg = OBC%segnum_u(I,j) + if (l_seg == OBC_NONE) cycle + + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then e_anom(i+1,j) = e_anom(i,j) - elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_W) then e_anom(i,j) = e_anom(i+1,j) endif enddo ; enddo @@ -2335,9 +2339,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, if (CS%BT_OBC%apply_v_OBCs) then ! copy back the value for v-points on the boundary. !GOMP parallel do default(shared) do J=js-1,je ; do I=is,ie - if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then + l_seg = OBC%segnum_v(i,J) + if (l_seg == OBC_NONE) cycle + + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then e_anom(i,j+1) = e_anom(i,j) - elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_S) then e_anom(i,j) = e_anom(i,j+1) endif enddo ; enddo diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index c594d31494..995827959d 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -263,6 +263,7 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & real :: du_lim ! The velocity change that give a relative CFL of 1 [L T-1 ~> m s-1]. real :: dx_E, dx_W ! Effective x-grid spacings to the east and west [L ~> m]. integer :: i, j, k, ish, ieh, jsh, jeh, n, nz + integer :: l_seg logical :: local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC logical :: local_Flather_OBC, local_open_BC, is_simple type(OBC_segment_type), pointer :: segment => NULL() @@ -303,7 +304,8 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & !$OMP uh,dt,US,G,GV,CS,local_specified_BC,OBC,uhbt,set_BT_cont, & !$OMP CFL_dt,I_dt,u_cor,BT_cont, local_Flather_OBC) & !$OMP private(do_I,duhdu,du,du_max_CFL,du_min_CFL,uh_tot_0,duhdu_tot_0, & -!$OMP is_simple,FAuI,visc_rem_max,I_vrm,du_lim,dx_E,dx_W,any_simple_OBC ) & +!$OMP is_simple,FAuI,visc_rem_max,I_vrm,du_lim,dx_E,dx_W, & +!$OMP any_simple_OBC,l_seg) & !$OMP firstprivate(visc_rem) do j=jsh,jeh do I=ish-1,ieh ; do_I(I) = .true. ; visc_rem_max(I) = 0.0 ; enddo @@ -318,8 +320,12 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & dt, G, US, j, ish, ieh, do_I, CS%vol_CFL, OBC) if (local_specified_BC) then do I=ish-1,ieh - if (OBC%segment(OBC%segnum_u(I,j))%specified) & - uh(I,j,k) = OBC%segment(OBC%segnum_u(I,j))%normal_trans(I,j,k) + l_seg = OBC%segnum_u(I,j) + + if (l_seg /= OBC_NONE) then + if (OBC%segment(l_seg)%specified) & + uh(I,j,k) = OBC%segment(l_seg)%normal_trans(I,j,k) + endif enddo endif enddo @@ -408,9 +414,13 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & any_simple_OBC = .false. if (present(uhbt) .or. set_BT_cont) then if (local_specified_BC .or. local_Flather_OBC) then ; do I=ish-1,ieh + l_seg = OBC%segnum_u(I,j) + ! Avoid reconciling barotropic/baroclinic transports if transport is specified - is_simple = OBC%segment(OBC%segnum_u(I,j))%specified - do_I(I) = .not.(OBC%segnum_u(I,j) /= OBC_NONE .and. is_simple) + is_simple = .false. + if (l_seg /= OBC_NONE) & + is_simple = OBC%segment(l_seg)%specified + do_I(I) = .not. (l_seg /= OBC_NONE .and. is_simple) any_simple_OBC = any_simple_OBC .or. is_simple enddo ; else ; do I=ish-1,ieh do_I(I) = .true. @@ -425,8 +435,12 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & if (present(u_cor)) then ; do k=1,nz do I=ish-1,ieh ; u_cor(I,j,k) = u(I,j,k) + du(I) * visc_rem(I,k) ; enddo if (local_specified_BC) then ; do I=ish-1,ieh - if (OBC%segment(OBC%segnum_u(I,j))%specified) & - u_cor(I,j,k) = OBC%segment(OBC%segnum_u(I,j))%normal_vel(I,j,k) + l_seg = OBC%segnum_u(I,j) + + if (l_seg /= OBC_NONE) then + if (OBC%segment(l_seg)%specified) & + u_cor(I,j,k) = OBC%segment(l_seg)%normal_vel(I,j,k) + endif enddo ; endif enddo ; endif ! u-corrected @@ -438,9 +452,15 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, US, CS, LB, uhbt, OBC, & visc_rem_max, j, ish, ieh, do_I) if (any_simple_OBC) then do I=ish-1,ieh - do_I(I) = OBC%segment(OBC%segnum_u(I,j))%specified + l_seg = OBC%segnum_u(I,j) + + do_I(I) = .false. + if (l_seg /= OBC_NONE) & + do_I(I) = OBC%segment(l_seg)%specified + if (do_I(I)) FAuI(I) = GV%H_subroundoff*G%dy_Cu(I,j) enddo + ! NOTE: do_I(I) should prevent access to segment OBC_NONE do k=1,nz ; do I=ish-1,ieh ; if (do_I(I)) then if ((abs(OBC%segment(OBC%segnum_u(I,j))%normal_vel(I,j,k)) > 0.0) .and. & (OBC%segment(OBC%segnum_u(I,j))%specified)) & @@ -529,6 +549,7 @@ subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, US, j, & ! with the same units as h_in. real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2]. integer :: i + integer :: l_seg logical :: local_open_BC local_open_BC = .false. @@ -561,13 +582,17 @@ subroutine zonal_flux_layer(u, h, h_L, h_R, uh, duhdu, visc_rem, dt, G, US, j, & if (local_open_BC) then do I=ish-1,ieh ; if (do_I(I)) then - if (OBC%segment(OBC%segnum_u(I,j))%open) then - if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then - uh(I) = G%dy_Cu(I,j) * u(I) * h(i) - duhdu(I) = G%dy_Cu(I,j) * h(i) * visc_rem(I) - else - uh(I) = G%dy_Cu(I,j) * u(I) * h(i+1) - duhdu(I) = G%dy_Cu(I,j) * h(i+1) * visc_rem(I) + l_seg = OBC%segnum_u(I,j) + + if (l_seg /= OBC_NONE) then + if (OBC%segment(l_seg)%open) then + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then + uh(I) = G%dy_Cu(I,j) * u(I) * h(i) + duhdu(I) = G%dy_Cu(I,j) * h(i) * visc_rem(I) + else + uh(I) = G%dy_Cu(I,j) * u(I) * h(i+1) + duhdu(I) = G%dy_Cu(I,j) * h(i+1) * visc_rem(I) + endif endif endif endif ; enddo @@ -1062,6 +1087,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, & real :: dv_lim ! The velocity change that give a relative CFL of 1 [L T-1 ~> m s-1]. real :: dy_N, dy_S ! Effective y-grid spacings to the north and south [L ~> m]. integer :: i, j, k, ish, ieh, jsh, jeh, n, nz + integer :: l_seg logical :: local_specified_BC, use_visc_rem, set_BT_cont, any_simple_OBC logical :: local_Flather_OBC, is_simple, local_open_BC type(OBC_segment_type), pointer :: segment => NULL() @@ -1103,7 +1129,7 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, & !$OMP set_BT_cont,CFL_dt,I_dt,v_cor,BT_cont, local_Flather_OBC ) & !$OMP private(do_I,dvhdv,dv,dv_max_CFL,dv_min_CFL,vh_tot_0, & !$OMP dvhdv_tot_0,visc_rem_max,I_vrm,dv_lim,dy_N, & -!$OMP is_simple,FAvi,dy_S,any_simple_OBC ) & +!$OMP is_simple,FAvi,dy_S,any_simple_OBC,l_seg) & !$OMP firstprivate(visc_rem) do J=jsh-1,jeh do i=ish,ieh ; do_I(i) = .true. ; visc_rem_max(I) = 0.0 ; enddo @@ -1118,8 +1144,12 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, & dt, G, US, J, ish, ieh, do_I, CS%vol_CFL, OBC) if (local_specified_BC) then do i=ish,ieh - if (OBC%segment(OBC%segnum_v(i,J))%specified) & - vh(i,J,k) = OBC%segment(OBC%segnum_v(i,J))%normal_trans(i,J,k) + l_seg = OBC%segnum_v(i,J) + + if (l_seg /= OBC_NONE) then + if (OBC%segment(l_seg)%specified) & + vh(i,J,k) = OBC%segment(l_seg)%normal_trans(i,J,k) + endif enddo endif enddo ! k-loop @@ -1204,9 +1234,13 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, & any_simple_OBC = .false. if (present(vhbt) .or. set_BT_cont) then if (local_specified_BC .or. local_Flather_OBC) then ; do i=ish,ieh + l_seg = OBC%segnum_v(i,J) + ! Avoid reconciling barotropic/baroclinic transports if transport is specified - is_simple = OBC%segment(OBC%segnum_v(i,J))%specified - do_I(i) = .not.(OBC%segnum_v(i,J) /= OBC_NONE .and. is_simple) + is_simple = .false. + if (l_seg /= OBC_NONE) & + is_simple = OBC%segment(l_seg)%specified + do_I(i) = .not.(l_seg /= OBC_NONE .and. is_simple) any_simple_OBC = any_simple_OBC .or. is_simple enddo ; else ; do i=ish,ieh do_I(i) = .true. @@ -1221,8 +1255,12 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, & if (present(v_cor)) then ; do k=1,nz do i=ish,ieh ; v_cor(i,J,k) = v(i,J,k) + dv(i) * visc_rem(i,k) ; enddo if (local_specified_BC) then ; do i=ish,ieh - if (OBC%segment(OBC%segnum_v(i,J))%specified) & - v_cor(i,J,k) = OBC%segment(OBC%segnum_v(i,J))%normal_vel(i,J,k) + l_seg = OBC%segnum_v(i,J) + + if (l_seg /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J))%specified) & + v_cor(i,J,k) = OBC%segment(OBC%segnum_v(i,J))%normal_vel(i,J,k) + endif enddo ; endif enddo ; endif ! v-corrected endif @@ -1233,9 +1271,15 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, US, CS, LB, vhbt, OBC, & visc_rem_max, J, ish, ieh, do_I) if (any_simple_OBC) then do i=ish,ieh - do_I(i) = (OBC%segment(OBC%segnum_v(i,J))%specified) + l_seg = OBC%segnum_v(i,J) + + do_I(I) = .false. + if(l_seg /= OBC_NONE) & + do_I(i) = (OBC%segment(l_seg)%specified) + if (do_I(i)) FAvi(i) = GV%H_subroundoff*G%dx_Cv(i,J) enddo + ! NOTE: do_I(I) should prevent access to segment OBC_NONE do k=1,nz ; do i=ish,ieh ; if (do_I(i)) then if ((abs(OBC%segment(OBC%segnum_v(i,J))%normal_vel(i,J,k)) > 0.0) .and. & (OBC%segment(OBC%segnum_v(i,J))%specified)) & @@ -1327,6 +1371,7 @@ subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, & ! with the same units as h, i.e. [H ~> m or kg m-2]. real :: h_marg ! The marginal thickness of a flux [H ~> m or kg m-2]. integer :: i + integer :: l_seg logical :: local_open_BC local_open_BC = .false. @@ -1360,13 +1405,17 @@ subroutine merid_flux_layer(v, h, h_L, h_R, vh, dvhdv, visc_rem, dt, G, US, J, & if (local_open_BC) then do i=ish,ieh ; if (do_I(i)) then - if (OBC%segment(OBC%segnum_v(i,J))%open) then - if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then - vh(i) = G%dx_Cv(i,J) * v(i) * h(i,j) - dvhdv(i) = G%dx_Cv(i,J) * h(i,j) * visc_rem(i) - else - vh(i) = G%dx_Cv(i,J) * v(i) * h(i,j+1) - dvhdv(i) = G%dx_Cv(i,J) * h(i,j+1) * visc_rem(i) + l_seg = OBC%segnum_v(i,J) + + if (l_seg /= OBC_NONE) then + if (OBC%segment(l_seg)%open) then + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then + vh(i) = G%dx_Cv(i,J) * v(i) * h(i,j) + dvhdv(i) = G%dx_Cv(i,J) * h(i,j) * visc_rem(i) + else + vh(i) = G%dx_Cv(i,J) * v(i) * h(i,j+1) + dvhdv(i) = G%dx_Cv(i,J) * h(i,j+1) * visc_rem(i) + endif endif endif endif ; enddo diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 7a33dc7d77..c134366cd0 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -8,7 +8,7 @@ module MOM_isopycnal_slopes use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_EOS, only : calculate_density_derivs -use MOM_open_boundary, only : ocean_OBC_type +use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S implicit none ; private @@ -105,6 +105,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & integer, dimension(2) :: EOSdom_u, EOSdom_v ! Domains for the equation of state calculations at u and v points integer :: is, ie, js, je, nz, IsdB integer :: i, j, k + integer :: l_seg logical :: local_open_u_BC, local_open_v_BC if (present(halo)) then @@ -183,7 +184,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, & !$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, & - !$OMP drdx,mag_grad2,Slope,slope2_Ratio) + !$OMP drdx,mag_grad2,Slope,slope2_Ratio,l_seg) do j=js,je ; do K=nz,2,-1 if (.not.(use_EOS)) then drdiA = 0.0 ; drdiB = 0.0 @@ -260,15 +261,18 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & slope_x(I,j,K) = (Z_to_L*(e(i,j,K)-e(i+1,j,K))) * G%IdxCu(I,j) endif if (local_open_u_BC) then - if (OBC%segment(OBC%segnum_u(I,j))%open) then - slope_x(I,j,K) = 0. - ! This and/or the masking code below is to make slopes match inside - ! land mask. Might not be necessary except for DEBUG output. -! if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then -! slope_x(I+1,j,K) = 0. -! else -! slope_x(I-1,j,K) = 0. -! endif + l_seg = OBC%segnum_u(I,j) + if (l_seg /= OBC_NONE) then + if (OBC%segment(l_seg)%open) then + slope_x(I,j,K) = 0. + ! This and/or the masking code below is to make slopes match inside + ! land mask. Might not be necessary except for DEBUG output. +! if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then +! slope_x(I+1,j,K) = 0. +! else +! slope_x(I-1,j,K) = 0. +! endif + endif endif slope_x(I,j,K) = slope_x(I,j,k) * max(g%mask2dT(i,j),g%mask2dT(i+1,j)) endif @@ -286,7 +290,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, & !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, & !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, & - !$OMP drdy,mag_grad2,Slope,slope2_Ratio) + !$OMP drdy,mag_grad2,Slope,slope2_Ratio,l_seg) do j=js-1,je ; do K=nz,2,-1 if (.not.(use_EOS)) then drdjA = 0.0 ; drdjB = 0.0 @@ -360,15 +364,18 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, & slope_y(i,J,K) = (Z_to_L*(e(i,j,K)-e(i,j+1,K))) * G%IdyCv(i,J) endif if (local_open_v_BC) then - if (OBC%segment(OBC%segnum_v(i,J))%open) then - slope_y(i,J,K) = 0. - ! This and/or the masking code below is to make slopes match inside - ! land mask. Might not be necessary except for DEBUG output. -! if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then -! slope_y(i,J+1,K) = 0. -! else -! slope_y(i,J-1,K) = 0. -! endif + l_seg = OBC%segnum_v(i,J) + if (l_seg /= OBC_NONE) then + if (OBC%segment(l_seg)%open) then + slope_y(i,J,K) = 0. + ! This and/or the masking code below is to make slopes match inside + ! land mask. Might not be necessary except for DEBUG output. +! if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then +! slope_y(i,J+1,K) = 0. +! else +! slope_y(i,J-1,K) = 0. +! endif + endif endif slope_y(i,J,K) = slope_y(i,J,k) * max(g%mask2dT(i,j),g%mask2dT(i,j+1)) endif diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 31f037c66e..b09cdc50c6 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -445,9 +445,8 @@ subroutine open_boundary_config(G, US, param_file, OBC) call get_param(param_file, mdl, "REENTRANT_Y", reentrant_y, default=.false.) ! Allocate everything - ! Note the 0-segment is needed when %segnum_u/v(:,:) = 0 - allocate(OBC%segment(0:OBC%number_of_segments)) - do l=0,OBC%number_of_segments + allocate(OBC%segment(1:OBC%number_of_segments)) + do l=1,OBC%number_of_segments OBC%segment(l)%Flather = .false. OBC%segment(l)%radiation = .false. OBC%segment(l)%radiation_tan = .false. @@ -4971,7 +4970,7 @@ subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) ! Segment rotation allocate(OBC%segment(0:OBC%number_of_segments)) - do l = 0, OBC%number_of_segments + do l = 1, OBC%number_of_segments call rotate_OBC_segment_config(OBC_in%segment(l), G_in, OBC%segment(l), G, turns) ! Data up to setup_[uv]_point_obc is needed for allocate_obc_segment_data! call allocate_OBC_segment_data(OBC, OBC%segment(l)) @@ -5168,7 +5167,7 @@ subroutine rotate_OBC_init(OBC_in, G, GV, US, param_file, tv, restart_CSp, OBC) "If true, Temperature and salinity are used as state "//& "variables.", default=.true., do_not_log=.true.) - do l = 0, OBC%number_of_segments + do l = 1, OBC%number_of_segments call rotate_OBC_segment_data(OBC_in%segment(l), OBC%segment(l), G%HI%turns) enddo diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index e0def91821..9a0a994450 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -17,7 +17,7 @@ module MOM_lateral_mixing_coeffs use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_speed, only : wave_speed, wave_speed_CS, wave_speed_init -use MOM_open_boundary, only : ocean_OBC_type +use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE implicit none ; private @@ -499,6 +499,7 @@ subroutine calc_Visbeck_coeffs(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, CS, O real :: H_geom ! The geometric mean of Hup*Hdn [H ~> m or kg m-2]. integer :: is, ie, js, je, nz integer :: i, j, k, kb_max + integer :: l_seg real :: S2max, wNE, wSE, wSW, wNW real :: H_u(SZIB_(G)), H_v(SZI_(G)) real :: S2_u(SZIB_(G), SZJ_(G)) @@ -568,8 +569,12 @@ subroutine calc_Visbeck_coeffs(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, CS, O CS%SN_u(I,j) = 0. endif if (local_open_u_BC) then - if (OBC%segment(OBC%segnum_u(I,j))%open) then - CS%SN_u(i,J) = 0. + l_seg = OBC%segnum_u(I,j) + + if (l_seg /= OBC_NONE) then + if (OBC%segment(l_seg)%open) then + CS%SN_u(i,J) = 0. + endif endif endif enddo @@ -609,8 +614,12 @@ subroutine calc_Visbeck_coeffs(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, CS, O CS%SN_v(i,J) = 0. endif if (local_open_v_BC) then - if (OBC%segment(OBC%segnum_v(i,J))%open) then - CS%SN_v(i,J) = 0. + l_seg = OBC%segnum_v(i,J) + + if (l_seg /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(i,J))%open) then + CS%SN_v(i,J) = 0. + endif endif endif enddo @@ -657,6 +666,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop real :: one_meter ! One meter in thickness units [H ~> m or kg m-2]. integer :: is, ie, js, je, nz integer :: i, j, k, kb_max + integer :: l_seg real :: S2N2_u_local(SZIB_(G), SZJ_(G),SZK_(G)) real :: S2N2_v_local(SZI_(G), SZJB_(G),SZK_(G)) logical :: local_open_u_BC, local_open_v_BC @@ -754,8 +764,12 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop CS%SN_u(I,j) = 0.0 endif if (local_open_u_BC) then - if (OBC%segment(OBC%segnum_u(I,j))%open) then - CS%SN_u(I,j) = 0. + l_seg = OBC%segnum_u(I,j) + + if (l_seg /= OBC_NONE) then + if (OBC%segment(l_seg)%open) then + CS%SN_u(I,j) = 0. + endif endif endif enddo @@ -776,8 +790,12 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e, calculate_slop CS%SN_v(I,j) = 0.0 endif if (local_open_v_BC) then - if (OBC%segment(OBC%segnum_v(I,j))%open) then - CS%SN_v(I,j) = 0. + l_seg = OBC%segnum_v(I,j) + + if (l_seg /= OBC_NONE) then + if (OBC%segment(OBC%segnum_v(I,j))%open) then + CS%SN_v(I,j) = 0. + endif endif endif enddo diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 86f828e5fa..42babae7d8 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -1677,7 +1677,9 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC) logical :: domore, do_i(SZI_(G)) integer :: i, j, k, is, ie, js, je, nz + integer :: l_seg logical :: local_open_u_BC, local_open_v_BC + logical :: has_obc local_open_u_BC = .false. local_open_v_BC = .false. @@ -1718,15 +1720,21 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC) do k=nz,1,-1 domore = .false. do i=is,ie ; if (do_i(i)) then + ! Determine if grid point is an OBC + has_obc = .false. if (local_open_v_BC) then - if (OBC%segment(OBC%segnum_v(i,J))%open) then - if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then - hvel = GV%H_to_Z*h(i,j,k) - else - hvel = GV%H_to_Z*h(i,j+1,k) - endif + l_seg = OBC%segnum_v(i,J) + if (l_seg /= OBC_NONE) then + has_obc = OBC%segment(l_seg)%open + endif + endif + + ! Compute h based on OBC state + if (has_obc) then + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then + hvel = GV%H_to_Z*h(i,j,k) else - hvel = 0.5*GV%H_to_Z*(h(i,j,k) + h(i,j+1,k)) + hvel = GV%H_to_Z*h(i,j+1,k) endif else hvel = 0.5*GV%H_to_Z*(h(i,j,k) + h(i,j+1,k)) @@ -1760,15 +1768,21 @@ subroutine set_BBL_TKE(u, v, h, fluxes, visc, G, GV, US, CS, OBC) endif ; enddo do k=nz,1,-1 ; domore = .false. do I=is-1,ie ; if (do_i(I)) then + ! Determine if grid point is an OBC + has_obc = .false. if (local_open_u_BC) then - if (OBC%segment(OBC%segnum_u(I,j))%open) then - if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then - hvel = GV%H_to_Z*h(i,j,k) - else - hvel = GV%H_to_Z*h(i+1,j,k) - endif - else - hvel = 0.5*GV%H_to_Z*(h(i,j,k) + h(i+1,j,k)) + l_seg = OBC%segnum_u(I,j) + if (l_seg /= OBC_NONE) then + has_obc = OBC%segment(l_seg)%open + endif + endif + + ! Compute h based on OBC state + if (has_obc) then + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then + hvel = GV%H_to_Z*h(i,j,k) + else ! OBC_DIRECTION_W + hvel = GV%H_to_Z*h(i+1,j,k) endif else hvel = 0.5*GV%H_to_Z*(h(i,j,k) + h(i+1,j,k)) From a243225b3e8a18b74b3ebbdc58ac80f231be2408 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 30 Jul 2020 15:14:26 -0400 Subject: [PATCH 4/4] OBC: Remove segment 0 refs in mask_outside_OBCs There were a few remaining segnum_[uv] references in MOM_open_boundary which could reference segment zero. This patch fixes those references. --- src/core/MOM_open_boundary.F90 | 29 +++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index b09cdc50c6..74cad3cf0c 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -4416,6 +4416,7 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) ! Local variables integer :: isd, ied, IsdB, IedB, jsd, jed, JsdB, JedB, n integer :: i, j + integer :: l_seg logical :: fatal_error = .False. real :: min_depth integer, parameter :: cin = 3, cout = 4, cland = -1, cedge = -2 @@ -4457,38 +4458,50 @@ subroutine mask_outside_OBCs(G, US, param_file, OBC) enddo do j=G%jsd,G%jed ; do i=G%IsdB+1,G%IedB-1 - if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + l_seg = OBC%segnum_u(I,j) + if (l_seg == OBC_NONE) cycle + + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_W) then if (color(i,j) == 0.0) color(i,j) = cout if (color(i+1,j) == 0.0) color(i+1,j) = cin - elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then + elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then if (color(i,j) == 0.0) color(i,j) = cin if (color(i+1,j) == 0.0) color(i+1,j) = cout endif enddo ; enddo do J=G%JsdB+1,G%JedB-1 ; do i=G%isd,G%ied - if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + l_seg = OBC%segnum_v(i,J) + if (l_seg == OBC_NONE) cycle + + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_S) then if (color(i,j) == 0.0) color(i,j) = cout if (color(i,j+1) == 0.0) color(i,j+1) = cin - elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then + elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then if (color(i,j) == 0.0) color(i,j) = cin if (color(i,j+1) == 0.0) color(i,j+1) = cout endif enddo ; enddo do J=G%JsdB+1,G%JedB-1 ; do i=G%isd,G%ied - if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + l_seg = OBC%segnum_v(i,J) + if (l_seg == OBC_NONE) cycle + + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_S) then if (color2(i,j) == 0.0) color2(i,j) = cout if (color2(i,j+1) == 0.0) color2(i,j+1) = cin - elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then + elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then if (color2(i,j) == 0.0) color2(i,j) = cin if (color2(i,j+1) == 0.0) color2(i,j+1) = cout endif enddo ; enddo do j=G%jsd,G%jed ; do i=G%IsdB+1,G%IedB-1 - if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + l_seg = OBC%segnum_u(I,j) + if (l_seg == OBC_NONE) cycle + + if (OBC%segment(l_seg)%direction == OBC_DIRECTION_W) then if (color2(i,j) == 0.0) color2(i,j) = cout if (color2(i+1,j) == 0.0) color2(i+1,j) = cin - elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then + elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then if (color2(i,j) == 0.0) color2(i,j) = cin if (color2(i+1,j) == 0.0) color2(i+1,j) = cout endif