Skip to content

Commit

Permalink
(+) porous topography implementation (#3)
Browse files Browse the repository at this point in the history 
* reads in porous topography parameters from CHANNEL_LIST_FILE

*new module to compute curve fit for porous topography

*porous constraints used to modify continuity_PPM, CoriolisAdv, and Rayleigh bottom channel drag
  • Loading branch information
@sditkovsky
sditkovsky authored Nov 17, 2021
1 parent 688bff9 commit 8b5e10a
Show file tree
Hide file tree
Showing 14 changed files with 545 additions and 110 deletions.
46 changes: 39 additions & 7 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,7 @@ module MOM
use MOM_unit_scaling, only : unit_scale_type, unit_scaling_init
use MOM_unit_scaling, only : unit_scaling_end, fix_restart_unit_scaling
use MOM_variables, only : surface, allocate_surface_state, deallocate_surface_state
use MOM_variables, only : thermo_var_ptrs, vertvisc_type
use MOM_variables, only : thermo_var_ptrs, vertvisc_type, porous_barrier_ptrs
use MOM_variables, only : accel_diag_ptrs, cont_diag_ptrs, ocean_internal_state
use MOM_variables, only : rotate_surface_state
use MOM_verticalGrid, only : verticalGrid_type, verticalGridInit, verticalGridEnd
Expand All @@ -136,6 +136,8 @@ module MOM
use MOM_wave_interface, only : wave_parameters_CS, waves_end
use MOM_wave_interface, only : Update_Stokes_Drift

use MOM_porous_barriers, only : porous_widths

! ODA modules
use MOM_oda_driver_mod, only : ODA_CS, oda, init_oda, oda_end
use MOM_oda_driver_mod, only : set_prior_tracer, set_analysis_time, apply_oda_tracer_increments
Expand Down Expand Up @@ -399,6 +401,15 @@ module MOM
type(ODA_CS), pointer :: odaCS => NULL() !< a pointer to the control structure for handling
!! ensemble model state vectors and data assimilation
!! increments and priors
type(porous_barrier_ptrs) :: pbv !< porous barrier fractional cell metrics
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) &
:: por_face_areaU !< fractional open area of U-faces [nondim]
real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) &
:: por_face_areaV !< fractional open area of V-faces [nondim]
real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NK_INTERFACE_) &
:: por_layer_widthU !< fractional open width of U-faces [nondim]
real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NK_INTERFACE_) &
:: por_layer_widthV !< fractional open width of V-faces [nondim]
type(particles), pointer :: particles => NULL() !<Lagrangian particles
end type MOM_control_struct

Expand Down Expand Up @@ -1019,6 +1030,8 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB

real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)+1) :: eta_por ! layer interface heights
!! for porous topo. [Z ~> m or 1/eta_to_m]
G => CS%G ; GV => CS%GV ; US => CS%US ; IDs => CS%IDs
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
Expand Down Expand Up @@ -1047,13 +1060,16 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
call diag_update_remap_grids(CS%diag)
endif

!update porous barrier fractional cell metrics
call porous_widths(h, CS%tv, G, GV, US, eta_por, CS%pbv)

! The bottom boundary layer properties need to be recalculated.
if (bbl_time_int > 0.0) then
call enable_averages(bbl_time_int, &
Time_local + real_to_time(US%T_to_s*(bbl_time_int-dt)), CS%diag)
! Calculate the BBL properties and store them inside visc (u,h).
call cpu_clock_begin(id_clock_BBL_visc)
call set_viscous_BBL(CS%u, CS%v, CS%h, CS%tv, CS%visc, G, GV, US, CS%set_visc_CSp)
call set_viscous_BBL(CS%u, CS%v, CS%h, CS%tv, CS%visc, G, GV, US, CS%set_visc_CSp, CS%pbv)
call cpu_clock_end(id_clock_BBL_visc)
if (showCallTree) call callTree_wayPoint("done with set_viscous_BBL (step_MOM)")
call disable_averaging(CS%diag)
Expand All @@ -1076,7 +1092,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
call step_MOM_dyn_split_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, &
p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, &
CS%eta_av_bc, G, GV, US, CS%dyn_split_RK2_CSp, calc_dtbt, CS%VarMix, &
CS%MEKE, CS%thickness_diffuse_CSp, waves=waves)
CS%MEKE, CS%thickness_diffuse_CSp, CS%pbv, waves=waves)
if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_split (step_MOM)")

elseif (CS%do_dynamics) then ! ------------------------------------ not SPLIT
Expand All @@ -1090,11 +1106,11 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, &
if (CS%use_RK2) then
call step_MOM_dyn_unsplit_RK2(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, &
p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, &
CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE)
CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_RK2_CSp, CS%VarMix, CS%MEKE, CS%pbv)
else
call step_MOM_dyn_unsplit(u, v, h, CS%tv, CS%visc, Time_local, dt, forces, &
p_surf_begin, p_surf_end, CS%uh, CS%vh, CS%uhtr, CS%vhtr, &
CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, Waves=Waves)
CS%eta_av_bc, G, GV, US, CS%dyn_unsplit_CSp, CS%VarMix, CS%MEKE, CS%pbv, Waves=Waves)
endif
if (showCallTree) call callTree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)")

