!===============================================================================
! 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