!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! 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_residency_meta_mod !{
!
! Richard D. Slater
!
!
! Stephen M. Griffies
!
!
!
! Ocean residency module
!
!
!
! This module contains the meta definitions, subroutines and functions
! to use in the ocean residency package. These routines are used by all
! of the ocean residency modules and need to be defined in a separate
! module so that there are no circular references between the modules.
!
! For an overview of the ocean residency modules, see ocean_residency.F90.
!
!
!
! $Id: ocean_residency_meta.F90,v 16.0.2.1 2009/01/21 22:48:04 smg Exp $
!
!
! modules
!
use field_manager_mod, only: fm_string_len, fm_field_name_len
use fm_util_mod, only: fm_util_default_caller
use mpp_mod, only: stdout, mpp_error, FATAL
!
! force all variables to be "typed"
!
implicit none
!
! Set all variables to be private by default
private
!
! Public routines
!
public :: ocean_residency_get_instances
public :: ocean_residency_set_region_geog
public :: ocean_residency_set_region_2d
public :: ocean_residency_set_region_3d
!
! Private routines
!
!
! Public parameters
!
real, public, parameter :: secs_in_year_r = 1.0 / (86400.0 * 365.25)
real, public, parameter :: sec_per_day = 86400.0
!
! Private parameters
!
character(len=fm_field_name_len), parameter :: package_name = 'ocean_residency_meta'
character(len=48), parameter :: mod_name = 'ocean_residency_meta_mod'
character(len=48), parameter :: diag_name = 'ocean_residency_meta'
!
! Public types
!
type, public :: ocean_residency_instance_type
character(len=fm_field_name_len) :: name = ' '
integer :: tracer_index = 0
logical :: some_restore = .false.
logical :: some_fix = .false.
integer :: num_regions = 1
type(ocean_residency_region_type), dimension(:), pointer :: region => NULL()
real :: integrate_region_value = secs_in_year_r
integer, dimension(:,:,:), pointer :: index => NULL()
real, dimension(:,:,:), pointer :: mask => NULL()
logical :: union = .true.
integer :: id_restore_region = -1
integer :: id_change = -1
logical :: int_found = .false.
character(len=fm_field_name_len) :: int_module_name = ' '
real, dimension(:), pointer :: int_params => NULL()
logical, dimension(:), pointer :: int_flags => NULL()
character(len=fm_string_len), dimension(:), pointer :: int_strings => NULL()
end type ocean_residency_instance_type
type, public :: ocean_residency_region_type
character(len=fm_field_name_len) :: name = ' '
real, dimension(:,:,:), pointer :: mask => NULL()
real :: restore = 0.0
real :: restore_region_value = 0.0
logical :: swap = .false.
logical :: swap_module = .false.
logical :: found = .false.
real, dimension(:), pointer :: east_bnd => NULL()
real, dimension(:), pointer :: north_bnd => NULL()
real, dimension(:), pointer :: south_bnd => NULL()
real, dimension(:), pointer :: west_bnd => NULL()
real, dimension(:), pointer :: top_bnd => NULL()
real, dimension(:), pointer :: bottom_bnd => NULL()
integer :: num_geog_regions = 0
integer :: id_restore_region = -1
character(len=fm_field_name_len) :: module_name = ' '
real, dimension(:), pointer :: params => NULL()
logical, dimension(:), pointer :: flags => NULL()
character(len=fm_string_len), dimension(:), pointer :: strings => NULL()
end type ocean_residency_region_type
!
! Public variables
!
type(ocean_residency_instance_type), dimension(:), allocatable, save, public :: instance
integer, save, public :: num_instances = 0
!
! Private variables
!
!
! Interface definitions for overloaded routines
!
contains
!#######################################################################
!
!
!
! Return an array of instances which have the given module_name
! This is used by modules, such as the mixed layer module, to obtain
! a list of instances which use that module. Then, the module
! needs only to loop over those instance to perform its required tasks.
!
!
function ocean_residency_get_instances(module_name, integrand) &
result(array) !{
implicit none
!
!-----------------------------------------------------------------------
! Arguments
!-----------------------------------------------------------------------
!
character(len=*), intent(in) :: module_name
logical, intent(in), optional :: integrand
!
!-----------------------------------------------------------------------
! Result
!-----------------------------------------------------------------------
!
type(ocean_residency_instance_type), dimension(:), pointer :: array
!
! local parameters
!
character(len=64), parameter :: sub_name = 'ocean_residency_get_instances'
character(len=256), parameter :: error_header = &
'==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):'
character(len=256), parameter :: warn_header = &
'==>Warning from ' // trim(mod_name) // '(' // trim(sub_name) // '):'
character(len=256), parameter :: note_header = &
'==>Note from ' // trim(mod_name) // '(' // trim(sub_name) // '):'
!
!-----------------------------------------------------------------------
! local variables
!-----------------------------------------------------------------------
!
logical :: found
integer :: n
integer :: ind
integer :: nn
integer :: num_elements
logical :: do_integrand
character(len=fm_field_name_len) :: module_name_check
if (module_name .eq. ' ') then !{
call mpp_error(FATAL, trim(error_header) // ' No module name given')
endif !}
!
! check whether we want to look at modules for defining the region
! (integrand = f or not given) or whether we are looking to set
! the field we're integrating over to something other than time
! (integrand = t)
!
if (present(integrand)) then !{
do_integrand = integrand
else !}{
do_integrand = .false.
endif !}
!
! return the mask array pointer for tracer "tracer_index"
!
num_elements = 0
do n = 1, num_instances !{
if (do_integrand) then !{
if (module_name .eq. instance(n)%int_module_name) then
num_elements = num_elements + 1
endif
else !}{
do nn = 1, instance(n)%num_regions !{
if (module_name .eq. instance(n)%region(nn)%module_name) then
num_elements = num_elements + 1
exit
endif
enddo !} nn
endif !}
enddo !} n
if (num_elements .eq. 0) then !{
nullify(array)
elseif (num_elements .lt. 0) then !}{
call mpp_error(FATAL, trim(error_header) // ' Should never get here for ' // trim(module_name))
else !}{
allocate( array(num_elements) )
ind = 0
do n = 1, num_instances !{
found = .false.
if (do_integrand) then !{
if (module_name .eq. instance(n)%int_module_name) then
found = .true.
instance(n)%int_found = .true.
endif
else !}{
do nn = 1, instance(n)%num_regions !{
if (module_name .eq. instance(n)%region(nn)%module_name) then
found = .true.
instance(n)%region(nn)%found = .true.
endif
enddo !} nn
endif !}
if (found) then !{
ind = ind + 1
allocate ( array(ind)%region(instance(n)%num_regions) )
array(ind)%num_regions = instance(n)%num_regions
array(ind)%region => instance(n)%region
array(ind)%name = instance(n)%name
array(ind)%tracer_index = instance(n)%tracer_index
array(ind)%integrate_region_value = instance(n)%integrate_region_value
array(ind)%index => instance(n)%index
array(ind)%mask => instance(n)%mask
array(ind)%union = instance(n)%union
array(ind)%id_change = instance(n)%id_change
array(ind)%int_module_name = instance(n)%int_module_name
array(ind)%int_params => instance(n)%int_params
array(ind)%int_flags => instance(n)%int_flags
array(ind)%int_strings => instance(n)%int_strings
endif !}
enddo !} n
if (num_elements .ne. ind) then !{
call mpp_error(FATAL, trim(error_header) // ' Number of elements do not match for ' // trim(module_name))
endif !}
endif !}
return
end function ocean_residency_get_instances !}
! NAME="ocean_residency_get_instances"
!#######################################################################
!
!
!
! Given a 2-d field of depths, determine the grid cells which fall inside (above) and
! outside of that range, then set the mask appropriately, as to whether
! we are interested in those points inside (swap=false) or outside (swap=true) of the region.
!
! Note: if the output array is not initialized, then it is assumed that it is already
! filled with a region, and the resulting region will be the intersection of the
! existing region and the newly specified region.
!
! Arguments:
!
! Input:
!
! isd: low dimension of first index
! ied: high dimension of first index
! jsd: low dimension of second index
! jed: high dimension of second index
! nk: dimension of third index
! control: array of depths that specifies the region
! depth_zwt: depth of bottom of grid cell
! restore_region_value: if supplied, the value to assign to array in the user-specified
! regions (default: 0.0)
! integrate_region_value: if supplied, the value to assign to array outside of the
! user-specified regions (default: secs_in_year_r)
! swap: if supplied and true then change the invert the defined region (default: true)
! initialize: if supplied and true, then initialize the region to the integrate_region_value (default: false)
! caller: if supplied, use for traceback of error messages (default: fm_default_caller)
!
! Input/Output:
! array: 3-d array which will contain the restore_region- and integrate_region-
! values.
!
!
!
subroutine ocean_residency_set_region_2d(isd, ied, jsd, jed, nk, control, array, depth_zwt, &
restore_region_value, integrate_region_value, swap, initialize, caller) !{
implicit none
!
! Local parameters
!
character(len=48), parameter :: sub_name = 'ocean_residency_set_region_2d'
!
! arguments
!
integer, intent(in) :: isd
integer, intent(in) :: ied
integer, intent(in) :: jsd
integer, intent(in) :: jed
integer, intent(in) :: nk
real, intent(in), dimension(isd:,jsd:) :: control
real, dimension(isd:,jsd:,:), intent(inout) :: array
real, dimension(isd:,jsd:,:), intent(in) :: depth_zwt
real, intent(in), optional :: restore_region_value
real, intent(in), optional :: integrate_region_value
logical, intent(in), optional :: swap
logical, intent(in), optional :: initialize
character(len=*), intent(in), optional :: caller
!
! Local variables
!
character(len=256) :: error_header
character(len=256) :: warn_header
character(len=256) :: note_header
character(len=128) :: caller_str
integer :: i
integer :: j
integer :: k
real :: restore_region_value_use
real :: integrate_region_value_use
logical :: swap_use
logical :: assume_restore_region
!
! set the caller string and headers
!
if (present(caller)) then !{
caller_str = '[' // trim(caller) // ']'
else !}{
caller_str = fm_util_default_caller
endif !}
error_header = '==>Error from ' // trim(mod_name) // &
'(' // trim(sub_name) // ')' // trim(caller_str) // ':'
warn_header = '==>Warning from ' // trim(mod_name) // &
'(' // trim(sub_name) // ')' // trim(caller_str) // ':'
note_header = '==>Note from ' // trim(mod_name) // &
'(' // trim(sub_name) // ')' // trim(caller_str) // ':'
!
! check that the array is as long as depth dimension
!
if (size(array,3) .lt. nk) then
call mpp_error(FATAL, trim(error_header) // ' Dimension 3 of output array too small')
endif
!
! set some default values in case they have not been set in the argument list
!
if (present(restore_region_value)) then !{
restore_region_value_use = restore_region_value
else !}{
restore_region_value_use = 0.0
endif !}
if (present(integrate_region_value)) then !{
integrate_region_value_use = integrate_region_value
else !}{
integrate_region_value_use = secs_in_year_r
endif !}
if (present(swap)) then !{
swap_use = swap
else !}{
swap_use = .false.
endif !}
!
! If initialize is true, then set all values to integrate_region_value,
! otherwise, array values are assumed to be initialized with the appropriate
! distribution of in- and out-region-values and the final value will be
! the intersection between this value and that determined below (i.e., both values
! must be "restore_region_value" to stay "restore_region_value", otherwise it becomes "integrate_region_value")
! If we initialize, we need to set assume_restore_region to true so that when
! calculating the region below, the intersection will have values in it.
! However, we don't want to initialize to restore_region_value, or else the behavior
! when no regions are found will be to have all locations restore_region.
!
if (present(initialize)) then !{
assume_restore_region = initialize
if (initialize) then !{
do k = 1, nk !{
do j = jsd, jed !{
do i = isd, ied !{
array(i,j,k) = integrate_region_value_use
enddo !} i
enddo !} j
enddo !} k
endif !}
else !}{
assume_restore_region = .false.
endif !}
do j = jsd, jed !{
do i = isd, ied !{
if (((control(i,j) .gt. 0.0) .eqv. (.not. swap_use)) .and. &
(assume_restore_region .or. array(i,j,1) .eq. restore_region_value_use)) then !{
array(i,j,1) = restore_region_value_use
else !}{
array(i,j,1) = integrate_region_value_use
endif !}
enddo !} i
enddo !} j
do k = 2, nk !{
do j = jsd, jed !{
do i = isd, ied !{
if (((control(i,j) .gt. depth_zwt(i,j,k-1)) .eqv. (.not. swap_use)) .and. &
(assume_restore_region .or. array(i,j,k) .eq. restore_region_value_use)) then !{
array(i,j,k) = restore_region_value_use
else !}{
array(i,j,k) = integrate_region_value_use
endif !}
enddo !} i
enddo !} j
enddo !} k
return
end subroutine ocean_residency_set_region_2d !}
! NAME="ocean_residency_set_region_2d"
!#######################################################################
!
!
!
! Set up an array covering the model domain with a user-specified
! value, in user-specified regions. There are a given number of
! 3-d regions specified by the values south_bnd, north_bnd, west_bnd,
! east_bnd, top_bnd and bottom_bnd.
! The longitudes are for a cyclic domain, and if west_bnd and east_bnd
! are on opposite sides of the cut, the correct thing will
! be done. east_bnd is considered to be east of west_bnd, so if east_bnd is
! less than west_bnd, then the region east of east_bnd to the cut will be
! filled, and the region from the cut to west_bnd will be filled.
!
! If the bottom bound is less than or equal to zero, then the top model box
! will be chosen. If the top bound is greater than or equal to the maximum
! model depth, then the bottom box will be chosed. Otherwise, if the grid
! cell center depth falls between top bound and bottom bound, then those cells
! shall be chosen.
!
! For longitude and latitude, if the grid cell center lies within the
! rectabgle defined by (west_bnd,south_bnd) and (east_bnd,north_bnd), then
! the whole grid cell is inside the region.
!
! Arrays of coordinates may be specified for irregular regions.
! The final region is the union of the multiple sets
! of coordinates. If swap is true, then the inverse of the defined
! region will be set.
!
! Note: if the output array is not initialized, then it is assumed that it is already
! filled with a region, and the resulting region will be the intersection of the
! existing region and the newly specified region.
!
! Arguments:
!
! Input:
!
! isd: low dimension of first index
! ied: high dimension of first index
! jsd: low dimension of second index
! jed: high dimension of second index
! nk: dimension of third index
! grid_xt: array of coordinates in the x-direction (typically longitude)
! grid_yt: array of coordinates in the y-direction (typically latitude)
! max_depth: maximum depth of the model
! depth_zt: depth of center of grid cell
! depth_zwt: depth of bottom of grid cell
! num_geog_regions: number of user-specified regions which will be filled
! west_bnd_in: 1-d array of western (starting) longitudes for the regions
! east_bnd_in: 1-d array of eastern (ending) longitudes for the regions
! south_bnd: 1-d array of southern (starting) latitudes for the regions
! north_bnd: 1-d array of northern (ending) latitudes for the regions
! top_bnd: 1-d array of southern (starting) depths for the regions
! bottom_bnd: 1-d array of northern (ending) depths for the regions
! Note: if south_bnd >= north_bnd, then nothing is done
! for that region
! kmt: array of indices for bottom grid cells
! name: character variable used in informative messages
! restore_region_value: if supplied, the value to assign to array in the user-specified
! regions (default: 0.0)
! integrate_region_value: if supplied, the value to assign to array outside of the
! user-specified regions (default: secs_in_year_r)
! swap: if supplied and true then change the invert the defined region (default: true)
!
! Input/Output:
! array: 3-d array which will contain the restore_region- and integrate_region-
! values.
!
!
!
subroutine ocean_residency_set_region_geog(isd, ied, jsd, jed, nk, array, &
grid_xt, grid_yt, max_depth, depth_zt, depth_zwt, &
num_geog_regions, &
west_bnd_in, east_bnd_in, south_bnd, north_bnd, &
top_bnd, bottom_bnd, kmt, name, &
restore_region_value, integrate_region_value, swap, &
initialize, caller) !{
implicit none
!
!-----------------------------------------------------------------------
! Arguments
!-----------------------------------------------------------------------
!
integer, intent(in) :: isd
integer, intent(in) :: ied
integer, intent(in) :: jsd
integer, intent(in) :: jed
integer, intent(in) :: nk
real, dimension(isd:,jsd:,:), intent(inout) :: array
real, dimension(isd:,jsd:), intent(in) :: grid_xt
real, dimension(isd:,jsd:), intent(in) :: grid_yt
real, dimension(isd:,jsd:,:), intent(in) :: depth_zt
real, dimension(isd:,jsd:,:), intent(in) :: depth_zwt
real, intent(in) :: max_depth
integer, intent(in) :: num_geog_regions
real, dimension(:), intent(in) :: west_bnd_in
real, dimension(:), intent(in) :: east_bnd_in
real, dimension(:), intent(in) :: south_bnd
real, dimension(:), intent(in) :: north_bnd
real, dimension(:), intent(in) :: top_bnd
real, dimension(:), intent(in) :: bottom_bnd
integer, dimension(isd:,jsd:), intent(in) :: kmt
real, intent(in), optional :: restore_region_value
real, intent(in), optional :: integrate_region_value
logical, intent(in), optional :: swap
character*(*), intent(in) :: name
logical, intent(in), optional :: initialize
character(len=*), intent(in), optional :: caller
!
! local parameters
!
character(len=64), parameter :: sub_name = 'ocean_residency_set_region_geog'
!
! local variables
!
integer :: i
integer :: j
integer :: k
integer :: n
real, dimension(:), allocatable :: west_bnd
real, dimension(:), allocatable :: east_bnd
character(len=256) :: error_header
character(len=256) :: warn_header
character(len=256) :: note_header
character(len=128) :: caller_str
real :: restore_region_value_use
real :: integrate_region_value_use
logical :: swap_use
logical :: restore_region
logical :: assume_restore_region
!
! set the caller string and headers
!
if (present(caller)) then !{
caller_str = '[' // trim(caller) // ']'
else !}{
caller_str = fm_util_default_caller
endif !}
error_header = '==>Error from ' // trim(mod_name) // &
'(' // trim(sub_name) // ')' // trim(caller_str) // ':'
warn_header = '==>Warning from ' // trim(mod_name) // &
'(' // trim(sub_name) // ')' // trim(caller_str) // ':'
note_header = '==>Note from ' // trim(mod_name) // &
'(' // trim(sub_name) // ')' // trim(caller_str) // ':'
!
! set some default values in case they have not been set in the argument list
!
if (present(restore_region_value)) then !{
restore_region_value_use = restore_region_value
else !}{
restore_region_value_use = 0.0
endif !}
if (present(integrate_region_value)) then !{
integrate_region_value_use = integrate_region_value
else !}{
integrate_region_value_use = secs_in_year_r
endif !}
if (present(swap)) then !{
swap_use = swap
else !}{
swap_use = .false.
endif !}
!
! check whether any regions have been specified
!
if (num_geog_regions .le. 0) then !{
!
! no regions specified, return all points out of region
! unless swap is true
!
if (swap_use) then !{
do k = 1, nk !{
do j = jsd, jed !{
do i = isd, ied !{
array(i,j,k) = restore_region_value_use
enddo !} i
enddo !} j
enddo !} k
else !}{
do k = 1, nk !{
do j = jsd, jed !{
do i = isd, ied !{
array(i,j,k) = integrate_region_value_use
enddo !} i
enddo !} j
enddo !} k
endif !}
else !}{
!
! regions are specified, check that all array large enough
!
if (size(west_bnd_in) .lt. num_geog_regions) then !{
call mpp_error(FATAL, trim(error_header) // ' west_bnd_in array too small for ' // trim(name))
elseif (size(east_bnd_in) .lt. num_geog_regions) then !}{
call mpp_error(FATAL, trim(error_header) // ' east_bnd_in array too small for ' // trim(name))
elseif (size(north_bnd) .lt. num_geog_regions) then !}{
call mpp_error(FATAL, trim(error_header) // ' north_bnd array too small for ' // trim(name))
elseif (size(south_bnd) .lt. num_geog_regions) then !}{
call mpp_error(FATAL, trim(error_header) // ' south_bnd array too small for ' // trim(name))
elseif (size(top_bnd) .lt. num_geog_regions) then !}{
call mpp_error(FATAL, trim(error_header) // ' top_bnd array too small for ' // trim(name))
elseif (size(bottom_bnd) .lt. num_geog_regions) then !}{
call mpp_error(FATAL, trim(error_header) // ' bottom_bnd array too small for ' // trim(name))
endif !}
!
! loop over the regions, checking for errors
!
do n = 1, num_geog_regions !{
if (north_bnd(n) .lt. south_bnd(n)) then !{
call mpp_error(FATAL, trim(error_header) // &
' North bound less than south bound for ' // trim(name))
elseif (top_bnd(n) .gt. bottom_bnd(n)) then !}{
call mpp_error(FATAL, trim(error_header) // &
' Top depth deeper than bottom depth for ' // trim(name))
endif !}
enddo !} n
!
! Find all boxes where the center lies in the
! 3-d rectangular region
!
!
! If initialize is true, then set all values to integrate_region_value,
! otherwise, array values are assumed to be initialized with the appropriate
! distribution of in- and out-region-values and the final value will be
! the intersection between this value and that determined below (i.e., both values
! must be "restore_region_value" to stay "restore_region_value", otherwise it becomes "integrate_region_value")
! If we initialize, we need to set assume_restore_region to true so that when
! calculating the region below, the intersection will have values in it.
! However, we don't want to initialize to restore_region_value, or else the behavior
! when no regions are found will be to have all locations restore_region.
!
if (present(initialize)) then !{
assume_restore_region = initialize
if (initialize) then !{
do k = 1, nk !{
do j = jsd, jed !{
do i = isd, ied !{
array(i,j,k) = integrate_region_value_use
enddo !} i
enddo !} j
enddo !} k
endif !}
else !}{
assume_restore_region = .false.
endif !}
!
! save the longitudes in case they need to be modified
!
allocate( west_bnd(num_geog_regions) )
allocate( east_bnd(num_geog_regions) )
do n = 1, num_geog_regions !{
west_bnd(n) = west_bnd_in(n)
east_bnd(n) = east_bnd_in(n)
!
! make sure that the longitudes are in the range [0,360]
!
do while (west_bnd(n) .gt. 360.0) !{
west_bnd(n) = west_bnd(n) - 360.0
enddo !}
do while (west_bnd(n) .lt. 0.0) !{
west_bnd(n) = west_bnd(n) + 360.0
enddo !}
do while (east_bnd(n) .gt. 360.0) !{
east_bnd(n) = east_bnd(n) - 360.0
enddo !}
do while (east_bnd(n) .lt. 0.0) !{
east_bnd(n) = east_bnd(n) + 360.0
enddo !}
enddo !} n
do k = 1, nk !{
do j = jsd, jed !{
do i = isd, ied !{
if (k .le. kmt(i,j)) then !{
restore_region = .false.
do n = 1, num_geog_regions !{
if ((((top_bnd(n) .lt. depth_zt(i,j,k) .and. &
depth_zt(i,j,k) .le. bottom_bnd(n)) .or. &
(bottom_bnd(n) .le. 0.0 .and. k .eq. 1) .or. &
(top_bnd(n) .ge. max_depth .and. k .eq. kmt(i,j)) .or. &
(top_bnd(n) .lt. 0.0 .and. &
depth_zt(i,j,k) .ge. &
depth_zwt(i,j,kmt(i,j)) + top_bnd(n))) .and. &
north_bnd(n) .gt. grid_yt(i,j) .and. &
grid_yt(i,j) .ge. south_bnd(n) .and. &
lon_between(grid_xt(i,j), west_bnd(n), east_bnd(n)))) then !{
restore_region = .true.
endif !}
enddo !} n
restore_region = restore_region .eqv. (.not. swap_use)
if ((assume_restore_region .or. array(i,j,k) .eq. restore_region_value_use) .and. restore_region) then !{
array(i,j,k) = restore_region_value_use
else !}{
array(i,j,k) = integrate_region_value_use
endif !}
endif !}
enddo !} i
enddo !} j
enddo !} k
!
! Clean up
!
deallocate( west_bnd )
deallocate( east_bnd )
endif !}
return
contains
!
! Return true if w <= x_in < e, taking into account the
! periodicity of longitude.
!
! x_in value to test
! w west longitude of boundary
! e east longitude of boundary
!
function lon_between(x_in, w, e) !{
implicit none
!
! function definition
!
logical :: lon_between
!
! arguments
!
real, intent(in) :: x_in
real, intent(in) :: w
real, intent(in) :: e
!
! local variables
!
real :: x
!
! Save input values so we may modify them safely
!
x = x_in
!
! make sure that all longitudes are in the range [0,360]
!
do while (x .gt. 360.0) !{
x = x - 360.0
enddo !}
do while (x .lt. 0.0) !{
x = x + 360.0
enddo !}
if (w .gt. e) then !{
lon_between = w .le. x .or. x .lt. e
else !}{
lon_between = w .le. x .and. x .lt. e
endif !}
return
end function lon_between !}
end subroutine ocean_residency_set_region_geog !}
! NAME="ocean_residency_set_region_geog"
!#######################################################################
!
!
!
! Set up an array where the a grid box is in the region if the value
! of the specified property (temperature, say) is within the given bounds.
! Multiple values for the range may be given, and the resulting mask
! will be the union of the multiple regions. If swap is true, then the
! inverse of the selected region will be set.
!
! Note: if the output array is not initialized, then it is assumed that it is already
! filled with a region, and the resulting region will be the intersection of the
! existing region and the newly specified region.
!
! Arguments:
!
! Input:
!
! isd: low dimension of first index
! ied: high dimension of first index
! jsd: low dimension of second index
! jed: high dimension of second index
! nk: dimension of third index
! num_geog_regions: number of user-specified regions which will be filled
! bounds: 1-d array of pairs of bounding values. The first value in
! the pair must be less than the second value
! kmt: array of indices for bottom grid cells
! name: character variable used in informative messages
! restore_region_value: if supplied, the value to assign to array in the user-specified
! regions (default: 0.0)
! integrate_region_value: if supplied, the value to assign to array outside of the
! user-specified regions (default: secs_in_year_r)
! swap: if supplied and true then change the invert the defined region (default: true)
! initialize: if supplied and true, then initialize the region to the integrate_region_value (default: false)
! caller: if supplied, use for traceback of error messages (default: fm_default_caller)
!
! Input/Output:
! array: 3-d array which will contain the restore_region- and integrate_region-
! values.
!
!
!
subroutine ocean_residency_set_region_3d(isd, ied, jsd, jed, nk, array, &
variable, bounds, kmt, name, &
restore_region_value, integrate_region_value, swap, &
initialize, caller) !{
implicit none
!
!-----------------------------------------------------------------------
! Arguments
!-----------------------------------------------------------------------
!
integer, intent(in) :: isd
integer, intent(in) :: ied
integer, intent(in) :: jsd
integer, intent(in) :: jed
integer, intent(in) :: nk
real, dimension(isd:,jsd:,:), intent(inout) :: array
real, dimension(isd:,jsd:,:), intent(in) :: variable
real, dimension(:), intent(in) :: bounds
integer, dimension(isd:,jsd:), intent(in) :: kmt
character*(*), intent(in) :: name
real, intent(in), optional :: restore_region_value
real, intent(in), optional :: integrate_region_value
logical, intent(in), optional :: swap
logical, intent(in), optional :: initialize
character(len=*), intent(in), optional :: caller
!
! local parameters
!
character(len=64), parameter :: sub_name = 'ocean_residency_set_region_3d'
!
! local variables
!
character(len=132) :: error_string
integer :: i
integer :: j
integer :: k
integer :: n
logical :: restore_region
character(len=256) :: error_header
character(len=256) :: warn_header
character(len=256) :: note_header
character(len=128) :: caller_str
real :: restore_region_value_use
real :: integrate_region_value_use
logical :: swap_use
logical :: assume_restore_region
!
! set the caller string and headers
!
if (present(caller)) then !{
caller_str = '[' // trim(caller) // ']'
else !}{
caller_str = fm_util_default_caller
endif !}
error_header = '==>Error from ' // trim(mod_name) // &
'(' // trim(sub_name) // ')' // trim(caller_str) // ':'
warn_header = '==>Warning from ' // trim(mod_name) // &
'(' // trim(sub_name) // ')' // trim(caller_str) // ':'
note_header = '==>Note from ' // trim(mod_name) // &
'(' // trim(sub_name) // ')' // trim(caller_str) // ':'
!
! check that an even number of bounds have been given
!
if (mod(size(bounds),2) .ne. 0) then
call mpp_error(FATAL, trim(error_header) // ' Odd number of bounds given')
endif
!
! check that the bounds are correctly set (i.e., lesser value first)
!
do n = 1, size(bounds), 2 !{
if (bounds(n+1) .le. bounds(n)) then !{
write (error_string,*) ' Bounds out of order: ', bounds(n+1), ' <= ', bounds(n)
call mpp_error(FATAL, trim(error_header) // trim(error_string))
endif !}
enddo !} n
!
! set some default values in case they have not been set in the argument list
!
if (present(restore_region_value)) then !{
restore_region_value_use = restore_region_value
else !}{
restore_region_value_use = 0.0
endif !}
if (present(integrate_region_value)) then !{
integrate_region_value_use = integrate_region_value
else !}{
integrate_region_value_use = secs_in_year_r
endif !}
if (present(swap)) then !{
swap_use = swap
else !}{
swap_use = .false.
endif !}
!
! If initialize is true, then set all values to integrate_region_value,
! otherwise, array values are assumed to be initialized with the appropriate
! distribution of in- and out-region-values and the final value will be
! the intersection between this value and that determined below (i.e., both values
! must be "restore_region_value" to stay "restore_region_value", otherwise it becomes "integrate_region_value")
! If we initialize, we need to set assume_restore_region to true so that when
! calculating the region below, the intersection will have values in it.
! However, we don't want to initialize to restore_region_value, or else the behavior
! when no regions are found will be to have all locations restore_region.
!
if (present(initialize)) then !{
assume_restore_region = initialize
if (initialize) then !{
do k = 1, nk !{
do j = jsd, jed !{
do i = isd, ied !{
array(i,j,k) = integrate_region_value_use
enddo !} i
enddo !} j
enddo !} k
endif !}
else !}{
assume_restore_region = .false.
endif !}
!
! find all boxes where the center lies in the region
!
do k = 1, nk !{
do j = jsd, jed !{
do i = isd, ied !{
if (k .le. kmt(i,j)) then !{
restore_region = .false.
do n = 1, size(bounds), 2 !{
if ((bounds(n+1) .ge. variable(i,j,k) .and. variable(i,j,k) .gt. bounds(n))) then !{
restore_region = .true.
exit
endif !}
enddo !} n
restore_region = restore_region .eqv. (.not. swap_use)
if ((assume_restore_region .or. array(i,j,k) .eq. restore_region_value_use) .and. restore_region) then !{
array(i,j,k) = restore_region_value_use
else !}{
array(i,j,k) = integrate_region_value_use
endif !}
endif !}
enddo !} i
enddo !} j
enddo !} k
return
end subroutine ocean_residency_set_region_3d !}
! NAME="ocean_residency_set_region_3d"
end module ocean_residency_meta_mod !}