module ocean_advection_velocity_mod
#define COMP isc:iec,jsc:jec
!
! S.M. Griffies
!
!
!
! Advection velocity components for tracer and momenta transport, with
! options for B-grid or C-grid.
!
!
!
! The module computes the horizontal and vertical components to the
! advection velocity on the face of tracer and velocity cells.
! Options are provided for either the B-grid or C-grid.
!
! The three velocity components are related by continuity.
!
! All components are density weighted, and the
! horizontal components are thickness weighted.
!
! Some diagnostics related to fluid mass transport
! classified according to both depth and density classes
! are also computed.
!
!
!
!
!
!
! 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: Elements of MOM (2012)
!
!
!
! The expressions for the horizontal components for tracer advection
! allow for a proper conversion between pressure work and buoyancy.
!
!
!
! The B-grid remapping operators are derived from considerations of linear
! interpolation and volume conservation. A "remapping error" is
! computed to determine consistency between the tracer and velocity
! grid advection velocities. This error is roundoff only for cases
! where the horizontal tracer and velocity grids are linearly related,
! as is the case for the spherical coordinate version of MOM. The
! tripolar version of MOM does not have tracer and velocity grids
! related linearly, and so the "remapping error" is nontrivial. The
! significance of this error is unclear. No adverse effects have been
! identified.
!
!
!
! The vertical velocity components for both the tracer and velocity cells
! are diagnosed via continuity (either volume or mass conservation depending
! on the use of the Boussinesq or non-Boussinesq versions of MOM).
!
!
!
!
!
!
!
! For debugging.
!
!
!
! For adding an inflow velocity from the northern boundary.
! Default is inflow_nboundary=.false.
!
!
!
! For reading in a file with specified zonal, meridional,
! and vertical components to the advective velocity.
! The file should have velocity at the east face of T-cell,
! north face, and bottom, just as on a C-grid. The units
! should be m/s for each component. MOM then multiplies
! but the appropriate thickness and density factors to
! generate transport for use in the model.
! Default read_advection_velocity=.false.
!
!
!
! For reading in a file with specified zonal, meridional,
! and vertical components to the advection transport.
! The file should have transport at the east face of T-cell,
! north face, and bottom, just as on a C-grid. The units
! should be (kg/m^3)*(m^2/sec) for horz components, and
! (kg/m^3)*(m/sec) for vertical component. MOM then uses
! these fields directly to initialize uhrho_et, vhrho_nt,
! and wrho_bt.
! Default read_advection_transport=.false.
!
!
!
! When reading in the advection velocity components, we
! may choose to keep them constant in time. This facilitates
! idealized tests of tracer advection.
! Default constant_advection_velocity=.false.
!
!
!
! This is a check value used to determine if the time steps will result in
! linearly stable advection. If set to a number < 0, then model will estimate the
! value as a function of maximum grid size.
! Note that this time step check is not rigorous, and it depends on the details of
! the advection scheme. Nonetheless, it provides some useful warning for setting the
! time steps in the model.
!
!
!
use axis_utils_mod, only: nearest_index
use constants_mod, only: epsln, c2dbars
use diag_manager_mod, only: send_data, register_diag_field, register_static_field
use fms_mod, only: write_version_number, open_namelist_file
use fms_mod, only: close_file, check_nml_error
use fms_mod, only: read_data, file_exist
use fms_mod, only: FATAL, WARNING
use fms_io_mod, only: register_restart_field, save_restart
use fms_io_mod, only: restore_state, restart_file_type
use mpp_domains_mod, only: mpp_update_domains
use mpp_mod, only: input_nml_file, mpp_error, stdout, stdlog
use mpp_mod, only: mpp_min, mpp_max, mpp_error, mpp_pe
use ocean_domains_mod, only: get_local_indices
use ocean_operators_mod, only: BAX, BAY, BDX_ET, BDY_NT, BDX_EU, BDY_NU
use ocean_operators_mod, only: REMAP_ET_TO_EU, REMAP_NT_TO_NU, REMAP_BT_TO_BU
use ocean_parameters_mod,only: missing_value, grav, rho0, rho0r
use ocean_parameters_mod,only: MOM_BGRID, MOM_CGRID, DEPTH_BASED
use ocean_types_mod, only: ocean_domain_type, ocean_grid_type
use ocean_types_mod, only: ocean_grid_type, ocean_thickness_type
use ocean_types_mod, only: ocean_adv_vel_type, ocean_velocity_type
use ocean_types_mod, only: ocean_time_type, ocean_time_steps_type
use ocean_types_mod, only: ocean_density_type, ocean_lagrangian_type
use ocean_workspace_mod, only: wrk1, wrk1_2d, wrk2_2d
use ocean_obc_mod, only: ocean_obc_adjust_advel
use ocean_util_mod, only: write_timestamp, diagnose_2d, diagnose_3d, diagnose_3d_u, diagnose_sum
use ocean_util_mod, only: write_chksum_3d, write_chksum_2d
implicit none
private
integer :: unit = 6
! for vertical coordinate class
integer :: vert_coordinate_class
! for Bgrid or Cgrid
integer :: horz_grid
! for diagnostics
integer :: id_wt =-1
integer :: id_wrho_bt =-1
integer :: id_uhrho_et =-1
integer :: id_vhrho_nt =-1
integer :: id_wrho_bu =-1
integer :: id_uhrho_eu =-1
integer :: id_vhrho_nu =-1
integer :: id_rhodz_vorticity_z=-1
integer :: id_mass_source =-1
integer :: id_horz_diverge_t =-1
integer :: id_horz_diverge_u =-1
integer :: id_courant_uet =-1
integer :: id_courant_vnt =-1
integer :: id_courant_wbt =-1
integer :: id_courant_ueu =-1
integer :: id_courant_vnu =-1
integer :: id_courant_wbu =-1
integer :: id_ue =-1
integer :: id_vn =-1
integer :: id_wb =-1
integer :: id_ue_rhodzt =-1
integer :: id_vn_rhodzt =-1
integer :: id_wb_rho =-1
integer :: id_eta_tend_press =-1
integer :: id_eta_tend_press_glob=-1
real :: dtts
real :: dtuv
real :: cellarea_r
logical :: used
!for restart file
type(restart_file_type), save :: adv_vel_restart
integer :: id_adv_vel_restart
#include
#ifdef MOM_STATIC_ARRAYS
real, dimension(isd:ied,jsd:jed) :: tmp
#else
real, dimension(:,:), allocatable :: tmp
#endif
real, dimension(:,:,:), allocatable :: vt_inflow ! m/s
real, dimension(:,:,:), allocatable :: mask_inflow
real, dimension(:,:,:), allocatable :: ue !m/s from file
real, dimension(:,:,:), allocatable :: vn !m/s from file
real, dimension(:,:,:), allocatable :: wb !m/s from file
real, dimension(:,:), allocatable :: cgrav2
type(ocean_domain_type), pointer :: Dom =>NULL()
type(ocean_grid_type), pointer :: Grd =>NULL()
public ocean_advection_velocity_init
public ocean_advection_velocity
public ocean_advection_velocity_end
private check_vert_cfl_blobs
private read_advection
private inflow_nboundary_init
character(len=128) :: version=&
'$Id: ocean_advection_velocity.F90,v 20.0.4.1 2014/01/22 17:19:10 smg Exp $'
character (len=128) :: tagname = &
'$Name: mom5_tikal_22jan2014_smg $'
logical :: have_obc = .false.
logical :: module_is_initialized = .FALSE.
logical :: prescribe_advection_velocity = .false.
real :: max_advection_velocity = -1.0
logical :: debug_this_module = .false.
logical :: inflow_nboundary = .false.
logical :: read_advection_velocity = .false.
logical :: read_advection_transport = .false.
logical :: constant_advection_velocity = .false.
namelist /ocean_advection_velocity_nml/ debug_this_module, max_advection_velocity, inflow_nboundary, &
read_advection_velocity, read_advection_transport, &
constant_advection_velocity
contains
!#######################################################################
!
!
!
! Initialize the advection velocity module
!
!
subroutine ocean_advection_velocity_init(Grid, Domain, Time, Time_steps, Thickness, Adv_vel, &
ver_coordinate_class, hor_grid, obc, use_blobs, &
introduce_blobs, debug)
type(ocean_grid_type), intent(in), target :: Grid
type(ocean_domain_type), intent(in), target :: Domain
type(ocean_time_type), intent(in), target :: Time
type(ocean_time_steps_type), intent(in) :: Time_steps
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_adv_vel_type), intent(inout) :: Adv_vel
integer, intent(in) :: ver_coordinate_class
integer, intent(in) :: hor_grid
logical, intent(in) :: obc
logical, intent(in) :: use_blobs
logical, intent(in) :: introduce_blobs
logical, intent(in), optional :: debug
character(len=128) :: filename
integer :: ioun, io_status, ierr
real :: gridmin, gridmax, dtadv, a, b
integer :: i, j, icg, jcg
real :: vertical_factor = 2.0
real :: max_dt_for_advection
real :: max_dt_for_advection0
real, dimension(7) :: res = (/0.0625, 0.125, 0.25, 0.50, 1.00, 2.0, 4.0/) ! max resolution (deg)
real, dimension(7) :: vel = (/2.0, 1.8, 1.6, 1.45, 1.25, 0.9, 0.4/) ! estimated max advection velocity (m/sec)
integer :: stdoutunit,stdlogunit
stdoutunit=stdout();stdlogunit=stdlog()
if ( module_is_initialized ) then
call mpp_error(FATAL,&
'==>Error from ocean_advection_velocity_mod: module already initialized.')
endif
call write_version_number(version, tagname)
module_is_initialized = .TRUE.
have_obc = obc
vert_coordinate_class = ver_coordinate_class
! provide for namelist over-ride of defaults
#ifdef INTERNAL_FILE_NML
read (input_nml_file, nml=ocean_advection_velocity_nml, iostat=io_status)
ierr = check_nml_error(io_status,'ocean_advection_velocity_nml')
#else
ioun = open_namelist_file ()
read (ioun, ocean_advection_velocity_nml,iostat=io_status)
ierr = check_nml_error(io_status,'ocean_advection_velocity_nml')
call close_file(ioun)
#endif
write (stdoutunit,'(/)')
write (stdoutunit, ocean_advection_velocity_nml)
write (stdlogunit, ocean_advection_velocity_nml)
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_advection_velocity_mod with debug_this_module=.true.'
endif
Dom => Domain
Grd => Grid
dtts = Time_steps%dtts
dtuv = Time_steps%dtuv
cellarea_r = 1.0/(epsln + Grd%tcellsurf)
horz_grid = hor_grid
#ifndef MOM_STATIC_ARRAYS
call get_local_indices(Dom, isd, ied, jsd, jed, isc, iec, jsc, jec)
nk = Grd%nk
! allocate variables
allocate ( Adv_vel%uhrho_et(isd:ied,jsd:jed,nk) )
allocate ( Adv_vel%vhrho_nt(isd:ied,jsd:jed,nk) )
allocate ( Adv_vel%uhrho_eu(isd:ied,jsd:jed,nk) )
allocate ( Adv_vel%vhrho_nu(isd:ied,jsd:jed,nk) )
allocate ( Adv_vel%wrho_bt(isd:ied,jsd:jed,0:nk) )
allocate ( Adv_vel%wrho_bu(isd:ied,jsd:jed,0:nk) )
allocate ( Adv_vel%diverge_t(isd:ied,jsd:jed,nk) )
allocate ( Adv_vel%diverge_u(isd:ied,jsd:jed,nk) )
allocate ( tmp(isd:ied,jsd:jed) )
#endif
allocate ( cgrav2(isd:ied,jsd:jed) )
Adv_vel%uhrho_et(:,:,:) = 0.0
Adv_vel%vhrho_nt(:,:,:) = 0.0
Adv_vel%wrho_bt(:,:,:) = 0.0
Adv_vel%uhrho_eu(:,:,:) = 0.0
Adv_vel%vhrho_nu(:,:,:) = 0.0
Adv_vel%wrho_bu(:,:,:) = 0.0
Adv_vel%diverge_t(:,:,:) = 0.0
Adv_vel%diverge_u(:,:,:) = 0.0
tmp(:,:) = 0.0
cgrav2(:,:) = grav*Grd%ht(:,:)*Grd%tmask(:,:,1)
if(inflow_nboundary) then
call inflow_nboundary_init
endif
if(read_advection_velocity .or. read_advection_transport) then
call read_advection(Time, Thickness, Adv_vel)
prescribe_advection_velocity = .true.
endif
if(constant_advection_velocity) then
write(stdoutunit,'(a)') &
'==>Note: Using ocean_advection_velocity_mod with constant_advection_velocity=.true.'
write(stdoutunit,'(a)') &
' Will hold the advection velocity components constant in time. '
prescribe_advection_velocity = .true.
endif
if(use_blobs .and. horz_grid == MOM_CGRID) then
write(stdoutunit,*) &
'==>Warning from advection velocity module: MOM C-grid with blobs is incomplete. Blobs have only been tested with Bgrid.'
endif
! initialise the restart if blobs are used
if (use_blobs) then
filename = 'ocean_adv_vel.res.nc'
id_adv_vel_restart = register_restart_field(adv_vel_restart, filename, 'wrho_bt', &
Adv_vel%wrho_bt(:,:,:), domain=Dom%domain2d)
if (.not. introduce_blobs) then
if(file_exist('INPUT/'//trim(filename))) then
write(stdoutunit, '(/a)') 'Reading in advective velocity data from '//trim(filename)
call restore_state( adv_vel_restart, id_adv_vel_restart )
else
if (.NOT.Time%init) then
call mpp_error(FATAL,&
'Expecting file '//trim(filename)//' to exist.&
&This file was not found and Time%init=.false.')
endif
endif
endif
endif
! estimate the max advective speed for the given resolution
if (max_advection_velocity < 0.0) then
gridmax = 0.0
do j=jsc,jec
do i=isc,iec
if (Grd%kmt(i,j) /= 0) then
gridmax = max(gridmax,Grd%dyt(i,j),Grd%dxt(i,j))
endif
enddo
enddo
call mpp_max (gridmax)
gridmax = gridmax/111.324e3 ! convert to degrees from metres
i = nearest_index(gridmax, res)
if (gridmax <= res(1)) then
max_advection_velocity = vel(1)
elseif (gridmax >= res(7)) then
max_advection_velocity = vel(7)
elseif (gridmax < res(i)) then
a = res(i)-gridmax
b = gridmax-res(i-1)
max_advection_velocity = (a*vel(i-1) + b*vel(i))/(a+b)
elseif (gridmax > res(i)) then
a = gridmax-res(i)
b = res(i+1)-gridmax
max_advection_velocity = (a*vel(i+1) + b*vel(i))/(a+b)
endif
write (stdoutunit,'(/a,f5.2,a/a/)') &
' Note: The max_advection_velocity is estimated to be ',max_advection_velocity,&
' m/s based on grid resolution.',' A more appropriate value can be specified via namelist.'
endif
icg = isc; jcg = jsc; gridmin = 1.0e20; max_dt_for_advection = 1.e6; gridmax = 0.0
do j=jsc,jec
do i=isc,iec
if (Grd%kmt(i,j) /= 0) then
gridmin = min(gridmin,Grd%dyt(i,j),Grd%dxt(i,j))
gridmax = max(gridmax,Grd%dyt(i,j),Grd%dxt(i,j))
dtadv = 0.5*gridmin/max_advection_velocity
if (dtadv < max_dt_for_advection) then
max_dt_for_advection = dtadv; icg = i; jcg = j
endif
endif
enddo
enddo
max_dt_for_advection = nint(max_dt_for_advection/vertical_factor) !
max_dt_for_advection = max_dt_for_advection + 1.e-12*mpp_pe() ! to separate redundancies
max_dt_for_advection0 = max_dt_for_advection
call mpp_min (max_dt_for_advection)
call mpp_max (gridmax)
if (max_dt_for_advection == max_dt_for_advection0) then
if (max(dtts,dtuv) <= max_dt_for_advection) then
write (unit,'(/a,i4,a,i4,a,f8.3,a,f8.3,a)') &
' Note: Advection stability most nearly violated at T-cell (i,j) = (',&
icg+Dom%ioff,',',jcg+Dom%joff,'), (lon,lat) = (',Grd%xt(icg,jcg),',',Grd%yt(icg,jcg),').'
else
write (unit,'(/a,i4,a,i4,a,f8.3,a,f8.3,a)') &
'=>Error: Advection stability violated at T-cell (i,j) = (',&
icg+Dom%ioff,',',jcg+Dom%joff,'), (lon,lat) = (',Grd%xt(icg,jcg),',',Grd%yt(icg,jcg),').'
call mpp_error(FATAL,&
'==>Error from ocean_advection_velocity_mod: tracer advection stability is violated.')
endif
write (unit,'(a,f7.2,a,/a,f12.2,a,/a,f12.2,a,f12.2,a)') &
' Assuming a maximum advection velocity of ',max_advection_velocity,' m/s,', &
' Linear stability requires max(dtuv,dtts) be less than ',max_dt_for_advection,' sec.', &
' Model is now using (dtuv,dtts) = (',dtuv,',',dtts,') sec.'
endif
max_dt_for_advection = nint(max_dt_for_advection)
! register advective velocity components for diagnostic output
id_horz_diverge_t = register_diag_field ('ocean_model', 'horz_diverge_t', &
Grd%vel_axes_wt(1:3), Time%model_time, &
'horizontal (on k=constant) divergence of (uhrho,vhrho) at T-points', '(kg/m^3)*m/sec', &
missing_value=missing_value, range=(/-10.e4,10.e4/))
id_wt = register_diag_field ('ocean_model', 'wt', Grd%vel_axes_wt(1:3), Time%model_time, &
'dia-surface velocity T-points', 'm/sec', missing_value=missing_value, range=(/-10.e4,10.e4/))
id_wrho_bt = register_diag_field ('ocean_model', 'wrhot', Grd%vel_axes_wt(1:3), Time%model_time, &
'rho*dia-surface velocity T-points', '(kg/m^3)*m/sec', missing_value=missing_value, range=(/-10.e4,10.e4/))
id_uhrho_et = register_diag_field ('ocean_model', 'uhrho_et', Grd%tracer_axes(1:3), Time%model_time, &
'uhrho_et on horz face of T-cells', '(kg/m^3)*m^2/sec', missing_value=missing_value, range=(/-100.e4,100.e4/))
id_vhrho_nt = register_diag_field ('ocean_model', 'vhrho_nt', Grd%tracer_axes(1:3), Time%model_time, &
'vhrho_nt on horz face of T-cells', '(kg/m^3)*m^2/sec', missing_value=missing_value, range=(/-100.e4,100.e4/))
id_horz_diverge_u = register_diag_field ('ocean_model', 'horz_diverge_u', &
Grd%vel_axes_wu(1:3), Time%model_time, &
'horizontal (on k=constant) divergence of (uhrho,vhrho) at U-points (only relevant for B-grid)',&
'(kg/m^3)*m/sec', missing_value=missing_value, range=(/-10.e4,10.e4/))
id_wrho_bu = register_diag_field ('ocean_model', 'wrhou', Grd%vel_axes_wu(1:3), Time%model_time, &
'rho*dia-surface velocity U-points (only relevant for B-grid)', '(kg/m^3)*m/sec', &
missing_value=missing_value, range=(/-10.e4,10.e4/))
id_uhrho_eu = register_diag_field ('ocean_model', 'uhrho_eu', Grd%vel_axes_uv(1:3), Time%model_time, &
'uhrho_eu on horz face of U-cells (only relevant for B-grid)', '(kg/m^3)*m^2/sec', &
missing_value=missing_value, range=(/-100.e4,100.e4/))
id_vhrho_nu = register_diag_field ('ocean_model', 'vhrho_nu', Grd%vel_axes_uv(1:3), Time%model_time, &
'vhrho_nu on horz face of U-cells (only relevant for B-grid)', '(kg/m^3)*m^2/sec', &
missing_value=missing_value, range=(/-100.e4,100.e4/))
! mass-based vorticity
id_rhodz_vorticity_z = register_diag_field ('ocean_model', 'rhodz_vorticity_z', Grd%vel_axes_uv(1:3), &
Time%model_time, 'vertical rhodz vorticity: (rho*dz*v)_x-(rho*dz*u)_y', &
'(kg*m^2)/sec', missing_value=missing_value, range=(/-1e8,1e8/))
! register Courant numbers for diagnostic output
id_courant_uet = register_diag_field ('ocean_model', 'courant_uet', Grd%tracer_axes(1:3), Time%model_time, &
'Courant number [uet*dtts/dxt]', 'none', missing_value=missing_value, range=(/-1e9,1e9/))
id_courant_vnt = register_diag_field ('ocean_model', 'courant_vnt', Grd%tracer_axes(1:3), Time%model_time, &
'Courant number [vnt*dtts/dyt]', 'none', missing_value=missing_value, range=(/-1e9,1e9/))
id_courant_wbt = register_diag_field ('ocean_model', 'courant_wbt', Grd%tracer_axes(1:3), Time%model_time, &
'Courant number [wbt*dtts/dzt]', 'none', missing_value=missing_value, range=(/-1e9,1e9/))
id_courant_ueu = register_diag_field ('ocean_model', 'courant_ueu', Grd%vel_axes_uv(1:3), Time%model_time, &
'Courant number [ueu*dtuv/dxu]', 'none', missing_value=missing_value, range=(/-1e9,1e9/))
id_courant_vnu = register_diag_field ('ocean_model', 'courant_vnu', Grd%vel_axes_uv(1:3), Time%model_time, &
'Courant number [vnu*dtuv/dyu]', 'none', missing_value=missing_value, range=(/-1e9,1e9/))
id_courant_wbu = register_diag_field ('ocean_model', 'courant_wbu', Grd%vel_axes_uv(1:3), Time%model_time, &
'Courant number [wnu*dtuv/dzu]', 'none', missing_value=missing_value, range=(/-1e9,1e9/))
! register mass source
id_mass_source = register_diag_field ('ocean_model', 'mass_source', Grd%tracer_axes(1:3), Time%model_time, &
'Mass source on T-grid', '(kg/m^3)*(m/s)', missing_value=missing_value, range=(/-1e9,1e9/))
! eta tendency from motion across pressure surfaces (when vert_coord=press)
id_eta_tend_press = register_diag_field ('ocean_model', 'eta_tend_press', &
Grd%tracer_axes(1:2),Time%model_time, &
'non-Bouss steric tendency from motion across pressure surfaces (when vert coord = press)', &
'm/s', missing_value=missing_value, range=(/-1e9,1e9/))
id_eta_tend_press_glob = register_diag_field ('ocean_model', 'eta_tend_press_glob', &
Time%model_time, &
'non-Bouss steric tendency from motion across pressure surfaces (when vert coord = press)', &
'm/s', missing_value=missing_value, range=(/-1e9,1e9/))
if(debug_this_module .and. use_blobs) then
write(stdoutunit,*) ' '
write(stdoutunit,*) 'From ocean_advection_velocity_mod: beginning chksums'
call write_timestamp(Time%model_time)
call write_chksum_3d('Adv_vel%wrho_bt', Adv_vel%wrho_bt(COMP,1:))
call write_chksum_2d('Adv_vel%wrho_bt(0)', Adv_vel%wrho_bt(COMP,0))
endif
end subroutine ocean_advection_velocity_init
! NAME="ocean_advection_velocity_init"
!#######################################################################
!
!
!
! Compute thickness weighted and density weighted advection velocity
! components for the B-grid on the T-cells and U-cells.
!
! U-cell advective components are not used for Cgrid, so we do not
! update the u-cell components. Instead, simply keep values at their
! initial setting of 0.0.
!
!
!
subroutine ocean_advection_velocity (Velocity, Time, Thickness, Dens, pme, river, Adv_vel, &
L_system, use_blobs)
type(ocean_velocity_type), intent(in) :: Velocity
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_density_type), intent(in) :: Dens
real, dimension(isd:,jsd:), intent(in) :: pme
real, dimension(isd:,jsd:), intent(in) :: river
type(ocean_adv_vel_type), intent(inout) :: Adv_vel
type(ocean_lagrangian_type),intent(in) :: L_system
logical, intent(in) :: use_blobs
integer :: i, j, k, tau
integer :: stdoutunit
stdoutunit=stdout()
if ( .not. module_is_initialized ) &
call mpp_error(FATAL, &
'==>Error from ocean_advection_velocity_mod (ocean_advection_velocity): module must be initialized')
tau = Time%tau
wrk1 = 0.0
tmp = 0.0
! compute thickness weighted and density weighted C-grid advection
! velocity components (uhrho_et, vhrho_nt) to advect tracer fields.
! Remap these to get the corresponding advection velocity components
! (uhrho_eu,vhrho_nu) to advect B-grid horizontal momentum.
if(.not. prescribe_advection_velocity .and. horz_grid == MOM_BGRID) then
do k=1,nk
Adv_vel%uhrho_et(:,:,k) = &
BAY(Velocity%u(:,:,k,1,tau)*Grd%dyu(:,:)*Thickness%rho_dzu(:,:,k,tau))/Grd%dyte(:,:)
Adv_vel%vhrho_nt(:,:,k) = &
BAX(Velocity%u(:,:,k,2,tau)*Grd%dxu(:,:)*Thickness%rho_dzu(:,:,k,tau))/Grd%dxtn(:,:)
Adv_vel%uhrho_eu(:,:,k) = REMAP_ET_TO_EU(Adv_vel%uhrho_et(:,:,k))
Adv_vel%vhrho_nu(:,:,k) = REMAP_NT_TO_NU(Adv_vel%vhrho_nt(:,:,k))
Adv_vel%diverge_t(:,:,k)= Grd%tmask(:,:,k) &
*(BDX_ET(Adv_vel%uhrho_et(:,:,k)) + BDY_NT(Adv_vel%vhrho_nt(:,:,k)))
Adv_vel%diverge_u(:,:,k)= Grd%umask(:,:,k) &
*(BDX_EU(Adv_vel%uhrho_eu(:,:,k)) + BDY_NU(Adv_vel%vhrho_nu(:,:,k)))
enddo
! dia-surface mass flux through ocean surface arises from water crossing ocean surface.
! minus sign arises from the model sign convention.
Adv_vel%wrho_bt(:,:,0) = -(pme(:,:) + river(:,:))
Adv_vel%wrho_bu(:,:,0) = Grd%umask(:,:,1)*REMAP_BT_TO_BU(Adv_vel%wrho_bt(:,:,0))
if (use_blobs) then
do k=1,nk
tmp(:,:) = Thickness%rho_dzt_tendency(:,:,k) - Thickness%mass_source(:,:,k) - L_system%conv_blob(:,:,k)
Adv_vel%wrho_bt(:,:,k) = (tmp(:,:) + Adv_vel%diverge_t(:,:,k) + Adv_vel%wrho_bt(:,:,k-1)) &
*Grd%tmask(:,:,k)
tmp(:,:) = Grd%umask(:,:,k)*REMAP_BT_TO_BU(tmp(:,:))
Adv_vel%wrho_bu(:,:,k) = tmp(:,:) + Adv_vel%diverge_u(:,:,k) + Adv_vel%wrho_bu(:,:,k-1)
enddo
if (debug_this_module) then
call check_vert_cfl_blobs(Time, Adv_vel, Thickness, L_system)
endif
else !use_blobs=.false.
do k=1,nk
tmp(:,:) = Thickness%rho_dzt_tendency(:,:,k) - Thickness%mass_source(:,:,k)
Adv_vel%wrho_bt(:,:,k) = (tmp(:,:) + Adv_vel%diverge_t(:,:,k) + Adv_vel%wrho_bt(:,:,k-1)) &
*Grd%tmask(:,:,k)
tmp(:,:) = Grd%umask(:,:,k)*REMAP_BT_TO_BU(tmp(:,:))
Adv_vel%wrho_bu(:,:,k) = tmp(:,:) + Adv_vel%diverge_u(:,:,k) + Adv_vel%wrho_bu(:,:,k-1)
enddo
endif
endif ! endif for .not. prescribe_advection_velocity and MOM_BGRID
! smg: June 2012: C-grid yet to be updated for blobs.
if(.not. prescribe_advection_velocity .and. horz_grid == MOM_CGRID) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
Adv_vel%uhrho_et(i,j,k) = Velocity%u(i,j,k,1,tau)*Thickness%rho_dzten(i,j,k,1) &
*Grd%dxt_dxter(i,j)*Grd%dyt_dyter(i,j)
Adv_vel%vhrho_nt(i,j,k) = Velocity%u(i,j,k,2,tau)*Thickness%rho_dzten(i,j,k,2) &
*Grd%dxt_dxtnr(i,j)*Grd%dyt_dytnr(i,j)
enddo
enddo
Adv_vel%diverge_t(:,:,k) = Grd%tmask(:,:,k)*(BDX_ET(Adv_vel%uhrho_et(:,:,k)) + BDY_NT(Adv_vel%vhrho_nt(:,:,k)))
enddo
! dia-surface mass flux through ocean surface arises from water crossing ocean surface.
! minus sign arises from the model sign convention.
Adv_vel%wrho_bt(:,:,0) = -(pme(:,:) + river(:,:))
do k=1,nk
tmp(:,:) = Thickness%rho_dzt_tendency(:,:,k) - Thickness%mass_source(:,:,k)
Adv_vel%wrho_bt(:,:,k) = (tmp(:,:) + Adv_vel%diverge_t(:,:,k) + Adv_vel%wrho_bt(:,:,k-1)) &
*Grd%tmask(:,:,k)
enddo
endif ! endif for .not. prescribe_advection_velocity and MOM_CGRID
if (have_obc) then
call ocean_obc_adjust_advel(Adv_vel)
endif
! specify a value for the inflow meridional transport
if(inflow_nboundary) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
Adv_vel%vhrho_nt(i,j,k) = Adv_vel%vhrho_nt(i,j,k) &
+ Grd%tmask(i,j,k)*mask_inflow(i,j,k)*rho0*Thickness%dzt(i,j,k)*vt_inflow(i,j,k)
enddo
enddo
enddo
endif
if(debug_this_module) then
write(stdoutunit,*) ' '
write(stdoutunit,*) 'From ocean_advection_velocity_mod: chksums'
call write_timestamp(Time%model_time)
call write_chksum_3d('Velocity%u(1)', Velocity%u(COMP,:,1,tau))
call write_chksum_3d('Velocity%u(2)', Velocity%u(COMP,:,2,tau))
! To make it equivalent to the checksum in the thickness module
call write_chksum_3d('rho_dzu', Thickness%rho_dzu(COMP,:,tau)*Grd%umask(COMP,:))
if (use_blobs) then
call write_chksum_3d('L_system%conv_blob', L_system%conv_blob(COMP,:))
endif
call write_chksum_3d('Adv_vel%uhrho_et', Adv_vel%uhrho_et(COMP,:))
call write_chksum_3d('Adv_vel%vhrho_nt', Adv_vel%vhrho_nt(COMP,:))
call write_chksum_3d('Adv_vel%uhrho_eu', Adv_vel%uhrho_eu(COMP,:))
call write_chksum_3d('Adv_vel%vhrho_nu', Adv_vel%vhrho_nu(COMP,:))
call write_chksum_3d('Adv_vel%wrho_bt', Adv_vel%wrho_bt(COMP,1:))
call write_chksum_3d('Adv_vel%wrho_bu', Adv_vel%wrho_bu(COMP,1:))
call write_chksum_2d('Adv_vel%wrho_bt(0)', Adv_vel%wrho_bt(COMP,0))
call write_chksum_2d('Adv_vel%wrho_bu(0)', Adv_vel%wrho_bu(COMP,0))
call write_chksum_3d('Adv_vel%diverge_t', Adv_vel%diverge_t(COMP,:))
call write_chksum_3d('Adv_vel%diverge_u', Adv_vel%diverge_u(COMP,:))
endif
! send advective velocity components to diagnostic manager
if (id_wt > 0) then
if(vert_coordinate_class==DEPTH_BASED) then
call diagnose_3d(Time, Grd, id_wt, rho0r*Adv_vel%wrho_bt(:,:,1:nk))
else
wrk1(:,:,:) = 0.0
do k=1,nk
do j=jsd,jed
do i=isd,ied
if(Grd%tmask(i,j,k) > 0.0) then
wrk1(i,j,k) = Adv_vel%wrho_bt(i,j,k)/Dens%rho(i,j,k,tau)
endif
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_wt, wrk1(:,:,:))
endif
endif
call diagnose_3d(Time, Grd, id_uhrho_et, Adv_vel%uhrho_et(:,:,:))
call diagnose_3d(Time, Grd, id_vhrho_nt, Adv_vel%vhrho_nt(:,:,:))
call diagnose_3d(Time, Grd, id_wrho_bt, Adv_vel%wrho_bt(:,:,1:nk))
call diagnose_3d_u(Time, Grd, id_uhrho_eu, Adv_vel%uhrho_eu(:,:,:))
call diagnose_3d_u(Time, Grd, id_vhrho_nu, Adv_vel%vhrho_nu(:,:,:))
call diagnose_3d_u(Time, Grd, id_wrho_bu, Adv_vel%wrho_bu(:,:,1:nk))
call diagnose_3d(Time, Grd, id_horz_diverge_t, Adv_vel%diverge_t(:,:,:))
call diagnose_3d_u(Time, Grd, id_horz_diverge_u, Adv_vel%diverge_u(:,:,:))
! send mass and thickness weighted vertical vorticity
if(id_rhodz_vorticity_z > 0) then
wrk1=0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = (Adv_vel%vhrho_nt(i+1,j,k)-Adv_vel%vhrho_nt(i,j,k))*Grd%dxur(i,j) &
-(Adv_vel%uhrho_et(i,j+1,k)-Adv_vel%uhrho_et(i,j,k))*Grd%dyur(i,j)
enddo
enddo
enddo
call diagnose_3d_u(Time, Grd, id_rhodz_vorticity_z, wrk1(:,:,:))
endif
! send Courant numbers to diagnostic manager
if (id_courant_uet > 0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = dtts*Adv_vel%uhrho_et(i,j,k)/(Thickness%rho_dzt(i,j,k,tau)*Grd%dxt(i,j))
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_courant_uet, wrk1(:,:,:))
endif
if (id_courant_vnt > 0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = dtts*Adv_vel%vhrho_nt(i,j,k)/(Thickness%rho_dzt(i,j,k,tau)*Grd%dyt(i,j))
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_courant_vnt, wrk1(:,:,:))
endif
if (id_courant_wbt > 0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = dtts*Adv_vel%wrho_bt(i,j,k)/Thickness%rho_dzt(i,j,k,tau)
enddo
enddo
enddo
call diagnose_3d(Time, Grd, id_courant_wbt, wrk1(:,:,:))
endif
if (id_courant_ueu > 0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = dtuv*Adv_vel%uhrho_eu(i,j,k)/(Thickness%rho_dzu(i,j,k,tau)*Grd%dxu(i,j))
enddo
enddo
enddo
call diagnose_3d_u(Time, Grd, id_courant_ueu, wrk1(:,:,:))
endif
if (id_courant_vnu > 0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = dtuv*Adv_vel%vhrho_nu(i,j,k)/(Thickness%rho_dzu(i,j,k,tau)*Grd%dyu(i,j))
enddo
enddo
enddo
call diagnose_3d_u(Time, Grd, id_courant_vnu, wrk1(:,:,:))
endif
if (id_courant_wbu > 0) then
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1(i,j,k) = dtuv*Adv_vel%wrho_bu(i,j,k)/Thickness%rho_dzu(i,j,k,tau)
enddo
enddo
enddo
call diagnose_3d_u(Time, Grd, id_courant_wbu, wrk1(:,:,:))
endif
if(id_eta_tend_press > 0 .or. id_eta_tend_press_glob > 0) then
! do not include k=0, as that is pme+river contribution,
! and that is not what is desired for this diagnostic.
! note that drhodP is in units of (kg/m^3)/dbars, so
! need to multiply by c2dbars to get squared sound speed
! in (m^2/s^2).
wrk1_2d(:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = wrk1_2d(i,j) + Dens%drhodP(i,j,k)*Adv_vel%wrho_bt(i,j,k) &
*Thickness%dzt(i,j,k)/(Dens%rho(i,j,k,tau)+epsln)
enddo
enddo
enddo
do j=jsc,jec
do i=isc,iec
wrk1_2d(i,j) = c2dbars*grav*wrk1_2d(i,j)
enddo
enddo
call diagnose_2d(Time, Grd, id_eta_tend_press, wrk1_2d(:,:))
call diagnose_sum(Time, Grd, Dom, id_eta_tend_press_glob, wrk1_2d, cellarea_r)
endif
end subroutine ocean_advection_velocity
! NAME="advection_velocity"
!#######################################################################
!
!
!
! This routine is useful to check for strong vertical motions in the
! presence of Lagrangian blobs. The CFL check is the same as in the
! module ocean_adv_vel_diag.F90, but we have more local control here
! so it is useful in the development stages of the blobs.
!
!
subroutine check_vert_cfl_blobs(Time, Adv_vel, Thickness, L_system)
type(ocean_time_type), intent(in) :: Time
type(ocean_adv_vel_type), intent(in) :: Adv_vel
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_lagrangian_type), intent(in) :: L_system
integer :: i,j,k,tau
integer :: stdoutunit
stdoutunit=stdout()
tau = Time%tau
do k=1,nk
do j=jsc,jec
do i=isc,iec
if (dtts*Adv_vel%wrho_bt(i,j,k)/Thickness%rho_dztT(i,j,k,tau) > 1.0) then
write (stdoutunit, '(/,a)') 'Vertical CFL condition violated'
call write_timestamp(Time%model_time)
print ('(3(a,i3),a,x,7(/,a20,x,es11.4))'), '(',i,',',j,',',k,')', &
'CFL number =', dtts*Adv_vel%wrho_bt(i,j,k)/Thickness%rho_dztT(i,j,k,tau), &
'wrho_bt =', Adv_vel%wrho_bt(i,j,k), &
'w =', rho0r*Adv_vel%wrho_bt(i,j,k), &
'rho_dzt_tendency =', Thickness%rho_dzt_tendency(i,j,k), &
'mass_source =', -Thickness%mass_source(i,j,k), &
'conv_blob =', -L_system%conv_blob(i,j,k), &
'diverge_t =', Adv_vel%diverge_t(i,j,k)
call mpp_error(WARNING, &
'==>Error from ocean_advection_velocity_mod: Vertical CFL condition violated in presence of blobs.')
endif
enddo
enddo
enddo
end subroutine check_vert_cfl_blobs
! NAME="check_vert_cfl_blobs"
!#######################################################################
!
!
!
! To read in advection velocity components or advection transport
! components.
!
! This routine assumes that the read-in components are masked.
!
! A/ read_advection_velocity=.true.
! Assume (x,y,z) advection velocity components read in and are in units
! of meter/sec, and placed on the T-cell faces, as in a C-grid ocean model.
! An application of this option is for developing idealized test cases
! for tracer advection.
!
! B/ read_advection_transport=.true.
! Assume (x,y,z) advection transport components read in from a file.
! The file should have transport at the east face of T-cell,
! north face, and bottom, just as on a C-grid. The units should be
! uhrho_et = (kg/m^3)*(m^2/sec)
! vhrho_nt = (kg/m^3)*(m^2/sec)
! wrho_bt = (kg/m^3)*(m/sec)
! MOM then uses these read fields to initialize uhrho_et, vhrho_nt,
! and wrho_bt.
! An application of this option is for diagnosing the advection
! operator and advection fluxes when reading in the time mean
! uhrho_et, vhrho_nt, and wrho_bt and initializing tracer concentration
! with a time mean value. Assuming we have also saved the time mean
! advection operator, we may then diagnose the eddy correlation
! contribution to tracer transport.
!
!
!
subroutine read_advection(Time, Thickness, Adv_vel)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_adv_vel_type), intent(inout) :: Adv_vel
integer :: i,j,k
character(len=128) :: filename
allocate (ue(isd:ied,jsd:jed,nk))
ue(:,:,:) = 0.0
allocate (vn(isd:ied,jsd:jed,nk))
vn(:,:,:) = 0.0
allocate (wb(isd:ied,jsd:jed,nk))
wb(:,:,:) = 0.0
filename = 'INPUT/ocean_advect_velocity.nc'
if (file_exist(trim(filename))) then
call read_data(filename,'ue',ue,Dom%domain2d,timelevel=1)
call read_data(filename,'vn',vn,Dom%domain2d,timelevel=1)
call read_data(filename,'wb',wb,Dom%domain2d,timelevel=1)
do k=1,nk
do j=jsd,jed
do i=isd,ied
ue(i,j,k) = ue(i,j,k)*Grd%tmask(i,j,k)
vn(i,j,k) = vn(i,j,k)*Grd%tmask(i,j,k)
wb(i,j,k) = wb(i,j,k)*Grd%tmask(i,j,k)
enddo
enddo
enddo
call mpp_update_domains(ue(:,:,:), Dom%domain2d)
call mpp_update_domains(vn(:,:,:), Dom%domain2d)
call mpp_update_domains(wb(:,:,:), Dom%domain2d)
else
call mpp_error(FATAL,&
'==>Error from ocean_advection_velocity_mod: read_advect_velocity cannot find ocean_advect_velocity.nc.')
endif
if(read_advection_velocity) then
if(horz_grid == MOM_BGRID) then
do k=1,nk
Adv_vel%uhrho_et(:,:,k) = rho0*ue(:,:,k)*BAY(Thickness%dzu(:,:,k))
Adv_vel%vhrho_nt(:,:,k) = rho0*vn(:,:,k)*BAX(Thickness%dzu(:,:,k))
Adv_vel%wrho_bt(:,:,k) = rho0*wb(:,:,k)
Adv_vel%uhrho_eu(:,:,k) = REMAP_ET_TO_EU(Adv_vel%uhrho_et(:,:,k))
Adv_vel%vhrho_nu(:,:,k) = REMAP_NT_TO_NU(Adv_vel%vhrho_nt(:,:,k))
Adv_vel%wrho_bu(:,:,k) = REMAP_BT_TO_BU(Adv_vel%wrho_bt(:,:,k))
enddo
else
! U-cell advective components are not used for Cgrid
do k=1,nk
Adv_vel%uhrho_et(:,:,k) = rho0*ue(:,:,k)*Thickness%dzten(:,:,k,1)
Adv_vel%vhrho_nt(:,:,k) = rho0*vn(:,:,k)*Thickness%dzten(:,:,k,2)
Adv_vel%wrho_bt(:,:,k) = rho0*wb(:,:,k)
enddo
endif
id_ue = register_static_field('ocean_model','ue', Grd%tracer_axes_flux_x(1:3), &
'zonal C-grid velocity from file','m/s', range=(/-1e3,1e3/))
id_vn = register_static_field('ocean_model','vn', Grd%tracer_axes_flux_y(1:3), &
'meridional C-grid velocity from file','m/s', range=(/-1e3,1e3/))
id_wb = register_static_field('ocean_model','wb', Grd%tracer_axes_wt(1:3), &
'vertical C-grid velocity from file','m/s', range=(/-1e3,1e3/))
endif
if(read_advection_transport) then
if(horz_grid == MOM_BGRID) then
do k=1,nk
Adv_vel%uhrho_et(:,:,k) = ue(:,:,k)
Adv_vel%vhrho_nt(:,:,k) = vn(:,:,k)
Adv_vel%wrho_bt(:,:,k) = wb(:,:,k)
Adv_vel%uhrho_eu(:,:,k) = REMAP_ET_TO_EU(Adv_vel%uhrho_et(:,:,k))
Adv_vel%vhrho_nu(:,:,k) = REMAP_NT_TO_NU(Adv_vel%vhrho_nt(:,:,k))
Adv_vel%wrho_bu(:,:,k) = REMAP_BT_TO_BU(Adv_vel%wrho_bt(:,:,k))
enddo
else
! U-cell advective components are not used for Cgrid
do k=1,nk
Adv_vel%uhrho_et(:,:,k) = ue(:,:,k)
Adv_vel%vhrho_nt(:,:,k) = vn(:,:,k)
Adv_vel%wrho_bt(:,:,k) = wb(:,:,k)
enddo
endif
id_ue = register_static_field('ocean_model','ue', Grd%tracer_axes_flux_x(1:3), &
'zonal C-grid transport from file','(kg/m^3)*(m^2/sec)', range=(/-1e12,1e12/))
id_vn = register_static_field('ocean_model','vn', Grd%tracer_axes_flux_y(1:3), &
'meridional C-grid transport from file','(kg/m^3)*(m^2/sec)', range=(/-1e12,1e12/))
id_wb = register_static_field('ocean_model','wb', Grd%tracer_axes_wt(1:3), &
'vertical C-grid transport from file','m/s', range=(/-1e12,1e12/))
endif
if (id_ue > 0) used = send_data(id_ue, ue(isc:iec,jsc:jec,:), Time%model_time)
if (id_vn > 0) used = send_data(id_vn, vn(isc:iec,jsc:jec,:), Time%model_time)
if (id_wb > 0) used = send_data(id_wb, wb(isc:iec,jsc:jec,:), Time%model_time)
id_ue_rhodzt = register_static_field('ocean_model','ue_rhodzt', Grd%tracer_axes_flux_x(1:3), &
'zonal C-grid velocity*rhodzt from file','(kg/m3)*m^2/s', range=(/-1e6,1e6/))
if (id_ue_rhodzt > 0) used = send_data(id_ue_rhodzt, Adv_vel%uhrho_et(isc:iec,jsc:jec,:), Time%model_time)
id_vn_rhodzt = register_static_field('ocean_model','vn_rhodzt', Grd%tracer_axes_flux_y(1:3), &
'merid C-grid velocity*rhodzt from file','(kg/m3)*m^2/s', range=(/-1e6,1e6/))
if (id_vn_rhodzt > 0) used = send_data(id_vn_rhodzt, Adv_vel%vhrho_nt(isc:iec,jsc:jec,:), Time%model_time)
id_wb_rho = register_static_field('ocean_model','wb_rho', Grd%tracer_axes_wt(1:3), &
'vertical C-grid velocity*rho from file','(kg/m3)*m/s', range=(/-1e6,1e6/))
if (id_wb_rho > 0) used = send_data(id_wb_rho, Adv_vel%wrho_bt(isc:iec,jsc:jec,1:nk), Time%model_time)
end subroutine read_advection
! NAME="read_advection"
!#######################################################################
!
!
!
! Initialize the advection velocity used for specifying a nonzero
! southward inflow introduced to the domain from the northern boundary.
!
!
!
subroutine inflow_nboundary_init
character(len=128) :: filename
allocate ( vt_inflow(isd:ied,jsd:jed,nk) )
allocate ( mask_inflow(isd:ied,jsd:jed,nk) )
vt_inflow(:,:,:) = 0.0
mask_inflow(:,:,:) = 0.0
filename = 'INPUT/mask_inflow.nc'
if (file_exist(trim(filename))) then
call read_data(filename,'mask_inflow',mask_inflow,Dom%domain2d,timelevel=1)
call mpp_update_domains(mask_inflow(:,:,:), Dom%domain2d)
else
call mpp_error(FATAL,&
'==>Error from ocean_advection_velocity_mod: inflow_nboundary_init cannot find mask_inflow file.')
endif
filename = 'INPUT/vt_inflow.nc'
if (file_exist(trim(filename))) then
call read_data(filename,'vt_inflow',vt_inflow,Dom%domain2d,timelevel=1)
call mpp_update_domains(vt_inflow(:,:,:), Dom%domain2d)
else
call mpp_error(FATAL,&
'==>Error from ocean_advection_velocity_mod: inflow_nboundary_init cannot find vt_inflow.')
endif
end subroutine inflow_nboundary_init
! NAME="inflow_nboundary_init"
!#######################################################################
!
!
!
! Write a restart if necessary
!
!
!
subroutine ocean_advection_velocity_end(Time, Adv_vel, use_blobs)
type(ocean_time_type), intent(in) :: Time
type(ocean_adv_vel_type), intent(in) :: Adv_vel
logical, intent(in) :: use_blobs
integer :: stdoutunit
stdoutunit = stdout()
if (use_blobs) call save_restart(adv_vel_restart)
if (debug_this_module) then
write(stdoutunit,*) ' '
write(stdoutunit,*) 'From ocean_advection_velocity_mod: ending chksums'
call write_timestamp(Time%model_time)
call write_chksum_3d('Adv_vel%uhrho_et', Adv_vel%uhrho_et(COMP,:))
call write_chksum_3d('Adv_vel%vhrho_nt', Adv_vel%vhrho_nt(COMP,:))
call write_chksum_3d('Adv_vel%uhrho_eu', Adv_vel%uhrho_eu(COMP,:))
call write_chksum_3d('Adv_vel%vhrho_nu', Adv_vel%vhrho_nu(COMP,:))
call write_chksum_3d('Adv_vel%wrho_bt', Adv_vel%wrho_bt(COMP,1:))
call write_chksum_3d('Adv_vel%wrho_bu', Adv_vel%wrho_bu(COMP,1:))
call write_chksum_2d('Adv_vel%wrho_bt(0)', Adv_vel%wrho_bt(COMP,0))
call write_chksum_2d('Adv_vel%wrho_bu(0)', Adv_vel%wrho_bu(COMP,0))
call write_chksum_3d('Adv_vel%diverge_t', Adv_vel%diverge_t(COMP,:))
call write_chksum_3d('Adv_vel%diverge_u', Adv_vel%diverge_u(COMP,:))
endif
end subroutine ocean_advection_velocity_end
! NAME="ocean_advection_velocity_end"
end module ocean_advection_velocity_mod