!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! 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_mod !{
!
! Richard D. Slater
!
!
! Stephen M. Griffies
!
!
!
! Ocean residency module
!
!
!
!
! This module is a superset of the ocean age tracer module. It may
! be used to reproduce the age tracer, but it can also do much
! more. Unlike the ocean age tracer module, here you may specify
! a 3-d field by specifying a series of "rectangular prisms".
! The grid cells which occupy this field may vary with time due
! to the variations in the thickness of the surface layer.
!
! You may also specify a 3-d field by choosing one of three mixed-layers
! (KPP, density-derived or buoyancy-derived). You may also specify a
! region based on a a range of a prognostic or diagnostic ocean
! tracer (such as all points with a temperature of 10-20 degrees).
!
! You may specify either the mixed-layer method or the tracer range along
! with the geographic specification. If multiple methods are used, then
! the resultant field is the intersection of the two methods. In fact, the
! geographic method is always in use, and defaults to the whole ocean.
!
! There is an option of using the inverse of any method. This is sometimes
! easier to use than explicitly specifying the inverse. For instance,
! to get the temperatures outside of the range 10-20 degrees, one could
! specify the 10-20 degree range and set to true, or, specify two
! different regions, one less than 10 degrees and one greater than 20 degrees.
!
! By default, the values inside the specified region are set to a
! specified value each time-step (default is 0), and outside of this region
! the field is integrated over time in units of years (integrand is
! 1/(365.25*86400) by default, but this value can be changed. The inner region
! can be forced to be 0 (or the user-specified inner value) at each time-step
! (the default), or it can be restored to this value at a user-specified
! rate (given in days), or, it can be left alone (not integrated or set to
! any value).
!
! Finally, since we are just integrating a specified field by a constant value
! (by default), it is simple to make that field a simple function of another
! tracer. This is the final option. You may specify a prognostic or diagnostic
! variable to be the integrand, and have it scaled by a constant. One case
! which has been done is to integrate irradiance in the mixed-layer.
!
! This module is split into several different modules. This particular module
! is the only one called from outside, and it makes appropriate calls to
! the other modules to implement its features. This was done so that it would
! be easier to expand the features of the ocean residency, without cluttering
! the code too much. Ideally, this could all be done with classes, but there
! is currently no support for classes in Fortran (pre-F2003). The current
! modules are:
!
! ocean_residency.F90 This module, the control center, and also does the
! integrations and resetting in the inner regions.
! ocean_residency_meta.F90 Has specifications of the ocean_residency type
! and other utility routines. This is separate
! to stop circular references
! ocean_residency_ml.F90 Handles setting mixed-layer masks
! ocean_residency_range.F90 Handles setting the masks for the tracer ranges
! ocean_residency_integrand.F90 Handles setting the integrand to a non-constant
! value
!
! field_table namelist inputs:
!
! The following 6 arrays specify the bounds of boxes, the intersection of which
! will define the region where the integrand is set to restore_region_value. ALl arrays
! must have the same number of elements. If none are specified, then no region will
! be used, unless swap is true, which is a quick way of selecting the entire ocean.
! This may be useful when one really wants to use one of the other methods of selecting,
! such as mixed layer depth, or tracer range.
!
! Longitudes will be shifted to lie in the range 0-360 degrees. If the eastern side is
! greater than the western, then the selected region will consist of those grid cells from
! the eastern value to 360, and 0 to the western value. If the northern value is less than
! the southern value, or if the top depth is greater than the bottom depth, then
! an error occurs, and the model stops.
!
! Three special cases exist for the depth (which currently must be in meters).
! If the bottom value is less than or equal to zero, then the top box is selected.
! If the top value is greater than the maximum depth, then the bottom box is selected.
! If the top value is negative, then grid cells within "the absolute value of the top value"
! from the bottom are selected.
!
! Note that the geographic specification is used for all residency tracers, so either
! values for these 6 arrays must be given so as to select a region, or else swap must
! be set to true and no arrays given.
!
! east_bnd: array of boundary points of the eastern side of the box,
! in degrees longitude (default: NULL)
! north_bnd: array of boundary points of the northern side of the box,
! in degrees latitude (default: NULL)
! south_bnd: array of boundary points of the southern side of the box,
! in degrees latitude (default: NULL)
! west_bnd: array of boundary points of the western side of the box,
! in degrees longitude (default: NULL)
! top_bnd: array of boundary points of the top side of the box,
! in meters (default: NULL)
! bottom_bnd: array of boundary points of the bottom side of the box,
! in meters (default: NULL)
!
! swap: if true, then select the inverse of the specified geographic region,
! otherwise, just use the specified region (default: false)
! restore: restoring value for values in the defined regions
! negative => do nothing, 0 => force to integrate_region_value,
! positive => time scale in days to force to integrate_region_value
! (default: 0.0)
! restore_region_value: value to set the mask to for grid cells within
! the specified region (default: 0.0)
! integrate_region_value: value to set the mask to for grid cells outside
! the specified region (default: secs_in_year_r)
! swap_module: if true, then select the inverse of the region from the specified module,
! otherwise, just use the specified region (default: false)
! module_name: if set, then it will be used to select the alternate
! method of selecting the region (default: ' ' -- only geographic
! selection is used)
!
! For the following arrays, see the different extra modules for required and
! possible values for the module selected.
!
! params: an array of real parameters which may be used for the method
! being used for selection (default: NULL)
! flags: an array of real parameters which may be used for the method
! being used for selection (default: NULL)
! strings: an array of real parameters which may be used for the method
! being used for selection (default: NULL)
!
! For the following arrays, see the different extra modules for required and
! possible values for the module selected to specify the integrand.
! int_module_name: if set, then the name of the module used to set
! the integrand (default: ' ' -- integrate time, in years)
! int_params: an array of real parameters which may be used for the integrand
! (default: NULL)
! int_flags: an array of real parameters which may be used for the integrand
! (default: NULL)
! int_strings: an array of real parameters which may be used for the integrand
! (default: NULL)
!
!----------------------------------------------------------------------------------------
!
! Sample field table entries:
! ---------------------------
!
! "tracer_packages","ocean_mod","ocean_residency"
! names = age_surface, age_bottom_inv, kppbl_nil, kppbl_14d, kppbl_frc, kppbl_irr_14d, temp_15_20
! horizontal-advection-scheme = mdppm
! vertical-advection-scheme = mdppm
! units = yr
! min_tracer_limit=0.0
! /
!
! This is the same as the old age tracer with all surface
! values forced to zero
!
! "namelists","ocean_mod","ocean_residency/age_surface"
! south_bnd = -90.0
! north_bnd = 90.0
! west_bnd = 0.0
! east_bnd = 360.0
! top_bnd = 0.0
! bottom_bnd = 0.0
! /
!
! This integrates the age in the bottom box and forces to
! zero everywhere else (note that swap is true)
!
! "namelists","ocean_mod","ocean_residency/age_bottom_inv"
! south_bnd = -90.0
! north_bnd = 90.0
! west_bnd = 0.0
! east_bnd = 360.0
! top_bnd = 10000.0
! bottom_bnd = 10000.0
! swap = t
! /
!
! This integrates the age of the water in the KPP
! boundary layer, but lets the age outside of this region
! keep its value until it again mixes with the boundary layer.
! The module_name is set, and strings is set to pick the
! type of mixed layer desired
! (note that the global geographic area is explicitly specified,
! and that the restoring time scale is negative)
!
! "namelists","ocean_mod","ocean_residency/kppbl_nil"
! south_bnd = -90.0
! north_bnd = 90.0
! west_bnd = 0.0
! east_bnd = 360.0
! top_bnd = 0.0
! bottom_bnd = 10000.0
! restore = -1.0
! module_name = ocean_residency_ml
! strings = kpp_bl
! swap_module = t
! /
!
! This is the same as above, but forces the age outside
! the boundary layer to 0 with a 14 day time scale
! (note that the global geographic region is specified by
! setting swap to true, also note the value of restore)
!
! "namelists","ocean_mod","ocean_residency/kppbl_14d"
! swap = t
! restore = 14.0
! module_name = ocean_residency_ml
! strings = kpp_bl
! swap_module = t
! /
!
! This is the same as above, but forces age to zero outside
! the boundary layer (note that restore did not need to
! be explicitly specified, as the default is zero)
!
! "namelists","ocean_mod","ocean_residency/kppbl_frc"
! swap = t
! restore = 0.0
! module_name = ocean_residency_ml
! strings = kpp_bl
! swap_module = t
! /
!
! The following integrates irradiance in the boundary
! layer (note that the units needed to be changed for
! netCDF output purposes)
!
! "prog_tracers","ocean_mod","residency_kppbl_irr_14d"
! units = W-yr/m^2
! /
!
! "namelists","ocean_mod","ocean_residency/kppbl_irr_14d"
! swap = t
! restore = 14.0
! module_name = ocean_residency_ml
! strings = kpp_bl
! swap_module = t
! int_module_name = ocean_residency_integrand
! int_strings = irr
! /
!
! This specifies the region as the area with
! a temperature range of between 15 and 20 degrees
! (note that the params holds the variable
! range)
!
! "namelists","ocean_mod","ocean_residency/temp_15_20"
! swap = t
! module_name = ocean_residency_range
! strings = tracer_range, temp
! params = 15.0, 20.0
! /
!
!
!
! $Id: ocean_residency.F90,v 16.0.2.2.82.1 2009/10/10 00:43:01 nnz Exp $
!
!
! modules
!
use time_manager_mod, only: time_type
use field_manager_mod, only: fm_string_len, fm_path_name_len, fm_field_name_len
use field_manager_mod, only: fm_get_length, fm_get_value, fm_new_value
use mpp_mod, only: stdout, mpp_error, FATAL, mpp_pe, mpp_root_pe
use diag_manager_mod, only: register_diag_field, send_data
use ocean_tpm_util_mod, only: otpm_set_tracer_package, otpm_set_prog_tracer, otpm_set_diag_tracer
use fm_util_mod, only: fm_util_set_value
use fm_util_mod, only: fm_util_start_namelist, fm_util_end_namelist
use fm_util_mod, only: fm_util_get_string, fm_util_get_string_array
use fm_util_mod, only: fm_util_get_logical, fm_util_get_logical_array
use fm_util_mod, only: fm_util_get_integer, fm_util_get_real, fm_util_get_real_array
use fm_util_mod, only: fm_util_check_for_bad_fields
use ocean_types_mod, only: ocean_prog_tracer_type, ocean_diag_tracer_type
use ocean_types_mod, only: ocean_thickness_type, ocean_time_type, ocean_density_type
use ocean_residency_meta_mod, only: ocean_residency_set_region_geog
use ocean_residency_meta_mod, only: instance, num_instances, secs_in_year_r, sec_per_day
use ocean_residency_ml_mod, only: do_ocean_residency_ml
use ocean_residency_ml_mod, only: ocean_residency_ml_start, ocean_residency_ml_source
use ocean_residency_range_mod, only: do_ocean_residency_range
use ocean_residency_range_mod, only: ocean_residency_range_start, ocean_residency_range_source
use ocean_residency_integrand_mod, only: do_ocean_residency_integrand
use ocean_residency_integrand_mod, only: ocean_residency_integrand_start, ocean_residency_integrand_source
!
! force all variables to be "typed"
!
implicit none
!
! Set all variables to be private by default
private
!
! Public routines
!
public :: ocean_residency_init
public :: ocean_residency_source
public :: ocean_residency_start
public :: ocean_residency_tracer
!
! Private routines
!
!
! Public parameters
!
!
! Private parameters
!
character(len=fm_field_name_len), parameter :: package_name = 'ocean_residency'
character(len=48), parameter :: mod_name = 'ocean_residency_mod'
character(len=48), parameter :: diag_name = 'ocean_residency'
character(len=fm_string_len), parameter :: default_restart_file = 'ocean_residency.res.nc'
!
! Public variables
!
logical, public :: do_ocean_residency
!
! Public types
!
!
! Private variables
!
integer :: package_index
contains
!#######################################################################
!
!
!
! Set up any extra fields needed by the tracer packages
!
!
subroutine ocean_residency_init !{
implicit none
!
!-----------------------------------------------------------------------
! Arguments
!-----------------------------------------------------------------------
!
!
! local parameters
!
character(len=64), parameter :: sub_name = 'ocean_residency_init'
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
!-----------------------------------------------------------------------
!
integer :: n
integer :: nn
integer :: num_regions
character(len=fm_field_name_len) :: name
character(len=fm_path_name_len) :: path_to_names
character(len=fm_field_name_len+1) :: suffix
character(len=fm_field_name_len+3) :: long_suffix
character(len=256) :: caller_str
character(len=fm_string_len), pointer, dimension(:) :: good_list
integer :: stdoutunit
stdoutunit=stdout()
!
! Initialize the ocean residency package
!
package_index = otpm_set_tracer_package(package_name, &
units = 'yr', flux_units = 'm', &
min_tracer_limit = 0.0, max_tracer_limit = 1.0e+20, &
restart_file = default_restart_file, &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
!
! Check whether to use this package
!
path_to_names = '/ocean_mod/tracer_packages/' // trim(package_name) // '/names'
num_instances = fm_get_length(path_to_names)
if (num_instances .lt. 0) then !{
call mpp_error(FATAL, trim(error_header) // ' Could not get number of instances')
endif !}
!
! Check some things
!
write (stdoutunit,*)
if (num_instances .eq. 0) then !{
write (stdoutunit,*) &
trim(note_header), ' No instances'
do_ocean_residency = .false.
else !}{
write (stdoutunit,*) trim(note_header), ' ', num_instances, ' instances'
do_ocean_residency = .true.
endif !}
!
! Return if we don't want to use this package
!
if (.not. do_ocean_residency) then !{
return
endif !}
!
! allocate the instance array
!
allocate ( instance(num_instances) )
!
! loop over the names, saving them into the instance array
!
do n = 1, num_instances !{
if (fm_get_value(path_to_names, name, index = n)) then !{
instance(n)%name = name
else !}{
write (name,*) n
call mpp_error(FATAL, trim(error_header) // &
' Bad field name for index ' // trim(name))
endif !}
enddo !} n
do n = 1, num_instances !{
!
! determine the tracer name for this instance
!
name = instance(n)%name
if (name(1:1) .eq. '_') then !{
suffix = ' '
long_suffix = ' '
else !}{
suffix = '_' // name
long_suffix = ' (' // trim(name) // ')'
endif !}
instance(n)%tracer_index = otpm_set_prog_tracer('residency' // trim(suffix), package_name, &
longname = 'Residency' // trim(long_suffix), units = 'yr', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
enddo !} n
!
! set up the density diagnostic tracers
!
n = otpm_set_diag_tracer('rho_in_situ', longname = 'In situ density', units = 'kg/m^3', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
n = otpm_set_diag_tracer('rho_neutral', longname = 'Neutral density', units = 'kg/m^3', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
n = otpm_set_diag_tracer('rho_pot_0', longname = 'Potential density-0', units = 'kg/m^3', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
n = otpm_set_diag_tracer('rho_pot_1', longname = 'Potential density-1', units = 'kg/m^3', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
n = otpm_set_diag_tracer('rho_pot_2', longname = 'Potential density-2', units = 'kg/m^3', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
n = otpm_set_diag_tracer('rho_pot_3', longname = 'Potential density-3', units = 'kg/m^3', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
n = otpm_set_diag_tracer('rho_pot_4', longname = 'Potential density-4', units = 'kg/m^3', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
!
!-----------------------------------------------------------------------
! Set up the instance residency namelists
!-----------------------------------------------------------------------
!
!
! Add the package name to the list of good namelists, to be used
! later for a consistency check
!
if (fm_new_value('/ocean_mod/GOOD/good_namelists', package_name, append = .true.) .le. 0) then !{
call mpp_error(FATAL, trim(error_header) // &
' Could not add ' // trim(package_name) // ' to "good_namelists" list')
endif !}
caller_str = trim(mod_name) // '(' // trim(sub_name) // ')'
do n = 1, num_instances !{
!
! create the instance namelist
!
call fm_util_start_namelist(package_name, instance(n)%name, caller = caller_str, no_overwrite = .true., &
check = .true.)
call fm_util_set_value('num_regions', 1)
call fm_util_set_value('integrate_region_value', secs_in_year_r)
call fm_util_set_value('union', .true.)
call fm_util_set_value('int_module_name', ' ')
call fm_util_set_value('int_params', 0.0, index = 0)
call fm_util_set_value('int_flags', .false., index = 0)
call fm_util_set_value('int_strings', ' ', index = 0)
!
! get the number of regions so that we may set up the extra namelists
!
num_regions = fm_util_get_integer('num_regions', scalar = .true.)
!
! create namelists for each region
!
if (num_regions .le. 0) then
call mpp_error(FATAL,trim(error_header) // ' num_regions is non-positive for instance ' // &
trim(instance(n)%name) // ' in package ' // trim(package_name))
endif
do nn = 1, num_regions !{
if (num_regions .eq. 1) then
suffix = ' '
else
write (long_suffix,*) nn
suffix = '_' // long_suffix(scan(long_suffix,'0123456789'):)
endif
call fm_util_set_value('restore' // suffix, 0.0)
call fm_util_set_value('restore_region_value' // suffix, 0.0)
call fm_util_set_value('swap' // suffix, .false.)
call fm_util_set_value('swap_module' // suffix, .false.)
call fm_util_set_value('west_bnd' // suffix, 0.0, index = 0)
call fm_util_set_value('east_bnd' // suffix, 0.0, index = 0)
call fm_util_set_value('south_bnd' // suffix, 0.0, index = 0)
call fm_util_set_value('north_bnd' // suffix, 0.0, index = 0)
call fm_util_set_value('top_bnd' // suffix, 0.0, index = 0)
call fm_util_set_value('bottom_bnd' // suffix, 0.0, index = 0)
call fm_util_set_value('module_name' // suffix, ' ')
call fm_util_set_value('params' // suffix, 0.0, index = 0)
call fm_util_set_value('flags' // suffix, .false., index = 0)
call fm_util_set_value('strings' // suffix, ' ', index = 0)
enddo !} nn
call fm_util_end_namelist(package_name, instance(n)%name, check = .true., caller = caller_str)
enddo !} n
!
! Check for any errors in the number of fields in the namelists for this package
!
good_list => fm_util_get_string_array('/ocean_mod/GOOD/namelists/' // trim(package_name) // '/good_values', &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
if (associated(good_list)) then !{
call fm_util_check_for_bad_fields('/ocean_mod/namelists/' // trim(package_name), good_list, &
caller = trim(mod_name) // '(' // trim(sub_name) // ')')
deallocate(good_list)
else !}{
call mpp_error(FATAL,trim(error_header) // ' Empty "' // trim(package_name) // '" list')
endif !}
return
end subroutine ocean_residency_init !}
! NAME="ocean_residency_init"
!#######################################################################
!
!
!
! Calculate the source arrays for the tracer packages
!
!
subroutine ocean_residency_source(isc, iec, jsc, jec, isd, ied, jsd, jed, nk, &
T_prog, T_diag, Time, Thickness, Dens, grid_xt, grid_yt, grid_zw, &
grid_tmask, grid_kmt, hblt_depth)
!
! modules
!
implicit none
!
!-----------------------------------------------------------------------
! Arguments
!-----------------------------------------------------------------------
!
integer, intent(in) :: isc
integer, intent(in) :: iec
integer, intent(in) :: jsc
integer, intent(in) :: jec
integer, intent(in) :: isd
integer, intent(in) :: ied
integer, intent(in) :: jsd
integer, intent(in) :: jed
integer, intent(in) :: nk
type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog
type(ocean_diag_tracer_type), dimension(:), intent(inout) :: T_diag
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) :: grid_xt
real, dimension(isd:,jsd:), intent(in) :: grid_yt
real, dimension(:), intent(in) :: grid_zw
real, dimension(isd:,jsd:,:), intent(in) :: grid_tmask
integer, dimension(isd:,jsd:), intent(in) :: grid_kmt
real, intent(in), dimension(isd:,jsd:) :: hblt_depth
!
! local parameters
!
character(len=64), parameter :: sub_name = 'ocean_residency_source'
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
!
integer :: i
integer :: ind
integer :: j
integer :: k
integer :: n
integer :: nn
logical :: used
logical :: good
integer :: stdoutunit
stdoutunit=stdout()
!
! set the source values for the residency tracers
!
do n = 1, num_instances !{
!
! set the values via the input values
!
do nn = 1, instance(n)%num_regions !{
call ocean_residency_set_region_geog(isd, ied, jsd, jed, nk, instance(n)%region(nn)%mask, &
grid_xt, grid_yt, grid_zw(nk), Thickness%depth_zt, Thickness%depth_zwt, &
instance(n)%region(nn)%num_geog_regions, &
instance(n)%region(nn)%west_bnd, instance(n)%region(nn)%east_bnd, &
instance(n)%region(nn)%south_bnd, instance(n)%region(nn)%north_bnd, &
instance(n)%region(nn)%top_bnd, instance(n)%region(nn)%bottom_bnd, &
grid_kmt, t_prog(instance(n)%tracer_index)%name, &
restore_region_value = 1.0, integrate_region_value = 0.0, &
swap = instance(n)%region(nn)%swap, &
initialize = .true.)
enddo !} nn
enddo !} n
!
! set the source values for the residency tracers which need to be set in other
! modules
!
if (do_ocean_residency_ml) then !{
call ocean_residency_ml_source(isd, ied, jsd, jed, nk, &
t_prog, time, thickness, dens, Thickness%depth_zwt, hblt_depth)
endif !}
if (do_ocean_residency_range) then !{
call ocean_residency_range_source(isd, ied, jsd, jed, nk, Time%taum1, t_prog, t_diag, grid_kmt)
endif !}
!
! save out the mask fields for each region
!
do n = 1, num_instances !{
do nn = 1, instance(n)%num_regions !{
do k = 1, nk
do j = jsd, jed
do i = isd, ied
instance(n)%region(nn)%mask(i,j,k) = instance(n)%region(nn)%mask(i,j,k) * grid_tmask(i,j,k)
enddo
enddo
enddo
if (instance(n)%region(nn)%id_restore_region .gt. 0) then !{
used = send_data(instance(n)%region(nn)%id_restore_region, &
instance(n)%region(nn)%mask(:,:,:), &
Time%model_time, rmask = grid_tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif !}
enddo !} nn
enddo !} n
!
! merge the mask fields into one field
!
do n = 1, num_instances !{
if (instance(n)%num_regions .eq. 1) then !{
do k = 1, nk !{
do j = jsd, jed !{
do i = isd, ied !{
if (grid_tmask(i,j,k) .lt. 0.5) then !{
instance(n)%mask(i,j,k) = 0.0
instance(n)%index(i,j,k) = -1
elseif (instance(n)%region(1)%mask(i,j,k) .eq. 1.0) then !}{
instance(n)%mask(i,j,k) = instance(n)%region(1)%restore_region_value
instance(n)%index(i,j,k) = 1
else !}{
instance(n)%mask(i,j,k) = instance(n)%integrate_region_value
instance(n)%index(i,j,k) = 0
endif !}
enddo !} i
enddo !} j
enddo !} k
else !}{
do k = 1, nk !{
do j = jsd, jed !{
do i = isd, ied !{
if (grid_tmask(i,j,k) .gt. 0.5) then !{
!
! we're looking to set the restore_region_value for any point that has at least
! one region set
!
if (instance(n)%union) then !{
!
! find the first region with a value set
!
ind = 0
do nn = 1, instance(n)%num_regions !{
if (instance(n)%region(nn)%mask(i,j,k) .eq. 1.0) then !{
ind = nn
exit
endif !}
enddo !} nn
!
! if ind == 0, then this point is out of all regions
!
if (ind .ne. 0) then !{
!
! Quit with an error if a
! region with a different restore_region_value also has a value set
!
do nn = ind + 1, instance(n)%num_regions !{
if (instance(n)%region(nn)%mask(i,j,k) .eq. 1.0) then !{
if (instance(n)%region(nn)%restore_region_value .ne. &
instance(n)%region(ind)%restore_region_value) then !{
write (stdoutunit,*) trim(error_header), ' Grid point ', i, ', ', j, ', ',k, &
' in two regions for "', trim(instance(n)%name), '"'
call mpp_error(FATAL,trim(error_header) // ' Grid point in two regions for "' // &
trim(instance(n)%name) // '"')
elseif (instance(n)%region(nn)%restore .ne. instance(n)%region(ind)%restore) then !}{
write (stdoutunit,*) trim(error_header), ' Grid point ', i, ', ', j, ', ',k, &
' has different restore values for "', trim(instance(n)%name), '"'
call mpp_error(FATAL,trim(error_header) // ' Grid point has different restore values for "' // &
trim(instance(n)%name) // '"')
endif !}
endif !}
enddo !} nn
endif !}
else !}{ ! not a union, therefore an intersection
!
! Find the first region with a value set
!
ind = 0
do nn = 1, instance(n)%num_regions !{
if (instance(n)%region(nn)%mask(i,j,k) .eq. 1.0) then !{
ind = nn
exit
endif !}
enddo !} nn
!
! If ind == 0, then this point is out of all regions
!
if (ind .gt. 0) then !{
!
! Otherwise, check whether this value is set for all regions with
! the same restore_region_value. If so, use that value, otherwise
! use the integrate_region_value. Also, quit with an error if a
! region with a different restore_region_value also has a value set
!
good = .true.
do nn = ind + 1, instance(n)%num_regions !{
if (instance(n)%region(nn)%mask(i,j,k) .eq. 1.0) then !{
if (instance(n)%region(nn)%restore_region_value .ne. &
instance(n)%region(ind)%restore_region_value) then !{
write (stdoutunit,*) trim(error_header), ' Grid point ', i, ', ', j, ', ',k, &
' in two regions for "', trim(instance(n)%name), '"'
call mpp_error(FATAL,trim(error_header) // ' Grid point in two regions for "' // &
trim(instance(n)%name) // '"')
elseif (instance(n)%region(nn)%restore_region_value .ne. &
instance(n)%region(ind)%restore_region_value) then !}{
write (stdoutunit,*) trim(error_header), ' Grid point ', i, ', ', j, ', ',k, &
' has different restore values for "', trim(instance(n)%name), '"'
call mpp_error(FATAL,trim(error_header) // ' Grid point has different restore values for "' // &
trim(instance(n)%name) // '"')
else !}{
good = good .and. instance(n)%region(nn)%mask(i,j,k) .eq. 1.0
endif !}
endif !}
enddo !} nn
if (.not. good) then !{
ind = 0
endif !}
endif !}
endif !}
instance(n)%index(i,j,k) = ind
if (ind .gt. 0) then
instance(n)%mask(i,j,k) = instance(n)%region(ind)%restore_region_value
else
instance(n)%mask(i,j,k) = instance(n)%integrate_region_value
endif
else !}{
instance(n)%index(i,j,k) = -1
instance(n)%mask(i,j,k) = 0.0
endif !}
enddo !} i
enddo !} j
enddo !} k
endif !}
enddo !} n
!
! set the integrate_region value if something other than time is to be integrated
!
if (do_ocean_residency_integrand) then !{
call ocean_residency_integrand_source(isd, ied, jsd, jed, nk, Time%taum1, t_prog, t_diag, grid_tmask)
endif !}
do n = 1, num_instances !{
do k = 1, nk !{
do j = jsd, jed !{
do i = isd, ied !{
t_prog(instance(n)%tracer_index)%source(i,j,k) = t_prog(instance(n)%tracer_index)%source(i,j,k) + &
instance(n)%mask(i,j,k)
enddo !} i
enddo !} j
enddo !} k
!
! Save the diagnostic for the merged mask array
!
if (instance(n)%id_restore_region .gt. 0) then !{
used = send_data(instance(n)%id_restore_region, &
instance(n)%mask(:,:,:), &
Time%model_time, rmask = grid_tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif !}
enddo !} n
!
! Restore to the restore_region values
!
do n = 1, num_instances !{
if (instance(n)%some_restore) then !{
do k = 1, nk
do j = jsd, jed !{
do i = isd, ied !{
ind = instance(n)%index(i,j,k)
!if (instance(n)%region(ind)%mask(i,j,k) .eq. 1.0) then !{}
if (ind .gt. 0) then !{
if (instance(n)%region(ind)%restore .gt. 0.0) then !{
instance(n)%mask(i,j,k) = &
(instance(n)%region(ind)%restore_region_value - &
t_prog(instance(n)%tracer_index)%field(i,j,k,Time%tau)) / &
(instance(n)%region(ind)%restore * sec_per_day) * grid_tmask(i,j,k)
t_prog(instance(n)%tracer_index)%source(i,j,k) = &
t_prog(instance(n)%tracer_index)%source(i,j,k) + instance(n)%mask(i,j,k)
endif !}
endif !}
enddo !}i
enddo !} j
enddo !} k
endif !}
enddo !} n
return
end subroutine ocean_residency_source
! NAME="ocean_residency_source"
!#######################################################################
!
!
!
! Start the ocean residency package
!
! Residency surface area specification
!
! west_bnd : western longitude of residency region
! east_bnd : eastern longitude of residency region
! south_bnd : southern latitude of residency region
! north_bnd : northern latitude of residency region
! top_bnd : top depth of residency region
! bottom_bnd : bottom depth of residency region
!
! To set the volumes, a number of namelists are read,
! each containing the above values. You may specify up to
! num_geog_region rectangular cubes bounded by
! (west_bnd, east_bnd, north_bnd, south_bnd, top_bnd, bottom_bnd).
! Any grid box whose center is in one of these volumes will
! be considered to be part of the volume where the
! residency is reset to zero every time-step.
!
! north_bnd may not equal south_bnd, and west_bnd may not equal east_bnd
!
! top_depth may equal bottom_depth. In that case, then whatever vertical
! box contains that depth will define the vertical range for the box
!
! If south_bnd > north_bnd, then nothing will be done for that rectangle
!
! The initial surface area is empty, with the default rectangle
! setting the surface area to be empty
!
! More than num_geog_regions rectanglar volumes may be used to specify
! the volume by using more than one namelist
!
!
subroutine ocean_residency_start(isd, ied, jsd, jed, nk, model_time, grid_tracer_axes) !{
!
! modules
!
implicit none
!
!-----------------------------------------------------------------------
! Arguments
!-----------------------------------------------------------------------
!
integer, intent(in) :: isd
integer, intent(in) :: ied
integer, intent(in) :: jsd
integer, intent(in) :: jed
integer, intent(in) :: nk
type(time_type), intent(in) :: model_time
integer, dimension(:), intent(in) :: grid_tracer_axes
!
! local parameters
!
character(len=64), parameter :: sub_name = 'ocean_residency_start'
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
!
integer :: done
integer :: i
integer :: j
integer :: k
integer :: l
integer :: n
integer :: nn
character(len=256) :: caller_str
integer :: len_w
integer :: len_e
integer :: len_s
integer :: len_n
integer :: len_t
integer :: len_b
character(len=fm_field_name_len+3) :: long_suffix
character(len=fm_field_name_len+1) :: suffix
character(len=24) :: num_suffix
character(len=24) :: num_long_suffix
character(len=24) :: number_str
logical :: found_error
integer :: stdoutunit
stdoutunit=stdout()
!
!-----------------------------------------------------------------------
! save the instance namelist values
!-----------------------------------------------------------------------
!
caller_str = trim(mod_name) // '(' // trim(sub_name) // ')'
do n = 1, num_instances !{
call fm_util_start_namelist(package_name, instance(n)%name, caller = caller_str)
instance(n)%num_regions = fm_util_get_integer ('num_regions', scalar = .true.)
instance(n)%integrate_region_value = fm_util_get_real ('integrate_region_value', scalar = .true.)
instance(n)%union = fm_util_get_logical ('union', scalar = .true.)
instance(n)%int_module_name = fm_util_get_string ('int_module_name', scalar = .true.)
instance(n)%int_params => fm_util_get_real_array ('int_params')
instance(n)%int_flags => fm_util_get_logical_array('int_flags')
instance(n)%int_strings => fm_util_get_string_array ('int_strings')
!
! create namelists for each region
!
if (instance(n)%num_regions .le. 0) then
call mpp_error(FATAL,trim(error_header) // ' num_regions is non-positive for instance ' // &
trim(instance(n)%name) // ' in package ' // trim(package_name))
endif
!
! allocate storage for this instance and
! set all of the values to the default
!
allocate( instance(n)%region(instance(n)%num_regions) )
do nn = 1, instance(n)%num_regions !{
allocate( instance(n)%region(nn)%mask(isd:ied,jsd:jed,nk) )
do k = 1, nk
do j = jsd, jed
do i = isd, ied
instance(n)%region(nn)%mask(i,j,k) = 0.0
enddo
enddo
enddo
enddo !} nn
allocate( instance(n)%mask(isd:ied,jsd:jed,nk) )
do k = 1, nk
do j = jsd, jed
do i = isd, ied
instance(n)%mask(i,j,k) = 0.0
enddo
enddo
enddo
allocate( instance(n)%index(isd:ied,jsd:jed,nk) )
do k = 1, nk
do j = jsd, jed
do i = isd, ied
instance(n)%index(i,j,k) = 0
enddo
enddo
enddo
!instance(n)%mask => instance(n)%region(1)%mask
instance(n)%some_restore = .false.
instance(n)%some_fix = .false.
do nn = 1, instance(n)%num_regions !{
if (instance(n)%num_regions .eq. 1) then
suffix = ' '
else
write (long_suffix,*) nn
suffix = '_' // long_suffix(scan(long_suffix,'0123456789'):)
endif
instance(n)%region(nn)%restore = fm_util_get_real ('restore' // suffix, scalar = .true.)
instance(n)%region(nn)%restore_region_value = fm_util_get_real ('restore_region_value' // suffix, scalar = .true.)
instance(n)%region(nn)%swap = fm_util_get_logical ('swap' // suffix, scalar = .true.)
instance(n)%region(nn)%swap_module = fm_util_get_logical ('swap_module' // suffix, scalar = .true.)
instance(n)%region(nn)%west_bnd => fm_util_get_real_array ('west_bnd' // suffix)
instance(n)%region(nn)%east_bnd => fm_util_get_real_array ('east_bnd' // suffix)
instance(n)%region(nn)%south_bnd => fm_util_get_real_array ('south_bnd' // suffix)
instance(n)%region(nn)%north_bnd => fm_util_get_real_array ('north_bnd' // suffix)
instance(n)%region(nn)%top_bnd => fm_util_get_real_array ('top_bnd' // suffix)
instance(n)%region(nn)%bottom_bnd => fm_util_get_real_array ('bottom_bnd' // suffix)
instance(n)%region(nn)%module_name = fm_util_get_string ('module_name' // suffix, scalar = .true.)
instance(n)%region(nn)%params => fm_util_get_real_array ('params' // suffix)
instance(n)%region(nn)%flags => fm_util_get_logical_array('flags' // suffix)
instance(n)%region(nn)%strings => fm_util_get_string_array ('strings' // suffix)
instance(n)%some_restore = instance(n)%some_restore .or. instance(n)%region(nn)%restore .gt. 0.0
instance(n)%some_fix = instance(n)%some_fix .or. instance(n)%region(nn)%restore .eq. 0.0
enddo !} nn
call fm_util_end_namelist(package_name, instance(n)%name, caller = caller_str)
!
! Check some things
!
do nn = 1, instance(n)%num_regions !{
if (.not. associated(instance(n)%region(nn)%west_bnd) .and. &
.not. associated(instance(n)%region(nn)%east_bnd) .and. &
.not. associated(instance(n)%region(nn)%south_bnd) .and. &
.not. associated(instance(n)%region(nn)%north_bnd) .and. &
.not. associated(instance(n)%region(nn)%top_bnd) .and. &
.not. associated(instance(n)%region(nn)%bottom_bnd)) then !{
write (stdoutunit,*) trim(note_header), &
' No region specified, assuming it will be specified externally: ', trim(instance(n)%name), ' region: ', nn
instance(n)%region(nn)%num_geog_regions = 0
elseif (.not. associated(instance(n)%region(nn)%west_bnd) .or. &
.not. associated(instance(n)%region(nn)%east_bnd) .or. &
.not. associated(instance(n)%region(nn)%south_bnd) .or. &
.not. associated(instance(n)%region(nn)%north_bnd) .or. &
.not. associated(instance(n)%region(nn)%top_bnd) .or. &
.not. associated(instance(n)%region(nn)%bottom_bnd)) then !}{
call mpp_error(FATAL, trim(error_header) // ' Some regions not specified: ' // trim(instance(n)%name))
else !}{
len_w = size(instance(n)%region(nn)%west_bnd)
len_e = size(instance(n)%region(nn)%east_bnd)
len_s = size(instance(n)%region(nn)%south_bnd)
len_n = size(instance(n)%region(nn)%north_bnd)
len_t = size(instance(n)%region(nn)%top_bnd)
len_b = size(instance(n)%region(nn)%bottom_bnd)
if (len_e .ne. len_w .or. len_e .ne. len_s .or. len_e .ne. len_n .or. &
len_e .ne. len_t .or. len_e .ne. len_b) then !{
call mpp_error(FATAL, trim(error_header) // ' Region sizes are not equal: ' // trim(instance(n)%name))
endif !}
instance(n)%region(nn)%num_geog_regions = len_w
endif !}
enddo !} nn
!
! register the fields
!
if (instance(n)%name(1:1) .eq. '_') then !{
suffix = ' '
long_suffix = ' '
else !}{
suffix = '_' // instance(n)%name
long_suffix = ' (' // trim(instance(n)%name) // ')'
endif !}
instance(n)%id_change = register_diag_field(trim(diag_name), &
'residency_change' // trim(suffix), grid_tracer_axes(1:3), &
model_time, 'Residency change' // trim(long_suffix), ' ', &
missing_value = -1.0e+10)
instance(n)%id_restore_region = register_diag_field(trim(diag_name), &
'residency_in_merged_regions' // trim(suffix), grid_tracer_axes(1:3), &
model_time, 'Residency in merged regions' // trim(long_suffix), ' ', &
missing_value = -1.0e+10)
do nn = 1, instance(n)%num_regions !{
if (instance(n)%num_regions .eq. 1) then
num_suffix = ' '
else
write (number_str,*) nn
num_suffix = '_' // number_str(scan(number_str,'0123456789'):)
num_long_suffix = ' ' // number_str(scan(number_str,'0123456789'):)
endif
instance(n)%region(nn)%id_restore_region = register_diag_field(trim(diag_name), &
'residency_restore_region' // trim(num_suffix) // trim(suffix), grid_tracer_axes(1:3), &
model_time, 'Residency in region' // trim(num_long_suffix) // trim(long_suffix), ' ', &
missing_value = -1.0e+10)
enddo !} nn
enddo !} n
!
! call the start routines for instances controlled by external modules
!
call ocean_residency_ml_start
call ocean_residency_range_start
call ocean_residency_integrand_start
!
! check that the external tracers are all accounted for
!
found_error = .false.
do n = 1, num_instances !{
if (instance(n)%int_module_name .ne. ' ' .and. .not. instance(n)%int_found) then !}{
found_error = .true.
write (stdoutunit,*) trim(error_header), ' Instance "', trim(instance(n)%name), &
'" not found with integrand module name "', trim(instance(n)%int_module_name), '"'
endif !}
do nn = 1, instance(n)%num_regions !{
if (instance(n)%region(nn)%module_name .ne. ' ' .and. .not. instance(n)%region(nn)%found) then !{
found_error = .true.
write (stdoutunit,*) trim(error_header), ' Instance "', trim(instance(n)%name), &
'" not found with module name "', trim(instance(n)%region(nn)%module_name), '"'
endif !}
enddo !} nn
enddo !} n
if (found_error) then !{
call mpp_error(FATAL, trim(error_header) // ' Some external instances not found')
endif !}
return
end subroutine ocean_residency_start !}
! NAME="ocean_residency_start"
!#######################################################################
!
!
!
! Subroutine to do calculations needed every time-step after
! the continuity equation has been integrated
!
!
subroutine ocean_residency_tracer(isc, iec, jsc, jec, &
isd, ied, jsd, jed, nk, T_prog, grid_tmask, taup1, model_time, dtts) !{
implicit none
!
!-----------------------------------------------------------------------
! Arguments
!-----------------------------------------------------------------------
!
type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog
integer, intent(in) :: isc
integer, intent(in) :: iec
integer, intent(in) :: jsc
integer, intent(in) :: jec
integer, intent(in) :: isd
integer, intent(in) :: ied
integer, intent(in) :: jsd
integer, intent(in) :: jed
integer, intent(in) :: nk
integer, intent(in) :: taup1
type(time_type), intent(in) :: model_time
real, dimension(isd:,jsd:,:), intent(in) :: grid_tmask
real, intent(in) :: dtts
!
! local parameters
!
character(len=64), parameter :: sub_name = 'ocean_residency_tracer'
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
!
integer :: i
integer :: j
integer :: k
integer :: n
integer :: ind
logical :: used
!
! fix the restore_region values
!
do n = 1, num_instances !{
if (instance(n)%some_fix) then !{
do k = 1, nk
do j = jsd, jed !{
do i = isd, ied !{
ind = instance(n)%index(i,j,k)
!if (instance(n)%region(ind)%mask(i,j,k) .eq. 1.0) then !{}
if (ind .gt. 0) then !{
if (instance(n)%region(ind)%restore .eq. 0.0) then !{
instance(n)%mask(i,j,k) = &
(instance(n)%region(ind)%restore_region_value - &
t_prog(instance(n)%tracer_index)%field(i,j,k,taup1)) / &
dtts * grid_tmask(i,j,k)
t_prog(instance(n)%tracer_index)%field(i,j,k,taup1) = instance(n)%region(ind)%restore_region_value
endif !}
endif !}
enddo !}i
enddo !} j
enddo !} k
endif !}
if (instance(n)%id_change .gt. 0) then
used = send_data(instance(n)%id_change, &
instance(n)%mask(:,:,:), &
model_time, rmask = grid_tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
endif
enddo !} n
return
end subroutine ocean_residency_tracer !}
! NAME="ocean_residency_tracer"
end module ocean_residency_mod !}