Skip to content

Commit

Permalink
Merge pull request #151 from Hallberg-NOAA/drag_as_body_force
Browse files Browse the repository at this point in the history 
+Add option to apply bottom drag as a body force
  • Loading branch information
@marshallward
marshallward authored Jul 7, 2022
2 parents 1c0e1f8 + e9e0f41 commit 5cadb72
Show file tree
Hide file tree
Showing 7 changed files with 212 additions and 156 deletions.
4 changes: 2 additions & 2 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3515,13 +3515,13 @@ subroutine extract_surface_state(CS, sfc_state_in)
enddo ! end of j loop
endif ! melt_potential

if (allocated(sfc_state%taux_shelf) .and. associated(CS%visc%taux_shelf)) then
if (allocated(sfc_state%taux_shelf) .and. allocated(CS%visc%taux_shelf)) then
!$OMP parallel do default(shared)
do j=js,je ; do I=is-1,ie
sfc_state%taux_shelf(I,j) = CS%visc%taux_shelf(I,j)
enddo ; enddo
endif
if (allocated(sfc_state%tauy_shelf) .and. associated(CS%visc%tauy_shelf)) then
if (allocated(sfc_state%tauy_shelf) .and. allocated(CS%visc%tauy_shelf)) then
!$OMP parallel do default(shared)
do J=js-1,je ; do i=is,ie
sfc_state%tauy_shelf(i,J) = CS%visc%tauy_shelf(i,J)
Expand Down
46 changes: 23 additions & 23 deletions src/core/MOM_variables.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -220,39 +220,39 @@ module MOM_variables
type, public :: vertvisc_type
real :: Prandtl_turb !< The Prandtl number for the turbulent diffusion
!! that is captured in Kd_shear [nondim].
real, pointer, dimension(:,:) :: &
bbl_thick_u => NULL(), & !< The bottom boundary layer thickness at the u-points [Z ~> m].
bbl_thick_v => NULL(), & !< The bottom boundary layer thickness at the v-points [Z ~> m].
kv_bbl_u => NULL(), & !< The bottom boundary layer viscosity at the u-points [Z2 T-1 ~> m2 s-1].
kv_bbl_v => NULL(), & !< The bottom boundary layer viscosity at the v-points [Z2 T-1 ~> m2 s-1].
ustar_BBL => NULL() !< The turbulence velocity in the bottom boundary layer at h points [Z T-1 ~> m s-1].
real, pointer, dimension(:,:) :: TKE_BBL => NULL()
!< A term related to the bottom boundary layer source of turbulent kinetic
!! energy, currently in [Z3 T-3 ~> m3 s-3], but may at some time be changed
!! to [R Z3 T-3 ~> W m-2].
real, pointer, dimension(:,:) :: &
taux_shelf => NULL(), & !< The zonal stresses on the ocean under shelves [R Z L T-2 ~> Pa].
tauy_shelf => NULL() !< The meridional stresses on the ocean under shelves [R Z L T-2 ~> Pa].
real, pointer, dimension(:,:) :: tbl_thick_shelf_u => NULL()
real, allocatable, dimension(:,:) :: &
bbl_thick_u, & !< The bottom boundary layer thickness at the u-points [Z ~> m].
bbl_thick_v, & !< The bottom boundary layer thickness at the v-points [Z ~> m].
kv_bbl_u, & !< The bottom boundary layer viscosity at the u-points [Z2 T-1 ~> m2 s-1].
kv_bbl_v, & !< The bottom boundary layer viscosity at the v-points [Z2 T-1 ~> m2 s-1].
ustar_BBL, & !< The turbulence velocity in the bottom boundary layer at h points [Z T-1 ~> m s-1].
TKE_BBL, & !< A term related to the bottom boundary layer source of turbulent kinetic
!! energy, currently in [Z3 T-3 ~> m3 s-3], but may at some time be changed
!! to [R Z3 T-3 ~> W m-2].
taux_shelf, & !< The zonal stresses on the ocean under shelves [R Z L T-2 ~> Pa].
tauy_shelf !< The meridional stresses on the ocean under shelves [R Z L T-2 ~> Pa].
real, allocatable, dimension(:,:) :: tbl_thick_shelf_u
!< Thickness of the viscous top boundary layer under ice shelves at u-points [Z ~> m].
real, pointer, dimension(:,:) :: tbl_thick_shelf_v => NULL()
real, allocatable, dimension(:,:) :: tbl_thick_shelf_v
!< Thickness of the viscous top boundary layer under ice shelves at v-points [Z ~> m].
real, pointer, dimension(:,:) :: kv_tbl_shelf_u => NULL()
real, allocatable, dimension(:,:) :: kv_tbl_shelf_u
!< Viscosity in the viscous top boundary layer under ice shelves at u-points [Z2 T-1 ~> m2 s-1].
real, pointer, dimension(:,:) :: kv_tbl_shelf_v => NULL()
real, allocatable, dimension(:,:) :: kv_tbl_shelf_v
!< Viscosity in the viscous top boundary layer under ice shelves at v-points [Z2 T-1 ~> m2 s-1].
real, pointer, dimension(:,:) :: nkml_visc_u => NULL()
real, allocatable, dimension(:,:) :: nkml_visc_u
!< The number of layers in the viscous surface mixed layer at u-points [nondim].
!! This is not an integer because there may be fractional layers, and it is stored in
!! terms of layers, not depth, to facilitate the movement of the viscous boundary layer
!! with the flow.
real, pointer, dimension(:,:) :: nkml_visc_v => NULL()
real, allocatable, dimension(:,:) :: nkml_visc_v
!< The number of layers in the viscous surface mixed layer at v-points [nondim].
real, allocatable, dimension(:,:,:) :: &
Ray_u, & !< The Rayleigh drag velocity to be applied to each layer at u-points [Z T-1 ~> m s-1].
Ray_v !< The Rayleigh drag velocity to be applied to each layer at v-points [Z T-1 ~> m s-1].

