!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! GNU General Public License !!
!! !!
!! This file is part of the Flexible Modeling System (FMS). !!
!! !!
!! FMS is free software; you can redistribute it and/or modify !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation. !!
!! !!
!! FMS is distributed in the hope that it will be useful, !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !!
!! GNU General Public License for more details. !!
!! !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to: !!
!! Free Software Foundation, Inc. !!
!! 59 Temple Place, Suite 330 !!
!! Boston, MA 02111-1307 USA !!
!! or see: !!
!! http://www.gnu.org/licenses/gpl.txt !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module ocean_mixdownslope_mod
!
! S.M. Griffies
!
!
!
! Mixing of tracer between dense shallow parcel and
! deeper parcels downslope.
!
!
!
! Mixing of tracer properties as dense shallow parcel is discharged
! into deeper water to approach the parcel's depth of neutral buoyancy.
! This module can be characterized as a mixture of the approach from
! Campin and Goosse (1999) and slope convection.
!
!
!
!
!
! Campin and Goosse (1999): Parameterization of density-driven downsloping flow
! for a coarse-resolution model in z-coordinate", Tellus 51A, pages 412-430
!
!
!
! S.M. Griffies, Elements of mom4p1 (2006)
! NOAA/Geophysical Fluid Dynamics Laboratory
!
!
!
!
!
!
! For using this module. Default use_this_module=.false.
!
!
! For debugging
!
!
! Set true to do bitwise exact global sum. When it is false, the global
! sum will be non-bitwise_exact, but will significantly increase efficiency.
! The default value is false.
!
!
!
! Number of horizontally distant points used in search downslope.
! Note: it is not possible to have
! mixdownslope_npts greater than or equal to the computational domain
! extents, as this would require updates across multiple processors.
! Default mixdownslope_npts=1.
!
!
! Fraction of the central cell that participates in downslope mixing
! in any particular direction. Default mixdownslope_frac_central=0.25
!
!
!
! To place more weight on points further from central point. This may
! be done to enhance properties getting downslope. Default is
! mixdownslope_weight_far=.false.
!
!
! Width of the re-weighting function used to emphasize points further
! along in the search for exchange points. Default mixdownslope_width=1.
!
!
!
! For reading in a mask that selects regions of the domain
! where mixdownslope is allowed to function (mask=1) or not
! to function (mask=0). Default read_mixdownslope_mask=.false.,
! whereby mixdownslope_mask is set to tmask(k=1).
!
!
! For modifying the mixdownslope mask based on reading in
! the GFDL regional mask. Default mixdownslope_mask_gfdl=.false.
!
!
!
use constants_mod, only: epsln
use diag_manager_mod, only: register_diag_field, register_static_field, send_data
use fms_mod, only: write_version_number, error_mesg, read_data, FATAL, NOTE
use fms_mod, only: open_namelist_file, check_nml_error, close_file, stdout, stdlog
use mpp_domains_mod, only: mpp_update_domains, CGRID_NE
use mpp_domains_mod, only: mpp_global_sum, BITWISE_EXACT_SUM, NON_BITWISE_EXACT_SUM
use mpp_mod, only: mpp_error, mpp_chksum
use ocean_density_mod, only: density
use ocean_domains_mod, only: get_local_indices, set_ocean_domain
use ocean_parameters_mod,only: missing_value
use ocean_parameters_mod,only: TERRAIN_FOLLOWING
use ocean_types_mod, only: ocean_domain_type, ocean_grid_type, ocean_time_type, ocean_options_type
use ocean_types_mod, only: ocean_prog_tracer_type, ocean_density_type, ocean_thickness_type
use ocean_util_mod, only: write_timestamp
use ocean_workspace_mod, only: wrk1, wrk1_v
implicit none
private
public ocean_mixdownslope_init
public mixdownslope
type(ocean_domain_type), pointer :: Dom => NULL()
type(ocean_grid_type), pointer :: Grd => NULL()
#include
integer, dimension(:,:), allocatable :: ip ! shifting in i-direction
integer, dimension(:,:), allocatable :: jq ! shifting in j-direction
real, dimension(:,:), allocatable :: data ! for reading in mask information
real, dimension(:,:), allocatable :: area_k ! area on k-level, including tmask
real, dimension(:,:), allocatable :: slope_x ! topog slope (dimensionless) in the i-direction
real, dimension(:,:), allocatable :: slope_y ! topog slope (dimensionless) in the j-direction
real, dimension(:,:), allocatable :: mixdownslope_mask ! mask=1 when apply mixdownslope, 0 otherwise
real, dimension(:,:,:), allocatable :: topog_step ! where downslope flow is topgraphically possible
real, dimension(:,:,:), allocatable :: tend_mix ! (rho*dzt*tracer) tendency
real, dimension(:,:,:,:), allocatable :: mixdownslope_frac ! fraction of mixing occurring between two cells
! fields with extended halos
integer, dimension(:,:), allocatable :: kmt_ex
integer, dimension(:,:,:), allocatable :: kup_ex
integer, dimension(:,:,:,:), allocatable :: kdw_ex
real, dimension(:,:,:), allocatable :: topog_slope_ex
real, dimension(:,:,:), allocatable :: topog_step_ex
real, dimension(:,:,:), allocatable :: tracer_ex
real, dimension(:,:,:), allocatable :: temp_ex
real, dimension(:,:,:), allocatable :: salt_ex
real, dimension(:,:,:), allocatable :: press_ex
real, dimension(:,:,:), allocatable :: mass_ex
real, dimension(:,:,:), allocatable :: rho_ex
! domain for mixdownslope
type(ocean_domain_type), save :: Mixdownslope_domain
! for diagnostic manager
logical :: used
integer :: id_slope_x=-1
integer :: id_slope_y=-1
integer :: id_topog_step_1 =-1
integer :: id_topog_step_2 =-1
integer :: id_topog_step_3 =-1
integer :: id_topog_step_4 =-1
integer :: id_mixdownslope_mask =-1
integer :: id_neutral_rho_mixdown=-1
integer :: id_wdian_rho_mixdown =-1
integer, dimension(:), allocatable :: id_mixdownslope
character(len=128) :: version=&
'=>Using: ocean_mixdownslope.f90 ($Id: ocean_mixdownslope.F90,v 16.0.108.1 2009/10/10 00:42:23 nnz Exp $)'
character (len=128) :: tagname=&
'$Name: mom4p1_pubrel_dec2009_nnz $'
! number of prognostic tracers
integer :: num_prog_tracers=0
! processor zero writes to unit 6
integer :: unit=6
! initialization flag
logical :: module_is_initialized=.false.
! flag for mpp_global_sum
integer :: global_sum_flag
! halo
integer :: ijhalo=1
! time step
real :: dtime
real :: dtimer
! set from nml
logical :: use_this_module = .false. ! must be set .true. in nml to enable this scheme
logical :: debug_this_module = .false. ! for debugging
logical :: do_bitwise_exact_sum = .false. ! set true to get slower sum that is same for all PEs.
logical :: mixdownslope_weight_far= .false. ! to place more weight on points further from central point
logical :: read_mixdownslope_mask = .false. ! for applying a mask where do not apply mixdownslope
logical :: mixdownslope_mask_gfdl = .false. ! for case when applying GFDL regional mask
integer :: mixdownslope_npts = 1 ! number horizontal points searched for the downslope mixing
integer :: mixdownslope_width = 1 ! width for exponential that damps points near to central point
real :: mixdownslope_frac_central =0.25 ! fraction of central cell participating in downslope mixing
namelist /ocean_mixdownslope_nml/ use_this_module, debug_this_module, mixdownslope_npts, &
mixdownslope_width, mixdownslope_weight_far, &
mixdownslope_frac_central, do_bitwise_exact_sum, &
read_mixdownslope_mask, mixdownslope_mask_gfdl
contains
!#######################################################################
!
!
!
! Initial set up for mixing of tracers into the abyss next to topography.
!
!
subroutine ocean_mixdownslope_init(Grid, Domain, Time, T_prog, Ocean_options, &
vert_coordinate_type, dtim, 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_prog_tracer_type), intent(inout) :: T_prog(:)
type(ocean_options_type), intent(inout) :: Ocean_options
integer, intent(in) :: vert_coordinate_type
real, intent(in) :: dtim
logical, intent(in), optional :: debug
integer :: io_status, ioun, ierr
integer :: i,j,m,n,kbot
integer :: stdoutunit,stdlogunit
stdoutunit=stdout();stdlogunit=stdlog()
if ( module_is_initialized ) then
call mpp_error(FATAL, &
'==>Error from ocean_mixdownslope_mod (ocean_mixdownslope_init): module already initialized')
endif
module_is_initialized = .TRUE.
call write_version_number( version, tagname )
#ifndef MOM4_STATIC_ARRAYS
call get_local_indices(Domain,isd,ied,jsd,jed,isc,iec,jsc,jec)
nk = Grid%nk
#endif
Dom => Domain
Grd => Grid
dtime = dtim
dtimer = 1.0/dtime
num_prog_tracers = size(T_prog(:))
! provide for namelist over-ride of default values
ioun = open_namelist_file()
read (ioun,ocean_mixdownslope_nml,IOSTAT=io_status)
write (stdlogunit,ocean_mixdownslope_nml)
write (stdoutunit,'(/)')
write (stdoutunit,ocean_mixdownslope_nml)
ierr = check_nml_error(io_status,'ocean_mixdownslope_nml')
call close_file (ioun)
if(do_bitwise_exact_sum) then
global_sum_flag = BITWISE_EXACT_SUM
else
global_sum_flag = NON_BITWISE_EXACT_SUM
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_mixdownslope_mod with debug_this_module=.true.'
endif
if(.not. use_this_module) then
call mpp_error(NOTE,&
'==>From ocean_mixdownslope_mod: NOT using downslope mixing scheme.')
Ocean_options%mixdownslope = 'Did NOT use downslope mixing scheme.'
return
else
if(vert_coordinate_type == TERRAIN_FOLLOWING) then
call mpp_error(FATAL, &
'==>ocean_mixdownslope_mod: this module is NOT for use with TERRAIN_FOLLOWING vert coodinates.')
endif
Ocean_options%mixdownslope = 'Used the downslope mixing scheme. This scheme is experimental!'
call mpp_error(NOTE,&
'==>From ocean_mixdownslope_mod: USING downslope mixing scheme.')
endif
if(debug_this_module) then
call mpp_error(NOTE,'==>From ocean_mixdownslope_mod: USING debug_this_module')
endif
! allocate for tracer tendency array
allocate(tend_mix(isd:ied,jsd:jed,nk))
tend_mix(:,:,:) = 0.0
if(mixdownslope_npts<1) then
call mpp_error(FATAL, &
'==>ocean_mixdownslope_mod: mixdownslope_npts < 1 not allowed. In nml, set mixdownslope_npts>=1.')
endif
write(stdoutunit,'(a,i3)')' In ocean_mixdownslope_mod: mixdownslope_npts = ',mixdownslope_npts
write(stdoutunit,'(a)') ' Be sure this number is smaller than dimensions of computational domain.'
! set mixdownslope_mask=0.0 in those regions where mixdownslope is NOT applied
! and mixdownslope_mask=1.0 in regions where mixdownslope is applied. Default is
! mixdownslope everywhere (mixdownslope_mask(:,:) = Grd%tmask(:,:,1)).
allocate(mixdownslope_mask(isd:ied,jsd:jed))
mixdownslope_mask(:,:) = Grd%tmask(:,:,1)
if(read_mixdownslope_mask) then
allocate(data(isd:ied,jsd:jed))
call read_data('INPUT/mixdownslope_mask','mixdownslope_mask',data,Domain%domain2d)
do j=jsc,jec
do i=isc,iec
mixdownslope_mask(i,j) = data(i,j)
enddo
enddo
! the following is specific to the mask used at GFDL for
! the OM3 ocean model configuration.
! here, remove the Black Sea (labelled with 7.0), as this
! region contains some odd water masses which cause instabilities
! with the mixdownslope scheme.
if(mixdownslope_mask_gfdl) then
do j=jsc,jec
do i=isc,iec
if(mixdownslope_mask(i,j)==0.0 .or. mixdownslope_mask(i,j)==7.0) then
mixdownslope_mask(i,j)=0.0
else
mixdownslope_mask(i,j)=1.0
endif
enddo
enddo
endif
call mpp_update_domains(mixdownslope_mask(:,:), Dom%domain2d)
endif
! halo size for extended domain and larger arrays
ijhalo=mixdownslope_npts
! k-level of central point in search
allocate(kup_ex(isc-ijhalo:iec+ijhalo,jsc-ijhalo:jec+ijhalo,4))
kup_ex(:,:,:) = 0
! k-level of deeper points in horizontally adjacent columns
allocate(kdw_ex(isc-ijhalo:iec+ijhalo,jsc-ijhalo:jec+ijhalo,ijhalo,4))
kdw_ex(:,:,:,:) = 0
! fraction of a cell participating in mixing
allocate(mixdownslope_frac(isc-ijhalo:iec+ijhalo,jsc-ijhalo:jec+ijhalo,ijhalo,4))
mixdownslope_frac(:,:,:,:) = 0.0
! define extended domain
call set_ocean_domain(Mixdownslope_domain,Grd,xhalo=ijhalo,yhalo=ijhalo,&
name='mixdownslope',maskmap=Dom%maskmap)
! time independent arrays defined over Mixdownslope_domain
allocate(kmt_ex (isc-ijhalo:iec+ijhalo,jsc-ijhalo:jec+ijhalo))
kmt_ex(:,:) = 0
kmt_ex (isc:iec,jsc:jec) = Grd%kmt (isc:iec,jsc:jec)
call mpp_update_domains (kmt_ex(:,:), Mixdownslope_domain%domain2d)
! time dependent arrays defined over Mixdownslope_domain
allocate(temp_ex (isc-ijhalo:iec+ijhalo,jsc-ijhalo:jec+ijhalo,nk))
allocate(salt_ex (isc-ijhalo:iec+ijhalo,jsc-ijhalo:jec+ijhalo,nk))
allocate(tracer_ex (isc-ijhalo:iec+ijhalo,jsc-ijhalo:jec+ijhalo,nk))
allocate(press_ex (isc-ijhalo:iec+ijhalo,jsc-ijhalo:jec+ijhalo,nk))
allocate(rho_ex (isc-ijhalo:iec+ijhalo,jsc-ijhalo:jec+ijhalo,nk))
allocate(mass_ex (isc-ijhalo:iec+ijhalo,jsc-ijhalo:jec+ijhalo,nk))
temp_ex(:,:,:) = 0.0
salt_ex(:,:,:) = 0.0
tracer_ex(:,:,:) = 0.0
press_ex(:,:,:) = 0.0
rho_ex(:,:,:) = 0.0
mass_ex(:,:,:) = 0.0
! area of cells on k-level, including tmask
allocate (area_k(isd:ied,jsd:jed))
area_k(:,:) = 0.0
! compute topographic slope arrays for the i-slope and j-slope
! slopes are centered on the i-face and j-face of tracer cells
allocate (slope_x(isd:ied,jsd:jed))
allocate (slope_y(isd:ied,jsd:jed))
slope_x = 0.0
slope_y = 0.0
do j=jsc,jec
do i=isc,iec
slope_x(i,j) = (Grd%ht(i+1,j)-Grd%ht(i,j))*Grd%dxter(i,j)*Grd%tmask(i,j,1)*Grd%tmask(i+1,j,1)
slope_y(i,j) = (Grd%ht(i,j+1)-Grd%ht(i,j))*Grd%dytnr(i,j)*Grd%tmask(i,j,1)*Grd%tmask(i,j+1,1)
slope_x(i,j) = abs(slope_x(i,j))
slope_y(i,j) = abs(slope_y(i,j))
enddo
enddo
call mpp_update_domains(slope_x(:,:),Dom%domain2d)
call mpp_update_domains(slope_y(:,:),Dom%domain2d)
! topographic slope for the four surrounding directions
allocate (topog_slope_ex(isc-ijhalo:iec+ijhalo,jsc-ijhalo:jec+ijhalo,4))
topog_slope_ex(:,:,:) = 0.0
do j=jsc,jec
do i=isc,iec
m=1 ; topog_slope_ex(i,j,m) = slope_x(i,j)
m=2 ; topog_slope_ex(i,j,m) = slope_y(i,j)
m=3 ; topog_slope_ex(i,j,m) = slope_x(i-1,j)
m=4 ; topog_slope_ex(i,j,m) = slope_y(i,j-1)
enddo
enddo
call mpp_update_domains (topog_slope_ex(:,:,:), Mixdownslope_domain%domain2d)
! compute directions from an (i,j) point where topography deepens.
! these directions may potentially have downslope mixing.
! insist that downslope mixing occurs only when there are more
! kmt cells in the adjacent column. also insist that downslope
! mixing does not involve k=1 cells.
allocate (topog_step(isd:ied,jsd:jed,4))
topog_step(:,:,:) = 0.0
do j=jsc,jec
do i=isc,iec
kbot = Grd%kmt(i,j)
if(kbot > 1) then
if(Grd%kmt(i+1,j) > kbot) topog_step(i,j,1)=1.0
if(Grd%kmt(i,j+1) > kbot) topog_step(i,j,2)=1.0
if(Grd%kmt(i-1,j) > kbot) topog_step(i,j,3)=1.0
if(Grd%kmt(i,j-1) > kbot) topog_step(i,j,4)=1.0
endif
enddo
enddo
! block out the bipolar fold in order to ensure tracer conservation.
! The reason we do so is related to how the algorithm reaches between
! adjacent columns of tracer points. When the column straddles
! the bipolar fold, the code logic is not general and so it actually
! attempts to reach to a non-adjacent column. Shy of generalizing
! the logic, we simply do not consider mixdownslope for points along fold.
if(jec+Dom%joff==Dom%jeg) then
topog_step(:,jec,:) = 0.0
endif
! fill the larger array
allocate (topog_step_ex(isc-ijhalo:iec+ijhalo,jsc-ijhalo:jec+ijhalo,4))
topog_step_ex(:,:,:) = 0.0
topog_step_ex(isc:iec,jsc:jec,:) = topog_step(isc:iec,jsc:jec,:)
call mpp_update_domains (topog_step_ex(:,:,:), Mixdownslope_domain%domain2d)
! labels for the four quadranats surrounding a
! central tracer cell (moving counter-clockwise)
allocate(ip(0:ijhalo,4))
allocate(jq(0:ijhalo,4))
do n=0,ijhalo
m=1
ip(n,m) = 1+(n-1)
jq(n,m) = 0
m=2
ip(n,m) = 0
jq(n,m) = 1+(n-1)
m=3
ip(n,m) = -1-(n-1)
jq(n,m) = 0
m=4
ip(n,m) = 0
jq(n,m) = -1-(n-1)
enddo
! register/send diagnostics
id_mixdownslope_mask = register_static_field ('ocean_model', 'mixdownslope_mask', &
Grd%tracer_axes(1:2), &
'mixdownslope mask', 'dimensionless', missing_value=missing_value, range=(/-1.0,1e4/))
if (id_mixdownslope_mask > 0) used = send_data (id_mixdownslope_mask, &
mixdownslope_mask(isc:iec,jsc:jec), Time%model_time, rmask=Grd%tmask(isc:iec,jsc:jec,1))
id_slope_x = register_static_field ('ocean_model', 'slope_x', Grd%tracer_axes_flux_x(1:2), &
'|d(ht)/dx| on T-cell face', 'm/m', missing_value=missing_value, range=(/-1.e9,1.e9/))
if (id_slope_x > 0) used = send_data (id_slope_x, slope_x(isc:iec,jsc:jec), &
Time%model_time, rmask=Grd%tmask(isc:iec,jsc:jec,1))
id_slope_y = register_static_field ('ocean_model', 'slope_y', Grd%tracer_axes_flux_y(1:2), &
'|d(ht)/dy| on T-cell face', 'm/m', missing_value=missing_value, range=(/-1.e9,1.e9/))
if (id_slope_y > 0) used = send_data (id_slope_y, slope_y(isc:iec,jsc:jec), &
Time%model_time, rmask=Grd%tmask(isc:iec,jsc:jec,1))
id_topog_step_1 = register_static_field ('ocean_model', 'topog_step_1', Grd%tracer_axes(1:2), &
'topog_step_1', 'dimensionless', missing_value=missing_value, range=(/-1.0,1.0/))
if (id_topog_step_1 > 0) used = send_data (id_topog_step_1, topog_step(1,isc:iec,jsc:jec), &
Time%model_time, rmask=Grd%tmask(isc:iec,jsc:jec,1))
id_topog_step_2 = register_static_field ('ocean_model', 'topog_step_2', Grd%tracer_axes(1:2), &
'topog_step_2', 'dimensionless', missing_value=missing_value, range=(/-1.0,1.0/))
if (id_topog_step_2 > 0) used = send_data (id_topog_step_2, topog_step(2,isc:iec,jsc:jec), &
Time%model_time, rmask=Grd%tmask(isc:iec,jsc:jec,1))
id_topog_step_3 = register_static_field ('ocean_model', 'topog_step_3', Grd%tracer_axes(1:2), &
'topog_step_3', 'dimensionless', missing_value=missing_value, range=(/-1.0,1.0/))
if (id_topog_step_3 > 0) used = send_data (id_topog_step_3, topog_step(3,isc:iec,jsc:jec), &
Time%model_time, rmask=Grd%tmask(isc:iec,jsc:jec,1))
id_topog_step_4 = register_static_field ('ocean_model', 'topog_step_4', Grd%tracer_axes(1:2), &
'topog_step_4', 'dimensionless', missing_value=missing_value, range=(/-1.0,1.0/))
if (id_topog_step_4 > 0) used = send_data (id_topog_step_4, topog_step(4,isc:iec,jsc:jec), &
Time%model_time, rmask=Grd%tmask(isc:iec,jsc:jec,1))
allocate (id_mixdownslope(num_prog_tracers))
id_mixdownslope = -1
do n=1,num_prog_tracers
if(T_prog(n)%name == 'temp') then
id_mixdownslope(n) = register_diag_field ('ocean_model', &
'mixdownslope_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:3), Time%model_time, &
'cp*mixdownslope*rho*dzt*temp', &
'Watt/m^2', missing_value=missing_value, range=(/-1.e9,1.e9/))
else
id_mixdownslope(n) = register_diag_field ('ocean_model', 'mixdownslope_'//trim(T_prog(n)%name), &
Grd%tracer_axes(1:3), Time%model_time, &
'mixdownslope*rho*dzt*tracer for '//trim(T_prog(n)%name), &
trim(T_prog(n)%flux_units), missing_value=missing_value, range=(/-1.e9,1.e9/))
endif
enddo
id_neutral_rho_mixdown = register_diag_field ('ocean_model', 'neutral_rho_mixdown', &
Grd%tracer_axes(1:3), Time%model_time, &
'update of neutral density from mixdownange scheme', &
'rho*rho_dz/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
id_wdian_rho_mixdown = register_diag_field ('ocean_model', 'wdian_rho_mixdown', &
Grd%tracer_axes(1:3), Time%model_time, &
'dianeutral velocity component due to mixdownange scheme', &
'm/sec', missing_value=missing_value, range=(/-1.e10,1.e10/))
if (debug_this_module) then
write(stdoutunit,*) ' '
write(stdoutunit,*) '==Global sums from initialization of ocean_mixdownslope_mod== '
call write_timestamp(Time%model_time)
write(stdoutunit,'(a,es24.17)') &
'slope_x = ',mpp_global_sum(Dom%domain2d,slope_x(:,:), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'slope_y = ',mpp_global_sum(Dom%domain2d,slope_y(:,:), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'topog_slope_ex(1)= ',mpp_global_sum(Mixdownslope_domain%domain2d,topog_slope_ex(:,:,1), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'topog_slope_ex(2)= ',mpp_global_sum(Mixdownslope_domain%domain2d,topog_slope_ex(:,:,2), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'topog_slope_ex(3)= ',mpp_global_sum(Mixdownslope_domain%domain2d,topog_slope_ex(:,:,3), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'topog_slope_ex(4)= ',mpp_global_sum(Mixdownslope_domain%domain2d,topog_slope_ex(:,:,4), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'topog_step(1) = ',mpp_global_sum(Dom%domain2d,topog_step(:,:,1), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'topog_step(2) = ',mpp_global_sum(Dom%domain2d,topog_step(:,:,2), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'topog_step(3) = ',mpp_global_sum(Dom%domain2d,topog_step(:,:,3), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'topog_step(4) = ',mpp_global_sum(Dom%domain2d,topog_step(:,:,4), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'topog_step_ex(1) = ',mpp_global_sum(Mixdownslope_domain%domain2d,topog_step_ex(:,:,1), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'topog_step_ex(2) = ',mpp_global_sum(Mixdownslope_domain%domain2d,topog_step_ex(:,:,2), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'topog_step_ex(3) = ',mpp_global_sum(Mixdownslope_domain%domain2d,topog_step_ex(:,:,3), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'topog_step_ex(4) = ',mpp_global_sum(Mixdownslope_domain%domain2d,topog_step_ex(:,:,4), global_sum_flag)
endif
end subroutine ocean_mixdownslope_init
! NAME="ocean_mixdownslope_init"
!#######################################################################
!
!
!
! Compute thickness and density weighted tracer tendency [tracer*rho*m/s]
! due to exchange of tracer properties in regions where density-driven
! downslope transport is favorable.
!
! Allow for exchanges to occur over horizontally
! distant points, so long as the dense shallow parcel finds that it
! will sit on the bottom of the horizontally adjacent columns. Doing
! so requires a search algorithm, which requires some if-test logic
! as well as extended halos. Note that the halos cannot be extended
! to larger than the size of the computational domain on a processor.
! This restriction limits the extent that we can search horizontally.
!
! The rates for the exchange are functions of the topographic slope
! and the density differences between parcels.
!
! This scheme can be characterized as a slope convection based on
! logic incorporated into the overflow and overexchange schemes.
!
!
!
subroutine mixdownslope (Time, Thickness, T_prog, Dens, index_temp, index_salt)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_prog_tracer_type), intent(inout) :: T_prog(:)
type(ocean_density_type), intent(in) :: Dens
integer, intent(in) :: index_temp
integer, intent(in) :: index_salt
integer :: tau, taum1
integer :: i, j, k, m, n, nt, npts
integer :: kmtij, ku, kd
integer :: iip0, iip1, jjq0, jjq1
integer :: iip1r, jjq1r
integer :: nm1
real :: temp_so, salt_so, press, density_check
real :: weight, arg
real :: delta, delta_rho(ijhalo)
real :: mass_sum, tmix
real :: mixdownslope_total, mixdownslope_total_r
real :: tendency, temporary
integer :: stdoutunit
stdoutunit=stdout()
if(.not. use_this_module) return
if(.not. module_is_initialized ) then
call mpp_error(FATAL, &
'==>Error from ocean_mixdownslope_mod (mixdownslope): module must be initialized')
endif
tau = Time%tau
taum1 = Time%taum1
delta_rho(:) = 0.0
! extend some fields to extended domain
rho_ex = 0.0
mass_ex = 0.0
press_ex = 0.0
temp_ex = 0.0
salt_ex = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
rho_ex(i,j,k) = Dens%rho(i,j,k,tau)
mass_ex(i,j,k) = Thickness%rho_dzt(i,j,k,tau)*Grd%dat(i,j)
press_ex(i,j,k) = Dens%pressure_at_depth(i,j,k)
temp_ex(i,j,k) = T_prog(index_temp)%field(i,j,k,tau)
salt_ex(i,j,k) = T_prog(index_salt)%field(i,j,k,tau)
enddo
enddo
enddo
call mpp_update_domains(rho_ex(:,:,:), Mixdownslope_domain%domain2d, complete=.false.)
call mpp_update_domains(mass_ex(:,:,:), Mixdownslope_domain%domain2d, complete=.false.)
call mpp_update_domains(press_ex(:,:,:), Mixdownslope_domain%domain2d, complete=.false.)
call mpp_update_domains(temp_ex(:,:,:), Mixdownslope_domain%domain2d, complete=.false.)
call mpp_update_domains(salt_ex(:,:,:), Mixdownslope_domain%domain2d, complete=.true.)
! compute details for cells participating in downslope mixing
! note: "so" = "shallow ocean" cell...the central point
! note: "do" = cells within "deep ocean" columns
! this part of the compuatation is independent of the tracer
kup_ex(:,:,:) = 0
kdw_ex(:,:,:,:) = 0
mixdownslope_frac(:,:,:,:) = 0.0
do j=jsc,jec
do i=isc,iec
kmtij = max(1,kmt_ex(i,j))
temp_so = temp_ex(i,j,kmtij)
salt_so = salt_ex(i,j,kmtij)
! 4-directions surrounding each cell
do m=1,4
! search for density favorable mixing
npts=0
nloop1_npts: do n=1,mixdownslope_npts
iip0 = i+ip(n-1,m)
jjq0 = j+jq(n-1,m)
iip1 = i+ip(n,m)
jjq1 = j+jq(n,m)
! check if downslope mixing is topographically possible
if(topog_step_ex(iip0,jjq0,m)==1.0) then
! check if density of shallow ocean cell > density of deep ocean cell
if(rho_ex(i,j,kmtij) > rho_ex(iip1,jjq1,kmtij)) then
! kdw = k-level in do-column where central cell is neutrally buoyant (or at bottom)
kdw_ex(i,j,n,m) = 0
kloop1 : do k=kmtij+1,kmt_ex(iip1,jjq1)
press = press_ex(iip1,jjq1,k)
density_check = density(salt_so,temp_so,press)
if(density_check > rho_ex(iip1,jjq1,k)) then
kdw_ex(i,j,n,m) = k
delta_rho(n) = density_check-rho_ex(iip1,jjq1,k)
else
exit kloop1
endif
enddo kloop1
! check that the n-parcel has a k-level greater than the n-1 parcel.
! if not, then do not mix with the n-column.
nm1 = max(n-1,1)
if(n > 1) then
if(kdw_ex(i,j,n,m) <= kdw_ex(i,j,nm1,m)) then
kdw_ex(i,j,n,m) = 0
endif
endif
! set strength of downslope mixing between central cell and n-parcel
if(kdw_ex(i,j,n,m) > 0) then
! add this cell to the number of cells participating in mixing
npts=npts+1
if(n==1) then
kup_ex(i,j,m) = kmtij
delta = rho_ex(i,j,kmtij)-rho_ex(iip1,jjq1,kmtij)
else
delta = delta_rho(n)
endif
! compute mixing weight as product of slope and density difference
mixdownslope_frac(i,j,n,m) = topog_slope_ex(iip0,jjq0,m)*delta
endif ! kdw_ex(i,j,n,m) > 0 if-test
endif ! rho_ex(i,j,k) > rho_ex(iip1,jjq1,k) if-test
endif ! topog_step_ex(iip0,jjq0,m)==1.0 if-test
! if kdw is not on the bottom, then exit n-loop since finished with search
if(kdw_ex(i,j,n,m) < kmt_ex(iip1,jjq1)) then
exit nloop1_npts
endif
enddo nloop1_npts ! n-loop for horizontal search reaching out from central point
! place more weight on the farther points to
! encourage tracer exchange further downslope.
if(mixdownslope_weight_far .and. npts > 0) then
weight = mixdownslope_frac(i,j,1,m)
do n=1,npts
arg = float((n-npts)/mixdownslope_width)
mixdownslope_frac(i,j,n,m) = weight*exp(arg)
enddo
endif
! normalize mixdownslope_frac to set fraction of
! a cell that is being mixed with central cell.
if(npts > 0) then
mixdownslope_total=0.0
do n=1,npts
mixdownslope_total = mixdownslope_total + mixdownslope_frac(i,j,n,m)
enddo
mixdownslope_total_r = 1.0/(mixdownslope_total+epsln)
do n=1,npts
mixdownslope_frac(i,j,n,m) = mixdownslope_frac(i,j,n,m)*mixdownslope_total_r
enddo
endif
! if no exchange points
if(npts==0) then
mixdownslope_frac(i,j,:,m) = 0.0
endif
enddo ! m-loop
enddo ! i-loop
enddo ! j-loop
! extend arrays to wide halos
call mpp_update_domains(kup_ex(:,:,:), Mixdownslope_domain%domain2d)
call mpp_update_domains(kdw_ex(:,:,:,:), Mixdownslope_domain%domain2d)
call mpp_update_domains(mixdownslope_frac(:,:,:,:), Mixdownslope_domain%domain2d)
! compute tracer tendency for cells participating in downslope mixing
do nt=1,num_prog_tracers
! place tracer concentration in the wider array tracer_ex
if(nt==index_temp) then
do k=1,nk
do j=jsc-ijhalo,jec+ijhalo
do i=isc-ijhalo,iec+ijhalo
tracer_ex(i,j,k) = temp_ex(i,j,k)
enddo
enddo
enddo
elseif(nt==index_salt) then
do k=1,nk
do j=jsc-ijhalo,jec+ijhalo
do i=isc-ijhalo,iec+ijhalo
tracer_ex(i,j,k) = salt_ex(i,j,k)
enddo
enddo
enddo
else
do k=1,nk
do j=jsc,jec
do i=isc,iec
tracer_ex(i,j,k) = T_prog(nt)%field(i,j,k,tau)
enddo
enddo
enddo
call mpp_update_domains (tracer_ex(:,:,:), Mixdownslope_domain%domain2d)
endif
! compute tendency, noting that each (i,j,k) cell can generally be part
! of mixing as either a central shallow cell (n=0), or as one of the deep cells (n>0).
! for the tendency, we compute a mixing tracer concentration between two cells,
! and then back-out a tendency which is placed in tend_mix.
tend_mix(:,:,:) = 0.0
do j=jsc,jec
do i=isc,iec
do m=1,4
do n=1,mixdownslope_npts
! i,j is central cell at k=ku
iip1 = i+ip(n,m)
jjq1 = j+jq(n,m)
ku = kup_ex(i,j,m)
kd = kdw_ex(i,j,n,m)
if(ku > 0 .and. kd > ku) then
mass_sum = mixdownslope_frac_central*mass_ex(i,j,ku) &
+mixdownslope_frac(i,j,n,m)*mass_ex(iip1,jjq1,kd)
tmix = (mixdownslope_frac_central*mass_ex(i,j,ku)*tracer_ex(i,j,ku) &
+mixdownslope_frac(i,j,n,m)*mass_ex(iip1,jjq1,kd)*tracer_ex(iip1,jjq1,kd)) &
/mass_sum
tendency = dtimer*mixdownslope_frac_central*mass_ex(i,j,ku)/Grd%dat(i,j) &
*(tmix-tracer_ex(i,j,ku))
tend_mix(i,j,ku) = tend_mix(i,j,ku) + tendency
endif
! i,j is deep cell at k=kd
iip1r = i-ip(n,m)
jjq1r = j-jq(n,m)
ku = kup_ex(iip1r,jjq1r,m)
kd = kdw_ex(iip1r,jjq1r,n,m)
if(ku > 0 .and. kd > ku) then
mass_sum = mixdownslope_frac_central*mass_ex(iip1r,jjq1r,ku) &
+mixdownslope_frac(iip1r,jjq1r,n,m)*mass_ex(i,j,kd)
tmix = (mixdownslope_frac_central &
*mass_ex(iip1r,jjq1r,ku)*tracer_ex(iip1r,jjq1r,ku) &
+mixdownslope_frac(iip1r,jjq1r,n,m) &
*mass_ex(i,j,kd)*tracer_ex(i,j,kd)) &
/mass_sum
tendency = dtimer*mixdownslope_frac(iip1r,jjq1r,n,m)*mass_ex(i,j,kd)/Grd%dat(i,j) &
*(tmix-tracer_ex(i,j,kd))
tend_mix(i,j,kd) = tend_mix(i,j,kd) + tendency
endif
enddo ! n-loop
enddo ! m-loop
enddo ! i-loop
enddo ! j-loop
! fill tracer tendency array
do k=1,nk
do j=jsc,jec
do i=isc,iec
T_prog(nt)%th_tendency(i,j,k) = T_prog(nt)%th_tendency(i,j,k) &
+tend_mix(i,j,k)*mixdownslope_mask(i,j)
enddo
enddo
enddo
if(id_mixdownslope(nt) > 0) then
used = send_data (id_mixdownslope(nt), T_prog(nt)%conversion*tend_mix(:,:,:), &
Time%model_time, rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if(id_neutral_rho_mixdown > 0 .or. id_wdian_rho_mixdown > 0) then
if(nt==index_temp) then
wrk1_v(:,:,:,1) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1_v(i,j,k,1) = tend_mix(i,j,k)*mixdownslope_mask(i,j)
enddo
enddo
enddo
endif
if(nt==index_salt) then
wrk1_v(:,:,:,2) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
wrk1_v(i,j,k,2) = tend_mix(i,j,k)*mixdownslope_mask(i,j)
enddo
enddo
enddo
endif
endif
if(debug_this_module) then
write(stdoutunit,*) ' '
write(stdoutunit,*) '==Global sums for tendency from ocean_mixdownslope_mod== '
call write_timestamp(Time%model_time)
do k=1,nk
area_k(:,:) = Grd%dat(:,:)*Grd%tmask(:,:,k)
tend_mix(:,:,k) = dtime*tend_mix(:,:,k)*area_k(:,:)*T_prog(nt)%conversion
write(stdoutunit,'(a,i2,a,i2,a,es24.17)') 'tend_mix(',nt,',',k,') = ',&
mpp_global_sum(Dom%domain2d,tend_mix(:,:,k))
enddo
write(stdoutunit,'(a,i2,a,es24.17)') &
'tend_mix(',nt,') = ',&
mpp_global_sum(Dom%domain2d,tend_mix(:,:,:), global_sum_flag)
endif
enddo ! nt-end for num_prog_tracers
! send diagnostics for update of neutral density from mixdownslope scheme
if (id_neutral_rho_mixdown > 0) then
used = send_data (id_neutral_rho_mixdown, &
Dens%drhodT(:,:,:)*wrk1_v(:,:,:,1)+Dens%drhodS(:,:,:)*wrk1_v(:,:,:,2) , &
Time%model_time,rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if (id_wdian_rho_mixdown > 0) then
wrk1(:,:,:) = 0.0
do k=1,nk
do j=jsc,jec
do i=isc,iec
temporary = Dens%drhodz_zt(i,j,k)*Thickness%rho_dzt(i,j,k,tau)
wrk1(i,j,k) = Grd%tmask(i,j,k) &
*(Dens%drhodT(i,j,k)*wrk1_v(i,j,k,1)+Dens%drhodS(i,j,k)*wrk1_v(i,j,k,2)) &
/(epsln+temporary)
enddo
enddo
enddo
used = send_data (id_wdian_rho_mixdown, wrk1(:,:,:), &
Time%model_time,rmask=Grd%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
if (debug_this_module) then
write(stdoutunit,*) ' '
write(stdoutunit,*) '==Global sums from ocean_mixdownslope_mod== '
call write_timestamp(Time%model_time)
write(stdoutunit,'(a,es24.17)') &
'rho_ex = ',&
mpp_global_sum(Mixdownslope_domain%domain2d,rho_ex(:,:,:), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'mass_ex = ',&
mpp_global_sum(Mixdownslope_domain%domain2d,mass_ex(:,:,:), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'press_ex = ',&
mpp_global_sum(Mixdownslope_domain%domain2d,press_ex(:,:,:), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'temp_ex = ',&
mpp_global_sum(Mixdownslope_domain%domain2d,temp_ex(:,:,:), global_sum_flag)
write(stdoutunit,'(a,es24.17)') &
'salt_ex = ',&
mpp_global_sum(Mixdownslope_domain%domain2d,salt_ex(:,:,:), global_sum_flag)
write(stdoutunit,'(a,i24)') &
'kmt_ex = ',&
mpp_global_sum(Mixdownslope_domain%domain2d,kmt_ex(:,:), global_sum_flag)
do m=1,4
write(stdoutunit,'(a,i24)') &
'kup_ex = ', &
mpp_global_sum(Mixdownslope_domain%domain2d,kup_ex(:,:,m), global_sum_flag)
do n=1,mixdownslope_npts
write(stdoutunit,'(a,es24.17)') &
'mixdownslope_frac = ', &
mpp_global_sum(Mixdownslope_domain%domain2d,mixdownslope_frac(:,:,n,m), global_sum_flag)
write(stdoutunit,'(a,i24)') &
'kdw_ex = ', &
mpp_global_sum(Mixdownslope_domain%domain2d,kdw_ex(:,:,n,m), global_sum_flag)
enddo
enddo
write(stdoutunit,*) ' '
write(stdoutunit,*) '==Global check sums from ocean_mixdownslope_mod== '
call write_timestamp(Time%model_time)
write (stdoutunit, *) 'chksum for rho_ex = ',mpp_chksum(rho_ex(isc:iec,jsc:jec,:))
write (stdoutunit, *) 'chksum for mass_ex = ',mpp_chksum(mass_ex(isc:iec,jsc:jec,:))
write (stdoutunit, *) 'chksum for press_ex = ',mpp_chksum(press_ex(isc:iec,jsc:jec,:))
write (stdoutunit, *) 'chksum for temp_ex = ',mpp_chksum(temp_ex(isc:iec,jsc:jec,:))
write (stdoutunit, *) 'chksum for salt_ex = ',mpp_chksum(salt_ex(isc:iec,jsc:jec,:))
write (stdoutunit, *) 'chksum for kmt_ex = ',mpp_chksum(kmt_ex(isc:iec,jsc:jec))
do m=1,4
write (stdoutunit, *) 'chksum for kup_ex = ',&
mpp_chksum(kup_ex(isc:iec,jsc:jec,m))
do n=1,mixdownslope_npts
write (stdoutunit, *) 'chksum for kdw_ex = ',&
mpp_chksum(kdw_ex(isc:iec,jsc:jec,n,m))
write (stdoutunit, *) 'chksum for mixdownslope_frac = ',&
mpp_chksum(mixdownslope_frac(isc:iec,jsc:jec,n,m))
enddo
enddo
endif
end subroutine mixdownslope
! NAME="mixdownslope"
end module ocean_mixdownslope_mod