diff --git a/src/core/MOM_CoriolisAdv.F90 b/src/core/MOM_CoriolisAdv.F90 index b12d3e37e7..baf0bfdf86 100644 --- a/src/core/MOM_CoriolisAdv.F90 +++ b/src/core/MOM_CoriolisAdv.F90 @@ -79,6 +79,8 @@ module MOM_CoriolisAdv integer :: id_intz_gKEu_2d = -1, id_intz_gKEv_2d = -1 ! integer :: id_hf_rvxu = -1, id_hf_rvxv = -1 integer :: id_hf_rvxu_2d = -1, id_hf_rvxv_2d = -1 + integer :: id_h_gKEu = -1, id_h_gKEv = -1 + integer :: id_h_rvxu = -1, id_h_rvxv = -1 integer :: id_intz_rvxu_2d = -1, id_intz_rvxv_2d = -1 !>@} end type CoriolisAdv_CS @@ -230,6 +232,10 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) ! 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. +! Diagnostics for thickness multiplied momentum budget terms + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: h_gKEu, h_rvxv ! [L2 T-2 ~> m2 s-2]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: h_gKEv, h_rvxu ! [L2 T-2 ~> m2 s-2]. + ! Diagnostics for depth-integrated momentum budget terms real, dimension(SZIB_(G),SZJ_(G)) :: intz_gKEu_2d, intz_rvxv_2d ! [L2 T-2 ~> m2 s-2]. real, dimension(SZI_(G),SZJB_(G)) :: intz_gKEv_2d, intz_rvxu_2d ! [L2 T-2 ~> m2 s-2]. @@ -955,6 +961,36 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS) deallocate(hf_rvxu_2d) endif + if (CS%id_h_gKEu > 0) then + h_gKEu(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + h_gKEu(I,j,k) = AD%gradKEu(I,j,k) * AD%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_gKEu, h_gKEu, CS%diag) + endif + if (CS%id_h_gKEv > 0) then + h_gKEv(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + h_gKEv(i,J,k) = AD%gradKEv(i,J,k) * AD%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_gKEv, h_gKEv, CS%diag) + endif + + if (CS%id_h_rvxv > 0) then + h_rvxv(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + h_rvxv(I,j,k) = AD%rv_x_v(I,j,k) * AD%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_rvxv, h_rvxv, CS%diag) + endif + if (CS%id_h_rvxu > 0) then + h_rvxu(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + h_rvxu(i,J,k) = AD%rv_x_u(i,J,k) * AD%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_rvxu, h_rvxu, CS%diag) + endif + if (CS%id_intz_rvxv_2d > 0) then intz_rvxv_2d(:,:) = 0.0 do k=1,nz ; do j=js,je ; do I=Isq,Ieq @@ -1262,6 +1298,14 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS) call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,JsdB,JedB,nz) endif + CS%id_h_gKEu = register_diag_field('ocean_model', 'h_gKEu', diag%axesCuL, Time, & + 'Thickness Multiplied Zonal Acceleration from Grad. Kinetic Energy', & + 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_h_gKEu > 0) then + call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz) + endif + CS%id_intz_gKEu_2d = register_diag_field('ocean_model', 'intz_gKEu_2d', diag%axesCu1, Time, & 'Depth-integral of Zonal Acceleration from Grad. Kinetic Energy', & 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) @@ -1270,6 +1314,14 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS) call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz) endif + CS%id_h_gKEv = register_diag_field('ocean_model', 'h_gKEv', diag%axesCvL, Time, & + 'Thickness Multiplied Meridional Acceleration from Grad. Kinetic Energy', & + 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_h_gKEv > 0) then + call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(AD%diag_hv,isd,ied,JsdB,JedB,nz) + endif + CS%id_intz_gKEv_2d = register_diag_field('ocean_model', 'intz_gKEv_2d', diag%axesCv1, Time, & 'Depth-integral of Meridional Acceleration from Grad. Kinetic Energy', & 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) @@ -1310,6 +1362,14 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS) call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz) endif + CS%id_h_rvxu = register_diag_field('ocean_model', 'h_rvxu', diag%axesCvL, Time, & + 'Thickness Multiplied Meridional Acceleration from Relative Vorticity', & + 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_h_rvxu > 0) then + call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(AD%diag_hv,isd,ied,JsdB,JedB,nz) + endif + CS%id_intz_rvxu_2d = register_diag_field('ocean_model', 'intz_rvxu_2d', diag%axesCv1, Time, & 'Depth-integral of Meridional Acceleration from Relative Vorticity', & 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) @@ -1318,6 +1378,14 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS) call safe_alloc_ptr(AD%diag_hv,isd,ied,JsdB,JedB,nz) endif + CS%id_h_rvxv = register_diag_field('ocean_model', 'h_rvxv', diag%axesCuL, Time, & + 'Thickness Multiplied Zonal Acceleration from Relative Vorticity', & + 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_h_rvxv > 0) then + call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(AD%diag_hu,IsdB,IedB,jsd,jed,nz) + endif + CS%id_intz_rvxv_2d = register_diag_field('ocean_model', 'intz_rvxv_2d', diag%axesCu1, Time, & 'Depth-integral of Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', & 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index ef7da5c291..63f7103349 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -160,9 +160,11 @@ module MOM_dynamics_split_RK2 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_h_PFu = -1, id_h_PFv = -1 integer :: id_hf_PFu_2d = -1, id_hf_PFv_2d = -1 integer :: id_intz_PFu_2d = -1, id_intz_PFv_2d = -1 ! integer :: id_hf_CAu = -1, id_hf_CAv = -1 + integer :: id_h_CAu = -1, id_h_CAv = -1 integer :: id_hf_CAu_2d = -1, id_hf_CAv_2d = -1 integer :: id_intz_CAu_2d = -1, id_intz_CAv_2d = -1 @@ -170,6 +172,7 @@ module MOM_dynamics_split_RK2 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_h_u_BT_accel = -1, id_h_v_BT_accel = -1 integer :: id_hf_u_BT_accel_2d = -1, id_hf_v_BT_accel_2d = -1 integer :: id_intz_u_BT_accel_2d = -1, id_intz_v_BT_accel_2d = -1 !>@} @@ -339,12 +342,17 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s 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, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: & + h_PFu, h_CAu, h_u_BT_accel ! [L2 T-2 ~> m2 s-2]. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: & + h_PFv, h_CAv, h_v_BT_accel ! [L2 T-2 ~> m2 s-2]. + real, dimension(SZIB_(G),SZJ_(G)) :: & intz_PFu_2d, intz_CAu_2d, intz_u_BT_accel_2d ! [L2 T-2 ~> m2 s-2]. - real, dimension(SZI_(G),SZJB_(G)) :: & intz_PFv_2d, intz_CAv_2d, intz_v_BT_accel_2d ! [L2 T-2 ~> m2 s-2]. + real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s]. logical :: dyn_p_surf @@ -940,6 +948,21 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s deallocate(hf_PFv_2d) endif + if (CS%id_h_PFu > 0) then + h_PFu(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + h_PFu(I,j,k) = CS%PFu(I,j,k) * CS%ADp%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_PFu, h_PFu, CS%diag) + endif + if (CS%id_h_PFv > 0) then + h_PFv(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + h_PFv(i,J,k) = CS%PFv(i,J,k) * CS%ADp%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_PFv, h_PFv, CS%diag) + endif + !if (CS%id_hf_CAu > 0) then ! allocate(hf_CAu(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq @@ -988,6 +1011,21 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s deallocate(hf_CAv_2d) endif + if (CS%id_h_CAu > 0) then + h_CAu(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + h_CAu(I,j,k) = CS%CAu(I,j,k) * CS%ADp%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_CAu, h_CAu, CS%diag) + endif + if (CS%id_h_CAv > 0) then + h_CAv(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + h_CAv(i,J,k) = CS%CAv(i,J,k) * CS%ADp%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_CAv, h_CAv, CS%diag) + endif + !if (CS%id_hf_u_BT_accel > 0) then ! allocate(hf_u_BT_accel(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke)) ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq @@ -1036,6 +1074,21 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s deallocate(hf_v_BT_accel_2d) endif + if (CS%id_h_u_BT_accel > 0) then + h_u_BT_accel(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + h_u_BT_accel(I,j,k) = CS%u_accel_bt(I,j,k) * CS%ADp%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_u_BT_accel, h_u_BT_accel, CS%diag) + endif + if (CS%id_h_v_BT_accel > 0) then + h_v_BT_accel(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + h_v_BT_accel(i,J,k) = CS%v_accel_bt(i,J,k) * CS%ADp%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_v_BT_accel, h_v_BT_accel, CS%diag) + 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) @@ -1442,6 +1495,16 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param 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,JsdB,JedB,nz) + CS%id_h_PFu = register_diag_field('ocean_model', 'h_PFu', diag%axesCuL, Time, & + 'Thickness Multiplied Zonal Pressure Force Acceleration', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_h_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz) + + CS%id_h_PFv = register_diag_field('ocean_model', 'h_PFv', diag%axesCvL, Time, & + 'Thickness Multiplied Meridional Pressure Force Acceleration', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_h_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz) + CS%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, Time, & 'Depth-integral of Zonal Pressure Force Acceleration', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) @@ -1462,6 +1525,16 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param 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,JsdB,JedB,nz) + CS%id_h_CAu = register_diag_field('ocean_model', 'h_CAu', diag%axesCuL, Time, & + 'Thickness Multiplied Zonal Coriolis and Advective Acceleration', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_h_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz) + + CS%id_h_CAv = register_diag_field('ocean_model', 'h_CAv', diag%axesCvL, Time, & + 'Thickness Multiplied Meridional Coriolis and Advective Acceleration', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_h_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz) + CS%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, Time, & 'Depth-integral of Zonal Coriolis and Advective Acceleration', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) @@ -1502,6 +1575,16 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param 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,JsdB,JedB,nz) + CS%id_h_u_BT_accel = register_diag_field('ocean_model', 'h_u_BT_accel', diag%axesCuL, Time, & + 'Thickness Multiplied Barotropic Anomaly Zonal Acceleration', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_h_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz) + + CS%id_h_v_BT_accel = register_diag_field('ocean_model', 'h_v_BT_accel', diag%axesCvL, Time, & + 'Thickness Multiplied Barotropic Anomaly Meridional Acceleration', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if(CS%id_h_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz) + CS%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, Time, & 'Depth-integral of Barotropic Anomaly Zonal Acceleration', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 5ac8777a19..c963852195 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -116,6 +116,7 @@ module MOM_diagnostics 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_h_du_dt = -1, id_h_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 @@ -248,6 +249,9 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & 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]. + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: h_du_dt ! h x dudt [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: h_dv_dt ! h x dvdt [L2 T-2 ~> m2 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] @@ -325,6 +329,21 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, & deallocate(hf_dv_dt_2d) endif + if (CS%id_h_du_dt > 0) then + h_du_dt(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + h_du_dt(I,j,k) = CS%du_dt(I,j,k) * ADp%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_du_dt, h_du_dt, CS%diag) + endif + if (CS%id_h_dv_dt > 0) then + h_dv_dt(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + h_dv_dt(i,J,k) = CS%dv_dt(i,J,k) * ADp%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_dv_dt, h_dv_dt, CS%diag) + endif + call diag_restore_grids(CS%diag) call calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS) @@ -1789,6 +1808,26 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) endif + CS%id_h_du_dt = register_diag_field('ocean_model', 'h_du_dt', diag%axesCuL, Time, & + 'Thickness Multiplied Zonal Acceleration', 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_h_du_dt > 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_hu,IsdB,IedB,jsd,jed,nz) + endif + + CS%id_h_dv_dt = register_diag_field('ocean_model', 'h_dv_dt', diag%axesCvL, Time, & + 'Thickness Multiplied Meridional Acceleration', 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_h_dv_dt > 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_hv,isd,ied,JsdB,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, & diff --git a/src/parameterizations/lateral/MOM_hor_visc.F90 b/src/parameterizations/lateral/MOM_hor_visc.F90 index 661fb715e7..75b26de700 100644 --- a/src/parameterizations/lateral/MOM_hor_visc.F90 +++ b/src/parameterizations/lateral/MOM_hor_visc.F90 @@ -187,6 +187,7 @@ module MOM_hor_visc integer :: id_grid_Re_Ah = -1, id_grid_Re_Kh = -1 integer :: id_diffu = -1, id_diffv = -1 ! integer :: id_hf_diffu = -1, id_hf_diffv = -1 + integer :: id_h_diffu = -1, id_h_diffv = -1 integer :: id_hf_diffu_2d = -1, id_hf_diffv_2d = -1 integer :: id_intz_diffu_2d = -1, id_intz_diffv_2d = -1 integer :: id_Ah_h = -1, id_Ah_q = -1 @@ -277,6 +278,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, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: h_diffu ! h x diffu [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: h_diffv ! h x diffv [L2 T-2 ~> m2 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] @@ -1681,6 +1685,21 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US, endif endif + if (present(ADp) .and. (CS%id_h_diffu > 0)) then + h_diffu(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + h_diffu(I,j,k) = diffu(I,j,k) * ADp%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_diffu, h_diffu, CS%diag) + endif + if (present(ADp) .and. (CS%id_h_diffv > 0)) then + h_diffv(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + h_diffv(i,J,k) = diffv(i,J,k) * ADp%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_diffv, h_diffv, CS%diag) + endif + end subroutine horizontal_viscosity !> Allocates space for and calculates static variables used by horizontal_viscosity(). @@ -2397,6 +2416,20 @@ subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, MEKE, ADp) call safe_alloc_ptr(ADp%diag_hfrac_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke) endif + CS%id_h_diffu = register_diag_field('ocean_model', 'h_diffu', diag%axesCuL, Time, & + 'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if ((CS%id_h_diffu > 0) .and. (present(ADp))) then + call safe_alloc_ptr(ADp%diag_hu,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke) + endif + + CS%id_h_diffv = register_diag_field('ocean_model', 'h_diffv', diag%axesCvL, Time, & + 'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if ((CS%id_h_diffv > 0) .and. (present(ADp))) then + call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke) + endif + CS%id_intz_diffu_2d = register_diag_field('ocean_model', 'intz_diffu_2d', diag%axesCu1, Time, & 'Depth-integral of Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', & conversion=GV%H_to_m*US%L_T2_to_m_s2) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index d448751137..8d7444017a 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -125,6 +125,7 @@ module MOM_vert_friction integer :: 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_h_du_dt_visc = -1, id_h_dv_dt_visc = -1 integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1 !>@} @@ -213,6 +214,9 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & 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] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: h_du_dt_visc ! h x du_dt_visc [L2 T-2 ~> m2 s-2] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: h_dv_dt_visc ! h x h_dv_dt_visc [L2 T-2 ~> m2 s-2] + logical :: do_i(SZIB_(G)) logical :: DoStokesMixing @@ -499,6 +503,21 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, & deallocate(hf_dv_dt_visc_2d) endif + if (CS%id_h_du_dt_visc > 0) then + h_du_dt_visc(:,:,:) = 0.0 + do k=1,nz ; do j=js,je ; do I=Isq,Ieq + h_du_dt_visc(I,j,k) = ADp%du_dt_visc(I,j,k) * ADp%diag_hu(I,j,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_du_dt_visc, h_du_dt_visc, CS%diag) + endif + if (CS%id_h_dv_dt_visc > 0) then + h_dv_dt_visc(:,:,:) = 0.0 + do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie + h_dv_dt_visc(i,J,k) = ADp%dv_dt_visc(i,J,k) * ADp%diag_hv(i,J,k) + enddo ; enddo ; enddo + call post_data(CS%id_h_dv_dt_visc, h_dv_dt_visc, CS%diag) + endif + end subroutine vertvisc !> Calculate the fraction of momentum originally in a layer that remains @@ -1838,6 +1857,22 @@ subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, & call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz) endif + CS%id_h_du_dt_visc = register_diag_field('ocean_model', 'h_du_dt_visc', diag%axesCuL, Time, & + 'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_h_du_dt_visc > 0) then + call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz) + call safe_alloc_ptr(ADp%diag_hu,IsdB,IedB,jsd,jed,nz) + endif + + CS%id_h_dv_dt_visc = register_diag_field('ocean_model', 'h_dv_dt_visc', diag%axesCvL, Time, & + 'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', 'm2 s-2', & + conversion=GV%H_to_m*US%L_T2_to_m_s2) + if (CS%id_h_dv_dt_visc > 0) then + call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz) + call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,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)