Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

3D thickness x momentum diagnostics #1398

Merged
merged 4 commits into from
May 18, 2021
Merged

3D thickness x momentum diagnostics #1398

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f8854e7
3D thickness x momentum diagnostics
hmkhatri May 11, 2021
f1b2320
Added description of diagnostics arrays and corrected units in comments
hmkhatri May 17, 2021
aa32127
Comments about new diagnostics
hmkhatri May 17, 2021
cf3b41f
Merge branch 'dev/gfdl' into add_diagnostic
marshallward May 18, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
80 changes: 78 additions & 2 deletions src/core/MOM_CoriolisAdv.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -230,9 +232,13 @@ 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, allocatable, dimension(:,:,:) :: h_gKEu, h_gKEv ! h x gKEu, h x gKEv [H L T-2 ~> m2 s-2].
real, allocatable, dimension(:,:,:) :: h_rvxv, h_rvxu ! h x rvxv, h x rvxu [H L 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].
real, dimension(SZIB_(G),SZJ_(G)) :: intz_gKEu_2d, intz_rvxv_2d ! [H L T-2 ~> m2 s-2].
real, dimension(SZI_(G),SZJB_(G)) :: intz_gKEv_2d, intz_rvxu_2d ! [H L T-2 ~> m2 s-2].

! To work, the following fields must be set outside of the usual
! is to ie range before this subroutine is called:
Expand Down Expand Up @@ -955,6 +961,44 @@ 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
allocate(h_gKEu(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
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)
deallocate(h_gKEu)
endif
if (CS%id_h_gKEv > 0) then
allocate(h_gKEv(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
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)
deallocate(h_gKEv)
endif

if (CS%id_h_rvxv > 0) then
allocate(h_rvxv(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
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)
deallocate(h_rvxv)
endif
if (CS%id_h_rvxu > 0) then
allocate(h_rvxu(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
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)
deallocate(h_rvxu)
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
Expand Down Expand Up @@ -1262,6 +1306,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)
Expand All @@ -1270,6 +1322,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)
Expand Down Expand Up @@ -1310,6 +1370,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)
Expand All @@ -1318,6 +1386,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)
Expand Down
102 changes: 99 additions & 3 deletions src/core/MOM_dynamics_split_RK2.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,16 +160,19 @@ 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

! 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_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
!>@}
Expand Down Expand Up @@ -339,11 +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)) :: &
intz_PFu_2d, intz_CAu_2d, intz_u_BT_accel_2d ! [L2 T-2 ~> m2 s-2].
! Diagnostics for thickness x momentum budget terms
real, allocatable, dimension(:,:,:) :: &
h_PFu, h_PFv, & ! Pressure force accel. x thickness [H L T-2 ~> m2 s-2].
h_CAu, h_CAv, & ! Coriolis force accel. x thickness [H L T-2 ~> m2 s-2].
h_u_BT_accel, h_v_BT_accel ! barotropic correction accel. x thickness [H L T-2 ~> m2 s-2].

! Dignostics for layer-sum of thickness x momentum budget terms
real, dimension(SZIB_(G),SZJ_(G)) :: &
intz_PFu_2d, intz_CAu_2d, intz_u_BT_accel_2d ! [H L 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].
intz_PFv_2d, intz_CAv_2d, intz_v_BT_accel_2d ! [H L T-2 ~> m2 s-2].

real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s].

Expand Down Expand Up @@ -940,6 +949,25 @@ 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
allocate(h_PFu(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
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)
deallocate(h_PFu)
endif
if (CS%id_h_PFv > 0) then
allocate(h_PFv(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
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)
deallocate(h_PFv)
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
Expand Down Expand Up @@ -988,6 +1016,25 @@ 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
allocate(h_CAu(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
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)
deallocate(h_CAu)
endif
if (CS%id_h_CAv > 0) then
allocate(h_CAv(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
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)
deallocate(h_CAv)
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
Expand Down Expand Up @@ -1036,6 +1083,25 @@ 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
allocate(h_u_BT_accel(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
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)
deallocate(h_u_BT_accel)
endif
if (CS%id_h_v_BT_accel > 0) then
allocate(h_v_BT_accel(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
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)
deallocate(h_v_BT_accel)
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)
Expand Down Expand Up @@ -1442,6 +1508,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)
Expand All @@ -1462,6 +1538,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)
Expand Down Expand Up @@ -1502,6 +1588,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)
Expand Down
43 changes: 43 additions & 0 deletions src/diagnostics/MOM_diagnostics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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, allocatable, dimension(:,:,:) :: h_du_dt ! h x dudt [H L T-2 ~> m2 s-2]
real, allocatable, dimension(:,:,:) :: h_dv_dt ! h x dvdt [H L 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]
Expand Down Expand Up @@ -325,6 +329,25 @@ 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
allocate(h_du_dt(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
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)
deallocate(h_du_dt)
endif
if (CS%id_h_dv_dt > 0) then
allocate(h_dv_dt(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
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)
deallocate(h_dv_dt)
endif

call diag_restore_grids(CS%diag)

call calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, GV, US, CS)
Expand Down Expand Up @@ -1789,6 +1812,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, &
Expand Down
Loading