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

3D Thickness multiplied momentum diagnostics #1397

Closed
wants to merge 27 commits into from
Closed

3D Thickness multiplied momentum diagnostics #1397

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
27 commits
Select commit Hold shift + click to select a range
3acd777
Draft. Compiles but FPE in advection
adcroft Dec 23, 2020
ca51707
Fix FPE problems with NW2 tracers
adcroft Jan 3, 2021
77c67fc
Allow NW2 tracesr to be added mid-run
adcroft Jan 3, 2021
b333739
Cleanup NW2 registration
adcroft Jan 3, 2021
adb3554
Fixed NW2 tracer restarts
adcroft Jan 3, 2021
27cd2ed
Merge pull request #6 from NOAA-GFDL/dev-gfdl-main-candidate-2021-02-19
adcroft Feb 22, 2021
2f8467a
Merge branch 'nw2_tracers' of https://github.com/neerajabhamidipati/M…
adcroft Feb 22, 2021
d20ecf7
Remove use of parameter for NTR_MAX in NW2
adcroft Feb 22, 2021
8f6ec58
Replace NW2 2x10 tracers with 3x3
adcroft Feb 22, 2021
280d8b1
Merge branch 'neerajabhamidipati-nw2_tracers' into dev/cpt
adcroft Feb 22, 2021
48239bb
MOM_variables merge
hmkhatri Mar 5, 2021
373da29
Return layer thicknesses from barotropic module
hmkhatri Mar 5, 2021
a6055e1
Some thickness multiplied 3D diagnostics
hmkhatri Mar 5, 2021
1616df4
Thickness multiplied diagnostics for stress terms
hmkhatri Mar 5, 2021
c987de7
typo correct
hmkhatri Mar 5, 2021
6779f46
More diagnostics
hmkhatri Mar 5, 2021
b8c1717
Thickness multiplied momentum budget terms
hmkhatri Mar 5, 2021
c01a3e7
resolve merge conflit and merge dev/gfdl
hmkhatri Apr 26, 2021
3dd66f5
remove duplicate initialization of array
hmkhatri Apr 26, 2021
3f5682f
remove conflicts
hmkhatri Apr 27, 2021
709165b
remove conflicts
hmkhatri Apr 28, 2021
f62e5b7
conflict resolve
hmkhatri Apr 28, 2021
b89cb21
Merge branch 'dev/gfdl' into cpt_TWA
hmkhatri Apr 28, 2021
8ec6202
Merge branch 'dev/gfdl' into cpt_TWA
hmkhatri May 3, 2021
26ce31d
submodule correct
hmkhatri May 4, 2021
0b94041
cvmix correct
hmkhatri May 4, 2021
e663258
Merge branch 'dev/gfdl' into cpt_TWA
hmkhatri May 10, 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
68 changes: 68 additions & 0 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,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].
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down
85 changes: 84 additions & 1 deletion 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,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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
39 changes: 39 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, 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]
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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, &
Expand Down
Loading