Hello,
At Matt Harrison's suggestion, I am providing the sponge code for
anyone who may need it. This will replace the existing code in
mom4/ocean_param/sources/sponge/ocean_sponges.F90
of the mom4_beta2 release of Oct. 2002. The comments at the top of the
file will explain its use.
Matt will merge this code in the next release of mom4 but till then
you are welcome to use it. If you have any problems/questions please
let me know.
Swathi
module ocean_sponge_mod
!-----------------------------------------------------------------------
! GNU General Public License
! This file is a part of MOM.
!
! MOM 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; either version 2 of
! the License, or (at your option) any later version.
!
! MOM 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.
!
! For the full text of the GNU General Public License,
! write to: Free Software Foundation, Inc.,
! 675 Mass Ave, Cambridge, MA 02139, USA.
! or see: http://www.gnu.org/licenses/gpl.html
!-----------------------------------------------------------------------
!
!<CONTACT EMAIL="[EMAIL PROTECTED]"> P.S.Swathi
!</CONTACT>
!
!<CONTACT EMAIL="[EMAIL PROTECTED]"> Bonnie Samuels
!</CONTACT>
!
!<CONTACT EMAIL="[EMAIL PROTECTED]"> R.W. Hallberg
!</CONTACT>
!
!<OVERVIEW>
! Thickness weighted tracer tendency [tracer*meter/sec] from sponges.
!</OVERVIEW>
!
!<DESCRIPTION>
! Based on the code of bls and rwh. Sponge time rates are fixed in
! time but sponge values can vary in time. Blindly apply sponge.
! Some redundant calculations are send_data calls to avoid
! problems in diag_manager..Hard coded for n=1 temp and n=2 salt.
!
! This module applies sponges to tracers. The sponges
! can occur at any location and with any distribution in the domain, and
! with any time step and damping rate. Sponges occur where positive
! inverse restore times occur in the field passed to sponge_init. An
! array of tracer tendencies due to the sponges is augmented through a
! call to sponge_tracer_source. The array of tracer tendencies must be
! reset to zero between calls.
!
! Different damping rates (but constatnt in time) can be specified for each field.
! No sponges are applied to fields for
! which uniformly zero.
!
! The user is responsible for providing time-dependent data on
! the model grid of values towards which the tracers are being driven.
!</DESCRIPTION>
!
!<NAMELIST NAME="ocean_sponge_nml">
! <DATA NAME="restoretime_temp_path" TYPE="character*128">
! Path where the restore time field for temp is located
! </DATA>
! <DATA NAME="restoretime_salt_path" TYPE="character*128">
! Path where the restore time field for salt is located
! </DATA>
! <DATA NAME="restoretime_temp_name" TYPE="character*128">
! Name of the restore time field for temp in the netcdf file
! </DATA>
! <DATA NAME="restoretime_salt_name" TYPE="character*128">
! Name of the restore time field for salt in the netcdf file
! </DATA>
! <DATA NAME="sponge_temp_path" TYPE="character*128">
! Path where the sponge potential temp field is located
! </DATA>
! <DATA NAME="sponge_temp_name" TYPE="character*128">
! Name of the sponge potential temp field within the netcdf file
! </DATA>
! <DATA NAME="sponge_salt_path" TYPE="character*128">
! Path where the sponge salinity field is located
! </DATA>
! <DATA NAME="sponge_salt_name" TYPE="character*128">
! Name of the sponge salinity field within the netcdf file
! </DATA>
! <DATA NAME="num_sponges" TYPE="integer">
! Number of tracers that are damped to the sponge
! </DATA>
! <DATA NAME="target_date" TYPE="time_type">
! Time that values within sponge are damped towards
! </DATA>
!</NAMELIST>
use fms_mod, only: write_version_number, error_mesg, FATAL
use time_manager_mod, only: time_type, set_date
use time_manager_mod, only: operator( + ), operator(-), operator( // ), operator(>),
operator(==), operator(<=)
use time_interp_external_mod, only: time_interp_external, init_external_field
use diag_manager_mod, only : register_diag_field, send_data
use field_manager_mod, only : method_type, default_method
!
use mpp_io_mod, only: mpp_open, mpp_close, mpp_read
use mpp_io_mod, only: MPP_RDONLY, MPP_NETCDF, MPP_MULTI, MPP_ASCII, MPP_SINGLE
use mpp_io_mod, only: fieldtype, mpp_get_atts, mpp_get_info, mpp_get_fields
use mpp_mod, only: mpp_chksum, mpp_error, WARNING
use fms_mod, only : lowercase
use ocean_domains_mod, only: iscomp, iecomp, jscomp, jecomp, halo
use ocean_tmngr_mod, only : tau, taum1, model_time
implicit none
type, public :: ocean_sponge_type
character(len=32) :: name, units, tracer_name
character(len=128) :: longname
integer :: sponge_id ! index to flux_manager
integer :: tracer_id ! index to tracer_manager
integer :: num_inputs
real :: convert
character(len=32) :: input_name
character(len=128) :: input_file
integer :: input_id ! index to time_interp_external_mod
integer :: diag_id ! diagnostics manager ids
real :: range(2)
integer :: num_methods
type(method_type) :: methods
end type ocean_sponge_type
type(ocean_sponge_type), public, dimension(:), allocatable :: tracer_sponge(:)
type(ocean_sponge_type),public :: default_sponge
!
private
integer,allocatable :: num_col(:) ! The number of columns which are subject to
damping within this domain.
integer, allocatable :: isp(:,:), jsp(:,:) ! The indices of each column.
real, allocatable :: sponge_val(:,:) ! The values that are being damped towards, in
the units
! of the respective fields. The 2 indices
are (in order):
! the column number and the vertical levels
real, allocatable :: inv_damp(:,:,:), damp_coef(:,:,:) ! The inverse damping rates for
each field, and
! the related coefficint
{(1-exp(-dt/rate))/dt} based on
! the previous time step.
real,allocatable :: sponge_mask(:,:,:,:) ! sponge_mask.
logical, allocatable :: sponge_anywhere(:) ! Denotes whether each field has nonzero
inverse damping
! rates anywhere in this domain.
logical :: first_call = .true.
public ocean_sponge_init
public sponge_tracer_source
character(len=256) :: version = &
'=>Using: ocean_sponges.F90 ($Id: ocean_sponges.F90,v 1.1.1.1 2002/12/23 11:35:30
swathi Exp $)'
character (len=128) :: tagname = &
'$Name: $'
logical :: module_is_initialized = .FALSE.
contains
!#######################################################################
! <SUBROUTINE NAME="ocean_sponge_init">
!
! <DESCRIPTION>
! This subroutine is intended to be used to initialize the sponges.
! Everything in this subroutine is a user prototype, and should be replacable.
! </DESCRIPTION>
!
! <IN NAME="dt" UNITS="seconds" TYPE="real">
! The tracer timestep (or the time step for the sponges), in seconds
! </IN>
subroutine ocean_sponge_init(dt)
use mpp_mod, only: stdout, stdlog
use ocean_domains_mod, only: iscomp, iecomp, jscomp, jecomp, halo, local_domain
use ocean_grids_mod, only: nk, tmask,tracer_axes
use ocean_tmngr_mod
real,intent(in) :: dt
real,allocatable :: inv_resttime(:,:,:,:) ! The inverse damping rates in s-1.
logical :: error=.false.
integer :: i, j, k, nt, kk, k_found, col, nf, ll
character*128 :: restoretime_temp_path='',restoretime_temp_name='coeff'
character*128 :: restoretime_salt_path='',restoretime_salt_name='coeff'
character*128 :: sponge_temp_path='',sponge_temp_name='temp'
character*128 :: sponge_salt_path='',sponge_salt_name='salt'
!
character(len=32) :: name
integer :: unit_t, unit_s, ioun, io_status
integer :: ndim, nvar, natt, ntime, llen
integer :: num_sponges = 2
type(fieldtype), allocatable, dimension(:) :: fields
logical :: verbose_flag = .false.
real,parameter :: missing = -1.0e10
namelist /ocean_sponge_nml/ restoretime_temp_path,restoretime_temp_name,&
restoretime_salt_path,restoretime_salt_name, &
sponge_temp_path,sponge_temp_name,sponge_salt_path,sponge_salt_name,&
num_sponges
if ( module_is_initialized ) call &
error_mesg( 'ocean_sponge_init', &
'attempted reinitialization', &
FATAL )
module_is_initialized = .TRUE.
call write_version_number( version, tagname )
! provide for namelist over-ride of defaults
call mpp_open (ioun, 'input.nml', action=MPP_RDONLY, form=MPP_ASCII)
read (ioun, ocean_sponge_nml,iostat=io_status)
! write (stdlog(), ocean_sponge_nml)
write (stdout(), ocean_sponge_nml)
if (io_status .gt. 0) then
call mpp_error(FATAL, '=>Error reading input.nml at ocean_sponge_nml')
endif
call mpp_close (ioun)
!
allocate (inv_resttime(ij_bounds,nk,num_sponges))
allocate (sponge_mask(ij_bounds,nk,num_sponges))
allocate (num_col(num_sponges))
allocate (inv_damp((iecomp-iscomp+1)*(jecomp-jscomp+1),num_sponges,nk))
allocate (damp_coef((iecomp-iscomp+1)*(jecomp-jscomp+1),num_sponges,nk))
allocate (isp((iecomp-iscomp+1)*(jecomp-jscomp+1),num_sponges))
allocate (jsp((iecomp-iscomp+1)*(jecomp-jscomp+1),num_sponges))
allocate (sponge_anywhere(num_sponges))
!
inv_resttime(:,:,:,:) = 0.0
sponge_mask(:,:,:,:) = 0.0
sponge_anywhere(:) = .false.
num_col(:) = 0
isp(:,:) = 0
jsp(:,:) = 0
inv_damp(:,:,:) = 0.0
damp_coef(:,:,:) = 0.0
!
! read the 3D restore time data
! restore time units are sec-1
write (stdout(),*) ' Note:Reading Sponge restoring times from ',
trim(restoretime_temp_path)
call mpp_open (unit_t, trim(restoretime_temp_path), action=MPP_RDONLY,
form=MPP_NETCDF,threading=MPP_MULTI,fileset=MPP_SINGLE)
call mpp_get_info(unit_t, ndim, nvar, natt, ntime)
allocate(fields(nvar))
call mpp_get_fields(unit_t,fields)
io_status = 0
do i=1,nvar
call mpp_get_atts(fields(i),name=name)
if (lowercase(trim(name)) == lowercase(trim(restoretime_temp_name))) then
call mpp_read(unit_t, fields(i),local_domain,inv_resttime(:,:,:,1))
io_status=io_status+1
endif
enddo
if (io_status .eq. 0) then
call mpp_error(FATAL, '=>Error reading sponge restoretime_temp')
endif
call mpp_close(unit_t)
deallocate(fields)
!
write (stdout(),*) ' Note:Reading Sponge restoring times from ',
trim(restoretime_salt_path)
call mpp_open (unit_s, trim(restoretime_salt_path), action=MPP_RDONLY,
form=MPP_NETCDF,threading=MPP_MULTI,fileset=MPP_SINGLE)
call mpp_get_info(unit_s, ndim, nvar, natt, ntime)
allocate(fields(nvar))
call mpp_get_fields(unit_s,fields)
io_status = 0
do i=1,nvar
call mpp_get_atts(fields(i),name=name)
if (lowercase(trim(name)) == lowercase(trim(restoretime_salt_name))) then
call mpp_read(unit_s, fields(i),local_domain,inv_resttime(:,:,:,2))
io_status=io_status+1
endif
enddo
if (io_status .eq. 0) then
call mpp_error(FATAL, '=>Error reading sponge restoretime_salt')
endif
call mpp_close(unit_s)
deallocate(fields)
!
do ll = 1,num_sponges
! reset land points to zero
do k=1,nk
do j=jscomp,jecomp
do i=iscomp,iecomp
inv_resttime(i,j,k,ll)=inv_resttime(i,j,k,ll)*tmask(i,j,k)
if (inv_resttime(i,j,k,ll) > 0) sponge_mask(i,j,k,ll)= 1.0
enddo
enddo
enddo
write (stdout(),*) 'll, Checksum: Original inv_resttime =', ll,mpp_chksum(
inv_resttime(iscomp:iecomp,jscomp:jecomp,1:nk,ll))
! Count the columns that are in sponges
! May have to change this
num_col(ll) = 0
do j=jscomp,jecomp
do i=iscomp,iecomp
k_found = 0
do k=1,nk
if ((inv_resttime(i,j,k,ll) > 0.0) .and. (k_found == 0)) then
num_col(ll) = num_col(ll) + 1
k_found = 1
endif
enddo
enddo
enddo
!
if (num_col(ll) == 0) then
write(stdout(),*)'There are no sponges in this domain'
sponge_anywhere(ll) = .false.
else
sponge_anywhere(ll) = .true.
col = 1
do j=jscomp,jecomp
do i=iscomp,iecomp
k_found = 0
do kk=1,nk
if ((inv_resttime(i,j,kk,ll) > 0.0) .and. (k_found == 0)) then
k_found = 1
isp(col,ll) = i
jsp(col,ll) = j
do k=1,nk
inv_damp(col,ll,k) = inv_resttime(i,j,k,ll)
enddo
col = col + 1
endif
enddo
enddo
enddo
endif
call calculate_damp_coef(dt, ll)
enddo
!
deallocate (inv_resttime)
! Get the sponge data fields
allocate(tracer_sponge(num_sponges))
!
! Some default values
!
default_sponge%name = 'none'
default_sponge%longname = 'none'
default_sponge%units = 'none'
default_sponge%convert = 1.0
default_sponge%input_name = 'none'
default_sponge%input_file = 'none'
default_sponge%range = (/-1.0e10,1.0e10/)
default_sponge%input_id = -1
default_sponge%num_methods = 0
default_sponge%methods%method_type = 'none'
default_sponge%diag_id = -1
!
tracer_sponge = default_sponge
!
tracer_sponge(1)%name = 'temp'
tracer_sponge(1)%longname = 'Sponge target temperature'
tracer_sponge(1)%units = 'deg C'
tracer_sponge(1)%convert = 1.0
tracer_sponge(1)%input_name = 'temp'
tracer_sponge(1)%input_file = trim(sponge_temp_path)
tracer_sponge(1)%range = (/-10,45/)
tracer_sponge(1)%input_id = init_external_field(tracer_sponge(1)%input_file,&
tracer_sponge(1)%input_name,&
MPP_NETCDF,MPP_MULTI,local_domain,'deg_C',verbose_flag)
if(tracer_sponge(1)%input_id == -1) then
call mpp_error(FATAL,'=>Error reading temp sponge field inside
ocean_sponges_init')
endif
tracer_sponge(1)%num_methods = 1
tracer_sponge(1)%methods%method_type = 'sponge_temp'
tracer_sponge(1)%diag_id =
register_diag_field('ocean_sponges',trim(tracer_sponge(1)%name)//'_sponge_data',&
tracer_axes(1:3),model_time,tracer_sponge(1)%longname,tracer_sponge(1)%units,&
missing_value=missing,range=tracer_sponge(1)%range)
!
tracer_sponge(2)%name = 'salt'
tracer_sponge(2)%longname = 'Sponge target salinity'
tracer_sponge(2)%units = 'g/kg'
tracer_sponge(2)%convert = 1.0
tracer_sponge(2)%input_name = 'salt'
tracer_sponge(2)%input_file = trim(sponge_salt_path)
tracer_sponge(2)%range = (/0,45/)
tracer_sponge(2)%input_id = init_external_field(tracer_sponge(2)%input_file,&
tracer_sponge(2)%input_name,&
MPP_NETCDF,MPP_MULTI,local_domain,'g/kg',verbose_flag)
if(tracer_sponge(2)%input_id == -1) then
call mpp_error(FATAL,'=>Error reading salt sponge field inside
ocean_sponges_init')
endif
tracer_sponge(2)%num_methods = 1
tracer_sponge(2)%methods%method_type = 'sponge_salt'
tracer_sponge(2)%diag_id =
register_diag_field('ocean_sponges',trim(tracer_sponge(2)%name)//'_sponge_data',&
tracer_axes(1:3),model_time,tracer_sponge(2)%longname,tracer_sponge(2)%units,&
missing_value=missing,range=tracer_sponge(2)%range)
end subroutine ocean_sponge_init
!#######################################################################
! <SUBROUTINE NAME="sponge_tracer_source">
!
! <DESCRIPTION>
! This subroutine calculates thickness weighted time tendencies of the
! tracers due to the sponges.
! </DESCRIPTION>
!
! <IN NAME="t" UNITS="tracer" TYPE="real" DIMS="(ij_bounds,nk)">
! The tracer array at the appropriate damping time level
! </IN>
! <INOUT NAME="source" TYPE="real" DIMS="(ij_bounds,nk)">
! The tracer source from sponges in [tracer * meter/sec].
! This array is assumed to be previously initialized to 0, and may be
! accumulated through several calls to different subroutines.
! </INOUT>
! <IN NAME="nf" TYPE="integer">
! The tracer number upon which to work
! </IN>
! <IN NAME="time" TYPE="(time_type)">
! The time at which the sponges are being applied
! </IN>
! <IN NAME="dt_in" TYPE="real">
! The time step in seconds
! </IN>
subroutine sponge_tracer_source(t, source, nf, time, dt_in)
use mpp_mod, only: stdout, stdlog
use ocean_grids_mod, only: nk, tmask,dht
real,intent(in), dimension(ij_bounds,nk) :: t
real,intent(inout), dimension(ij_bounds,nk) :: source
integer, intent(in) :: nf
type(time_type),intent(in) :: time
real,intent(in) :: dt_in
real, dimension(ij_bounds,nk) :: data
logical :: verbose_flag = .false.
real,parameter :: missing = -1.0e10
logical :: used
!
integer :: k, col
! Sponges currently allowed for temp and salt
if (nf > 2) then
write(stdout(),*)'Sponges only for temp and salt, nf = ',nf
return
endif
data = 0.0
!
! We allow writing even when there are no sponges because of
! problems with the diag_manager
if (sponge_anywhere(nf)) then
allocate(sponge_val(num_col(nf),nk))
call time_interp_external(tracer_sponge(nf)%input_id, model_time, data, verbose
= verbose_flag)
if (tracer_sponge(nf)%diag_id > 0) then
used = send_data(tracer_sponge(nf)%diag_id, &
data(___,1:nk)*sponge_mask(___,1:nk,nf), &
model_time,rmask=tmask(___,1:nk))
endif
do k=1,nk
do col=1,num_col(nf)
sponge_val(col,k) = data(isp(col,nf),jsp(col,nf),k)
enddo
enddo
do col=1,num_col(nf)
do k=1,nk
source(isp(col,nf),jsp(col,nf),k) = source(isp(col,nf),jsp(col,nf),k) &
+ dht(isp(col,nf),jsp(col,nf),k)*damp_coef(col,nf,k)*(sponge_val(col,k) -
t(isp(col,nf),jsp(col,nf),k))
enddo
enddo
deallocate(sponge_val)
else
if (tracer_sponge(nf)%diag_id > 0) then
used = send_data(tracer_sponge(nf)%diag_id, &
data(___,1:nk), &
model_time,rmask=tmask(___,1:nk))
endif
endif
end subroutine sponge_tracer_source
! </SUBROUTINE> NAME="sponge_tracer_source"
!#######################################################################
! <SUBROUTINE NAME="calculate_damp_coef">
!
! <DESCRIPTION>
!
! This subroutine calculates the damping coefficent {damp_coef =
! (1-exp(-dt/rate))/dt} with appropriate checks to handle the limits of
! huge or small damping rates gracefully.
! </DESCRIPTION>
!
! <IN NAME="dt" UNITS="seconds" TYPE="real">
! The time step in seconds
! </IN>
! <IN NAME="nf" TYPE="integer">
! The number of the field whose damping coefficient is being set
! </IN>
subroutine calculate_damp_coef(dt, nf)
use ocean_grids_mod, only: nk
real,intent(in) :: dt
integer, intent(in) :: nf
real :: inv_dt ! The inverse of the time step, in s-1.
integer :: col, kk
inv_dt = 1.0 / dt
do col=1,num_col(nf)
do kk=1,nk
if (dt*inv_damp(col,nf,kk) > 4.0e-8) then
if (dt*inv_damp(col,nf,kk) > 37.0) then
damp_coef(col,nf,kk) = inv_dt
else
damp_coef(col,nf,kk) = (1.0 - exp(-dt*inv_damp(col,nf,kk))) * inv_dt
endif
else if (dt*inv_damp(col,nf,kk) > 0.0) then
damp_coef(col,nf,kk) = inv_damp(col,nf,kk)
else
damp_coef(col,nf,kk) = 0.0
endif
enddo
enddo
end subroutine calculate_damp_coef
! </SUBROUTINE> NAME="calculate_damp_coef"
end module ocean_sponge_mod
