!===============================================================================
! SVN $Id: shr_tInterp_mod.F90 22636 2010-04-28 22:48:43Z tcraig $
! SVN $URL: https://svn-ccsm-models.cgd.ucar.edu/csm_share/branch_tags/cesm1_0_rel_tags/cesm1_0_rel01_share3_100616/shr/shr_tInterp_mod.F90 $
!===============================================================================
!BOP ===========================================================================
!
! !MODULE: shr_tInterp_mod -- time interpolation routines
!
! !DESCRIPTION:
! shared time interpolation factor routines
!
! !REVISION HISTORY:
! 2004-Dec-10 - J. Schramm - first version
! 2005-Apr-10 - T. Craig - updated for shr bundles
!
! !INTERFACE: ------------------------------------------------------------------
module shr_tInterp_mod 1,7
! !USES:
use shr_sys_mod
! shared system calls
use shr_cal_mod
! shared calendar type and methods
use shr_kind_mod
! kinds for strong typing
use shr_const_mod
! shared constants
use shr_log_mod
, only: s_loglev => shr_log_Level
use shr_log_mod
, only: s_logunit => shr_log_Unit
use shr_orb_mod
, only: shr_orb_cosz, shr_orb_decl
implicit none
private ! except
! !PUBLIC TYPES:
! no public types
! !PUBLIC MEMBER FUNCTIONS:
public :: shr_tInterp_getFactors ! get time-interp factors
public :: shr_tInterp_getAvgCosz ! get cosz, time avg of
public :: shr_tInterp_getCosz ! get cosz
public :: shr_tInterp_setAbort ! set abort on error
public :: shr_tInterp_setDebug ! set debug level
public :: shr_tInterp_getDebug ! get debug level
! !PUBLIC DATA MEMBERS:
! no public data
!EOP
real(SHR_KIND_R8),parameter :: c0 = 0.0_SHR_KIND_R8
real(SHR_KIND_R8),parameter :: c1 = 1.0_SHR_KIND_R8
real(SHR_KIND_R8),parameter :: eps = 1.0E-12_SHR_KIND_R8
logical ,save :: doabort = .true.
integer ,save :: debug = 0
!===============================================================================
contains
!===============================================================================
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_tInterp_getFactors -- calculate time interpolation factors
!
! !DESCRIPTION:
! Returns two interpolation factors
! Legal algorithms are (algo):
! lower - sets factors to data 1 (lower-bound), f1=1, f2=0
! upper - sets factors to data 2 (upper-bound), f1=0, f2=1
! nearest - sets factors to nearest data in time
! linear - sets factors to linear interpolation between lb and ub
! \newline
! call shr\_tInterp\_getFactors(D1,s1,D2,s2,D,s,f1,f2,'linear',rc)
! \newline
! time of 2 >= time of 1 for all algos
! time of 2 >= time of model data >= time of 1 for linear
!
! !REVISION HISTORY:
! 2005-Apr-10 - T. Craig - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_tInterp_getFactors(D1,S1,D2,S2,Din,Sin,f1,f2,algo,rc) 1,7
implicit none
! !USES:
! !INPUT/OUTPUT PARAMETERS:
integer(SHR_KIND_IN),intent(in) :: D1,S1 ! LB date & sec (20010115,3600)
integer(SHR_KIND_IN),intent(in) :: D2,S2 ! UB date & sec
integer(SHR_KIND_IN),intent(in) :: Din,Sin ! desired/model date & sec
real(SHR_KIND_R8) ,intent(out) :: f1 ! wgt for 1
real(SHR_KIND_R8) ,intent(out) :: f2 ! wgt for 2
character(*) ,intent(in) ,optional :: algo ! algorithm
integer(SHR_KIND_IN),intent(out),optional :: rc ! return code
!EOP
!----- local ------
real(SHR_KIND_R8) :: etime1 ! elapsed days for 1 (lower-bound)
real(SHR_KIND_R8) :: etime2 ! elapsed days for 2 (upper-bound)
real(SHR_KIND_R8) :: etimeIn ! elapsed days for model date
integer(SHR_KIND_IN) :: eday ! elapsed days
character(SHR_KIND_CS) :: lalgo ! local algo variable
integer(SHR_KIND_IN) :: lrc ! local rc
!----- formats -----
character(*),parameter :: subName = "(shr_tInterp_getFactors) "
character(*),parameter :: F00 = "('(shr_tInterp_getFactors) ',8a)"
character(*),parameter :: F01 = "('(shr_tInterp_getFactors) ',a,2f17.8)"
character(*),parameter :: F02 = "('(shr_tInterp_getFactors) ',a,3f17.8)"
character(*),parameter :: F03 = "('(shr_tInterp_getFactors) ',2a,3(i9.8,i6))"
!-------------------------------------------------------------------------------
! Computes time interpolation factors
!-------------------------------------------------------------------------------
lrc = 0
if (present(algo)) then
lalgo = algo
else
lalgo = 'linear'
endif
!--- compute elapsed time ---
call shr_cal_date2eday
(D1, eday)
etime1 = real(eday,SHR_KIND_R8) + real(s1,SHR_KIND_R8)/shr_const_cDay
call shr_cal_date2eday
(D2, eday)
etime2 = real(eday,SHR_KIND_R8) + real(s2,SHR_KIND_R8)/shr_const_cDay
call shr_cal_date2eday
(Din, eday)
etimein = real(eday,SHR_KIND_R8) + real(sin,SHR_KIND_R8)/shr_const_cDay
!DML
! write(s_logunit,*) subName,' ETIME ',etime1,etime2,etimein
! --- always check that 1 <= 2, although we could relax this requirement ---
if (etime2 < etime1) then
if (s_loglev > 0) write(s_logunit,F01) ' ERROR: etime2 < etime1 ',etime2,etime1
lrc = 1
call shr_tInterp_abort
(subName//' etime2 < etime1 ')
endif
f1 = c0
f2 = c0
! --- set interpolation factors ---
if (trim(lalgo) == 'lower') then
if (etime1 < etime2) then
f1 = c1
else
f2 = c1
endif
elseif (trim(lalgo) == 'upper') then
if (etime1 < etime2) then
f2 = c1
else
f1 = c1
endif
elseif (trim(lalgo) == 'nearest') then
if (abs(etimein-etime1) <= abs(etime2-etimein)) then
f1 = c1
f2 = c0
else
f1 = c0
f2 = c1
endif
elseif (trim(lalgo) == 'linear') then
!--- check that etimein is between etime1 and etime2 ---
if (etime2 < etimein .or. etime1 > etimein) then
write(s_logunit,F02) ' ERROR illegal linear times: ',etime1,etimein,etime2
lrc = 1
call shr_tInterp_abort
(subName//' illegal etimes ')
endif
if (etime2 == etime1) then
f1 = 0.5_SHR_KIND_R8
f2 = 0.5_SHR_KIND_R8
else
f1 = (etime2-etimein)/(etime2-etime1)
f2 = c1 - f1
endif
else
if (s_loglev > 0) write(s_logunit,F00) 'ERROR: illegal lalgo option: ',trim(lalgo)
lrc = 1
call shr_tInterp_abort
(subName//' illegal algo option '//trim(lalgo))
endif
!--- check that f1 and f2 are OK, each between 0 and 1 and they sum to 1 ---
if (f1 < c0-eps .or. f1 > c1+eps .or. &
f2 < c0-eps .or. f2 > c1+eps .or. &
abs(f1+f2-c1) > eps) then
if (s_loglev > 0) write(s_logunit,F01) 'ERROR: illegal tInterp values ',f1,f2
lrc = 1
call shr_tInterp_abort
(subName//' illegal tInterp values ')
endif
if (debug > 0 .and. s_loglev > 0) then
write(s_logunit,F03) 'DEBUG: algo,D1,S1,Din,Sin,D2,S2=',trim(lAlgo),D1,S1,Din,Sin,D2,S2
write(s_logunit,F01) 'DEBUG: algo,f1,f2= '//trim(lAlgo),f1,f2
endif
if (present(rc)) rc = lrc
end subroutine shr_tInterp_getFactors
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_tInterp_getAvgCosz -- returns avg cos(z) between LB and UB
!
! !DESCRIPTION:
!
! Returns a time-average of cos(z) over the time interval [LB,UB].
! The time-avg is calculated via a right-sum Riemann sum where the partitian
! width is the model dt, and the left-most partitian starts at the LB.
!
! !REVISION HISTORY:
! 2010-Apr - B. Kauffman - change to t-avg cosz computation, uses model dt
! 2009-Oct - T. Craig - migrated from dshr code
! 2008-Jun - E. Kluzek - first version
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_tInterp_getAvgCosz(tavCosz,lonr,latr,ymd1,tod1,ymd2,tod2,eccen,mvelpp,lambm0,obliqr,dt) 1,7
implicit none
! !USES:
! !INPUT/OUTPUT PARAMETERS:
real(SHR_KIND_R8) ,intent(out) :: tavCosz(:) ! t-avg of cosz over [LB,UB]
real(SHR_KIND_R8) ,intent(in) :: latr(:) ! latitude
real(SHR_KIND_R8) ,intent(in) :: lonr(:) ! longitude
integer(SHR_KIND_IN),intent(in) :: ymd1,tod1 ! date of lb
integer(SHR_KIND_IN),intent(in) :: ymd2,tod2 ! date of ub
real(SHR_KIND_R8) ,intent(in) :: eccen ! orb param
real(SHR_KIND_R8) ,intent(in) :: mvelpp ! orb param
real(SHR_KIND_R8) ,intent(in) :: lambm0 ! orb param
real(SHR_KIND_R8) ,intent(in) :: obliqr ! orb param
integer(SHR_KIND_IN),intent(in) :: dt ! model time step in secs
!EOP
!----- local ------
integer(SHR_KIND_IN) :: lsize ! size of local data
integer(SHR_KIND_IN) :: eday1, eday2 ! LB,UB elapsed days, integer
real(SHR_KIND_R8) :: reday1,reday2 ! LB,UB elapsed days, including fraction
integer(SHR_KIND_IN) :: n ! number of time-samples in t-avg
integer(SHR_KIND_IN) :: eday ! time sample's elapsed days, integer
real(SHR_KIND_R8) :: reday ! time sample's elapsed days, including fraction
integer(SHR_KIND_IN) :: ymd,tod,ymd0,tod0 ! used to compute time of time-sample
real(SHR_KIND_R8),pointer :: cosz(:) ! cos(zenith angle)
!----- formats -----
character(*),parameter :: subName = "(shr_tInterp_getAvgCosz) "
character(*),parameter :: F00 = "('(shr_tInterp_getAvgCosz) ',8a)"
character(*),parameter :: F01 = "('(shr_tInterp_getAvgCosz) ',a,2f17.8)"
character(*),parameter :: F02 = "('(shr_tInterp_getAvgCosz) ',a,3f17.8)"
character(*),parameter :: F03 = "('(shr_tInterp_getAvgCosz) ',2a,3(i9.8,i6))"
!-------------------------------------------------------------------------------
! Computes time avg cosz over interval [LB,UB]
!-------------------------------------------------------------------------------
lsize = size(lonr)
allocate(cosz(lsize))
if (lsize < 1 .or. size(latr) /= lsize .or. size(tavCosz) /= lsize) then
call shr_sys_abort
(subname//' ERROR: lon lat tavCosz sizes disagree')
endif
!--- get LB & UB dates ---
call shr_cal_date2eday
(ymd1,eday1)
call shr_cal_date2eday
(ymd2,eday2)
reday1 = eday1 + tod1/SHR_CONST_CDAY
reday2 = eday2 + tod2/SHR_CONST_CDAY
if (reday1 > reday2) call shr_sys_abort
(subname//'ERROR: lower-bound > upper-bound')
!--- compute time average ---
tavCosz = 0.0_SHR_KIND_R8 ! initialize partial sum
n = 0 ! number of samples in t-avg
reday = reday1 ! mid [LB,UB] interval t-step starts at LB
ymd = ymd1
tod = tod1
do while( reday < reday2) ! mid-interval t-steps thru interval [LB,UB]
!--- advance to next time in [LB,UB] ---
n = n + 1
ymd0 = ymd
tod0 = tod
call shr_cal_advDateInt
(dt,'seconds',ymd0,tod0,ymd,tod)
call shr_cal_date2eday
(ymd ,eday ) ! convert to elapsed days
reday = eday + tod/SHR_CONST_CDAY ! add in partial day
!--- get next cosz value for t-avg ---
call shr_tInterp_getCosz
(cosz,lonr,latr,ymd,tod,eccen,mvelpp,lambm0,obliqr)
tavCosz = tavCosz + cosz ! add to partial sum
end do
tavCosz = tavCosz/real(n,SHR_KIND_R8) ! form t-avg
deallocate( cosz )
end subroutine shr_tInterp_getAvgCosz
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_tInterp_getCosz -- calculate cosine(solar-zenith angle).
!
! !DESCRIPTION:
!
! Calculate the cos(solar-zenith angle).
!
! !REVISION HISTORY:
! 2010-Apr - B. Kauffman - returns cosz
! 2009-Oct - T. Craig - added
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_tInterp_getCosz(cosz,lonr,latr,ymd,tod,eccen,mvelpp,lambm0,obliqr) 2,4
implicit none
! !USES:
! !INPUT/OUTPUT PARAMETERS:
real(SHR_KIND_R8) ,intent(out) :: cosz(:) ! cos(zenith angle)
real(SHR_KIND_R8) ,intent(in) :: latr(:) ! latitude
real(SHR_KIND_R8) ,intent(in) :: lonr(:) ! longitude
integer(SHR_KIND_IN),intent(in) :: ymd,tod ! date of interest
real(SHR_KIND_R8) ,intent(in) :: eccen ! orb param
real(SHR_KIND_R8) ,intent(in) :: mvelpp ! orb param
real(SHR_KIND_R8) ,intent(in) :: lambm0 ! orb param
real(SHR_KIND_R8) ,intent(in) :: obliqr ! orb param
!EOP
!----- local ------
integer(SHR_KIND_IN) :: n
integer(SHR_KIND_IN) :: lsize
real(SHR_KIND_R8) :: calday ! julian days
real(SHR_KIND_R8) :: declin,eccf ! orb params
real(SHR_KIND_R8),parameter :: solZenMin = 0.001_SHR_KIND_R8 ! min solar zenith angle
!----- formats -----
character(*),parameter :: subName = "(shr_tInterp_getCosz) "
character(*),parameter :: F00 = "('(shr_tInterp_getCosz) ',8a)"
character(*),parameter :: F01 = "('(shr_tInterp_getCosz) ',a,2f17.8)"
character(*),parameter :: F02 = "('(shr_tInterp_getCosz) ',a,3f17.8)"
character(*),parameter :: F03 = "('(shr_tInterp_getCosz) ',2a,3(i9.8,i6))"
!-------------------------------------------------------------------------------
! Returns cos(zenith angle)
!-------------------------------------------------------------------------------
lsize = size(lonr)
if (lsize < 1 .or. size(latr) /= lsize .or. size(cosz) /= lsize) then
call shr_sys_abort
(subname//' ERROR: lon lat cosz sizes disagree')
endif
call shr_cal_date2julian
(ymd,tod,calday)
call shr_orb_decl
( calday, eccen, mvelpp, lambm0, obliqr, declin, eccf)
do n = 1,lsize
cosz(n) = max(solZenMin, shr_orb_cosz
( calday, latr(n), lonr(n), declin ))
end do
end subroutine shr_tInterp_getCosz
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_tInterp_setAbort -- Set local shr_tInterp abort flag
!
! !DESCRIPTION:
! Set local shr_tInterp abort flag, true = abort, false = print and continue
! \newline
! call shr\_tInterp\_setAbort(.false.)
!
! !REVISION HISTORY:
! 2005-Apr-10 - T. Craig - first prototype
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_tInterp_setAbort(flag)
implicit none
! !INPUT/OUTPUT PARAMETERS:
logical,intent(in) :: flag
!EOP
!--- formats ---
character(*),parameter :: subName = "(shr_tInterp_setAbort) "
character(*),parameter :: F00 = "('(shr_tInterp_setAbort) ',a) "
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
doabort = flag
end subroutine shr_tInterp_setAbort
!===============================================================================
!XXBOP ===========================================================================
!
! !IROUTINE: shr_tInterp_abort -- local interface for abort
!
! !DESCRIPTION:
! Local interface for shr\_tInterp abort calls
! \newline
! call shr\_tInterp\_abort(subName//' ERROR illegal option')
!
! !REVISION HISTORY:
! 2005-Apr-10 - T. Craig - first prototype
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_tInterp_abort(string) 4,1
implicit none
! !INPUT/OUTPUT PARAMETERS:
character(*),optional,intent(in) :: string
!XXEOP
!--- formats ---
character(SHR_KIND_CL) :: lstring
character(*),parameter :: subName = "(shr_tInterp_abort) "
character(*),parameter :: F00 = "('(shr_tInterp_abort) ',a) "
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
lstring = ''
if (present(string)) lstring = string
if (doabort) then
call shr_sys_abort
(lstring)
else
write(s_logunit,F00) trim(lstring)
endif
end subroutine shr_tInterp_abort
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_tInterp_getDebug -- Get local shr_tInterp debug level
!
! !DESCRIPTION:
! Get local shr_tInterp debug level, 0 = production
! \newline
! call shr\_tInterp\_getDebug(level)
!
! !REVISION HISTORY:
! 2005-Jun-14 - B. Kauffman
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_tInterp_getDebug(level)
implicit none
! !INPUT/OUTPUT PARAMETERS:
integer,intent(out) :: level
!EOP
!--- formats ---
character(*),parameter :: subName = "(shr_tInterp_getDebug) "
character(*),parameter :: F00 = "('(shr_tInterp_getDebug) ',a) "
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
level = debug
end subroutine shr_tInterp_getDebug
!===============================================================================
!BOP ===========================================================================
!
! !IROUTINE: shr_tInterp_setDebug -- Set local shr_tInterp debug level
!
! !DESCRIPTION:
! Set local shr_tInterp debug level, 0 = production
! \newline
! call shr\_tInterp\_setDebug(2)
!
! !REVISION HISTORY:
! 2005-Apr-10 - T. Craig - first prototype
!
! !INTERFACE: ------------------------------------------------------------------
subroutine shr_tInterp_setDebug(iflag)
implicit none
! !INPUT/OUTPUT PARAMETERS:
integer,intent(in) :: iflag
!EOP
!--- formats ---
character(*),parameter :: subName = "(shr_tInterp_setDebug) "
character(*),parameter :: F01 = "('(shr_tInterp_setDebug) ',a,i3) "
!-------------------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
debug = iflag
if (debug>0 .and. s_loglev > 0) write(s_logunit,F01) "DEBUG: level changed to ",debug
end subroutine shr_tInterp_setDebug
!===============================================================================
!===============================================================================
end module shr_tInterp_mod