Skip to content

Commit

Permalink
Momentum budget terms multiplied by fractional layer-thicknesses (#1163)
Browse files Browse the repository at this point in the history 
Squash merge:

* New diagnostics for barotropic momentum budget calculations are created. In these, different budget terms multiplied by fractional layer thicknesses are saved. Thus, these terms can be added over the whole to obtain the barotropic momentum budget. 'btstep' subroutine in MOM_barotropic module is modified to return fractional thicknesses. Pressure force acceleration multiplied by fractional thickness as diagnostics are added in this commit.

* More thickness weighted diagnostics

* Implemented Andrew Shao's suggestions

* More fractional thickness multiplied diagnostics

* Some barotropic diagnostics obtained from fractional thickness multiplied momentum budget terms

* Removed trailing spaces

* All fractional-thickness multiplied diagnostics implemented

* define diagnostic variables as pointers

* Shorter initialization of 2D arrays

* Modifications in vertical friction diagnostics

* Diagnostics initialization as pointers to save memory allocation

* Vertical friction module

* Removed commented lines and trailing spaces

* Diagnostic description change

* Fractional thickness-weighted diagnostics for acceleration due to relative vorticity and gradient of kinetic energy

* Modifications in vertical friction diagnostics (fractional thickness-weighted)

* Modifications in diag_hfrac_u and diag_hfrac_v calls

* Made allocation of hfrac at u/v points optional. These arrays are only allocated if any of the relevant diagnostics are called.

* Fractional-thickness weighted 3D diagnostics are now commented out as we do not the proper grid remapping option.

* Initialization of 2D diagnostic arrays changed from pointer type to allocatable array style.

* Trailing spaces removed and line lengths reduced to be under 120 characters

* Comment modified

Reviewed by Robert Hallberg <Robert.Hallberg@noaa.gov>
  • Loading branch information
@hmkhatri
hmkhatri authored Jul 29, 2020
1 parent df33724 commit 4c030e6
Show file tree
Hide file tree
Showing 7 changed files with 614 additions and 10 deletions.
154 changes: 154 additions & 0 deletions src/core/MOM_CoriolisAdv.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,10 @@ module MOM_CoriolisAdv
!>@{ Diagnostic IDs
integer :: id_rv = -1, id_PV = -1, id_gKEu = -1, id_gKEv = -1
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_hf_rvxu = -1, id_hf_rvxv = -1
integer :: id_hf_rvxu_2d = -1, id_hf_rvxv_2d = -1
!>@}
end type CoriolisAdv_CS

Expand Down Expand Up @@ -211,6 +215,16 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
real :: QUHeff,QVHeff ! More temporary variables [H L2 T-1 s-1 ~> m3 s-2 or kg s-2].
integer :: i, j, k, n, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz

! 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].
!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].
! 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.

! To work, the following fields must be set outside of the usual
! is to ie range before this subroutine is called:
! v(is-1:ie+2,js-1:je+1), u(is-1:ie+1,js-1:je+2), h(is-1:ie+2,js-1:je+2),
Expand Down Expand Up @@ -828,6 +842,82 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
if (CS%id_gKEv>0) call post_data(CS%id_gKEv, AD%gradKEv, CS%diag)
if (CS%id_rvxu > 0) call post_data(CS%id_rvxu, AD%rv_x_u, CS%diag)
if (CS%id_rvxv > 0) call post_data(CS%id_rvxv, AD%rv_x_v, CS%diag)

! Diagnostics for terms multiplied by fractional thicknesses

! 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.
!if (CS%id_hf_gKEu > 0) then
! allocate(hf_gKEu(G%IsdB:G%IedB,G%jsd:G%jed,G%ke))
! do k=1,nz ; do j=js,je ; do I=Isq,Ieq
! hf_gKEu(I,j,k) = AD%gradKEu(I,j,k) * AD%diag_hfrac_u(I,j,k)
! enddo ; enddo ; enddo
! call post_data(CS%id_hf_gKEu, hf_gKEu, CS%diag)
!endif

!if (CS%id_hf_gKEv > 0) then
! allocate(hf_gKEv(G%isd:G%ied,G%JsdB:G%JedB,G%ke))
! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
! hf_gKEv(i,J,k) = AD%gradKEv(i,J,k) * AD%diag_hfrac_v(i,J,k)
! enddo ; enddo ; enddo
! call post_data(CS%id_hf_gKEv, hf_gKEv, CS%diag)
!endif

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

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

