#include <misc.h>
#include <preproc.h>
module accumulMod 5,7
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: accumulMod
!
! !DESCRIPTION:
! This module contains generic subroutines that can be used to
! define, accumulate and extract user-specified fields over
! user-defined intervals. Each interval and accumulation type is
! unique to each field processed.
! Subroutine [init_accumulator] defines the values of the accumulated
! field data structure. Subroutine [update_accum_field] does
! the actual accumulation for a given field.
! Four types of accumulations are possible:
! - Average over time interval. Time average fields are only
! valid at the end of the averaging interval.
! - Running mean over time interval. Running means are valid once the
! length of the simulation exceeds the
! - Running accumulation over time interval. Accumulated fields are
! continuously accumulated. The trigger value "-99999." resets
! the accumulation to zero.
!
! !USES:
use shr_kind_mod, only: r8 => shr_kind_r8
use abortutils, only: endrun
use clm_varctl, only: iulog
!
! !PUBLIC TYPES:
implicit none
save
!
! !PUBLIC MEMBER FUNCTIONS:
public :: accumulRest ! Write/read restart of accumulation fields
public :: init_accum_field ! Initialize an accumulator field
public :: print_accum_fields ! Print info about accumulator fields
public :: extract_accum_field ! Extracts the current value of an accumulator field
public :: update_accum_field ! Update the current value of an accumulator field
interface extract_accum_field 31
module procedure extract_accum_field_sl! Extract current val of single-level accumulator field
module procedure extract_accum_field_ml! Extract current val of multi-level accumulator field
end interface
interface update_accum_field ! Updates the current value of an accumulator field 17
module procedure update_accum_field_sl! Update single-level accumulator field
module procedure update_accum_field_ml! Update multi-level accumulator field
end interface
!
! !REVISION HISTORY:
! Created by Sam Levis
! Updated to clm2.1 data structures by Mariana Vertenstein
! Updated to include all subgrid type and multilevel fields, M. Vertenstein 03/2003
!
!
private
!
! !PRIVATE TYPES:
!EOP
!
type accum_field
character(len= 8) :: name !field name
character(len=128) :: desc !field description
character(len= 8) :: units !field units
character(len= 8) :: acctype !accumulation type: ["timeavg","runmean","runaccum"]
character(len= 8) :: type1d !subgrid type: ["gridcell","landunit","column" or "pft"]
character(len= 8) :: type2d !type2d ('','levsoi','numrad',..etc. )
integer :: beg1d !subgrid type beginning index
integer :: end1d !subgrid type ending index
integer :: num1d !total subgrid points
integer :: numlev !number of vertical levels in field
real(r8):: initval !initial value of accumulated field
real(r8), pointer :: val(:,:) !accumulated field
integer :: period !field accumulation period (in model time steps)
end type accum_field
integer, parameter :: max_accum = 100 !maximum number of accumulated fields
type (accum_field) :: accum(max_accum) !array accumulated fields
integer :: naccflds = 0 !accumulator field counter
!------------------------------------------------------------------------
contains
!------------------------------------------------------------------------
!BOP
!
! !IROUTINE: init_accum_field
!
! !INTERFACE:
subroutine init_accum_field (name, units, desc, & 17,8
accum_type, accum_period, numlev, subgrid_type, init_value, type2d)
!
! !DESCRIPTION:
! Initialize accumulation fields. This subroutine sets:
! o name of accumulated field
! o units of accumulated field
! o accumulation type of accumulated field
! o description of accumulated fields: accdes
! o accumulation period for accumulated field (in iterations)
! o initial value of accumulated field
!
! !USES:
use shr_const_mod, only: SHR_CONST_CDAY
use clm_time_manager, only : get_step_size
use decompMod, only : get_proc_bounds, get_proc_global
!
! !ARGUMENTS:
implicit none
character(len=*), intent(in) :: name !field name
character(len=*), intent(in) :: units !field units
character(len=*), intent(in) :: desc !field description
character(len=*), intent(in) :: accum_type !field type: tavg, runm, runa, ins
integer , intent(in) :: accum_period !field accumulation period
character(len=*), intent(in) :: subgrid_type !["gridcell","landunit","column" or "pft"]
integer , intent(in) :: numlev !number of vertical levels
real(r8), intent(in) :: init_value !field initial or reset value
character(len=*), intent(in), optional :: type2d !level type (optional) - needed if numlev > 1
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 03/2003
!
!
! !LOCAL VARIABLES:
!EOP
integer :: nf ! field index
integer :: beg1d,end1d ! beggining and end subgrid indices
integer :: num1d ! total number subgrid indices
integer :: begp, endp ! per-proc beginning and ending pft indices
integer :: begc, endc ! per-proc beginning and ending column indices
integer :: begl, endl ! per-proc beginning and ending landunit indices
integer :: begg, endg ! per-proc gridcell ending gridcell indices
integer :: numg ! total number of gridcells across all processors
integer :: numl ! total number of landunits across all processors
integer :: numc ! total number of columns across all processors
integer :: nump ! total number of pfts across all processors
!------------------------------------------------------------------------
! Determine necessary indices
call get_proc_bounds(begg, endg, begl, endl, begc, endc, begp, endp)
call get_proc_global(numg, numl, numc, nump)
! update field index
! Consistency check that number of accumulated does not exceed maximum.
naccflds = naccflds + 1
if (naccflds > max_accum) then
write(iulog,*) 'ACCUMULINIT error: user-defined accumulation fields ', &
'equal to ',naccflds,' exceeds max_accum'
call endrun
end if
nf = naccflds
! Note accumulation period must be converted from days
! to number of iterations
accum(nf)%name = trim(name)
accum(nf)%units = trim(units)
accum(nf)%desc = trim(desc)
accum(nf)%acctype = trim(accum_type)
accum(nf)%initval = init_value
accum(nf)%period = accum_period
if (accum(nf)%period < 0) then
accum(nf)%period = -accum(nf)%period * nint(SHR_CONST_CDAY) / get_step_size()
end if
select case (trim(subgrid_type))
case ('gridcell')
beg1d = begg
end1d = endg
num1d = numg
case ('landunit')
beg1d = begl
end1d = endl
num1d = numl
case ('column')
beg1d = begc
end1d = endc
num1d = numc
case ('pft')
beg1d = begp
end1d = endp
num1d = nump
case default
write(iulog,*)'ACCUMULINIT: unknown subgrid type ',subgrid_type
call endrun ()
end select
accum(nf)%type1d = trim(subgrid_type)
accum(nf)%beg1d = beg1d
accum(nf)%end1d = end1d
accum(nf)%num1d = num1d
accum(nf)%numlev = numlev
if (present(type2d)) then
accum(nf)%type2d = type2d
else
accum(nf)%type2d = ' '
end if
! Allocate and initialize accumulation field
allocate(accum(nf)%val(beg1d:end1d,numlev))
accum(nf)%val(beg1d:end1d,numlev) = init_value
end subroutine init_accum_field
!------------------------------------------------------------------------
!BOP
!
! !IROUTINE: print_accum_fields
!
! !INTERFACE:
subroutine print_accum_fields() 1,2
!
! !DESCRIPTION:
! Diagnostic printout of accumulated fields
!
! !USES:
use spmdMod, only : masterproc
use nanMod, only : bigint
!
! !ARGUMENTS:
implicit none
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein 03/2003
!
!EOP
!
integer :: i,nf !indices
!------------------------------------------------------------------------
if (masterproc) then
write(iulog,*)
write(iulog,*) 'Initializing variables for time accumulation .....'
write(iulog,'(72a1)') ("-",i=1,60)
write(iulog,*) 'Accumulated fields'
write(iulog,1002)
write(iulog,'(72a1)') ("_",i=1,71)
do nf = 1, naccflds
if (accum(nf)%period /= bigint) then
write(iulog,1003) nf,accum(nf)%name,accum(nf)%units,&
accum(nf)%acctype, accum(nf)%period, accum(nf)%initval, &
accum(nf)%desc
else
write(iulog,1004) nf,accum(nf)%name,accum(nf)%units,&
accum(nf)%acctype, accum(nf)%initval, accum(nf)%desc
endif
end do
write(iulog,'(72a1)') ("_",i=1,71)
write(iulog,*)
write(iulog,'(72a1)') ("-",i=1,60)
write(iulog,*) 'Successfully initialized variables for accumulation'
write(iulog,*)
endif
1002 format(' No',' Name ',' Units ',' Type ','Period',' Inival',' Description')
1003 format((1x,i2),(1x,a8),(1x,a8),(1x,a8), (1x,i5),(1x,f4.0),(1x,a40))
1004 format((1x,i2),(1x,a8),(1x,a8),(1x,a8),' N.A.',(1x,f4.0),(1x,a40))
end subroutine print_accum_fields
!------------------------------------------------------------------------
!BOP
!
! !IROUTINE: extract_accum_field_sl
!
! !INTERFACE:
subroutine extract_accum_field_sl (name, field, nstep) 1,3
!
! !DESCRIPTION:
! Extract single-level accumulated field.
! This routine extracts the field values from the multi-level
! accumulation field. It extracts the current value except if
! the field type is a time average. In this case, an absurd value
! is assigned to indicate the time average is not yet valid.
!
! !USES:
use clm_varcon, only : spval
!
! !ARGUMENTS:
implicit none
character(len=*), intent(in) :: name !field name
real(r8), pointer, dimension(:) :: field !field values for current time step
integer , intent(in) :: nstep !timestep index
!
! !REVISION HISTORY:
! Created by Sam Levis
! Updated to clm2.1 data structures by Mariana Vertenstein
! Updated to include all subgrid type and multilevel fields, Mariana Vertenstein 03-2003
!
!
! !LOCAL VARIABLES:
!EOP
integer :: i,k,nf !indices
integer :: beg,end !subgrid beginning,ending indices
!------------------------------------------------------------------------
! find field index. return if "name" is not on list
nf = 0
do i = 1, naccflds
if (name == accum(i)%name) nf = i
end do
if (nf == 0) then
write(iulog,*) 'EXTRACT_ACCUM_FIELD_SL error: field name ',name,' not found'
call endrun
endif
! error check
beg = accum(nf)%beg1d
end = accum(nf)%end1d
if (size(field,dim=1) /= end-beg+1) then
write(iulog,*)'ERROR in extract_accum_field for field ',accum(nf)%name
write(iulog,*)'size of first dimension of field is ',&
size(field,dim=1),' and should be ',end-beg+1
call endrun
endif
! extract field
if (accum(nf)%acctype == 'timeavg' .and. &
mod(nstep,accum(nf)%period) /= 0) then
!dir$ concurrent
!cdir nodep
do k = beg,end
field(k) = spval !assign absurd value when avg not ready
end do
else
!dir$ concurrent
!cdir nodep
do k = beg,end
field(k) = accum(nf)%val(k,1)
end do
end if
end subroutine extract_accum_field_sl
!------------------------------------------------------------------------
!BOP
!
! !IROUTINE: extract_accum_field_ml
!
! !INTERFACE:
subroutine extract_accum_field_ml (name, field, nstep) 1,4
!
! !DESCRIPTION:
! Extract mutli-level accumulated field.
! This routine extracts the field values from the multi-level
! accumulation field. It extracts the current value except if
! the field type is a time average. In this case, an absurd value
! is assigned to indicate the time average is not yet valid.
!
! !USES:
use clm_varcon, only : spval
!
! !ARGUMENTS:
implicit none
character(len=*), intent(in) :: name !field name
real(r8), pointer, dimension(:,:) :: field !field values for current time step
integer, intent(in) :: nstep !timestep index
!
! !REVISION HISTORY:
! Created by Sam Levis
! Updated to clm2.1 data structures by Mariana Vertenstein
! Updated to include all subgrid type and multilevel fields, M. Vertenstein 03/2003
!
!
! !LOCAL VARIABLES:
!EOP
integer :: i,j,k,nf !indices
integer :: beg,end !subgrid beginning,ending indices
integer :: numlev !number of vertical levels
!------------------------------------------------------------------------
! find field index. return if "name" is not on list
nf = 0
do i = 1, naccflds
if (name == accum(i)%name) nf = i
end do
if (nf == 0) then
write(iulog,*) 'EXTRACT_ACCUM_FIELD_ML error: field name ',name,' not found'
call endrun
endif
! error check
numlev = accum(nf)%numlev
beg = accum(nf)%beg1d
end = accum(nf)%end1d
if (size(field,dim=1) /= end-beg+1) then
write(iulog,*)'ERROR in extract_accum_field for field ',accum(nf)%name
write(iulog,*)'size of first dimension of field is ',&
size(field,dim=1),' and should be ',end-beg+1
call endrun
else if (size(field,dim=2) /= numlev) then
write(iulog,*)'ERROR in extract_accum_field for field ',accum(nf)%name
write(iulog,*)'size of second dimension of field iis ',&
size(field,dim=2),' and should be ',numlev
call endrun
endif
!extract field
if (accum(nf)%acctype == 'timeavg' .and. &
mod(nstep,accum(nf)%period) /= 0) then
do j = 1,numlev
!dir$ concurrent
!cdir nodep
do k = beg,end
field(k,j) = spval !assign absurd value when avg not ready
end do
end do
else
do j = 1,numlev
!dir$ concurrent
!cdir nodep
do k = beg,end
field(k,j) = accum(nf)%val(k,j)
end do
end do
end if
end subroutine extract_accum_field_ml
!------------------------------------------------------------------------
!BOP
!
! !IROUTINE: update_accum_field_sl
!
! !INTERFACE:
subroutine update_accum_field_sl (name, field, nstep) 1,3
!
! !DESCRIPTION:
! Accumulate single level field over specified time interval.
! The appropriate field is accumulated in the array [accval].
!
! !ARGUMENTS:
implicit none
character(len=*), intent(in) :: name !field name
real(r8), pointer, dimension(:) :: field !field values for current time step
integer , intent(in) :: nstep !time step index
!
! !REVISION HISTORY:
! Created by Sam Levis
! Updated to clm2.1 data structures by Mariana Vertenstein
! Updated to include all subgrid type and multilevel fields by M. Vertenstein 03/2003
!
!
! !LOCAL VARIABLES:
!EOP
integer :: i,k,nf !indices
integer :: accper !temporary accumulation period
integer :: beg,end !subgrid beginning,ending indices
!------------------------------------------------------------------------
! find field index. return if "name" is not on list
nf = 0
do i = 1, naccflds
if (name == accum(i)%name) nf = i
end do
if (nf == 0) then
write(iulog,*) 'UPDATE_ACCUM_FIELD_SL error: field name ',name,' not found'
call endrun
endif
! error check
beg = accum(nf)%beg1d
end = accum(nf)%end1d
if (size(field,dim=1) /= end-beg+1) then
write(iulog,*)'ERROR in UPDATE_ACCUM_FIELD_SL for field ',accum(nf)%name
write(iulog,*)'size of first dimension of field is ',size(field,dim=1),&
' and should be ',end-beg+1
call endrun
endif
! accumulate field
if (accum(nf)%acctype /= 'timeavg' .AND. &
accum(nf)%acctype /= 'runmean' .AND. &
accum(nf)%acctype /= 'runaccum') then
write(iulog,*) 'UPDATE_ACCUM_FIELD_SL error: incorrect accumulation type'
write(iulog,*) ' was specified for field ',name
write(iulog,*)' accumulation type specified is ',accum(nf)%acctype
write(iulog,*)' only [timeavg, runmean, runaccum] are currently acceptable'
call endrun()
end if
! reset accumulated field value if necessary and update
! accumulation field
! running mean never reset
if (accum(nf)%acctype == 'timeavg') then
!time average field reset every accumulation period
!normalize at end of accumulation period
if ((mod(nstep,accum(nf)%period) == 1) .and. (nstep /= 0)) then
accum(nf)%val(beg:end,1) = 0._r8
end if
accum(nf)%val(beg:end,1) = accum(nf)%val(beg:end,1) + field(beg:end)
if (mod(nstep,accum(nf)%period) == 0) then
accum(nf)%val(beg:end,1) = accum(nf)%val(beg:end,1) / accum(nf)%period
endif
else if (accum(nf)%acctype == 'runmean') then
!running mean - reset accumulation period until greater than nstep
accper = min (nstep,accum(nf)%period)
accum(nf)%val(beg:end,1) = ((accper-1)*accum(nf)%val(beg:end,1) + field(beg:end)) / accper
else if (accum(nf)%acctype == 'runaccum') then
!running accumulation field reset at trigger -99999
!dir$ concurrent
!cdir nodep
do k = beg,end
if (nint(field(k)) == -99999) then
accum(nf)%val(k,1) = 0._r8
end if
end do
accum(nf)%val(beg:end,1) = min(max(accum(nf)%val(beg:end,1) + field(beg:end), 0._r8), 99999._r8)
end if
end subroutine update_accum_field_sl
!------------------------------------------------------------------------
!BOP
!
! !IROUTINE: update_accum_field_ml
!
! !INTERFACE:
subroutine update_accum_field_ml (name, field, nstep) 1,4
!
! !DESCRIPTION:
! Accumulate multi level field over specified time interval.
!
! !ARGUMENTS:
implicit none
character(len=*), intent(in) :: name !field name
real(r8), pointer, dimension(:,:) :: field !field values for current time step
integer , intent(in) :: nstep !time step index
!
! !REVISION HISTORY:
! Created by Sam Levis
! Updated to clm2.1 data structures by Mariana Vertenstein
! Updated to include all subgrid type and multilevel fields by M. Vertenstein 03/2003
!
!
! !LOCAL VARIABLES:
!EOP
integer :: i,j,k,nf !indices
integer :: accper !temporary accumulation period
integer :: beg,end !subgrid beginning,ending indices
integer :: numlev !number of vertical levels
!------------------------------------------------------------------------
! find field index. return if "name" is not on list
nf = 0
do i = 1, naccflds
if (name == accum(i)%name) nf = i
end do
if (nf == 0) then
write(iulog,*) 'UPDATE_ACCUM_FIELD_ML error: field name ',name,' not found'
call endrun
endif
! error check
numlev = accum(nf)%numlev
beg = accum(nf)%beg1d
end = accum(nf)%end1d
if (size(field,dim=1) /= end-beg+1) then
write(iulog,*)'ERROR in UPDATE_ACCUM_FIELD_ML for field ',accum(nf)%name
write(iulog,*)'size of first dimension of field is ',size(field,dim=1),&
' and should be ',end-beg+1
call endrun
else if (size(field,dim=2) /= numlev) then
write(iulog,*)'ERROR in UPDATE_ACCUM_FIELD_ML for field ',accum(nf)%name
write(iulog,*)'size of second dimension of field is ',size(field,dim=2),&
' and should be ',numlev
call endrun
endif
! accumulate field
if (accum(nf)%acctype /= 'timeavg' .AND. &
accum(nf)%acctype /= 'runmean' .AND. &
accum(nf)%acctype /= 'runaccum') then
write(iulog,*) 'UPDATE_ACCUM_FIELD_ML error: incorrect accumulation type'
write(iulog,*) ' was specified for field ',name
write(iulog,*)' accumulation type specified is ',accum(nf)%acctype
write(iulog,*)' only [timeavg, runmean, runaccum] are currently acceptable'
call endrun()
end if
! accumulate field
! reset accumulated field value if necessary and update
! accumulation field
! running mean never reset
if (accum(nf)%acctype == 'timeavg') then
!time average field reset every accumulation period
!normalize at end of accumulation period
if ((mod(nstep,accum(nf)%period) == 1) .and. (nstep /= 0)) then
accum(nf)%val(beg:end,1:numlev) = 0._r8
endif
accum(nf)%val(beg:end,1:numlev) = accum(nf)%val(beg:end,1:numlev) + field(beg:end,1:numlev)
if (mod(nstep,accum(nf)%period) == 0) then
accum(nf)%val(beg:end,1:numlev) = accum(nf)%val(beg:end,1:numlev) / accum(nf)%period
endif
else if (accum(nf)%acctype == 'runmean') then
!running mean - reset accumulation period until greater than nstep
accper = min (nstep,accum(nf)%period)
accum(nf)%val(beg:end,1:numlev) = &
((accper-1)*accum(nf)%val(beg:end,1:numlev) + field(beg:end,1:numlev)) / accper
else if (accum(nf)%acctype == 'runaccum') then
!running accumulation field reset at trigger -99999
do j = 1,numlev
!dir$ concurrent
!cdir nodep
do k = beg,end
if (nint(field(k,j)) == -99999) then
accum(nf)%val(k,j) = 0._r8
end if
end do
end do
accum(nf)%val(beg:end,1:numlev) = &
min(max(accum(nf)%val(beg:end,1:numlev) + field(beg:end,1:numlev), 0._r8), 99999._r8)
end if
end subroutine update_accum_field_ml
!------------------------------------------------------------------------
!BOP
!
! !IROUTINE: accumulRest
!
! !INTERFACE:
subroutine accumulRest( ncid, flag ) 3,11
!
! !DESCRIPTION:
! Read/write accumulation restart data
!
! !USES:
use ncdio
!
! !ARGUMENTS:
implicit none
integer, intent(in) :: ncid !netcdf unit
character(len=*), intent(in) :: flag !'define','read', or 'write'
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
!
!
! !LOCAL VARIABLES:
!EOP
integer :: nf,k,j ! indices
integer :: beg1d, end1d ! buffer bounds
integer :: ier ! error status
real(r8), pointer :: rbuf1d(:) ! temporary 1d buffer
logical :: readvar ! determine if variable is on initial file
character(len=128) :: varname ! temporary
character(len= 32) :: subname='AccumRest' ! subroutine name
!------------------------------------------------------------------------
do nf = 1,naccflds
! Note = below need to allocate rbuf for single level variables, since
! accum(nf)%val is always 2d
varname = trim(accum(nf)%name) // '_VALUE'
if (flag == 'define') then
if (accum(nf)%numlev == 1) then
call ncd_defvar(ncid=ncid, varname=varname, xtype=nf_double, &
dim1name=accum(nf)%type1d, &
long_name=accum(nf)%desc, units=accum(nf)%units)
else
call ncd_defvar(ncid=ncid, varname=varname, xtype=nf_double, &
dim1name=accum(nf)%type1d, dim2name=accum(nf)%type2d, &
long_name=accum(nf)%desc, units=accum(nf)%units)
end if
else if (flag == 'read' .or. flag == 'write') then
if (accum(nf)%numlev == 1) then
beg1d = accum(nf)%beg1d
end1d = accum(nf)%end1d
allocate(rbuf1d(beg1d:end1d))
if (flag == 'write') then
rbuf1d(beg1d:end1d) = accum(nf)%val(beg1d:end1d,1)
end if
call ncd_iolocal(varname=varname, data=rbuf1d, &
dim1name=accum(nf)%type1d, &
ncid=ncid, flag=flag, readvar=readvar)
if (flag == 'read' .and. readvar) then
accum(nf)%val(beg1d:end1d,1) = rbuf1d(beg1d:end1d)
end if
deallocate(rbuf1d)
else
call ncd_iolocal(varname=varname, data=accum(nf)%val, &
dim1name=accum(nf)%type1d, dim2name=accum(nf)%type2d, &
ncid=ncid, flag=flag, readvar=readvar)
end if
if (flag=='read' .and. .not. readvar) then
if (is_restart()) call endrun()
end if
end if
varname = trim(accum(nf)%name) // '_PERIOD'
if (flag == 'define') then
call ncd_defvar(ncid=ncid, varname=varname, xtype=nf_int, &
long_name='', units='time steps')
else if (flag == 'read' .or. flag == 'write') then
call ncd_ioglobal(varname=varname, data=accum(nf)%period, &
ncid=ncid, flag=flag, readvar=readvar, bcast=.true.)
if (flag=='read' .and. .not. readvar) then
if (is_restart()) call endrun()
end if
end if
end do
end subroutine accumulRest
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: is_restart
!
! !INTERFACE:
logical function is_restart( ) 307,6
!
! !DESCRIPTION:
! Determine if restart run
!
! !USES:
use clm_varctl, only : nsrest
!
! !ARGUMENTS:
implicit none
!
! !CALLED FROM:
! subroutine initialize in this module
!
! !REVISION HISTORY:
! Created by Mariana Vertenstein
!
!EOP
!-----------------------------------------------------------------------
if (nsrest == 1) then
is_restart = .true.
else
is_restart = .false.
end if
end function is_restart
end module accumulMod