Skip to content

Commit

Permalink
Bugfix in MOM_porous_barriers
Browse files Browse the repository at this point in the history 
Fix a bug that layer/interface weights may not be reset to 1.0 if no
calculations are needed. The bug has negligible impact for barotropic
applications but may affect baroclinic runs which are yet to performed.
  • Loading branch information
@herrwang0 @marshallward
herrwang0 authored and marshallward committed Oct 11, 2024
1 parent f90b071 commit 80d8b5f
Show file tree
Hide file tree
Showing 4 changed files with 45 additions and 29 deletions.
8 changes: 4 additions & 4 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2799,10 +2799,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
CS%time_in_cycle = 0.0 ; CS%time_in_thermo_cycle = 0.0

!allocate porous topography variables
allocate(CS%pbv%por_face_areaU(IsdB:IedB,jsd:jed,nz)) ; CS%pbv%por_face_areaU(:,:,:) = 1.0
allocate(CS%pbv%por_face_areaV(isd:ied,JsdB:JedB,nz)) ; CS%pbv%por_face_areaV(:,:,:) = 1.0
allocate(CS%pbv%por_layer_widthU(IsdB:IedB,jsd:jed,nz+1)) ; CS%pbv%por_layer_widthU(:,:,:) = 1.0
allocate(CS%pbv%por_layer_widthV(isd:ied,JsdB:JedB,nz+1)) ; CS%pbv%por_layer_widthV(:,:,:) = 1.0
allocate(CS%pbv%por_face_areaU(IsdB:IedB,jsd:jed,nz), source=1.0)
allocate(CS%pbv%por_face_areaV(isd:ied,JsdB:JedB,nz), source=1.0)
allocate(CS%pbv%por_layer_widthU(IsdB:IedB,jsd:jed,nz+1), source=1.0)
allocate(CS%pbv%por_layer_widthV(isd:ied,JsdB:JedB,nz+1), source=1.0)

! Use the Wright equation of state by default, unless otherwise specified
! Note: this line and the following block ought to be in a separate
Expand Down
64 changes: 40 additions & 24 deletions src/core/MOM_porous_barriers.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,16 +122,20 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt)
A_layer_prev(I,j) = A_layer
endif ; enddo ; enddo ; enddo
else
do k=nk,1,-1 ; do j=js,je ; do I=Isq,Ieq ; if (do_I(I,j)) then
call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), &
eta_u(I,j,K), A_layer, do_I(I,j))
if (eta_u(I,j,K) - (eta_u(I,j,K+1)+dz_min) > 0.0) then
pbv%por_face_areaU(I,j,k) = min(1.0, (A_layer - A_layer_prev(I,j)) / (eta_u(I,j,K) - eta_u(I,j,K+1)))
do k=nk,1,-1 ; do j=js,je ; do I=Isq,Ieq
if (do_I(I,j)) then
call calc_por_layer(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), &
eta_u(I,j,K), A_layer, do_I(I,j))
if (eta_u(I,j,K) - (eta_u(I,j,K+1)+dz_min) > 0.0) then
pbv%por_face_areaU(I,j,k) = min(1.0, (A_layer - A_layer_prev(I,j)) / (eta_u(I,j,K) - eta_u(I,j,K+1)))
else
pbv%por_face_areaU(I,j,k) = 0.0 ! use calc_por_interface() might be a better choice
endif
A_layer_prev(I,j) = A_layer
else
pbv%por_face_areaU(I,j,k) = 0.0 ! use calc_por_interface() might be a better choice
pbv%por_face_areaU(I,j,k) = 1.0
endif
A_layer_prev(I,j) = A_layer
endif ; enddo ; enddo ; enddo
enddo ; enddo ; enddo
endif

