Skip to content

Commit

Permalink
Merge pull request #1398 from hmkhatri/add_diagnostic
Browse files Browse the repository at this point in the history 
3D thickness x momentum diagnostics
  • Loading branch information
@marshallward
marshallward authored May 18, 2021
2 parents 75a8549 + cf3b41f commit 61330e4
Show file tree
Hide file tree
Showing 5 changed files with 298 additions and 7 deletions.
80 changes: 78 additions & 2 deletions src/core/MOM_CoriolisAdv.F90
View file
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
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
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

0 comments on commit 61330e4

Please sign in to comment.