! The following elements are pointers so they can be used as targets for pointers in the restart registry.
real, pointer, dimension(:,:) :: &
MLD => NULL() !< Instantaneous active mixing layer depth [Z ~> m].
real, pointer, dimension(:,:,:) :: &
Ray_u => NULL(), & !< The Rayleigh drag velocity to be applied to each layer at u-points [Z T-1 ~> m s-1].
Ray_v => NULL() !< The Rayleigh drag velocity to be applied to each layer at v-points [Z T-1 ~> m s-1].
MLD => NULL() !< Instantaneous active mixing layer depth [Z ~> m].
real, pointer, dimension(:,:,:) :: Kd_shear => NULL()
!< The shear-driven turbulent diapycnal diffusivity at the interfaces between layers
!! in tracer columns [Z2 T-1 ~> m2 s-1].
Expand Down
52 changes: 33 additions & 19 deletions src/diagnostics/MOM_PointAccel.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,8 @@ module MOM_PointAccel
!! written by this PE during the current run.
integer :: max_writes !< The maximum number of times any PE can write out
!! a column's worth of accelerations during a run.
logical :: full_column !< If true, write out the accelerations in all massive layers,
!! otherwise just document the ones with large velocities.
type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock.
type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
!! regulate the timing of diagnostic output.
Expand Down Expand Up @@ -80,11 +82,12 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
!! call to PointAccel_init.
real, intent(in) :: vel_rpt !< The velocity magnitude that triggers a report [L T-1 ~> m s-1]
real, optional, intent(in) :: str !< The surface wind stress [R L Z T-2 ~> Pa]
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), &
optional, intent(in) :: a !< The layer coupling coefficients from vertvisc [Z T-1 ~> m s-1].
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
optional, intent(in) :: hv !< The layer thicknesses at velocity grid points,
!! from vertvisc [H ~> m or kg m-2].