! v-points
Expand All @@ -154,16 +158,20 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt)
A_layer_prev(i,J) = A_layer
endif ; enddo ; enddo ; enddo
else
do k=nk,1,-1 ; do J=Jsq,Jeq ; do i=is,ie ; if (do_I(i,J)) then
call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), &
eta_v(i,J,K), A_layer, do_I(i,J))
if (eta_v(i,J,K) - (eta_v(i,J,K+1)+dz_min) > 0.0) then
pbv%por_face_areaV(i,J,k) = min(1.0, (A_layer - A_layer_prev(i,J)) / (eta_v(i,J,K) - eta_v(i,J,K+1)))
do k=nk,1,-1 ; do J=Jsq,Jeq ; do i=is,ie
if (do_I(i,J)) then
call calc_por_layer(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), &
eta_v(i,J,K), A_layer, do_I(i,J))
if (eta_v(i,J,K) - (eta_v(i,J,K+1)+dz_min) > 0.0) then
pbv%por_face_areaV(i,J,k) = min(1.0, (A_layer - A_layer_prev(i,J)) / (eta_v(i,J,K) - eta_v(i,J,K+1)))
else
pbv%por_face_areaV(i,J,k) = 0.0 ! use calc_por_interface() might be a better choice
endif
A_layer_prev(i,J) = A_layer
else
pbv%por_face_areaV(i,J,k) = 0.0 ! use calc_por_interface() might be a better choice
pbv%por_face_areaV(i,J,k) = 1.0
endif
A_layer_prev(i,J) = A_layer
endif ; enddo ; enddo ; enddo
enddo ; enddo ; enddo
endif

if (CS%debug) then
Expand Down Expand Up @@ -231,10 +239,14 @@ subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt)
eta_u(I,j,K), pbv%por_layer_widthU(I,j,K), do_I(I,j))
endif ; enddo ; enddo ; enddo
else
do K=1,nk+1 ; do j=js,je ; do I=Isq,Ieq ; if (do_I(I,j)) then
call calc_por_interface(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), &
eta_u(I,j,K), pbv%por_layer_widthU(I,j,K), do_I(I,j))
endif ; enddo ; enddo ; enddo
do K=1,nk+1 ; do j=js,je ; do I=Isq,Ieq
if (do_I(I,j)) then
call calc_por_interface(G%porous_DminU(I,j), G%porous_DmaxU(I,j), G%porous_DavgU(I,j), &
eta_u(I,j,K), pbv%por_layer_widthU(I,j,K), do_I(I,j))
else
pbv%por_layer_widthU(I,j,K) = 1.0
endif
enddo ; enddo ; enddo
endif

! v-points
Expand All @@ -249,10 +261,14 @@ subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt)
eta_v(i,J,K), pbv%por_layer_widthV(i,J,K), do_I(i,J))
endif ; enddo ; enddo ; enddo
else
do K=1,nk+1 ; do J=Jsq,Jeq ; do i=is,ie ; if (do_I(i,J)) then
call calc_por_interface(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), &
eta_v(i,J,K), pbv%por_layer_widthV(i,J,K), do_I(i,J))
endif ; enddo ; enddo ; enddo
do K=1,nk+1 ; do J=Jsq,Jeq ; do i=is,ie
if (do_I(i,J)) then
call calc_por_interface(G%porous_DminV(i,J), G%porous_DmaxV(i,J), G%porous_DavgV(i,J), &
eta_v(i,J,K), pbv%por_layer_widthV(i,J,K), do_I(i,J))
else
pbv%por_layer_widthV(i,J,K) = 1.0
endif
enddo ; enddo ; enddo
endif

if (CS%debug) then
Expand Down
1 change: 0 additions & 1 deletion src/core/MOM_variables.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -323,7 +323,6 @@ module MOM_variables
end type BT_cont_type

!> Container for grids modifying cell metric at porous barriers
! TODO: rename porous_barrier_type to porous_barrier_type
type, public :: porous_barrier_type
! Each of the following fields has nz layers.
real, allocatable :: por_face_areaU(:,:,:) !< fractional open area of U-faces [nondim]
Expand Down
1 change: 1 addition & 0 deletions src/initialization/MOM_fixed_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,6 +146,7 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF, write_geom, output_dir)
endif

! Read sub-grid scale topography parameters at velocity points used for porous barrier calculation
! TODO: The following routine call may eventually be merged as one of the CHANNEL_CONFIG options
call get_param(PF, mdl, "SUBGRID_TOPO_AT_VEL", read_porous_file, &
"If true, use variables from TOPO_AT_VEL_FILE as parameters for porous barrier.", &
default=.False.)
Expand Down

0 comments on commit 80d8b5f

Please sign in to comment.