Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

(+) porous topography implementation #3

Merged
merged 9 commits into from
Nov 30, 2021
Merged

(+) porous topography implementation #3

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
8b5e10a
(+) porous topography implementation (#3)
sditkovsky Nov 17, 2021
bc82609
small fixes to porous topo code
sditkovsky Nov 19, 2021
e717e46
preserve dimensions in porous topo interpolation
sditkovsky Nov 24, 2021
04d932e
Remove Z_to_m rescaling in porous module
sditkovsky Nov 24, 2021
a98c3ef
small bug fix
sditkovsky Nov 24, 2021
66e1493
small bug fix 2
sditkovsky Nov 24, 2021
62786c4
fix merge issue
sditkovsky Nov 24, 2021
ec84be5
change conversion H_to_MKS to H_to_m for u/veffA
sditkovsky Nov 29, 2021
3f38239
Merge branch 'dev/gfdl' into porous_topo_sjd
sditkovsky Nov 29, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 39 additions & 7 deletions src/core/MOM.F90
View file
Open in desktop
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
Hallberg-NOAA marked this conversation as resolved.
Show resolved Hide resolved
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
Open in desktop
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
Open in desktop
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