! Local variables
real :: CFL ! The local velocity-based CFL number [nondim]
real :: Angstrom ! A negligibly small thickness [H ~> m or kg m-2]
Expand Down Expand Up @@ -145,6 +148,9 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
if (ke < ks) then
ks = 1; ke = nz; write(file,'("U: Unable to set ks & ke.")')
endif
if (CS%full_column) then
ks = 1 ; ke = nz
endif

call get_date(CS%Time, yr, mo, day, hr, minute, sec)
call get_time((CS%Time - set_date(yr, 1, 1, 0, 0, 0)), sec, yearday)
Expand Down Expand Up @@ -217,21 +223,23 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
(vel_scale*ADp%du_other(I,j,k)) ; enddo
endif
if (present(a)) then
write(file,'(/,"a: ")', advance='no')
do k=ks,ke+1 ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (US%Z_to_m*a(I,j,k)*dt) ; enddo
write(file,'(/,"a: ",ES10.3," ")', advance='no') US%Z_to_m*a(I,j,ks)*dt
do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') (US%Z_to_m*a(I,j,K)*dt) ; enddo
endif
if (present(hv)) then
write(file,'(/,"hvel: ")', advance='no')
do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hv(I,j,k) ; enddo
endif
write(file,'(/,"Stress: ",ES10.3)') vel_scale*US%Z_to_m * (str*dt / GV%Rho0)
if (present(str)) then
write(file,'(/,"Stress: ",ES10.3)', advance='no') vel_scale*US%Z_to_m * (str*dt / GV%Rho0)
endif

if (associated(CS%u_accel_bt)) then
write(file,'("dubt: ")', advance='no')
write(file,'(/,"dubt: ")', advance='no')
do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') &
(vel_scale*dt*CS%u_accel_bt(I,j,k)) ; enddo
write(file,'(/)')
endif
write(file,'(/)')

write(file,'(/,"h--: ")', advance='no')
do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (h_scale*hin(i,j-1,k)) ; enddo
Expand All @@ -249,14 +257,12 @@ subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st

e(nz+1) = -US%Z_to_m*(G%bathyT(i,j) + G%Z_ref)
do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i,j,k) ; enddo
write(file,'(/,"e-: ")', advance='no')
write(file,'(ES10.3," ")', advance='no') e(ks)
write(file,'(/,"e-: ",ES10.3," ")', advance='no') e(ks)
do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo

e(nz+1) = -US%Z_to_m*(G%bathyT(i+1,j) + G%Z_ref)
do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i+1,j,k) ; enddo
write(file,'(/,"e+: ")', advance='no')
write(file,'(ES10.3," ")', advance='no') e(ks)
write(file,'(/,"e+: ",ES10.3," ")', advance='no') e(ks)
do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo
if (associated(CS%T)) then
write(file,'(/,"T-: ")', advance='no')
Expand Down Expand Up @@ -415,11 +421,12 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
!! call to PointAccel_init.
real, intent(in) :: vel_rpt !< The velocity magnitude that triggers a report [L T-1 ~> m s-1]
real, optional, intent(in) :: str !< The surface wind stress [R L Z T-2 ~> Pa]
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), &
optional, intent(in) :: a !< The layer coupling coefficients from vertvisc [Z T-1 ~> m s-1].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
optional, intent(in) :: hv !< The layer thicknesses at velocity grid points,
!! from vertvisc [H ~> m or kg m-2].

! Local variables
real :: CFL ! The local velocity-based CFL number [nondim]
real :: Angstrom ! A negligibly small thickness [H ~> m or kg m-2]
Expand Down Expand Up @@ -479,6 +486,9 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
if (ke < ks) then
ks = 1; ke = nz; write(file,'("V: Unable to set ks & ke.")')
endif
if (CS%full_column) then
ks = 1 ; ke = nz
endif