Expand Down Expand Up @@ -1296,6 +1312,9 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
integer :: halo_sz ! The size of a halo where data must be valid.
integer :: i, j, k, is, ie, js, je, nz

real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)+1) :: eta_por ! layer interface heights
!! for porous topo. [Z ~> m or 1/eta_to_m]

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
showCallTree = callTree_showQuery()
if (showCallTree) call callTree_enter("step_MOM_thermo(), MOM.F90")
Expand Down Expand Up @@ -1331,7 +1350,9 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
! DIABATIC_FIRST=True. Otherwise diabatic() is called after the dynamics
! and set_viscous_BBL is called as a part of the dynamic stepping.
call cpu_clock_begin(id_clock_BBL_visc)
call set_viscous_BBL(u, v, h, tv, CS%visc, G, GV, US, CS%set_visc_CSp)
!update porous barrier fractional cell metrics
call porous_widths(h, CS%tv, G, GV, US, eta_por, CS%pbv)
call set_viscous_BBL(u, v, h, tv, CS%visc, G, GV, US, CS%set_visc_CSp, CS%pbv)
call cpu_clock_end(id_clock_BBL_visc)
if (showCallTree) call callTree_wayPoint("done with set_viscous_BBL (step_MOM_thermo)")
endif
Expand Down Expand Up @@ -2330,6 +2351,13 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
ALLOC_(CS%eta_av_bc(isd:ied,jsd:jed)) ; CS%eta_av_bc(:,:) = 0.0 ! -G%Z_ref
CS%time_in_cycle = 0.0 ; CS%time_in_thermo_cycle = 0.0

!allocate porous topography variables
ALLOC_(CS%por_face_areaU(IsdB:IedB,jsd:jed,nz)) ; CS%por_face_areaU(:,:,:) = 1.0
ALLOC_(CS%por_face_areaV(isd:ied,JsdB:JedB,nz)) ; CS%por_face_areaV(:,:,:) = 1.0
ALLOC_(CS%por_layer_widthU(IsdB:IedB,jsd:jed,nz+1)) ; CS%por_layer_widthU(:,:,:) = 1.0
ALLOC_(CS%por_layer_widthV(isd:ied,JsdB:JedB,nz+1)) ; CS%por_layer_widthV(:,:,:) = 1.0
CS%pbv%por_face_areaU => CS%por_face_areaU; CS%pbv%por_face_areaV=> CS%por_face_areaV
CS%pbv%por_layer_widthU => CS%por_layer_widthU; CS%pbv%por_layer_widthV => CS%por_layer_widthV
! Use the Wright equation of state by default, unless otherwise specified
! Note: this line and the following block ought to be in a separate
! initialization routine for tv.
Expand Down Expand Up @@ -2647,7 +2675,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
CS%dt, CS%ADp, CS%CDp, MOM_internal_state, CS%VarMix, CS%MEKE, &
CS%thickness_diffuse_CSp, &
CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, &
CS%visc, dirs, CS%ntrunc, calc_dtbt=calc_dtbt, cont_stencil=CS%cont_stencil)
CS%visc, dirs, CS%ntrunc, CS%pbv, calc_dtbt=calc_dtbt, cont_stencil=CS%cont_stencil)
if (CS%dtbt_reset_period > 0.0) then
CS%dtbt_reset_interval = real_to_time(CS%dtbt_reset_period)
! Set dtbt_reset_time to be the next even multiple of dtbt_reset_interval.
Expand Down Expand Up @@ -3581,6 +3609,10 @@ subroutine MOM_end(CS)

if (CS%use_ALE_algorithm) call ALE_end(CS%ALE_CSp)

!deallocate porous topography variables
DEALLOC_(CS%por_face_areaU) ; DEALLOC_(CS%por_face_areaV)
DEALLOC_(CS%por_layer_widthU) ; DEALLOC_(CS%por_layer_widthV)

! NOTE: Allocated in PressureForce_FV_Bouss
if (associated(CS%tv%varT)) deallocate(CS%tv%varT)

Expand Down
20 changes: 11 additions & 9 deletions src/core/MOM_CoriolisAdv.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ module MOM_CoriolisAdv
use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S
use MOM_string_functions, only : uppercase
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : accel_diag_ptrs
use MOM_variables, only : accel_diag_ptrs, porous_barrier_ptrs
use MOM_verticalGrid, only : verticalGrid_type

implicit none ; private
Expand Down Expand Up @@ -117,7 +117,7 @@ module MOM_CoriolisAdv
contains

!> Calculates the Coriolis and momentum advection contributions to the acceleration.
subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS, pbv)
type(ocean_grid_type), intent(in) :: G !< Ocen grid structure
type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]
Expand All @@ -135,6 +135,7 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
type(accel_diag_ptrs), intent(inout) :: AD !< Storage for acceleration diagnostics
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(CoriolisAdv_CS), pointer :: CS !< Control structure for MOM_CoriolisAdv
type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics

! Local variables
real, dimension(SZIB_(G),SZJB_(G)) :: &
Expand Down Expand Up @@ -285,7 +286,8 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
enddo ; enddo

!$OMP parallel do default(private) shared(u,v,h,uh,vh,CAu,CAv,G,GV,CS,AD,Area_h,Area_q,&
!$OMP RV,PV,is,ie,js,je,Isq,Ieq,Jsq,Jeq,nz,vol_neglect,h_tiny,OBC,eps_vel)
!$OMP RV,PV,is,ie,js,je,Isq,Ieq,Jsq,Jeq,nz,vol_neglect,h_tiny,OBC,eps_vel, &
!$OMP pbv)
do k=1,nz

! Here the second order accurate layer potential vorticities, q,
Expand All @@ -306,10 +308,10 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)

if (CS%Coriolis_En_Dis) then
do j=Jsq,Jeq+1 ; do I=is-1,ie
uh_center(I,j) = 0.5 * (G%dy_Cu(I,j) * u(I,j,k)) * (h(i,j,k) + h(i+1,j,k))
uh_center(I,j) = 0.5 * ((G%dy_Cu(I,j)*pbv%por_face_areaU(I,j,k)) * u(I,j,k)) * (h(i,j,k) + h(i+1,j,k))
enddo ; enddo
do J=js-1,je ; do i=Isq,Ieq+1
vh_center(i,J) = 0.5 * (G%dx_Cv(i,J) * v(i,J,k)) * (h(i,j,k) + h(i,j+1,k))
vh_center(i,J) = 0.5 * ((G%dx_Cv(i,J)*pbv%por_face_areaV(i,J,k)) * v(i,J,k)) * (h(i,j,k) + h(i,j+1,k))
enddo ; enddo
endif

Expand Down Expand Up @@ -352,9 +354,9 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
if (CS%Coriolis_En_Dis) then
do i = max(Isq-1,OBC%segment(n)%HI%isd), min(Ieq+2,OBC%segment(n)%HI%ied)
if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
vh_center(i,J) = G%dx_Cv(i,J) * v(i,J,k) * h(i,j,k)
vh_center(i,J) = (G%dx_Cv(i,J)*pbv%por_face_areaV(i,J,k)) * v(i,J,k) * h(i,j,k)
else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
vh_center(i,J) = G%dx_Cv(i,J) * v(i,J,k) * h(i,j+1,k)
vh_center(i,J) = (G%dx_Cv(i,J)*pbv%por_face_areaV(i,J,k)) * v(i,J,k) * h(i,j+1,k)
endif
enddo
endif
Expand Down Expand Up @@ -391,9 +393,9 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
if (CS%Coriolis_En_Dis) then
do j = max(Jsq-1,OBC%segment(n)%HI%jsd), min(Jeq+2,OBC%segment(n)%HI%jed)
if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
uh_center(I,j) = G%dy_Cu(I,j) * u(I,j,k) * h(i,j,k)
uh_center(I,j) = (G%dy_Cu(I,j)*pbv%por_face_areaU(I,j,k)) * u(I,j,k) * h(i,j,k)
else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
uh_center(I,j) = G%dy_Cu(I,j) * u(I,j,k) * h(i+1,j,k)
uh_center(I,j) = (G%dy_Cu(I,j)*pbv%por_face_areaU(I,j,k)) * u(I,j,k) * h(i+1,j,k)
endif
enddo
endif
Expand Down
7 changes: 4 additions & 3 deletions src/core/MOM_continuity.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ module MOM_continuity
use MOM_grid, only : ocean_grid_type
use MOM_open_boundary, only : ocean_OBC_type
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : BT_cont_type
use MOM_variables, only : BT_cont_type, porous_barrier_ptrs
use MOM_verticalGrid, only : verticalGrid_type

implicit none ; private
Expand All @@ -39,7 +39,7 @@ module MOM_continuity

!> Time steps the layer thicknesses, using a monotonically limited, directionally split PPM scheme,
!! based on Lin (1994).
subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vhbt, &
subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhbt, vhbt, &
visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont)
type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure.
type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
Expand All @@ -61,6 +61,7 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vhbt,
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(continuity_CS), pointer :: CS !< Control structure for mom_continuity.
type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
type(porous_barrier_ptrs), intent(in) :: pbv !< porous barrier fractional cell metrics
real, dimension(SZIB_(G),SZJ_(G)), &
optional, intent(in) :: uhbt !< The vertically summed volume
!! flux through zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1].
Expand Down Expand Up @@ -95,7 +96,7 @@ subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vhbt,
" one must be present in call to continuity.")

if (CS%continuity_scheme == PPM_SCHEME) then
call continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS%PPM_CSp, OBC, uhbt, vhbt, &
call continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS%PPM_CSp, OBC, pbv, uhbt, vhbt, &
visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont=BT_cont)
else
call MOM_error(FATAL, "continuity: Unrecognized value of continuity_scheme")
Expand Down
Loading

0 comments on commit 8b5e10a

Please sign in to comment.