!if (CS%id_hf_rvxv > 0) then
! allocate(hf_rvxv(G%IsdB:G%IedB,G%jsd:G%jed,G%ke))
! do k=1,nz ; do j=js,je ; do I=Isq,Ieq
! hf_rvxv(I,j,k) = AD%rv_x_v(I,j,k) * AD%diag_hfrac_u(I,j,k)
! enddo ; enddo ; enddo
! call post_data(CS%id_hf_rvxv, hf_rvxv, CS%diag)
!endif

!if (CS%id_hf_rvxu > 0) then
! allocate(hf_rvxu(G%isd:G%ied,G%JsdB:G%JedB,G%ke))
! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
! hf_rvxu(i,J,k) = AD%rv_x_u(i,J,k) * AD%diag_hfrac_v(i,J,k)
! enddo ; enddo ; enddo
! call post_data(CS%id_hf_rvxu, hf_rvxu, CS%diag)
!endif

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

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

end subroutine CorAdCalc
Expand Down Expand Up @@ -1087,6 +1177,70 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
'Zonal Acceleration from Relative Vorticity', 'm-1 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)
!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)
!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%axesCuL, Time, &
'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', &
'm-1 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%axesCvL, Time, &
'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Grad. Kinetic Energy', &
'm-1 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_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)
!if (CS%id_hf_rvxu > 0) then
! call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
! call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
!endif

!CS%id_hf_rvxv = register_diag_field('ocean_model', 'hf_rvxv', diag%axesCuL, Time, &
! 'Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
! 'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
!if (CS%id_hf_rvxv > 0) then
! call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz)
! call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
!endif

CS%id_hf_rvxu_2d = register_diag_field('ocean_model', 'hf_rvxu_2d', diag%axesCvL, Time, &
'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)
call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
endif

CS%id_hf_rvxv_2d = register_diag_field('ocean_model', 'hf_rvxv_2d', diag%axesCuL, Time, &
'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
'm-1 s-2', conversion=US%L_T2_to_m_s2)
if (CS%id_hf_rvxv_2d > 0) then
call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz)
call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
endif

end subroutine CoriolisAdv_init

!> Destructor for coriolisadv_cs
Expand Down
15 changes: 14 additions & 1 deletion src/core/MOM_barotropic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ module MOM_barotropic
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : BT_cont_type, alloc_bt_cont_type
use MOM_verticalGrid, only : verticalGrid_type
use MOM_variables, only : accel_diag_ptrs

implicit none ; private

Expand Down Expand Up @@ -405,7 +406,7 @@ module MOM_barotropic
subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, &
eta_PF_in, U_Cor, V_Cor, accel_layer_u, accel_layer_v, &
eta_out, uhbtav, vhbtav, G, GV, US, CS, &
visc_rem_u, visc_rem_v, etaav, OBC, BT_cont, eta_PF_start, &
visc_rem_u, visc_rem_v, etaav, ADp, OBC, BT_cont, eta_PF_start, &
taux_bot, tauy_bot, uh0, vh0, u_uh0, v_vh0)
type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
Expand Down Expand Up @@ -458,6 +459,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: visc_rem_v !< Ditto for meridional direction [nondim].
real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: etaav !< The free surface height or column mass
!! averaged over the barotropic integration [H ~> m or kg m-2].
type(accel_diag_ptrs), optional, pointer :: ADp !< Acceleration diagnostic pointers
type(ocean_OBC_type), optional, pointer :: OBC !< The open boundary condition structure.
type(BT_cont_type), optional, pointer :: BT_cont !< A structure with elements that describe
!! the effective open face areas as a function of barotropic
Expand Down Expand Up @@ -2583,6 +2585,17 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
if (CS%id_frhatv1 > 0) CS%frhatv1(:,:,:) = CS%frhatv(:,:,:)
endif

if ((present(ADp)) .and. (associated(ADp%diag_hfrac_u))) then
do k=1,nz ; do j=js,je ; do I=is-1,ie
ADp%diag_hfrac_u(I,j,k) = CS%frhatu(I,j,k)
enddo ; enddo ; enddo
endif
if ((present(ADp)) .and. (associated(ADp%diag_hfrac_v))) then
do k=1,nz ; do J=js-1,je ; do i=is,ie
ADp%diag_hfrac_v(i,J,k) = CS%frhatv(i,J,k)
enddo ; enddo ; enddo
endif

if (G%nonblocking_updates) then
if (find_etaav) call complete_group_pass(CS%pass_etaav, G%Domain)
call complete_group_pass(CS%pass_ubta_uhbta, G%Domain)
Expand Down
Loading

0 comments on commit 4c030e6

Please sign in to comment.