Skip to content

Commit

Permalink
Depth-integrated momentum budget diagnostics
Browse files Browse the repository at this point in the history
  • Loading branch information
@hmkhatri
hmkhatri committed Feb 25, 2021
1 parent 0393546 commit 73addc4
Show file tree
Hide file tree
Showing 4 changed files with 135 additions and 25 deletions.
96 changes: 86 additions & 10 deletions src/core/MOM_CoriolisAdv.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,8 +76,10 @@ module MOM_CoriolisAdv
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_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_intz_rvxu_2d = -1, id_intz_rvxv_2d = -1
!>@}
end type CoriolisAdv_CS

Expand Down Expand Up @@ -218,7 +220,9 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
! 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].
hf_rvxu_2d, hf_rvxv_2d, & ! Depth sum of hf_rvxu, hf_rvxv [L T-2 ~> m s-2].
intz_gKEu_2d, intz_gKEv_2d, & ! Depth integral of gKEu, gKEv [L T-2 ~> m s-2].
intz_rvxu_2d, intz_rvxv_2d ! Depth integral of rvxu, 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].
Expand Down Expand Up @@ -883,6 +887,26 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
deallocate(hf_gKEv_2d)
endif

if (CS%id_intz_gKEu_2d > 0) then
allocate(intz_gKEu_2d(G%IsdB:G%IedB,G%jsd:G%jed))
intz_gKEu_2d(:,:) = 0.0
do k=1,nz ; do j=js,je ; do I=Isq,Ieq
intz_gKEu_2d(I,j) = intz_gKEu_2d(I,j) + AD%gradKEu(I,j,k) * AD%diag_hu(I,j,k)
enddo ; enddo ; enddo
call post_data(CS%id_intz_gKEu_2d, intz_gKEu_2d, CS%diag)
deallocate(intz_gKEu_2d)
endif

if (CS%id_intz_gKEv_2d > 0) then
allocate(intz_gKEv_2d(G%isd:G%ied,G%JsdB:G%JedB))
intz_gKEv_2d(:,:) = 0.0
do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
intz_gKEv_2d(i,J) = intz_gKEv_2d(i,J) + AD%gradKEv(i,J,k) * AD%diag_hv(i,J,k)
enddo ; enddo ; enddo
call post_data(CS%id_intz_gKEv_2d, intz_gKEv_2d, CS%diag)
deallocate(intz_gKEv_2d)
endif

!if (CS%id_hf_rvxv > 0) then
! allocate(hf_rvxv(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
! do k=1,nz ; do j=js,je ; do I=Isq,Ieq
Expand Down Expand Up @@ -917,6 +941,26 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
enddo ; enddo ; enddo
call post_data(CS%id_hf_rvxu_2d, hf_rvxu_2d, CS%diag)
deallocate(hf_rvxu_2d)
endif

if (CS%id_intz_rvxv_2d > 0) then
allocate(intz_rvxv_2d(G%IsdB:G%IedB,G%jsd:G%jed))
intz_rvxv_2d(:,:) = 0.0
do k=1,nz ; do j=js,je ; do I=Isq,Ieq
intz_rvxv_2d(I,j) = intz_rvxv_2d(I,j) + AD%rv_x_v(I,j,k) * AD%diag_hu(I,j,k)
enddo ; enddo ; enddo
call post_data(CS%id_intz_rvxv_2d, intz_rvxv_2d, CS%diag)
deallocate(intz_rvxv_2d)
endif

if (CS%id_intz_rvxu_2d > 0) then
allocate(intz_rvxu_2d(G%isd:G%ied,G%JsdB:G%JedB))
intz_rvxu_2d(:,:) = 0.0
do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
intz_rvxu_2d(i,J) = intz_rvxu_2d(i,J) + AD%rv_x_u(i,J,k) * AD%diag_hv(i,J,k)
enddo ; enddo ; enddo
call post_data(CS%id_intz_rvxu_2d, intz_rvxu_2d, CS%diag)
deallocate(intz_rvxu_2d)
endif
endif

Expand Down Expand Up @@ -1163,53 +1207,69 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
'Potential Vorticity', 'm-1 s-1', conversion=GV%m_to_H*US%s_to_T)

CS%id_gKEu = register_diag_field('ocean_model', 'gKEu', diag%axesCuL, Time, &
'Zonal Acceleration from Grad. Kinetic Energy', 'm-1 s-2', conversion=US%L_T2_to_m_s2)
'Zonal Acceleration from Grad. Kinetic Energy', 'm s-2', conversion=US%L_T2_to_m_s2)
if (CS%id_gKEu > 0) call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz)