call get_date(CS%Time, yr, mo, day, hr, minute, sec)
call get_time((CS%Time - set_date(yr, 1, 1, 0, 0, 0)), sec, yearday)
Expand Down Expand Up @@ -556,21 +566,23 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st
(vel_scale*ADp%dv_other(i,J,k)) ; enddo
endif
if (present(a)) then
write(file,'(/,"a: ")', advance='no')
do k=ks,ke+1 ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') (US%Z_to_m*a(i,j,k)*dt) ; enddo
write(file,'(/,"a: ",ES10.3," ")', advance='no') US%Z_to_m*a(i,J,ks)*dt
do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') (US%Z_to_m*a(i,J,K)*dt) ; enddo
endif
if (present(hv)) then
write(file,'(/,"hvel: ")', advance='no')
do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hv(i,J,k) ; enddo
endif
write(file,'(/,"Stress: ",ES10.3)') vel_scale*US%Z_to_m * (str*dt / GV%Rho0)
if (present(str)) then
write(file,'(/,"Stress: ",ES10.3)', advance='no') vel_scale*US%Z_to_m * (str*dt / GV%Rho0)
endif

if (associated(CS%v_accel_bt)) then
write(file,'("dvbt: ")', advance='no')
do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') &
(vel_scale*dt*CS%v_accel_bt(i,J,k)) ; enddo
write(file,'(/)')
endif
write(file,'(/)')

write(file,'("h--: ")', advance='no')
do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ")', advance='no') h_scale*hin(i-1,j,k) ; enddo
Expand All @@ -587,14 +599,12 @@ subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt, G, GV, US, CS, vel_rpt, st

e(nz+1) = -US%Z_to_m*(G%bathyT(i,j) + G%Z_ref)
do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i,j,k) ; enddo
write(file,'(/,"e-: ")', advance='no')
write(file,'(ES10.3," ")', advance='no') e(ks)
write(file,'(/,"e-: ",ES10.3," ")', advance='no') e(ks)
do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo

e(nz+1) = -US%Z_to_m*(G%bathyT(i,j+1) + G%Z_ref)
do k=nz,1,-1 ; e(K) = e(K+1) + h_scale*hin(i,j+1,k) ; enddo
write(file,'(/,"e+: ")', advance='no')
write(file,'(ES10.3," ")', advance='no') e(ks)
write(file,'(/,"e+: ",ES10.3," ")', advance='no') e(ks)
do K=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ")', advance='no') e(K) ; enddo
if (associated(CS%T)) then
write(file,'(/,"T-: ")', advance='no')
Expand Down Expand Up @@ -773,6 +783,10 @@ subroutine PointAccel_init(MIS, Time, G, param_file, diag, dirs, CS)
call get_param(param_file, mdl, "MAX_TRUNC_FILE_SIZE_PER_PE", CS%max_writes, &
"The maximum number of columns of truncations that any PE "//&
"will write out during a run.", default=50, debuggingParam=.true.)
call get_param(param_file, mdl, "DEBUG_FULL_COLUMN", CS%full_column, &
"If true, write out the accelerations in all massive layers; otherwise "//&
"just document the ones with large velocities.", &
default=.false., debuggingParam=.true.)

if (len_trim(dirs%output_directory) > 0) then
if (len_trim(CS%u_trunc_file) > 0) &
Expand Down
2 changes: 1 addition & 1 deletion src/parameterizations/lateral/MOM_MEKE.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,7 +248,7 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h


! Calculate drag_rate_visc(i,j) which accounts for the model bottom mean flow
if (CS%visc_drag) then
if (CS%visc_drag .and. allocated(visc%Kv_bbl_u) .and. allocated(visc%Kv_bbl_v)) then
!$OMP parallel do default(shared)
do j=js,je ; do I=is-1,ie
drag_vel_u(I,j) = 0.0
Expand Down
Loading

0 comments on commit 5cadb72

Please sign in to comment.