Skip to content

Commit

Permalink
Merge pull request #64 from NOAA-GFDL/dev/gfdl
Browse files Browse the repository at this point in the history 
Sync with NOAA-GFDL dev/gfdl
  • Loading branch information
@wrongkindofdoctor
wrongkindofdoctor authored Jul 31, 2020
2 parents de7f95a + 3d557d9 commit 0cf3cb9
Show file tree
Hide file tree
Showing 12 changed files with 824 additions and 110 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
30 changes: 25 additions & 5 deletions 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 @@ -687,6 +689,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
integer :: is, ie, js, je, nz, Isq, Ieq, Jsq, Jeq
integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
integer :: ioff, joff
integer :: l_seg

if (.not.associated(CS)) call MOM_error(FATAL, &
"btstep: Module MOM_barotropic must be initialized before it is used.")
Expand Down Expand Up @@ -2324,9 +2327,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
if (CS%BT_OBC%apply_u_OBCs) then ! copy back the value for u-points on the boundary.
!GOMP parallel do default(shared)
do j=js,je ; do I=is-1,ie
if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then
l_seg = OBC%segnum_u(I,j)
if (l_seg == OBC_NONE) cycle

if (OBC%segment(l_seg)%direction == OBC_DIRECTION_E) then
e_anom(i+1,j) = e_anom(i,j)
elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then
elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_W) then
e_anom(i,j) = e_anom(i+1,j)
endif
enddo ; enddo
Expand All @@ -2335,9 +2341,12 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
if (CS%BT_OBC%apply_v_OBCs) then ! copy back the value for v-points on the boundary.
!GOMP parallel do default(shared)
do J=js-1,je ; do I=is,ie
if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then
l_seg = OBC%segnum_v(i,J)
if (l_seg == OBC_NONE) cycle

if (OBC%segment(l_seg)%direction == OBC_DIRECTION_N) then
e_anom(i,j+1) = e_anom(i,j)
elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then
elseif (OBC%segment(l_seg)%direction == OBC_DIRECTION_S) then
e_anom(i,j) = e_anom(i,j+1)
endif
enddo ; enddo
Expand Down Expand Up @@ -2583,6 +2592,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 0cf3cb9

Please sign in to comment.