CS%id_gKEv = register_diag_field('ocean_model', 'gKEv', diag%axesCvL, Time, &
'Meridional Acceleration from Grad. Kinetic Energy', 'm-1 s-2', conversion=US%L_T2_to_m_s2)
'Meridional Acceleration from Grad. Kinetic Energy', 'm s-2', conversion=US%L_T2_to_m_s2)
if (CS%id_gKEv > 0) call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)

CS%id_rvxu = register_diag_field('ocean_model', 'rvxu', diag%axesCvL, Time, &
'Meridional Acceleration from Relative Vorticity', 'm-1 s-2', conversion=US%L_T2_to_m_s2)
'Meridional Acceleration from Relative Vorticity', 'm s-2', conversion=US%L_T2_to_m_s2)
if (CS%id_rvxu > 0) call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)

CS%id_rvxv = register_diag_field('ocean_model', 'rvxv', diag%axesCuL, Time, &
'Zonal Acceleration from Relative Vorticity', 'm-1 s-2', conversion=US%L_T2_to_m_s2)
'Zonal Acceleration from Relative Vorticity', 'm 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)
! 'm 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)
! 'm 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%axesCu1, Time, &
'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', &
'm-1 s-2', conversion=US%L_T2_to_m_s2)
'm 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%axesCv1, Time, &
'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Grad. Kinetic Energy', &
'm-1 s-2', conversion=US%L_T2_to_m_s2)
'm 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_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=US%L_T_to_m_s**2)
if (CS%id_intz_gKEu_2d > 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_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=US%L_T_to_m_s**2)
if (CS%id_intz_gKEv_2d > 0) then
call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
call safe_alloc_ptr(AD%diag_hv,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)
Expand All @@ -1227,7 +1287,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
!endif

CS%id_hf_rvxu_2d = register_diag_field('ocean_model', 'hf_rvxu_2d', diag%axesCv1, Time, &
'Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', &
'Depth-sum 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)
Expand All @@ -1242,6 +1302,22 @@ 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_intz_rvxu_2d = register_diag_field('ocean_model', 'intz_rvxu_2d', diag%axesCv1, Time, &
'Depth-integral of Meridional Acceleration from Relative Vorticity', &
'm-2 s-2', conversion=US%L_T_to_m_s**2)
if (CS%id_intz_rvxu_2d > 0) then
call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
call safe_alloc_ptr(AD%diag_hv,isd,ied,Jsd,JedB,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', &
'm-2 s-2', conversion=US%L_T_to_m_s**2)
if (CS%id_intz_rvxv_2d > 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

end subroutine CoriolisAdv_init

!> Destructor for coriolisadv_cs
Expand Down
24 changes: 12 additions & 12 deletions src/core/MOM_dynamics_split_RK2.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1452,13 +1452,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
if(CS%id_hf_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,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', 'm s-2', &
conversion=US%L_T2_to_m_s2)
'Depth-integral of Zonal Pressure Force Acceleration', 'm2 s-2', &
conversion=US%L_T_to_m_s**2)
if(CS%id_intz_PFu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)

CS%id_intz_PFv_2d = register_diag_field('ocean_model', 'intz_PFv_2d', diag%axesCv1, Time, &
'Depth-integral of Meridional Pressure Force Acceleration', 'm s-2', &
conversion=US%L_T2_to_m_s2)
'Depth-integral of Meridional Pressure Force Acceleration', 'm2 s-2', &
conversion=US%L_T_to_m_s**2)
if(CS%id_intz_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)

CS%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, Time, &
Expand All @@ -1472,13 +1472,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
if(CS%id_hf_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,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', 'm s-2', &
conversion=US%L_T2_to_m_s2)
'Depth-integral of Zonal Coriolis and Advective Acceleration', 'm2 s-2', &
conversion=US%L_T_to_m_s**2)
if(CS%id_intz_CAu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)

