module ocean_bbc_mod
!
! Matthew Harrison
!
!
! S.M. Griffies
!
!
! Martin Schmidt
!
!
! Hyun-Chul Lee
!
!
!
! Set bottom boundary conditions
!
!
!
! Set bottom boundary conditions
!
!
!
!
!
! R.C. Pacanowski and S.M. Griffies
! The MOM3 Manual (1999)
!
!
!
! S.M. Griffies, M.J. Harrison, R.C. Pacanowski, and A. Rosati
! A Technical Guide to MOM4 (2003)
!
!
!
! S.M. Griffies, 2012: Elements of MOM
!
!
!
!
!
!
!
! Dimensionless coefficient for quadratic bottom drag.
!
!
!
! For incorporating the bottom momentum drag implicitly in time.
! Default is bmf_implicit=.false.
!
!
!
! For determining bottom drag coefficient using a constant roughness length.
! Will take maximum between cdbot and the computed value using law of
! wall log-profile. This option is of use when we have very very
! refined vertical resolution (say on order of meters) near the bottom.
! Terrain following coordinates should use this option since they generally
! have very refined vertical grid spacing on topography.
! Default is cdbot_law_of_wall=.false.
!
!
! Bottom roughness length. Default is law_of_wall_rough_length=0.01m, following
! the default used in the Princeton Ocean Model (POM). This value
! corresponds to "Law of Wall" physics.
!
!
!
! For determining bottom drag coefficient using a map of the roughness length.
! This approach is more relevant for coarse models
! than the constant roughness length used in the cdbot_law_of_wall option.
! cdbot_lo <= cdbot(i,j) <= cdbot_hi.
! Default is cdbot_roughness_length=.false.
!
!
!
! For determining bottom drag coefficient using a map of the roughness length
! and tidal velocity amplitude.
! This approach is more relevant for coarse models
! than the constant roughness length used in the cdbot_law_of_wall option.
! cdbot_lo <= cdbot(i,j) <= cdbot_hi.
! Default is cdbot_roughness_uamp=.false.
!
!
!
! H0 in a parameterization of cdbot_roughness_uamp.
! Default is cdbot_HH=1100.0.
!
!
!
! U0 in a parameterization of cdbot_roughness_uamp.
! Default is cdbot_UU=1.0.
!
!
!
! For determining bottom drag coefficient using a map of the roughness length
! and the surface wind wave field. The modified drag coefficient is calculated
! following Grant and Mattsen. Likewise this method can be improved using
! more sophisticated wave models including swell.
! Default is cdbot_wave=.false.
!
!
!
! Residual bottom velocity due to unresolved fluctuations (e.g., waves and tides)
! that contribute to bottom dissipation. Should be set to zero when running
! with explicit representation of tidal forcing and when waves are well resolved.
! Default is uresidual=.05.
!
!
!
! Maximum magnitude of the bottom velocity used to compute the bottom
! momentum drag. Default is uvmag_max=10.0.
!
!
!
! Maximum magnitude of the bottom momentum drag.
! Default is bmf_max=1.0.
!
!
!
! For debugging purposes.
!
!
!
use constants_mod, only: epsln, pi
use diag_manager_mod, only: register_diag_field, register_static_field, send_data
use fms_mod, only: open_namelist_file, check_nml_error, write_version_number
use fms_mod, only: close_file, read_data
use mpp_mod, only: input_nml_file, mpp_error, NOTE, FATAL, WARNING, stdout, stdlog
use mpp_domains_mod, only: mpp_update_domains, BGRID_NE
use ocean_domains_mod, only: get_local_indices
use ocean_parameters_mod, only: missing_value, TERRAIN_FOLLOWING
use ocean_parameters_mod, only: onefourth, rho0, rho0r, cp_ocean, von_karman
use ocean_parameters_mod, only: MOM_BGRID, MOM_CGRID
use ocean_types_mod, only: ocean_velocity_type, ocean_domain_type
use ocean_types_mod, only: ocean_grid_type, ocean_prog_tracer_type
use ocean_types_mod, only: ocean_time_type, ocean_thickness_type
use ocean_types_mod, only: ocean_density_type, ocean_options_type
use ocean_types_mod, only: ocean_external_mode_type
use ocean_workspace_mod, only: wrk1_2d, wrk2_2d
use wave_types_mod, only: ocean_wave_type
use ocean_wave_mod, only: wave_model_is_initialised
use ocean_util_mod, only: diagnose_2d, diagnose_2d_u, diagnose_sum
implicit none
private
#include
integer :: num_prog_tracers
! for Bgrid or Cgrid
integer :: horz_grid
! diagnostics
integer :: id_eta_tend_geoheat =-1
integer :: id_eta_tend_geoheat_glob =-1
integer :: id_geo_heat =-1
integer :: id_gamma_bmf =-1
integer :: id_drag_coeff =-1
integer :: id_cdbot =-1
integer :: id_bmf_u =-1
integer :: id_bmf_v =-1
integer :: id_cur_wav_dr =-1
integer :: id_wave_s =-1
integer :: id_wave_u =-1
integer :: id_iter =-1
logical :: used
type(ocean_grid_type), pointer :: Grd =>NULL()
type(ocean_domain_type), pointer :: Dom =>NULL()
character(len=128) :: version=&
'$Id: ocean_bbc.F90,v 20.0 2013/12/14 00:10:36 fms Exp $'
character (len=128) :: tagname = &
'$Name: tikal $'
public :: ocean_bbc_init
public :: get_ocean_bbc
private :: current_wave_drag_diag
private :: wave_u_diag
logical :: module_is_initialized = .false.
logical :: first_call = .true.
real :: vonkarman2
real :: law_of_wall_rough_length_r
real :: cellarea_r
real :: cp_r
real :: uresidual2
real :: small=1.e-20
integer :: index_temp=-1
integer :: index_salt=-1
integer :: index_redist_heat=-1
integer :: bbc_geothermal
integer :: stdoutunit,stdlogunit
real, allocatable, dimension(:,:) :: cdbot_array ! dimensionless drag coefficient
real, allocatable, dimension(:,:) :: cdbot_lowall ! dimensionless drag coefficient from law of wall
real, allocatable, dimension(:,:) :: drag_coeff ! rho0 * dimensionless drag coefficient
real, allocatable, dimension(:,:) :: roughness_length ! for computing bottom drag coefficient
real, allocatable, dimension(:,:) :: tidal_uamp ! for computing bottom drag coefficient
real, allocatable, dimension(:,:) :: geo_heat ! geothermal heating
real, allocatable, dimension(:,:) :: wave_u ! wave action at the sea floor, wave boundary layer
real, allocatable, dimension(:,:) :: wave_s ! wave action at the sea floor, skin layer
real, allocatable, dimension(:,:) :: data
integer, allocatable, dimension(:,:) :: grd_kbot ! bottom k-level
real, allocatable, dimension(:, :) :: uresidual2_2d
real, allocatable, dimension(:, :) :: tide_speed_t
! nml variables
logical :: bmf_implicit = .false.
logical :: debug_this_module = .false.
logical :: cdbot_law_of_wall = .false.
logical :: cdbot_roughness_length = .false.
logical :: cdbot_wave = .false.
logical :: cdbot_roughness_uamp = .false.
logical :: use_geothermal_heating = .false.
real :: convert_geothermal = .001 ! for converting units from mW to W
real :: law_of_wall_rough_length = 0.01 ! metre
real :: cdbot = 2.5e-3 ! dimensionless
real :: uresidual = .05 ! m/sec
real :: cdbot_hi = 3.0e-3 ! hi-end of cdbot for cdbot_roughness_length and cdbot_roughness_uamp
real :: cdbot_lo = 1.0e-3 ! lo-end of cdbot for cdbot_roughness_length and cdbot_roughness_uamp
real :: cdbot_gamma = 40.0 ! for setting exp decay for cdbot w/ cdbot_roughness_length
real :: uvmag_max = 10.0 ! max velocity scale (m/s) for computing bmf
real :: bmf_max = 1.0 ! max bmf (N/m^2)
real :: cdbot_HH = 1100.0 ! H0 of cdbot_roughness_uamp (m)
real :: cdbot_UU = 1.0 ! U0 of cdbot_roughness_uamp (m/s)
logical :: read_tide_speed = .false.
real :: uresidual2_max = 0.05
namelist /ocean_bbc_nml/ bmf_implicit, cdbot, uresidual, cdbot_law_of_wall, law_of_wall_rough_length, &
cdbot_roughness_length, use_geothermal_heating, convert_geothermal, &
cdbot_hi, cdbot_lo, cdbot_gamma, uvmag_max, bmf_max, debug_this_module, &
cdbot_roughness_uamp, cdbot_HH, cdbot_UU, cdbot_wave
namelist /ocean_bbc_ofam_nml/ read_tide_speed, uresidual2_max
contains
!#######################################################################
!
!
!
! Initialize the bottom boundary condition module
!
!
subroutine ocean_bbc_init(Grid, Domain, Time, T_prog, Velocity, &
Ocean_options, vert_coordinate_type, hor_grid)
type(ocean_grid_type), target, intent(in) :: Grid
type(ocean_domain_type), target, intent(in) :: Domain
type(ocean_time_type), intent(in) :: Time
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
type(ocean_velocity_type), intent(inout) :: Velocity
type(ocean_options_type), intent(inout) :: Ocean_options
integer, intent(in) :: vert_coordinate_type
integer, intent(in) :: hor_grid
logical :: read_roughness=.false.
integer :: i,j,n, ioun, io_status, ierr
integer :: grd_axes(3)
real :: depth, cdbot_term1, cdbot_term2
stdoutunit=stdout()
stdlogunit=stdlog()
module_is_initialized = .TRUE.
call write_version_number(version, tagname)
#ifdef INTERNAL_FILE_NML
read (input_nml_file, nml=ocean_bbc_nml, iostat=io_status)
ierr = check_nml_error(io_status,'ocean_bbc_nml')
read (input_nml_file, nml=ocean_bbc_ofam_nml, iostat=io_status)
ierr = check_nml_error(io_status,'ocean_bbc_ofam_nml')
#else
ioun = open_namelist_file()
read(ioun, ocean_bbc_nml, iostat=io_status)
ierr = check_nml_error(io_status,'ocean_bbc_nml')
rewind(ioun)
read(ioun, ocean_bbc_ofam_nml, iostat=io_status)
ierr = check_nml_error(io_status,'ocean_bbc_ofam_nml')
call close_file(ioun)
#endif
write (stdoutunit,'(/)')
write (stdoutunit, ocean_bbc_nml)
write (stdlogunit, ocean_bbc_nml)
write (stdoutunit,'(/)')
write (stdoutunit, ocean_bbc_ofam_nml)
write (stdlogunit, ocean_bbc_ofam_nml)
#ifndef MOM_STATIC_ARRAYS
call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec)
nk = Grid%nk
#endif
Dom => Domain
Grd => Grid
horz_grid = hor_grid
num_prog_tracers = size(T_prog(:))
index_temp=-1;index_salt=-1;index_redist_heat=-1
do n=1,num_prog_tracers
if (T_prog(n)%name == 'temp') index_temp = n
if (T_prog(n)%name == 'salt') index_salt = n
if (T_prog(n)%name == 'redist_heat') index_redist_heat = n
enddo
do n=1, num_prog_tracers
#ifndef MOM_STATIC_ARRAYS
allocate(T_prog(n)%btf(isd:ied,jsd:jed))
#endif
T_prog(n)%btf = 0.0
enddo
allocate (cdbot_array(isd:ied,jsd:jed))
allocate (cdbot_lowall(isd:ied,jsd:jed))
allocate (drag_coeff(isd:ied,jsd:jed))
allocate (roughness_length(isd:ied,jsd:jed))
allocate (tidal_uamp(isd:ied,jsd:jed))
allocate (geo_heat(isd:ied,jsd:jed))
allocate (data(isd:ied,jsd:jed))
allocate (wave_u(isd:ied,jsd:jed))
allocate (wave_s(isd:ied,jsd:jed))
allocate (grd_kbot(isd:ied,jsd:jed))
vonkarman2 = von_karman*von_karman
geo_heat(:,:) = 0.0
data(:,:) = 0.0
roughness_length(:,:) = 0.0
tidal_uamp(:,:) = 0.0
wave_u(:,:) = 0.0
wave_s(:,:) = 0.0
if(horz_grid == MOM_BGRID) then
grd_kbot(:,:) = Grd%kmu(:,:)
grd_axes(:) = Grd%vel_axes_uv(:)
else
grd_kbot(:,:) = Grd%kmt(:,:)
grd_axes(:) = Grd%tracer_axes(:)
endif
if(cdbot_wave .or. cdbot_roughness_length .or. cdbot_roughness_uamp) then
read_roughness = .true.
endif
drag_coeff(:,:) = rho0*cdbot*Grd%mask(:,:,1)
cdbot_array(:,:) = cdbot*Grd%mask(:,:,1)
cdbot_lowall(:,:) = cdbot*Grd%mask(:,:,1)
cellarea_r = 1.0/(epsln + Grd%tcellsurf)
cp_r = 1.0/cp_ocean
uresidual2 = uresidual**2
if (read_tide_speed) then
allocate(uresidual2_2d(isd:ied, jsd:jed))
allocate(tide_speed_t(isd:ied, jsd:jed))
call read_data('INPUT/tideamp.nc', 'tideamp', tide_speed_t, Domain%domain2d)
write (stdout(),*) '==>ocean_bbc_mod: Completed read of tide_speed.'
uresidual2_2d(:,:) = min(tide_speed_t(:,:)**2, uresidual2_max)
call mpp_update_domains(uresidual2_2d(:,:), Dom%domain2d)
deallocate(tide_speed_t)
else
write(stdout(),'(a)') &
'==>Note: NOT reading tide_speed for ocean_vert_tidal_mod.'
call mpp_error(NOTE, &
'==>ocean_vert_tidal_mod: Setting tide_speed to default value.')
endif
law_of_wall_rough_length_r = 1.0/law_of_wall_rough_length
if(cdbot_law_of_wall) then
write(stdoutunit,'(a)') &
'==>Note from ocean_bbc_mod: Computing bottom drag coefficient from roughness length.'
Ocean_options%bottom_roughness = 'Used bottom drag coefficient from Law of Wall roughness length.'
elseif (cdbot_roughness_length) then
Ocean_options%bottom_roughness = 'Used bottom drag from map of bottom roughness using law of wall.'
elseif (cdbot_roughness_uamp) then
Ocean_options%bottom_roughness = 'Used bottom drag from map of bottom roughness w/ tidal amplitude from law of wall.'
elseif (cdbot_wave) then
Ocean_options%bottom_roughness = 'Used bottom drag from map of bottom roughness using law of wall.'&
//' Also considered influence of the wind wave field.'
else
Ocean_options%bottom_roughness = 'Used bottom drag from specified drag coefficient.'
endif
if(bmf_implicit) then
write(stdoutunit,'(a)') &
'==>Note from ocean_bbc_mod: Computing bottom drag implicitly in time.'
Ocean_options%bmf_implicit = 'Used bottom drag computed implicitly in time.'
Velocity%bmf_implicit = .true.
else
Ocean_options%bmf_implicit = 'Used bottom drag computed explicitly in time.'
Velocity%bmf_implicit = .false.
endif
! compute drag coefficient using law of wall with a roughness length
! read in from a map. This approach is more relevant for large-scale
! coarse models than the constant roughness length used in the
! cdbot_law_of_wall option.
! read in topographic roughness_length.
! assume roughness_length on B-grid velocity point if using MOM_BGRID.
! assume roughness_length on tracer point if using MOM_CGRID.
data = 0.0
if(cdbot_wave) then
data = law_of_wall_rough_length
endif
if(read_roughness) then
call read_data('INPUT/roughness_cdbot.nc','roughness_cdbot', data, Domain%domain2d)
if(horz_grid == MOM_BGRID) then
write (stdoutunit,*) '==>ocean_bbc_mod: read in topographic roughness length, assumed to be on B-grid velocity point.'
else
write (stdoutunit,*) '==>ocean_bbc_mod: read in topographic roughness length, assumed to be on tracer point.'
endif
do j=jsc,jec
do i=isc,iec
roughness_length(i,j) = Grd%mask(i,j,1)*max(0.0,data(i,j))
enddo
enddo
call mpp_update_domains(roughness_length(:,:), Dom%domain2d)
endif
! bottom drag from law of wall using roughness_length
if(cdbot_roughness_length) then
cdbot_lowall=0.0
do j=jsd,jed
do i=isd,ied
depth = Grd%ht(i,j)
if(depth > 0.0) then
cdbot_lowall(i,j)= (von_karman/( alog(depth/(roughness_length(i,j)+small)) + small ) )**2
endif
enddo
enddo
! bottom drag coefficient; the range of cdbot is
! cdbot_lo <= cdbot(i,j) <= cdbot_hi
do j=jsd,jed
do i=isd,ied
cdbot_array(i,j) = Grd%mask(i,j,1)*cdbot_hi*(1.0-exp(-cdbot_gamma*cdbot_lowall(i,j)))
cdbot_array(i,j) = max(cdbot_lo,cdbot_array(i,j))
enddo
enddo
endif
! scheme coded by H.C. Lee (May 2012)
if(cdbot_roughness_uamp) then
! read in tidal velocity amplitude; assumed on U-grid
data=0.0
call read_data('INPUT/tideamp.nc','TIDEAMP', data, Domain%domain2d)
if(horz_grid == MOM_BGRID) then
write (stdoutunit,*) '==>ocean_bbc_mod: completed read of tidal velocity amplitude, assumed to be on B-grid velocity point.'
else
write (stdoutunit,*) '==>ocean_bbc_mod: completed read of tidal velocity amplitude, assumed to be on tracer point.'
endif
do j=jsc,jec
do i=isc,iec
tidal_uamp(i,j) = Grd%mask(i,j,1)*max(0.0,data(i,j))
enddo
enddo
call mpp_update_domains(tidal_uamp(:,:), Dom%domain2d)
! bottom drag from law of wall using roughness_length and tidal velocity amplitude
do j=jsd,jed
do i=isd,ied
depth = Grd%ht(i,j)
cdbot_term1= cdbot_HH/(roughness_length(i,j)+small)
cdbot_term2= cdbot_UU/(tidal_uamp(i,j)+small)
if(depth > 0.0) then
cdbot_array(i,j)= (von_karman/( alog(cdbot_term1*cdbot_term2) + small ) )**2
endif
enddo
enddo
! bottom drag coefficient; the range of cdbot is
! cdbot_lo <= cdbot(i,j) <= cdbot_hi
do j=jsd,jed
do i=isd,ied
cdbot_array(i,j) = min(cdbot_hi,cdbot_array(i,j))
cdbot_array(i,j) = max(cdbot_lo,cdbot_array(i,j))
enddo
enddo
endif
! save cdbot into Velocity type
do j=jsd,jed
do i=isd,ied
Velocity%cdbot_array(i,j) = cdbot_array(i,j)
enddo
enddo
if(use_geothermal_heating) then
write(stdoutunit,'(a)') &
'==>Note from ocean_bbc_mod: Geothermal heating introduced at ocean bottom.'
Ocean_options%geothermal_heating = 'Used geothermal heating introduced at ocean bottom.'
! fill geothermal heating array
geo_heat = 0.0
call read_data('INPUT/geothermal_heating.nc', 'geo_heat', data(:,:), Dom%domain2d, timelevel=1)
! make sure there are no spurious negative values.
! also convert to units appropriate for temperature btf.
do j=jsc,jec
do i=isc,iec
geo_heat(i,j) = Grd%tmask(i,j,1)*max(0.0,data(i,j))*cp_r*convert_geothermal
enddo
enddo
else
Ocean_options%geothermal_heating = 'Did NOT introduce geothermal heating.'
endif
if(vert_coordinate_type == TERRAIN_FOLLOWING .and. .not. cdbot_law_of_wall) then
call mpp_error(WARNING,&
'==>Warning from ocean_bbc_init: recommend law_of_wall=.true. for TERRAIN_FOLLOWING coordinates.')
endif
id_drag_coeff = register_diag_field ('ocean_model', 'drag_coeff', &
grd_axes(1:2), Time%model_time, &
'Dimensionless bottom drag coefficient', 'dimensionless', &
missing_value=missing_value, range=(/-1.0,1.e3/))
id_gamma_bmf = register_diag_field ('ocean_model', 'gamma_bmf', &
grd_axes(1:2), Time%model_time, &
'Bottom drag factor rho0*cdbot*uvmag', 'kg/(m^2 sec)', &
missing_value=missing_value, range=(/-1.0,1.e12/))
id_bmf_u = register_diag_field ('ocean_model', 'bmf_u', &
grd_axes(1:2), Time%model_time, &
'Bottom u-stress via bottom drag', 'N/m^2', &
missing_value=missing_value, range=(/-1.0,1.e2/))
id_bmf_v = register_diag_field ('ocean_model', 'bmf_v', &
grd_axes(1:2), Time%model_time, &
'Bottom v-stress via bottom drag', 'N/m^2', &
missing_value=missing_value, range=(/-1.0,1.e2/))
id_cdbot = register_static_field ('ocean_model', 'cdbot', &
grd_axes(1:2), 'Bottom drag coefficient', 'dimensionless',&
missing_value=missing_value, range=(/-1.0,1e6/))
if (id_cdbot > 0) then
used = send_data (id_cdbot, cdbot_array(:,:), Time%model_time, &
rmask=Grd%mask(:,:,1), is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
id_geo_heat = register_static_field ('ocean_model', 'geo_heat', &
Grd%tracer_axes(1:2), 'Geothermal heating', 'W/m^2',&
missing_value=missing_value, range=(/-10.0,1e6/), &
standard_name='upward_geothermal_heat_flux_at_sea_floor')
call diagnose_2d(Time, Grd, id_geo_heat, geo_heat(:,:)*cp_ocean)
id_eta_tend_geoheat = register_diag_field('ocean_model', &
'eta_tend_geoheat', Grd%tracer_axes(1:2),Time%model_time , &
'non-Bouss steric sea level tendency from geothermal heating',&
'm/s',missing_value=missing_value,range=(/-1.e10,1.e10/))
id_eta_tend_geoheat_glob = register_diag_field('ocean_model', &
'eta_tend_geoheat_glob', Time%model_time, &
'non-Bouss steric sea level tendency from geothermal heating',&
'm/s',missing_value=missing_value,range=(/-1.e10,1.e10/))
id_cur_wav_dr = register_diag_field ('ocean_model', 'current_wave_stress',&
Grd%tracer_axes(1:2), Time%model_time, &
'combined current wave bottom stress', 'N/m^2', &
missing_value=missing_value, range=(/-1.0,1.e3/))
id_wave_s = register_diag_field ('ocean_model', 'wave_s', &
Grd%tracer_axes(1:2), Time%model_time, &
'wave skin friction velocity', 'm/s', &
missing_value=missing_value, range=(/-1.0,1.e3/))
id_wave_u = register_diag_field ('ocean_model', 'wave_u',&
Grd%tracer_axes(1:2), Time%model_time, &
'wave friction velocity', 'm/s', &
missing_value=missing_value, range=(/-1.0,1.e3/))
id_iter = register_diag_field ('ocean_model', 'iter',&
grd_axes(1:2), Time%model_time, &
'number of ustar iterations', 'm/s', &
missing_value=missing_value, range=(/0.0,1.e3/))
end subroutine ocean_bbc_init
! NAME="ocean_bbc_init"
!#######################################################################
!
!
!
! Set bottom boundary conditions for velocity and tracer.
!
! Dimensions of bottom momentum flux are
! N/m^2 = (kg/m^3)*(m^2/s^2).
!
! Note the use of rho0 for the conversion from m^2/s^2 to
! (kg/m^3)*(m^2/s^2). We do not know the precise value
! of cdbot, so the rho0 approximate value is well within
! our level of uncertainty. No reason therefore to
! use in situ rho for this conversion, even when using
! non-Boussinesq version of MOM.
!
!
! Note that bmf needs to be computed on the data domain since the
! halo values are accessed in ocean_vert_gotm.F90.
!
!
!
subroutine get_ocean_bbc(Time, Thickness, Dens, Velocity, T_prog, Waves)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_density_type), intent(in) :: Dens
type(ocean_velocity_type), intent(inout) :: Velocity
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
type(ocean_wave_type), intent(inout) :: Waves
integer :: i, j, kbot, n
integer :: taum1, tau, tstep
real :: uvmag, argument, distance
real :: umag, vmag
real :: rhobot_inv, alphabot
if (.not. module_is_initialized) then
call mpp_error(FATAL,'==>Error from ocean_bbc_mod (get_ocean_bbc): module must be initialized')
endif
taum1 = Time%taum1
tau = Time%tau
if(bmf_implicit) then
tstep=tau
else
tstep=taum1
endif
! Bottom tracer flux defaulted to zero.
do n=1,num_prog_tracers
do j=jsd,jed
do i=isd,ied
T_prog(n)%btf(i,j) = 0.0
enddo
enddo
enddo
! T_prog(index_temp)%btf<0 means heat enters ocean; hence the minus sign.
if (use_geothermal_heating) then
do j=jsc,jec
do i=isc,iec
T_prog(index_temp)%btf(i,j) = -geo_heat(i,j)
enddo
enddo
if(index_redist_heat > 0) then
do j=jsc,jec
do i=isc,iec
T_prog(index_redist_heat)%btf(i,j) = -geo_heat(i,j)
enddo
enddo
endif
endif
! set bottom drag coefficient
drag_coeff(:,:) = 0.0
if(horz_grid == MOM_BGRID) then
if(cdbot_law_of_wall) then
do j=jsd,jed
do i=isd,ied
kbot = grd_kbot(i,j)
if(kbot > 0) then
distance = Thickness%depth_zwu(i,j,kbot) - Thickness%depth_zu(i,j,kbot)
argument = law_of_wall_rough_length_r*(distance + law_of_wall_rough_length)
drag_coeff(i,j) = rho0*max(cdbot,vonkarman2/(log(argument))**2)
endif
enddo
enddo
else
do j=jsd,jed
do i=isd,ied
drag_coeff(i,j) = rho0*cdbot_array(i,j)*Grd%umask(i,j,1)
enddo
enddo
endif
else
if(cdbot_law_of_wall) then
do j=jsd,jed
do i=isd,ied
kbot = grd_kbot(i,j)
if(kbot > 0) then
distance = Thickness%depth_zwt(i,j,kbot) - Thickness%depth_zt(i,j,kbot)
argument = law_of_wall_rough_length_r*(distance + law_of_wall_rough_length)
drag_coeff(i,j) = rho0*max(cdbot,vonkarman2/(log(argument))**2)
endif
enddo
enddo
else
do j=jsd,jed
do i=isd,ied
drag_coeff(i,j) = rho0*cdbot_array(i,j)*Grd%tmask(i,j,1)
enddo
enddo
endif
endif ! horz_grid iftest
! modified stress if use idealized wave model
if(cdbot_wave) then
call current_wave_drag_diag(Thickness, Velocity, drag_coeff, Waves, Time, tstep)
endif
! set quadratic bottom momentum drag
! units of bmf are N/m^2
! introduce ceilings for uvmag and bmf
! to minimize potential for instability.
Velocity%bmf(:,:,:) = 0.0
! if cgrid, ignore offset in velocity components for uvmag calculation
do j=jsd,jed
do i=isd,ied
kbot = grd_kbot(i,j)
Velocity%gamma(i,j) = 0.0
if (kbot > 0) then
if (read_tide_speed) then
uvmag = sqrt(uresidual2_2d(i, j) &
+ Velocity%u(i, j, kbot, 1, taum1)**2 &
+ Velocity%u(i, j, kbot, 2, taum1)**2)
else
uvmag = sqrt(uresidual2 + Velocity%u(i, j, kbot, 1, tstep)**2 &
+ Velocity%u(i, j, kbot, 2, tstep)**2)
end if
uvmag = min(uvmag,uvmag_max)
Velocity%gamma(i,j) = drag_coeff(i,j)*uvmag
umag = abs(Velocity%u(i,j,kbot,1,tstep))
vmag = abs(Velocity%u(i,j,kbot,2,tstep))
Velocity%bmf(i,j,1) = min(bmf_max, Velocity%gamma(i,j)*umag)*sign(1.0,Velocity%u(i,j,kbot,1,tstep))
Velocity%bmf(i,j,2) = min(bmf_max, Velocity%gamma(i,j)*vmag)*sign(1.0,Velocity%u(i,j,kbot,2,tstep))
endif
enddo
enddo
if (read_tide_speed) then
deallocate(uresidual2_2d)
end if
! send diagnostics
if (id_drag_coeff > 0) then
used = send_data(id_drag_coeff, rho0r*drag_coeff(:,:), Time%model_time, &
rmask=Grd%mask(:,:,1), is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
if (id_gamma_bmf > 0) then
used = send_data(id_gamma_bmf, Velocity%gamma(:,:), Time%model_time, &
rmask=Grd%mask(:,:,1), is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
if (id_bmf_u > 0) then
used = send_data(id_bmf_u, Velocity%bmf(:,:,1), Time%model_time, &
rmask=Grd%mask(:,:,1), is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
if (id_bmf_v > 0) then
used = send_data(id_bmf_v, Velocity%bmf(:,:,2), Time%model_time, &
rmask=Grd%mask(:,:,1), is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
! some sea level diagnostics from geothermal heating
if (id_eta_tend_geoheat > 0 .or. id_eta_tend_geoheat_glob > 0) then
wrk1_2d = 0.0
do j=jsc,jec
do i=isc,iec
kbot=Grd%kmt(i,j)
if(kbot > 0) then
rhobot_inv = 1.0/(epsln+Dens%rho(i,j,kbot,tau))
alphabot = -rhobot_inv*Dens%drhodT(i,j,kbot)
wrk1_2d(i,j) = -alphabot*rhobot_inv*T_prog(index_temp)%btf(i,j)
endif
enddo
enddo
call diagnose_2d(Time, Grd, id_eta_tend_geoheat, wrk1_2d(:,:))
call diagnose_sum(Time, Grd, Dom, id_eta_tend_geoheat_glob, wrk1_2d, cellarea_r)
endif
end subroutine get_ocean_bbc
! NAME="get_ocean_bbc"
!#######################################################################
!
!
!
! calculates wave-current bottom shear stress
! using model of Grand and Madsen(1979) J. Geophys. Res. 84, C4, 1797
! see Signell et. al (1990) J. Geophys. Res. 95, C6, 9671
!
! input bot_vel: current velocity at u points
! Note: assumed that this is the velocity just above the bottom boundary layer.
!
! A relation between grain size, ripples steepness and and roughness is
! assumed. More general relations are possible but not needed, since
! the output is used only to parameterise resuspension of organic matter
! in the ecosystem model ocean_shared/generic_tracers/generic_ERGOM.F90.
!
! output: effective drag coefficient drag_coeff.
! It is valid for momentum flux from currents into the bottom, but
! from combined waves+current action.
! Velocity%current_wave_stress is the stress from waves and currents
! to the sediment.
!
! Note:
!
! 1/ drag_coeff in this routine arises from both (waves+current);
! That is, ustar**2/uref**2 1 meter above the bottom.
!
! 2/ bottom velocity is taken just above the bottom boundary layer, and
! assumed here to be at lowest u point.
!
! 3/ to understand the calculation of shear stress acting on grains
! at sediment surface ("wrk1_2d(i,j) = (ustar2/ucomb)*0.3152"),
! see Kuhrts et al. (2004) Eqs. (4,5,6,7).
! Assume a thin skin friction layer according to Smith, McLean (1977)
! thickness of the skin friction layer zskin scales with roughness length ruff
! grained sediments are characterised by median diameter d50[m]
! ripples at sea bottom have spacing lambda and height eta with steepness eta/lambda=0.1
! common approx. for grain roughness length = d50/12
! Nielsen (1983) form drag roughness length = (8/30)*eta*(eta/lambda) ==> (8/3000)*lambda
! Yalin (1977) lambda = 1000*d50 ==> form drag = (8/3)*d50
! ruff = grain + form drag = (1/12+8/3)*d50 = (33/12)*d50 = 33*grain
! Smith, McLean (1977) zskin = 0.09*grain*(lambda/grain)**0.8 = 165*grain = 5*ruff
! matching at zskin leads to log(zskin/ruff)/log(zskin/grain)) = log(5)/log(165) = 0.3152
! the current induced skin friction velocity is derived from matching skin friction
! to wave boundary layer.
!
! 4/ Algorithm has yet to be updated for Cgrid.
!
! April 2012
! martin.schmidt@io-warnemuende.de
!
!
!
subroutine current_wave_drag_diag(Thickness, Velocity, drag_coeff, Waves, Time, tstep)
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_velocity_type), intent(inout) :: Velocity
real, dimension(isd:,jsd:), intent(out) :: drag_coeff
type(ocean_wave_type), intent(inout) :: Waves
type(ocean_time_type), intent(in) :: Time
integer, intent(in) :: tstep
real,parameter:: twopi=2.*pi
real :: omega, cdrag, cold, cnew, ruff, rauh, dist, wave_p_u, wave_u_u
real :: bot_vel, u_vel, v_vel ! Velocity in the bottommost layer.
real :: ustar, ustar2, ucomb, ucskin ! friction velocities in wave & skin boundary layers
integer:: i, j, icount, kmu
integer, parameter:: itmax=10
if (first_call) then
if ( .not. wave_model_is_initialised() ) then
call mpp_error(FATAL,&
'==>Error from ocean_bbc_mod (get_ocean_bbc): enable a wave model for option cdbot_wave ')
endif
first_call = .false.
endif
! get the wave induced bottom friction velocities
! wave_u applies inside the wave-boundary layer, wave_s inside the skin friction layer
if (id_iter > 0) wrk2_2d = 0.0
wrk1_2d = 0.0
drag_coeff = 0.0
call wave_u_diag(wave_u, wave_s, Waves)
do j=jsc,jec
do i=isc,iec
kmu=Grd%kmu(i,j)
if (kmu >= 1) then
ruff = roughness_length(i,j) ! bottom roughness length (grain+form drag) [m]
dist = Thickness%depth_zwu(i,j,kmu) - Thickness%depth_zu(i,j,kmu)
cdrag = (von_karman/log(dist/ruff))**2 ! drag coefficient from currents at height=dist above sea bed
u_vel = Velocity%u(i,j,kmu,1,tstep)
v_vel = Velocity%u(i,j,kmu,2,tstep)
bot_vel=sqrt(u_vel*u_vel + v_vel*v_vel)
! find the effective bottom roughness rauh seen by bot_vel above the wave-boundary layer
! determined by matching log-profiles of outer and wave-boundary layer, see Kuhrts et al. (2004) Eq. (1,2)
! Waves%-quantities are defined at the t-grid.
! put wave_u onto the u-grid
wave_u_u = max(wave_u(i,j), wave_u(i+1,j), &
wave_u(i,j+1), wave_u(i+1,j+1))
if (wave_u_u > epsln) then
! put wave_p onto the u-grid
wave_p_u = max(Waves%wave_p(i,j), Waves%wave_p(i+1,j), &
Waves%wave_p(i,j+1), Waves%wave_p(i+1,j+1))
! wave orbital frequency derived from peak frequency [1/s]
omega=twopi*wave_p_u
! iterate now the bottom drag above the wave-boundary layer (Grand and Madsen)
cnew = cdrag ! use cdrag as first guess ! [dimensionless]
do icount=1,itmax
cold = cnew ! [dimensionless]
! calculate current friction velocity with effective bottom drag
ustar = bot_vel*sqrt(cold) ! [m/s]
! add wave friction to current friction velocity assuming same direction
ucomb = sqrt(ustar*ustar + wave_u_u*wave_u_u) ! [m/s]
! determine effective bottom roughness rauh from physical bottom roughness ruff
rauh = ruff*(2*von_karman*ucomb/(omega*ruff+epsln))**& ! [m]
(1.0-ustar/(ucomb+epsln))
! calculate bottom drag above the wave-boundary layer
cnew = (von_karman/log(dist/rauh))**2 ! [dimensionless]
if (id_iter > 0) wrk2_2d(i,j) = icount
if (abs(cold-cnew) .le. 0.01*abs(cnew)) exit
!TS iteration cnew(iter) better than cnew=cdrag
!!! in case of bad convergence, but this must never happen
!! if (icount.eq.itmax) cnew = cdrag
enddo
else ! no waves
cnew=cdrag
endif
! the effective drag coefficient seen by model velocities above the wave boundary layer
drag_coeff(i,j) = cnew * rho0
! the combined current-wave induced stress in wave boundary layer
ustar2 = bot_vel*bot_vel*cnew ! [m/s]
ucomb = sqrt(ustar2 + wave_u_u*wave_u_u) ! [m/s]
! calculate shear stress acting on grains at sediment surface
wrk1_2d(i,j) = (ustar2/ucomb)*0.3152
! the wave induced friction velocity wave_s is derived applying Nielsen (1992) for grain roughness
! store the total wave-current bottom stress in skin friction layer for later use in sediment dynamics
endif ! kmu.ge.1
enddo
enddo
call mpp_update_domains(wrk1_2d(:,:), drag_coeff(:,:), Dom%domain2d)
! place onto the T-grid
do j=jsc,jec
do i=isc,iec
ucskin = max(wrk1_2d(i,j ), wrk1_2d(i-1,j ), &
wrk1_2d(i,j-1), wrk1_2d(i-1,j-1))
Velocity%current_wave_stress(i,j) = rho0*(ucskin*ucskin + wave_s(i,j)*wave_s(i,j))
enddo
enddo
! diagnostics
call diagnose_2d(Time, Grd, id_cur_wav_dr, Velocity%current_wave_stress(:,:))
call diagnose_2d(Time, Grd, id_wave_s, wave_s(:,:))
call diagnose_2d(Time, Grd, id_wave_u, wave_u(:,:))
call diagnose_2d_u(Time, Grd, id_iter, wrk2_2d(:,:))
if(debug_this_module) then
write(stdoutunit,*) 'ocean_wave_model: end wave_drag_diag'
endif
end subroutine current_wave_drag_diag
! NAME="current_wave_drag_diag"
!#######################################################################
!
!
!
! calculates wave bottom shear stress velocity
! wave friction factor is parametrized by approximation of
! Nielsen (1992), Coastal bottom boundary layers and sediment transport!
!
! April 2012
! martin.schmidt@io-warnemuende.de
!
!
!
subroutine wave_u_diag(wave_u, wave_s, Waves)
real, dimension(isd:,jsd:), intent(inout) :: wave_u
real, dimension(isd:,jsd:), intent(inout) :: wave_s
type(ocean_wave_type), intent(inout) :: Waves
real :: omega, U_m, ampli, f_w, wk_times_depth, ruff, grain
integer :: i, j
real,parameter:: twopi=2.*pi
! these fields have units m/s
wave_u=0.0
wave_s=0.0
do j=jsd,jed
do i=isd,ied
! water_depth*wave number [dimensionless]
wk_times_depth = Grd%ht(i,j) * Waves%wave_k(i,j)
!TS uncleasr why wk_times_depth .gt. 0.3 ??? possibly to exclude land ht=0
!! if (wk_times_depth .gt. 0.3 .and. wk_times_depth.le.50.) then
if (wk_times_depth.le.50.) then
ruff = roughness_length(i,j) ! bottom roughness length (grain+form drag) [m]
omega=twopi*Waves%wave_p(i,j) ! wave orbital frequency [1/s]
ampli=0.5*Waves%height(i,j)/(sinh(wk_times_depth)+epsln) ! wave amplitude at sea bottom [m]
U_m=ampli*omega ! wave orbital velocity at sea bottom [m/s]
!Kuhrts et al. (2004) Eq. (3) fw >= 0.3 if ampli >= 43.9315 * ruff
if (ampli.ge.(43.9315*ruff)) then
f_w= exp(5.5*(30.*ruff/(ampli+epsln))**0.2 - 6.3) ! [dimensionless]
else
f_w= 0.3 ! [dimensionless]
endif
wave_u(i,j)=sqrt(0.5*f_w)*U_m ! [m/s]
!Kuhrts et al. (2004) Eq. (7) applied for grain roughness
grain = ruff/33.0
if (ampli.ge.(43.9315*grain)) then
f_w= exp(5.5*(30.*grain/(ampli+epsln))**0.2 - 6.3) ! [dimensionless]
else
f_w= 0.3 ! [dimensionless]
endif
wave_s(i,j)=sqrt(0.5*f_w)*U_m ! [m/s]
endif
enddo
enddo
if(debug_this_module) then
write(stdoutunit,*) 'ocean_wave_model: end wave_u_diag'
endif
end subroutine wave_u_diag
! NAME="wave_u_diag"
end module ocean_bbc_mod