module ocean_form_drag_mod
!
! S.M. Griffies
!
!
!
! Compute form drag as per Greatbatch and Lamb (1990) and/or
! Aiki etal (2004).
!
!
!
! Compute thickness weighted and density weighted tendency
! [velocity*(kg/m^3)*meter/sec] for velocity associated with
! parameterized form drag. Code employs the methods described
! in Greatbatch and Lamb (1990) as well as Aiki etal (2004).
!
! This scheme has not been updated for C-grid specific layout.
! However, there is little reason to expect that shifting the
! viscosity by 1/2 grid point will make sense physically for this
! scheme. Indeed, one does not expect this scheme to be relevant
! near boundaries anyhow. So use for both Bgrid and Cgrid should
! be reasonable, though with the caveat that differences will appear
! near boundaries.
!
!
!
!
!
!
! Greatbatch and Lamb, 1990: On parameterizing vertical mixing of
! momentum in non eddy-resolving ocean models. Journal of
! Physical Oceanography. vol. 20, pages 1634-1637.
!
!
!
! Greatbatch, 1998: Exploring the relationship betwen eddy-induced
! transport velocity, vertical momentum transfer, and the isopycnal
! flux of potential vorticity. Journal of Physical Oceanography,
! vol. 28, pages 422-432.
!
!
!
! Aiki, Jacobson, and Yamagata, 2004: Parameterizing ocean eddy transports
! from surface to bottom. Geophysical Research Letters, vol. 31.
!
!
!
! Ferreira and Marshall, 2006: Formulation and implementation of
! a residual-mean ocean circulation model. Ocean Modelling,
! vol. 13, pages 86--107.
!
!
!
! Danabosaglu and Marshall, 2007: Effects of vertical variations of thickness
! diffusivity in an ocean general circulation model.
! Ocean Modelling, vol. 18, pages 122-141.
!
!
!
!
!
!
! Needs to be true in order to use this scheme. Default is false.
!
!
! For debugging.
!
!
! For verbose initialization information.
!
!
!
! For using the Greatbatch form drag approach, which places
! the contributions from "GM" into the momentum equation.
!
! use_form_drag_gbatch, is the traditional transformed
! Eulerian mean approach as per Greatbatch and Lamb (1990)
! and Greatbatch (1998). In this approach, we modify
! visc_cbu so that the the "GM-effects" occur through
! vertical viscosity.
!
! Since thise scheme is experimental, it is not recommended
! for general use, so we set default
! use_form_drag_gbatch=.false.
!
!
!
! Logical to enable the use of a diabatic layer over which the
! eddy-induce velocity is assumed to be constant with depth.
! Default form_drag_gbatch_surf_layer=.false.
!
!
! Minimum number of vertical grid points in the surface turbulent
! boundary layer for use with form_drag_gbatch_surf_layer=.true.
! Default ksurf_blayer_min=3.
!
!
!
! For use of a vertical viscosity with the form_drag_gbatch that is
! equal to visc = form_drag_gbatch_alpha * f^2. This form of the
! vertical viscosity is used by Ferreira and Marshall, 2006.
! Default form_drag_gbatch_alpha_f2=.false.
!
!
! For use of a vertical viscosity with the form_drag_gbatch that is
! equal to visc = form_drag_gbatch_alpha * f^2.
! Default form_drag_gbatch_alpha=3e8
!
!
!
! To compute vertical viscosity according to visc=kappa*(f/N)**2
! with kappa=gm-diffusivity, f=Coriolis, and N=buoyancy frequency.
! This is the form suggested by Greatbatch and Lamb (1990).
! Default form_drag_gbatch_f2overN2=.false.
!
!
!
! To compute vertical viscosity according to visc=kappa*(f/Nb)**2
! with kappa=gm-diffusivity, f=Coriolis, and Nb=buoyancy frequency
! just below the diabatic mixed layer. This is the form suggested
! by Danabasoglu and Marshall (2007).
! Default form_drag_gbatch_f2overNb2=.false.
!
!
!
! To compute vertical viscosity according to visc=kappa*(f/No)**2
! with kappa=gm-diffusivity, f=Coriolis, and No=constant buoyancy
! frequency set via namelist.
! Default form_drag_gbatch_f2overNo2=.false.
!
!
! To compute vertical viscosity according to visc=kappa*(f/No)**2
! with kappa=gm-diffusivity, f=Coriolis, and No=constant buoyancy
! frequency. Default form_drag_gbatch_No=5e-3
!
!
!
! For vertically smoothing the squared buoyancy frequency for
! use in computing the vertical viscosity in the form drag
! approach. Default form_drag_gbatch_smooth_N2=.false.
!
!
! Number of 1-2-1 passes for vertically smoothing the squared
! buoyancy frequency. Default num_121_passes=1.
!
!
!
! Maximum vertical viscosity used for the form drag contribution
! to vertical friction from the Greatbatch TEM approach.
! Default visc_cbu_form_drag_max=1m^2/sec.
!
!
! For diagnostic purpuses, maximum form drag eddy induced velocity.
! Default vel_form_drag_max=1m/sec.
!
!
! Minimum sequared buoyancy frequency (N^2) for use in computing
! the vertical viscosity from the Greatbatch form drag scheme.
! Default N_squared_min=1e-10.
!
!
!
! For using the Aiki form drag approach.
!
!
!
! Dimensionless parameters from the Aiki etal scheme.
! Default cprime_aiki=0.3.
!
!
! Compute a dimensionless scale function proportional to
! the GM-diffusivity, for use with the Aiki etal form
! drag scheme. Default form_drag_aiki_scale_by_gm=.false.
!
!
!
! For implementing the Aiki form drag just in a selected number of
! bottom layers. Will still insist that the scheme conserves
! momentum, as a form drag scheme should do.
! Default form_drag_aiki_bottom_layer=.false.
!
!
! Number of klevels above the bottom that we choose to close the
! Aiki form drag scheme. Default form_drag_aiki_bottom_klevels=3.
! This should ideally be less than the minimum number of klevels
! in the model.
!
!
! For scaling the coefficient used by the form drag scheme with
! the bottom slope.
!
!
! For scaling the coefficient used by the form drag scheme with
! the bottom slope raised to the power "form_drag_aiki_gradH_power".
!
!
! For scaling setting the maximum of the bottom slope
! for use with scaling the form drag coefficient.
! Default form_drag_aiki_gradH_max=.05
!
!
!
use constants_mod, only: epsln
use diag_manager_mod, only: register_diag_field, register_static_field
use fms_mod, only: stdout, stdlog, FATAL, WARNING, NOTE
use fms_mod, only: write_version_number, open_namelist_file, check_nml_error, close_file
use mpp_mod, only: input_nml_file, mpp_error
use mpp_domains_mod, only: mpp_global_max, mpp_update_domains
use ocean_domains_mod, only: get_local_indices
use ocean_parameters_mod, only: missing_value, onefourth
use ocean_parameters_mod, only: rho0r, onehalf, grav
use ocean_types_mod, only: ocean_domain_type, ocean_grid_type
use ocean_types_mod, only: ocean_thickness_type, ocean_time_type
use ocean_types_mod, only: ocean_velocity_type, ocean_options_type
use ocean_types_mod, only: ocean_density_type, ocean_time_steps_type
use ocean_workspace_mod, only: wrk1, wrk2, wrk3
use ocean_workspace_mod, only: wrk1_v, wrk2_v, wrk3_v
use ocean_workspace_mod, only: wrk1_v2d, wrk2_v2d
use ocean_workspace_mod, only: wrk1_2d, wrk2_2d, wrk3_2d, wrk4_2d
use ocean_util_mod, only: diagnose_2d_u, diagnose_3d_u
implicit none
private
public compute_visc_form_drag
public form_drag_accel
public ocean_form_drag_init
type(ocean_domain_type), pointer :: Dom => NULL()
type(ocean_grid_type), pointer :: Grd => NULL()
real, dimension(:,:), allocatable :: gradH_scale ! (dimensionless) scaling according to grad of topog
real, dimension(:,:,:), allocatable :: agm_diffusivity ! (m^2/sec) GM diffusivity
real, dimension(:,:,:), allocatable :: N_squared ! (m^2/sec) GM diffusivity
! for diagnostics
logical :: used
integer :: id_form_drag_aiki_gradH_scale=-1
integer :: id_form_drag_aiki_u =-1
integer :: id_form_drag_aiki_v =-1
integer :: id_form_drag_aiki_power =-1
integer :: id_form_drag_aiki_coeff =-1
integer :: id_visc_cbu_form_drag_u =-1
integer :: id_visc_cbu_form_drag_v =-1
integer :: id_form_drag_velocity_u =-1
integer :: id_form_drag_velocity_v =-1
integer :: id_form_drag_velocity_fu =-1
integer :: id_form_drag_velocity_fv =-1
integer :: id_N_squared =-1
integer :: id_f2overN2 =-1
integer :: id_f2overNb2 =-1
integer :: id_ksurf_blayer_form_drag =-1
integer :: id_surface_blayer_form_drag=-1
character(len=128) :: version=&
'$Id: ocean_form_drag.F90,v 20.0 2013/12/14 00:16:34 fms Exp $'
character (len=128) :: tagname = &
'$Name: tikal $'
integer :: isd, ied, jsd, jed, isc, iec, jsc, jec, nk
real :: onepcprime_aiki2_r
real :: cprime_aiki2
! nml for Aiki form drag scheme
logical :: use_form_drag_aiki =.false.
logical :: form_drag_aiki_bottom_layer =.false.
logical :: form_drag_aiki_scale_by_gradH =.false.
real :: form_drag_aiki_gradH_max = 0.05
real :: form_drag_aiki_gradH_power = 1.0
integer :: form_drag_aiki_bottom_klevels = 3
real :: cprime_aiki = 0.3
! nml for Gbatch form drag scheme
logical :: use_form_drag_gbatch =.false.
logical :: form_drag_gbatch_surf_layer =.false.
logical :: form_drag_gbatch_alpha_f2 =.false.
logical :: form_drag_gbatch_f2overN2 =.false.
logical :: form_drag_gbatch_f2overNb2 =.false.
logical :: form_drag_gbatch_f2overNo2 =.false.
logical :: form_drag_gbatch_smooth_N2 =.false.
logical :: form_drag_aiki_scale_by_gm =.false.
integer :: num_121_passes = 1
integer :: ksurf_blayer_min = 3
real :: form_drag_gbatch_alpha = 3e8
real :: form_drag_gbatch_No = 5e-3
real :: visc_cbu_form_drag_max = 1.0
real :: vel_form_drag_max = 1.0
real :: N_squared_min = 1e-10
real :: agm_form_drag = 600.0
real :: agm_r = 0.0016666667
logical :: verbose_init = .true.
logical :: use_this_module = .false.
logical :: debug_this_module = .false.
logical :: module_is_initialized = .FALSE.
namelist /ocean_form_drag_nml/ verbose_init, use_this_module, debug_this_module, &
use_form_drag_aiki, cprime_aiki, form_drag_aiki_bottom_layer, &
form_drag_aiki_bottom_klevels, form_drag_aiki_scale_by_gradH, &
form_drag_aiki_gradH_max, form_drag_aiki_gradH_power, &
form_drag_aiki_scale_by_gm, &
use_form_drag_gbatch, visc_cbu_form_drag_max, vel_form_drag_max, &
N_squared_min, agm_form_drag, &
form_drag_gbatch_surf_layer, ksurf_blayer_min, &
form_drag_gbatch_alpha_f2, form_drag_gbatch_alpha, &
form_drag_gbatch_f2overN2, form_drag_gbatch_f2overNb2, &
form_drag_gbatch_smooth_N2, num_121_passes, &
form_drag_gbatch_f2overNo2, form_drag_gbatch_No
contains
!#######################################################################
!
!
!
! Initial set up for parameterized pressure form drag.
!
!
subroutine ocean_form_drag_init(Grid, Domain, Time, Time_steps, Ocean_options, debug)
type(ocean_grid_type), intent(in), target :: Grid
type(ocean_domain_type), intent(in), target :: Domain
type(ocean_time_type), intent(in) :: Time
type(ocean_time_steps_type), intent(in) :: Time_steps
type(ocean_options_type), intent(inout) :: Ocean_options
logical, intent(in), optional :: debug
integer :: io_status, ioun, ierr
integer :: i,j
integer :: num_visc_schemes
integer :: stdoutunit,stdlogunit
stdoutunit=stdout();stdlogunit=stdlog()
if (module_is_initialized ) then
call mpp_error(FATAL, &
'==>Error in ocean_form_drag_mod (ocean_form_drag_init): module already initialized')
endif
module_is_initialized = .TRUE.
call write_version_number(version, tagname)
call get_local_indices(Domain,isd,ied,jsd,jed,isc,iec,jsc,jec)
nk = Grid%nk
Dom => Domain
Grd => Grid
! provide for namelist over-ride of default values
#ifdef INTERNAL_FILE_NML
read (input_nml_file, nml=ocean_form_drag_nml, iostat=io_status)
ierr = check_nml_error(io_status,'ocean_form_drag_nml')
#else
ioun = open_namelist_file()
read (ioun,ocean_form_drag_nml,IOSTAT=io_status)
ierr = check_nml_error(io_status,'ocean_form_drag_nml')
call close_file (ioun)
#endif
write (stdlogunit,ocean_form_drag_nml)
write (stdoutunit,'(/)')
write (stdoutunit,ocean_form_drag_nml)
! this array needs allocation even if not using module
allocate (agm_diffusivity(isd:ied,jsd:jed,nk))
agm_diffusivity=0.0
if(use_this_module) then
if(use_form_drag_gbatch .and. use_form_drag_aiki) then
Ocean_options%form_drag = 'Used form drag from BOTH Gbatch and Aiki: experimental, so beware.'
write(stdoutunit,'(a)') &
'==>Note: ocean_form_drag_mod using BOTH Gbatch and Aiki schemes: experimental, so beware!'
write(stdoutunit,'(a)') &
'==>Note: May choose to set "use_gm_skew=.false." in ocean_neutral_physics to avoid double counting.'
endif
if(use_form_drag_gbatch .and. .not. use_form_drag_aiki) then
Ocean_options%form_drag = 'Used form drag from Greatbatch: experimental, so beware.'
write(stdoutunit,'(a)') &
'==>Note: ocean_form_drag_mod using Greatbatch scheme: experimental, so beware'
write(stdoutunit,'(a)') &
'==>Note: May choose to set "use_gm_skew=.false." in ocean_neutral_physics to avoid double counting.'
endif
if(.not. use_form_drag_gbatch .and. use_form_drag_aiki) then
Ocean_options%form_drag = 'Used form drag from Aiki etal: experimental, so beware.'
write(stdoutunit,'(a)') &
'==>Note: ocean_form_drag_mod using Aiki scheme: experimental, so beware'
write(stdoutunit,'(a)') &
'==>Note: May choose to set "use_gm_skew=.false." in ocean_neutral_physics to avoid double counting.'
endif
if(.not. use_form_drag_gbatch .and. .not. use_form_drag_aiki) then
write(stdoutunit,'(a)') &
'==>Note: ocean_form_drag_mod: use_this_module=.true. but use_form_drag_gbatch=.false. &use_form_drag_aiki=.false.'
Ocean_options%form_drag = 'Did NOT use parameterized form drag.'
endif
else
call mpp_error(NOTE, '==>Note from ocean_form_drag_mod: NOT USING this module')
Ocean_options%form_drag = 'Did NOT use parameterized form drag.'
return
endif
if (PRESENT(debug) .and. .not. debug_this_module) then
debug_this_module = debug
endif
if(debug_this_module) then
write(stdoutunit,'(a)') '==>Note: running ocean_form_drag_mod with debug_this_module=.true.'
endif
if(Time_steps%aidif /= 1.0) then
call mpp_error(FATAL, &
'==>Error from ocean_form_drag_mod: aidif=1.0 required to handle form drag implicitly.')
endif
if(use_form_drag_gbatch) then
num_visc_schemes = 0
if(form_drag_gbatch_alpha_f2) then
write(stdoutunit,'(1x,a)') &
'==> Note from ocean_form_drag_mod: Vertical viscosity for TEM method is set to constant*f^2'
num_visc_schemes = 1 + num_visc_schemes
endif
if(form_drag_gbatch_f2overN2) then
write(stdoutunit,'(1x,a)') &
'==> Note from ocean_form_drag_mod: Vertical viscosity for TEM method is set to kappa*(f/N)^2'
num_visc_schemes = 1 + num_visc_schemes
endif
if(form_drag_gbatch_f2overNb2) then
write(stdoutunit,'(1x,a)') &
'==> Note from ocean_form_drag_mod: Vertical viscosity for TEM method is set to kappa*(f/Nb)^2'
num_visc_schemes = 1 + num_visc_schemes
endif
if(form_drag_gbatch_f2overNo2) then
write(stdoutunit,'(1x,a)') &
'==> Note from ocean_form_drag_mod: Vertical viscosity for TEM method is set to kappa*(f/No)^2'
num_visc_schemes = 1 + num_visc_schemes
endif
if (num_visc_schemes == 0) then
call mpp_error(WARNING, &
'==>Warning from ocean_form_drag_init: no method chosen for vertical viscosity. Gbatch form drag will do nothing.')
endif
if (num_visc_schemes > 1) then
call mpp_error(FATAL, &
'==>Error from ocean_form_drag_init: can choose just one method to compute vertical viscosity')
endif
if(form_drag_gbatch_surf_layer) then
write(stdoutunit,'(1x,a)') &
'==> Note from ocean_form_drag_mod: Vertical viscosity for TEM method modified in surf turb layer'
endif
endif
! commonly used constants for Aiki scheme
onepcprime_aiki2_r = 1.0/(1.0 + cprime_aiki*cprime_aiki)
cprime_aiki2 = cprime_aiki*cprime_aiki
! squared buoyancy frequency
allocate (N_squared(isd:ied,jsd:jed,nk))
N_squared = 0.0
! scaling proportional to bottom topog gradient
allocate (gradH_scale(isd:ied,jsd:jed))
gradH_scale(:,:) = 1.0
if(form_drag_aiki_scale_by_gradH) then
do j=jsd,jed
do i=isd,ied
gradH_scale(i,j) = &
(Grd%gradH(i,j)/(epsln+form_drag_aiki_gradH_max))**form_drag_aiki_gradH_power
gradH_scale(i,j) = min(gradH_scale(i,j),1.0)
enddo
enddo
endif
! Aiki related diagnostics
id_form_drag_aiki_gradH_scale = register_static_field ('ocean_model', 'form_drag_aiki_gradH_scale', &
Grd%vel_axes_uv(1:2), 'Topography scaling for Aiki form drag', 'dimensionless', &
missing_value=missing_value, range=(/-1e6,1e6/))
call diagnose_2d_u(Time, Grd, id_form_drag_aiki_gradH_scale, gradH_scale(:,:))
id_form_drag_aiki_u = register_diag_field ('ocean_model', 'form_drag_aiki_u', &
Grd%vel_axes_u(1:3), Time%model_time, &
'Aiki i-component form stress', 'N/m2', &
missing_value=missing_value, range=(/-1e6,1e6/))
id_form_drag_aiki_v = register_diag_field ('ocean_model', 'form_drag_aiki_v', &
Grd%vel_axes_v(1:3), Time%model_time, &
'Aiki j-component form stress', 'N/m2', &
missing_value=missing_value, range=(/-1e6,1e6/))
id_form_drag_aiki_power = register_diag_field ('ocean_model', 'form_drag_aiki_power', &
Grd%vel_axes_uv(1:3), Time%model_time, &
'Power sink from Aiki form drag', 'Watt', &
missing_value=missing_value, range=(/-1e16,1e16/))
id_form_drag_aiki_coeff= register_diag_field ('ocean_model', 'form_drag_aiki_coeff', &
Grd%vel_axes_wu(1:3), Time%model_time, &
'dimensionless form scaling coefficient for Aiki scheme', &
'dimensionless', missing_value=missing_value, range=(/-10.0,1e6/))
! Greatbatch related diagnostics
id_N_squared = register_diag_field ('ocean_model', 'N_sqrd_gbatch', Grd%vel_axes_wu(1:3), &
Time%model_time, 'N^2 on T-cell bottom for Gbatch form drag', 'sec^-2', &
missing_value=missing_value, range=(/-1.e2,1e8/))
id_f2overN2 = register_diag_field ('ocean_model', 'f2overN2', Grd%vel_axes_wu(1:3), &
Time%model_time, '(f/N)**2 on U-cell bottom', 'dimensionless', &
missing_value=missing_value, range=(/-1.0,1e6/))
id_f2overNb2= register_diag_field ('ocean_model', 'f2overNb2', Grd%vel_axes_wu(1:2), &
Time%model_time, '(f/N)**2 on U-cell at base of surf mixed layer', 'sec^2',&
missing_value=missing_value, range=(/-1.0,1e6/))
id_visc_cbu_form_drag_u= register_diag_field ('ocean_model', 'visc_cbu_form_drag_u', &
Grd%vel_axes_u(1:3), Time%model_time, &
'vertical viscosity from form drag on u-velocity component', 'm^2/s', &
missing_value=missing_value, range=(/-10.0,1e6/), &
standard_name='ocean_vertical_momentum_diffusivity_due_to_form_drag')
id_visc_cbu_form_drag_v= register_diag_field ('ocean_model', 'visc_cbu_form_drag_v', &
Grd%vel_axes_v(1:3), Time%model_time, &
'vertical viscosity from form drag on v-velocity component', 'm^2/s', &
missing_value=missing_value, range=(/-10.0,1e6/))
id_form_drag_velocity_u= register_diag_field ('ocean_model', 'form_drag_velocity_u', &
Grd%vel_axes_u(1:3), Time%model_time, &
'eddy induced i-velocity from Gbatch form drag', 'm/s', &
missing_value=missing_value, range=(/-10.0,10.0/))
id_form_drag_velocity_v= register_diag_field ('ocean_model', 'form_drag_velocity_v', &
Grd%vel_axes_v(1:3), Time%model_time, &
'eddy induced j-velocity from Gbatch form drag', 'm/s', &
missing_value=missing_value, range=(/-10.0,10.0/))
id_form_drag_velocity_fu= register_diag_field ('ocean_model', 'form_drag_velocity_fu', &
Grd%vel_axes_uv(1:3), Time%model_time, &
'f * eddy induced i-velocity from Gbatch form drag', 'm/s', &
missing_value=missing_value, range=(/-10.0,10.0/))
id_form_drag_velocity_fv= register_diag_field ('ocean_model', 'form_drag_velocity_fv', &
Grd%vel_axes_uv(1:3), Time%model_time, &
'f * eddy induced j-velocity from Gbatch form drag', 'm/s', &
missing_value=missing_value, range=(/-10.0,10.0/))
id_ksurf_blayer_form_drag= register_diag_field ('ocean_model', 'ksurf_blayer_form_drag', &
Grd%vel_axes_uv(1:2), Time%model_time, &
'k-index at the base of the boundary layer for GBatch form drag', &
'unitless', missing_value=missing_value, range=(/-10.0,1e10/))
id_surface_blayer_form_drag= register_diag_field ('ocean_model', 'surface_blayer_form_drag', &
Grd%vel_axes_uv(1:2), Time%model_time, &
'U-cell depth of surface boundary layer used in Gbatch form drag', &
'm', missing_value=missing_value, range=(/-10.0,1e10/))
end subroutine ocean_form_drag_init
! NAME="ocean_form_drag_init"
!#######################################################################
!
!
! Compute vertical viscosity arising from parameterized
! form drag according to the methods from the Greatbatch TEM approach.
!
! Follow the method of Ferreira and Marshall (2006)
! for transitioning through the diabatic surface mixed layer.
! We take the surf_blthick from mixed layer schemes (e.g., KPP), and use
! this to define the region where the eddy induced velocity has
! zero vertical shear. Note that we do not introduce an additional
! transitional layer.
!
! There are four methods to compute the vertical viscosity within
! the interior, with one required to be chosen:
! 1/ form_drag_gbatch_f2overN2
! visc = kappa*(f/N)^2, with kappa=gm-diffusivity
! 2/ form_drag_gbatch_f2overNb2
! visc = kappa*(f/Nb)^2, with Nb=buoyancy freq just below surface mixed layer
! 3/ form_drag_gbatch_f2overNo2
! visc = kappa*(f/No)^2, with No=constant buoyancy freq
! 4/ form_drag_gbatch_alpha_f2
! visc = alpha*f^2, with alpha = constant (units of m^2*sec)
!
!
!
!
subroutine compute_visc_form_drag(Time, Thickness, Velocity, Dens, &
gm_diffusivity, surface_blt, visc_cbu_form_drag)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_velocity_type), intent(in) :: Velocity
type(ocean_density_type), intent(in) :: Dens
real, dimension(isd:,jsd:,:), intent(in) :: gm_diffusivity
real, dimension(isd:,jsd:), intent(in) :: surface_blt
real, dimension(isd:,jsd:,:,:), intent(inout) :: visc_cbu_form_drag
integer :: i,j,k,kb,m,n,tau,kbot
integer :: ksurf_blayer(isd:ied,jsd:jed)
real :: surface_blayerT(isd:ied,jsd:jed)
real :: surface_blayerU(isd:ied,jsd:jed)
real :: N_squaredU, sign_switch
real :: tmp,N2_prev, No_r
logical :: below_surf_blayer
if(.not. use_this_module .or. .not. use_form_drag_gbatch) then
visc_cbu_form_drag(:,:,:,:) = 0.0
agm_diffusivity(:,:,:) = gm_diffusivity(:,:,:)
return
endif
tau = Time%tau
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3_v(:,:,:,:) = 0.0
visc_cbu_form_drag(:,:,:,:) = 0.0
No_r = 1.0/(epsln+form_drag_gbatch_No)
agm_diffusivity(:,:,:) = gm_diffusivity(:,:,:)
! compute N^2
do k=1,nk
do j=jsd,jed
do i=isd,ied
N_squared(i,j,k) = -grav*rho0r*Dens%drhodz_zt(i,j,k)
enddo
enddo
enddo
if(form_drag_gbatch_smooth_N2) then
do m=1,num_121_passes
do j=jsd,jed
do i=isd,ied
N2_prev = onefourth*N_squared(i,j,1)
kbot=Grd%kmt(i,j)
if(kbot>1) then
do k=2,kbot-2
tmp = N_squared(i,j,k)
N_squared(i,j,k) = N2_prev + onehalf*N_squared(i,j,k) &
+ onefourth*N_squared(i,j,k+1)
N2_prev = onefourth*tmp
enddo
endif
enddo
enddo
enddo
endif
! squared buoyancy frequency on U-cell
do k=1,nk-1
do j=jsc,jec
do i=isc,iec
if(Grd%umask(i,j,k) > 0.0) then
N_squaredU = onefourth*(N_squared(i,j,k) + N_squared(i+1,j,k) &
+ N_squared(i,j+1,k) + N_squared(i+1,j+1,k))
wrk2(i,j,k) = abs(N_squaredU)
wrk1(i,j,k) = Grd%f(i,j)**2/(wrk2(i,j,k) + N_squared_min)
endif
enddo
enddo
enddo
! viscosity from Greatbatch and Lamb (1990)
if(form_drag_gbatch_f2overN2) then
do k=1,nk-1
do j=jsc,jec
do i=isc,iec
if(Grd%umask(i,j,k) > 0.0) then
wrk3_v(i,j,k,1) = Grd%umask(i,j,k+1) &
*min(visc_cbu_form_drag_max,wrk1(i,j,k)*agm_diffusivity(i,j,k))
wrk3_v(i,j,k,2) = wrk3_v(i,j,k,1)
visc_cbu_form_drag(i,j,k,1) = wrk3_v(i,j,k,1)
visc_cbu_form_drag(i,j,k,2) = wrk3_v(i,j,k,2)
endif
enddo
enddo
enddo
endif
! viscosity from Greatbatch and Lamb (1990) with constant N
if(form_drag_gbatch_f2overNo2) then
do k=1,nk-1
do j=jsc,jec
do i=isc,iec
if(Grd%umask(i,j,k) > 0.0) then
wrk3_v(i,j,k,1) = Grd%umask(i,j,k+1) &
*min(visc_cbu_form_drag_max,agm_diffusivity(i,j,k)*(No_r*Grd%f(i,j))**2)
wrk3_v(i,j,k,2) = wrk3_v(i,j,k,1)
visc_cbu_form_drag(i,j,k,1) = wrk3_v(i,j,k,1)
visc_cbu_form_drag(i,j,k,2) = wrk3_v(i,j,k,2)
endif
enddo
enddo
enddo
endif
! viscosity from Ferreira and Marshall (2006)
if(form_drag_gbatch_alpha_f2) then
do k=1,nk-1
do j=jsc,jec
do i=isc,iec
if(Grd%umask(i,j,k) > 0.0) then
visc_cbu_form_drag(i,j,k,1) = form_drag_gbatch_alpha*Grd%f(i,j)**2
visc_cbu_form_drag(i,j,k,2) = form_drag_gbatch_alpha*Grd%f(i,j)**2
endif
enddo
enddo
enddo
endif
! determine some properties of the surface diabatic layer
if(form_drag_gbatch_surf_layer .or. form_drag_gbatch_f2overNb2) then
ksurf_blayer(:,:) = 1
surface_blayerT(:,:) = 0.0
surface_blayerU(:,:) = 0.0
wrk4_2d(:,:) = 0.0
! compute boundary layer depth on U-cell
do j=jsc,jec
do i=isc,iec
surface_blayerT(i,j) = surface_blt(i,j)
enddo
enddo
call mpp_update_domains (surface_blayerT(:,:), Dom%domain2d)
surface_blayerU(:,:) = 0.0
do j=jsc,jec
do i=isc,iec
surface_blayerU(i,j) = onefourth*Grd%umask(i,j,k) &
*(surface_blayerT(i,j) +surface_blayerT(i+1,j) &
+surface_blayerT(i,j+1)+surface_blayerT(i+1,j+1))
enddo
enddo
! compute k-level at or just below surface boundary layer.
! initialize ksurf_blayer to 1 rather than 0, to avoid
! accessing the zeroth element of an array.
wrk3_2d(:,:) = 0.0
do j=jsc,jec
do i=isc,iec
kb = Grd%kmt(i,j)
below_surf_blayer= .false.
ksurf_blayer(i,j)= ksurf_blayer_min
ksurf_blayer_loop: do k=ksurf_blayer_min,kb-1
if(Thickness%depth_zwu(i,j,k) >= surface_blayerU(i,j) &
.and. .not. below_surf_blayer) then
ksurf_blayer(i,j) = k
below_surf_blayer = .true.
exit ksurf_blayer_loop
endif
enddo ksurf_blayer_loop
! for cases where blayer reaches to the bottom.
! need to give ksurf_blayer one above the bottom so
! to be able to compute a nonzero vertical shear.
if(.not. below_surf_blayer) then
ksurf_blayer(i,j) = max(kb-1,ksurf_blayer_min)
endif
! need a real field to send to diagnostic manager
wrk3_2d(i,j) = float(ksurf_blayer(i,j))
! for the form_drag_gbatch_f2overNb2 scheme, save
! (f/Nb)**2, where Nb=buoyancy frequency one level beneath
! diabatic mixed layer.
k = min(nk,1+ksurf_blayer(i,j))
wrk4_2d(i,j) = Grd%f(i,j)**2/(wrk2(i,j,k) + N_squared_min)
enddo
enddo
endif
! viscosity from Danabasoglu and Marshall (2007)
if(form_drag_gbatch_f2overNb2) then
do k=1,nk-1
do j=jsc,jec
do i=isc,iec
if(Grd%umask(i,j,k) > 0.0) then
wrk3_v(i,j,k,1) = Grd%umask(i,j,k+1) &
*min(visc_cbu_form_drag_max,wrk4_2d(i,j)*agm_diffusivity(i,j,k))
wrk3_v(i,j,k,2) = wrk3_v(i,j,k,1)
visc_cbu_form_drag(i,j,k,1) = wrk3_v(i,j,k,1)
visc_cbu_form_drag(i,j,k,2) = wrk3_v(i,j,k,2)
endif
enddo
enddo
enddo
endif
! modify viscosity in surface boundary layer
if(form_drag_gbatch_surf_layer) then
do n=1,2
do j=jsc,jec
do i=isc,iec
! compute vertical shear at blayer base
k = ksurf_blayer(i,j)
wrk1_v2d(i,j,n) = Grd%umask(i,j,k+1) &
*(Velocity%u(i,j,k,n,tau)-Velocity%u(i,j,k+1,n,tau))/Thickness%dzwu(i,j,k)
! modify the vertical viscosity inside blayer.
! use linear tapering function, going to zero as approach ocean top.
do k=1,ksurf_blayer(i,j)
if(surface_blayerU(i,j) > 0.0) then
wrk3_v(i,j,k,n) = wrk3_v(i,j,ksurf_blayer(i,j),n) &
*(Thickness%depth_zwu(i,j,k)/surface_blayerU(i,j))*wrk1_v2d(i,j,n) &
*Thickness%dzwu(i,j,k)/(Velocity%u(i,j,k,n,tau)-Velocity%u(i,j,k+1,n,tau)+epsln)
wrk3_v(i,j,k,n) = Grd%umask(i,j,k+1)*min(visc_cbu_form_drag_max,abs(wrk3_v(i,j,k,n)))
visc_cbu_form_drag(i,j,k,n) = wrk3_v(i,j,k,n)
endif
enddo
enddo
enddo
enddo
endif
! send some diagnostics
call diagnose_3d_u(Time, Grd, id_N_squared, N_squared(:,:,:))
call diagnose_2d_u(Time, Grd, id_ksurf_blayer_form_drag, wrk3_2d(:,:))
call diagnose_2d_u(Time, Grd, id_surface_blayer_form_drag, surface_blayerU(:,:))
call diagnose_3d_u(Time, Grd, id_f2overN2, wrk1(:,:,:))
call diagnose_2d_u(Time, Grd, id_f2overNb2, wrk4_2d(:,:))
call diagnose_3d_u(Time, Grd, id_visc_cbu_form_drag_u, visc_cbu_form_drag(:,:,:,1))
call diagnose_3d_u(Time, Grd, id_visc_cbu_form_drag_v, visc_cbu_form_drag(:,:,:,2))
! diagnose the eddy induced velocity
if(id_form_drag_velocity_u > 0 .or. id_form_drag_velocity_v > 0 .or. &
id_form_drag_velocity_fu > 0 .or. id_form_drag_velocity_fv > 0) then
wrk1_v(:,:,:,:) = 0.0
wrk2_v(:,:,:,:) = 0.0
wrk1_2d(:,:) = 0.0
wrk2_2d(:,:) = 0.0
do n=1,2
sign_switch = (-1.0)**(n-1)
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = 0.0
enddo
enddo
do k=1,nk-1
do j=jsc,jec
do i=isc,iec
wrk2_2d(i,j) = Grd%umask(i,j,k+1)*wrk3_v(i,j,k,n) &
*(Velocity%u(i,j,k,n,tau)-Velocity%u(i,j,k+1,n,tau)) &
/Thickness%dzwu(i,j,k)
wrk2_v(i,j,k,n) = sign_switch*Grd%umask(i,j,k)*(wrk1_2d(i,j)-wrk2_2d(i,j)) &
/Thickness%dzu(i,j,k)
wrk1_v(i,j,k,n) = wrk2_v(i,j,k,n)/(Grd%f(i,j)+epsln)
wrk1_2d(i,j) = wrk2_2d(i,j)
enddo
enddo
enddo
do j=jsc,jec
do i=isc,iec
wrk2_v(i,j,nk,n) = sign_switch*Grd%umask(i,j,nk)*wrk1_2d(i,j) &
/Thickness%dzu(i,j,k)
wrk1_v(i,j,nk,n) = wrk2_v(i,j,nk,n)/(Grd%f(i,j)+epsln)
enddo
enddo
do k=1,nk
do j=jsc,jec
do i=isc,iec
if(abs(wrk1_v(i,j,k,n)) > vel_form_drag_max) then
wrk1_v(i,j,k,n) = vel_form_drag_max*sign(1.0,wrk1_v(i,j,k,n))
endif
enddo
enddo
enddo
enddo
call diagnose_3d_u(Time, Grd, id_form_drag_velocity_fu, wrk2_v(:,:,:,1))
call diagnose_3d_u(Time, Grd, id_form_drag_velocity_fv, wrk2_v(:,:,:,2))
call diagnose_3d_u(Time, Grd, id_form_drag_velocity_u, wrk1_v(:,:,:,1))
call diagnose_3d_u(Time, Grd, id_form_drag_velocity_v, wrk1_v(:,:,:,2))
endif
end subroutine compute_visc_form_drag
! NAME="compute_visc_form_drag"
!#######################################################################
!
!
!
! Compute the rho*dz weighted acceleration from the Aiki etal
! form drag scheme.
!
!
subroutine form_drag_accel(Time, Thickness, Velocity, energy_analysis_step)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_velocity_type), intent(inout) :: Velocity
logical, intent(in) :: energy_analysis_step
integer :: i,j,k,n
integer :: tau, taum1
integer :: kmu
real :: agm_ucell, agm_max, agm_max_r
if(.not. use_this_module .or. .not. use_form_drag_aiki) then
if(energy_analysis_step) then
do n=1,2
do k=1,nk
do j=jsc,jec
do i=isc,iec
Velocity%wrkv(i,j,k,n) = 0.0
enddo
enddo
enddo
enddo
endif
return
endif
if ( .not. module_is_initialized ) then
call mpp_error(FATAL, &
'==>Error in ocean_form_drag_mod (form_drag_accel): module needs to be initialized')
endif
taum1 = Time%taum1
tau = Time%tau
wrk1 = 1.0
! Compute a dimensionless scale function proportional to GM-diffusivity
if(form_drag_aiki_scale_by_gm) then
agm_max = mpp_global_max(Dom%domain2d, agm_diffusivity(:,:,1))
if(agm_max > 0.0) then
agm_max_r = 1.0/agm_max
do k=1,nk
do j=jsc,jec
do i=isc,iec
agm_ucell = onefourth*Grd%umask(i,j,k) &
*(agm_diffusivity(i,j,k) +agm_diffusivity(i+1,j,k) &
+agm_diffusivity(i,j+1,k)+agm_diffusivity(i+1,j+1,k))
wrk1(i,j,k) = agm_ucell*agm_r
enddo
enddo
enddo
endif
endif
wrk1_2d(:,:) = 0.0
wrk2_2d(:,:) = 0.0
wrk1_v2d(:,:,:) = 0.0
wrk2_v2d(:,:,:) = 0.0
wrk1_v(:,:,:,:) = 0.0
! pre-calculate some 2d coefficients and the inverse mass per area of a column
do j=jsc,jec
do i=isc,iec
wrk1_v2d(i,j,1) = abs(Grd%f(i,j))*cprime_aiki*onepcprime_aiki2_r
wrk1_v2d(i,j,2) = Grd%f(i,j)*cprime_aiki2*onepcprime_aiki2_r
if(Grd%umask(i,j,1) > 0.0) then
wrk1_2d(i,j) = 1.0/Thickness%mass_u(i,j,tau)
endif
enddo
enddo
! vertical integral of equation (11) in Aiki etal paper.
! NOTE: wrk1 is a dimensionless coefficient which
! scales the drag according to GM-diffusivity (from ocean_nphysics_mod)
! or slope of the topography (default)
do k=1,nk
do j=jsc,jec
do i=isc,iec
kmu = Grd%kmu(i,j)
if(k <= kmu) then
wrk1_v(i,j,k,1) = -wrk1_v2d(i,j,1)*(Velocity%u(i,j,k,1,tau)-Velocity%u(i,j,kmu,1,tau)) &
-wrk1_v2d(i,j,2)*(Velocity%u(i,j,k,2,tau)-Velocity%u(i,j,kmu,2,tau))
wrk1_v(i,j,k,2) = -wrk1_v2d(i,j,1)*(Velocity%u(i,j,k,2,tau)-Velocity%u(i,j,kmu,2,tau)) &
+wrk1_v2d(i,j,2)*(Velocity%u(i,j,k,1,tau)-Velocity%u(i,j,kmu,1,tau))
wrk1_v(i,j,k,1) = wrk1_v(i,j,k,1)*wrk1(i,j,k)*gradH_scale(i,j)
wrk1_v(i,j,k,2) = wrk1_v(i,j,k,2)*wrk1(i,j,k)*gradH_scale(i,j)
endif
enddo
enddo
enddo
if(form_drag_aiki_bottom_layer) then
! zero those levels not used
do j=jsc,jec
do i=isc,iec
do k=1,Grd%kmu(i,j)-(form_drag_aiki_bottom_klevels-1)
wrk1_v(i,j,k,1) = 0.0
wrk1_v(i,j,k,2) = 0.0
enddo
enddo
enddo
! mass per area in the retained levels
wrk1_2d(:,:) = 0.0
do j=jsc,jec
do i=isc,iec
do k=(Grd%kmu(i,j)-form_drag_aiki_bottom_klevels),Grd%kmu(i,j)
wrk1_2d(i,j) = wrk1_2d(i,j) + Thickness%rho_dzu(i,j,k,tau)
enddo
enddo
enddo
do j=jsc,jec
do i=isc,iec
if(wrk1_2d(i,j) > 0.0) then
wrk1_2d(i,j) = 1.0/wrk1_2d(i,j)
endif
enddo
enddo
endif
do n=1,2
! vertical integral
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk2_v2d(i,j,n) = wrk2_v2d(i,j,n) + Thickness%rho_dzu(i,j,k,tau)*wrk1_v(i,j,k,n)
enddo
enddo
enddo
! divide by mass per area of the column
do j=jsc,jec
do i=isc,iec
wrk2_v2d(i,j,n) = wrk2_v2d(i,j,n)*wrk1_2d(i,j)
enddo
enddo
enddo
! ensure SGS form drag has zero vertical integral
if(form_drag_aiki_bottom_layer) then
! compute form drag which has zero vertical integral over retained levels
do n=1,2
do j=jsc,jec
do i=isc,iec
do k=(Grd%kmu(i,j)-form_drag_aiki_bottom_klevels),Grd%kmu(i,j)
wrk1_v(i,j,k,n) = Thickness%rho_dzu(i,j,k,tau)*(wrk1_v(i,j,k,n) - wrk2_v2d(i,j,n))
enddo
enddo
enddo
enddo
else
! compute form drag which has zero vertical integral over full column
do n=1,2
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1_v(i,j,k,n) = Thickness%rho_dzu(i,j,k,tau)*(wrk1_v(i,j,k,n) - wrk2_v2d(i,j,n))
enddo
enddo
enddo
enddo
endif
! energetic analysis
if(energy_analysis_step) then
do n=1,2
do k=1,nk
do j=jsc,jec
do i=isc,iec
Velocity%wrkv(i,j,k,n) = wrk1_v(i,j,k,n)
enddo
enddo
enddo
enddo
! add to acceleration
else
do n=1,2
do k=1,nk
do j=jsc,jec
do i=isc,iec
Velocity%accel(i,j,k,n) = Velocity%accel(i,j,k,n) + wrk1_v(i,j,k,n)
enddo
enddo
enddo
enddo
call diagnose_3d_u(Time, Grd, id_form_drag_aiki_u, wrk1_v(isc:iec,jsc:jec,:,1))
call diagnose_3d_u(Time, Grd, id_form_drag_aiki_v, wrk1_v(isc:iec,jsc:jec,:,2))
if (id_form_drag_aiki_power > 0) then
wrk2(:,:,:) = 0.0
do n=1,2
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk2(i,j,k) = wrk2(i,j,k) + Velocity%u(i,j,k,n,tau)*wrk1_v(i,j,k,n)*Grd%dau(i,j)
enddo
enddo
enddo
enddo
call diagnose_3d_u(Time, Grd, id_form_drag_aiki_power, wrk2(isc:iec,jsc:jec,:))
endif
if (id_form_drag_aiki_coeff > 0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk2(i,j,k) = wrk1(i,j,k)*gradH_scale(i,j)
enddo
enddo
enddo
call diagnose_3d_u(Time, Grd, id_form_drag_aiki_coeff, wrk2(:,:,:))
endif
endif ! endif for energy_analysis_step
end subroutine form_drag_accel
! NAME="form_drag_accel"
end module ocean_form_drag_mod