CS%id_intz_CAv_2d = register_diag_field('ocean_model', 'intz_CAv_2d', diag%axesCv1, Time, &
'Depth-integral of Meridional Coriolis and Advective Acceleration', 'm s-2', &
conversion=US%L_T2_to_m_s2)
'Depth-integral of Meridional Coriolis and Advective Acceleration', 'm2 s-2', &
conversion=US%L_T_to_m_s**2)
if(CS%id_intz_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)

CS%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, Time, &
Expand Down Expand Up @@ -1512,13 +1512,13 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
if(CS%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,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', 'm s-2', &
conversion=US%L_T2_to_m_s2)
'Depth-integral of Barotropic Anomaly Zonal Acceleration', 'm2 s-2', &
conversion=US%L_T_to_m_s**2)
if(CS%id_intz_u_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)

CS%id_intz_v_BT_accel_2d = register_diag_field('ocean_model', 'intz_v_BT_accel_2d', diag%axesCv1, Time, &
'Depth-integral of Barotropic Anomaly Meridional Acceleration', 'm s-2', &
conversion=US%L_T2_to_m_s2)
'Depth-integral of Barotropic Anomaly Meridional Acceleration', 'm2 s-2', &
conversion=US%L_T_to_m_s**2)
if(CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,Jsd,JedB,nz)

id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
Expand Down
2 changes: 1 addition & 1 deletion src/core/MOM_variables.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -190,7 +190,7 @@ module MOM_variables
real, pointer :: diag_hfrac_v(:,:,:) => NULL() !< Fractional layer thickness at v points
real, pointer :: diag_hu(:,:,:) => NULL() !< layer thickness at u points
real, pointer :: diag_hv(:,:,:) => NULL() !< 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.
Expand Down
38 changes: 36 additions & 2 deletions src/parameterizations/lateral/MOM_hor_visc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -188,6 +188,7 @@ module MOM_hor_visc
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_intz_diffu_2d = -1, id_intz_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
Expand Down Expand Up @@ -276,8 +277,8 @@ 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, allocatable, dimension(:,:) :: hf_diffu_2d, hf_diffv_2d ! Depth sum of hf_diffu, hf_diffv [L T-2 ~> m s-2]
real, allocatable, dimension(:,:) :: intz_diffu_2d, intz_diffv_2d ! Integral of diffu, diffv [L2 T-2 ~> m2 s-2]

real, dimension(SZIB_(G),SZJB_(G)) :: &
dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1]
Expand Down Expand Up @@ -1662,6 +1663,25 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
deallocate(hf_diffv_2d)
endif

if (present(ADp) .and. (CS%id_intz_diffu_2d > 0)) then
allocate(intz_diffu_2d(G%IsdB:G%IedB,G%jsd:G%jed))
intz_diffu_2d(:,:) = 0.0
do k=1,nz ; do j=js,je ; do I=Isq,Ieq
intz_diffu_2d(I,j) = intz_diffu_2d(I,j) + diffu(I,j,k) * ADp%diag_hu(I,j,k)
enddo ; enddo ; enddo
call post_data(CS%id_intz_diffu_2d, intz_diffu_2d, CS%diag)
deallocate(intz_diffu_2d)
endif
if (present(ADp) .and. (CS%id_intz_diffv_2d > 0)) then
allocate(intz_diffv_2d(G%isd:G%ied,G%JsdB:G%JedB))
intz_diffv_2d(:,:) = 0.0
do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
intz_diffv_2d(i,J) = intz_diffv_2d(i,J) + diffv(i,J,k) * ADp%diag_hv(i,J,k)
enddo ; enddo ; enddo
call post_data(CS%id_intz_diffv_2d, intz_diffv_2d, CS%diag)
deallocate(intz_diffv_2d)
endif

end subroutine horizontal_viscosity

!> Allocates space for and calculates static variables used by horizontal_viscosity().
Expand Down Expand Up @@ -2378,6 +2398,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_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=US%L_T_to_m_s**2)
if ((CS%id_intz_diffu_2d > 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_intz_diffv_2d = register_diag_field('ocean_model', 'intz_diffv_2d', diag%axesCv1, Time, &
'Depth-integral of Meridional Acceleration from Horizontal Viscosity', 'm2 s-2', &
conversion=US%L_T_to_m_s**2)
if ((CS%id_intz_diffv_2d > 0) .and. (present(ADp))) then
call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%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, &
Expand Down

0 comments on commit 73addc4

Please sign in to comment.