CESM_Driver: shr_tInterp_mod.F90 Source File

00001 !===============================================================================
00002 ! SVN $Id: shr_tInterp_mod.F90 22636 2010-04-28 22:48:43Z tcraig $
00003 ! 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 $
00004 !===============================================================================
00005 !BOP ===========================================================================
00006 !
00007 ! !MODULE: shr_tInterp_mod -- time interpolation routines
00008 !
00009 ! !DESCRIPTION:
00010 !     shared time interpolation factor routines
00011 !
00012 ! !REVISION HISTORY:
00013 !     2004-Dec-10 - J. Schramm - first version
00014 !     2005-Apr-10 - T. Craig - updated for shr bundles
00015 !
00016 ! !INTERFACE: ------------------------------------------------------------------
00017 
00018 module shr_tInterp_mod
00019  
00020 ! !USES:
00021 
00022    use shr_sys_mod   ! shared system calls
00023    use shr_cal_mod   ! shared calendar type and methods
00024    use shr_kind_mod  ! kinds for strong typing
00025    use shr_const_mod ! shared constants
00026    use shr_log_mod, only: s_loglev  => shr_log_Level
00027    use shr_log_mod, only: s_logunit => shr_log_Unit
00028    use shr_orb_mod, only: shr_orb_cosz, shr_orb_decl
00029 
00030    implicit none
00031 
00032    private ! except
00033 
00034 ! !PUBLIC TYPES:
00035 
00036    ! no public types
00037 
00038 ! !PUBLIC MEMBER FUNCTIONS:
00039 
00040    public :: shr_tInterp_getFactors  ! get time-interp factors
00041    public :: shr_tInterp_getAvgCosz  ! get cosz, time avg of
00042    public :: shr_tInterp_getCosz     ! get cosz 
00043    public :: shr_tInterp_setAbort    ! set abort on error
00044    public :: shr_tInterp_setDebug    ! set debug level
00045    public :: shr_tInterp_getDebug    ! get debug level
00046 
00047 ! !PUBLIC DATA MEMBERS:
00048 
00049    ! no public data
00050 
00051 !EOP
00052 
00053    real(SHR_KIND_R8),parameter :: c0 = 0.0_SHR_KIND_R8
00054    real(SHR_KIND_R8),parameter :: c1 = 1.0_SHR_KIND_R8
00055    real(SHR_KIND_R8),parameter :: eps = 1.0E-12_SHR_KIND_R8
00056 
00057    logical ,save               :: doabort = .true.
00058    integer ,save               :: debug = 0
00059 
00060 !===============================================================================
00061 contains
00062 !===============================================================================
00063 
00064 !===============================================================================
00065 !BOP ===========================================================================
00066 !
00067 ! !IROUTINE: shr_tInterp_getFactors -- calculate time interpolation factors
00068 !
00069 ! !DESCRIPTION:
00070 !     Returns two interpolation factors
00071 !     Legal algorithms are (algo):
00072 !       lower   - sets factors to data 1 (lower-bound), f1=1, f2=0
00073 !       upper   - sets factors to data 2 (upper-bound), f1=0, f2=1
00074 !       nearest - sets factors to nearest data in time
00075 !       linear  - sets factors to linear interpolation between lb and ub
00076 !     \newline
00077 !     call shr\_tInterp\_getFactors(D1,s1,D2,s2,D,s,f1,f2,'linear',rc)
00078 !     \newline
00079 !     time of 2 >= time of 1  for all algos
00080 !     time of 2 >= time of model data >= time of 1  for linear
00081 !
00082 ! !REVISION HISTORY:
00083 !     2005-Apr-10 - T. Craig - first version
00084 !
00085 ! !INTERFACE: ------------------------------------------------------------------
00086 
00087 subroutine shr_tInterp_getFactors(D1,S1,D2,S2,Din,Sin,f1,f2,algo,rc)
00088 
00089    implicit none
00090 
00091 ! !USES:
00092 
00093 ! !INPUT/OUTPUT PARAMETERS:
00094 
00095    integer(SHR_KIND_IN),intent(in)           :: D1,S1   ! LB date & sec (20010115,3600)
00096    integer(SHR_KIND_IN),intent(in)           :: D2,S2   ! UB date & sec
00097    integer(SHR_KIND_IN),intent(in)           :: Din,Sin ! desired/model date & sec
00098    real(SHR_KIND_R8)   ,intent(out)          :: f1      ! wgt for 1
00099    real(SHR_KIND_R8)   ,intent(out)          :: f2      ! wgt for 2
00100    character(*)        ,intent(in) ,optional :: algo    ! algorithm
00101    integer(SHR_KIND_IN),intent(out),optional :: rc      ! return code
00102 
00103 !EOP
00104 
00105    !----- local  ------
00106    real(SHR_KIND_R8)      :: etime1      ! elapsed days for 1 (lower-bound)
00107    real(SHR_KIND_R8)      :: etime2      ! elapsed days for 2 (upper-bound)
00108    real(SHR_KIND_R8)      :: etimeIn     ! elapsed days for model date
00109    integer(SHR_KIND_IN)   :: eday        ! elapsed days
00110    character(SHR_KIND_CS) :: lalgo       ! local algo variable
00111    integer(SHR_KIND_IN)   :: lrc         ! local rc
00112 
00113    !----- formats -----
00114    character(*),parameter :: subName = "(shr_tInterp_getFactors) "
00115    character(*),parameter :: F00   = "('(shr_tInterp_getFactors) ',8a)" 
00116    character(*),parameter :: F01   = "('(shr_tInterp_getFactors) ',a,2f17.8)" 
00117    character(*),parameter :: F02   = "('(shr_tInterp_getFactors) ',a,3f17.8)" 
00118    character(*),parameter :: F03   = "('(shr_tInterp_getFactors) ',2a,3(i9.8,i6))" 
00119 
00120 !-------------------------------------------------------------------------------
00121 ! Computes time interpolation factors
00122 !-------------------------------------------------------------------------------
00123 
00124    lrc = 0
00125 
00126    if (present(algo)) then
00127      lalgo = algo
00128    else
00129      lalgo = 'linear'
00130    endif
00131 
00132    !--- compute elapsed time ---
00133    call shr_cal_date2eday(D1, eday)
00134    etime1 = real(eday,SHR_KIND_R8) + real(s1,SHR_KIND_R8)/shr_const_cDay
00135 
00136    call shr_cal_date2eday(D2, eday)
00137    etime2 = real(eday,SHR_KIND_R8) + real(s2,SHR_KIND_R8)/shr_const_cDay
00138 
00139    call shr_cal_date2eday(Din, eday)
00140    etimein = real(eday,SHR_KIND_R8) + real(sin,SHR_KIND_R8)/shr_const_cDay
00141 
00142    !DML
00143 !  write(s_logunit,*) subName,' ETIME ',etime1,etime2,etimein
00144 
00145    ! --- always check that 1 <= 2, although we could relax this requirement ---
00146    if (etime2 < etime1) then
00147      if (s_loglev > 0) write(s_logunit,F01) ' ERROR: etime2 < etime1 ',etime2,etime1
00148      lrc = 1
00149      call shr_tInterp_abort(subName//' etime2 < etime1 ')
00150    endif
00151 
00152    f1 = c0
00153    f2 = c0
00154    ! --- set interpolation factors ---
00155    if (trim(lalgo) == 'lower') then
00156      if (etime1 < etime2) then
00157        f1 = c1
00158      else
00159        f2 = c1
00160      endif
00161    elseif (trim(lalgo) == 'upper') then
00162      if (etime1 < etime2) then
00163        f2 = c1
00164      else
00165        f1 = c1
00166      endif
00167    elseif (trim(lalgo) == 'nearest') then
00168      if (abs(etimein-etime1) <= abs(etime2-etimein)) then
00169        f1 = c1
00170        f2 = c0
00171      else
00172        f1 = c0
00173        f2 = c1
00174      endif
00175    elseif (trim(lalgo) == 'linear') then
00176      !--- check that etimein is between etime1 and etime2 ---
00177      if (etime2 < etimein .or. etime1 > etimein) then
00178        write(s_logunit,F02) ' ERROR illegal linear times: ',etime1,etimein,etime2
00179        lrc = 1
00180        call shr_tInterp_abort(subName//' illegal etimes ')
00181      endif
00182      if (etime2 == etime1) then
00183        f1 = 0.5_SHR_KIND_R8
00184        f2 = 0.5_SHR_KIND_R8
00185      else
00186        f1 = (etime2-etimein)/(etime2-etime1)
00187        f2 = c1 - f1
00188      endif
00189    else
00190      if (s_loglev > 0) write(s_logunit,F00) 'ERROR: illegal lalgo option: ',trim(lalgo)
00191      lrc = 1
00192      call shr_tInterp_abort(subName//' illegal algo option '//trim(lalgo))
00193    endif
00194 
00195    !--- check that f1 and f2 are OK, each between 0 and 1 and they sum to 1 ---
00196    if (f1 < c0-eps .or. f1 > c1+eps .or. &
00197        f2 < c0-eps .or. f2 > c1+eps .or. &
00198        abs(f1+f2-c1) > eps) then
00199      if (s_loglev > 0) write(s_logunit,F01) 'ERROR: illegal tInterp values ',f1,f2
00200      lrc = 1
00201      call shr_tInterp_abort(subName//' illegal tInterp values ')
00202    endif
00203 
00204    if (debug > 0 .and. s_loglev > 0) then
00205      write(s_logunit,F03) 'DEBUG: algo,D1,S1,Din,Sin,D2,S2=',trim(lAlgo),D1,S1,Din,Sin,D2,S2
00206      write(s_logunit,F01) 'DEBUG: algo,f1,f2= '//trim(lAlgo),f1,f2
00207    endif
00208 
00209    if (present(rc)) rc = lrc
00210 
00211 end subroutine shr_tInterp_getFactors
00212 
00213 !===============================================================================
00214 !BOP ===========================================================================
00215 !
00216 ! !IROUTINE: shr_tInterp_getAvgCosz -- returns avg cos(z) between LB and UB
00217 !
00218 ! !DESCRIPTION:
00219 !
00220 !   Returns a time-average of cos(z) over the time interval [LB,UB].  
00221 !   The time-avg is calculated via a right-sum Riemann sum where the partitian 
00222 !   width is the model dt, and the left-most partitian starts at the LB.
00223 !
00224 ! !REVISION HISTORY:
00225 !     2010-Apr - B. Kauffman - change to t-avg cosz computation, uses model dt
00226 !     2009-Oct - T. Craig - migrated from dshr code
00227 !     2008-Jun - E. Kluzek - first version
00228 !
00229 ! !INTERFACE: ------------------------------------------------------------------
00230 
00231 subroutine shr_tInterp_getAvgCosz(tavCosz,lonr,latr,ymd1,tod1,ymd2,tod2,eccen,mvelpp,lambm0,obliqr,dt)
00232 
00233    implicit none
00234 
00235 ! !USES:
00236 
00237 ! !INPUT/OUTPUT PARAMETERS:
00238 
00239    real(SHR_KIND_R8)   ,intent(out)  :: tavCosz(:) ! t-avg of cosz over [LB,UB]
00240    real(SHR_KIND_R8)   ,intent(in)   :: latr(:)    ! latitude
00241    real(SHR_KIND_R8)   ,intent(in)   :: lonr(:)    ! longitude
00242    integer(SHR_KIND_IN),intent(in)   :: ymd1,tod1  ! date of lb
00243    integer(SHR_KIND_IN),intent(in)   :: ymd2,tod2  ! date of ub
00244    real(SHR_KIND_R8)   ,intent(in)   :: eccen      ! orb param
00245    real(SHR_KIND_R8)   ,intent(in)   :: mvelpp     ! orb param
00246    real(SHR_KIND_R8)   ,intent(in)   :: lambm0     ! orb param
00247    real(SHR_KIND_R8)   ,intent(in)   :: obliqr     ! orb param
00248    integer(SHR_KIND_IN),intent(in)   :: dt         ! model time step in secs
00249 
00250 !EOP
00251 
00252    !----- local  ------
00253    integer(SHR_KIND_IN) :: lsize              ! size of local data
00254    integer(SHR_KIND_IN) ::  eday1, eday2      ! LB,UB elapsed days, integer
00255    real(SHR_KIND_R8)    :: reday1,reday2      ! LB,UB elapsed days, including fraction
00256    integer(SHR_KIND_IN) :: n                  ! number of time-samples in t-avg
00257    integer(SHR_KIND_IN) ::  eday              ! time sample's elapsed days, integer
00258    real(SHR_KIND_R8)    :: reday              ! time sample's elapsed days, including fraction
00259    integer(SHR_KIND_IN) :: ymd,tod,ymd0,tod0  ! used to compute time of time-sample
00260    real(SHR_KIND_R8),pointer :: cosz(:)       ! cos(zenith angle)
00261 
00262    !----- formats -----
00263    character(*),parameter :: subName = "(shr_tInterp_getAvgCosz) "
00264    character(*),parameter :: F00   = "('(shr_tInterp_getAvgCosz) ',8a)" 
00265    character(*),parameter :: F01   = "('(shr_tInterp_getAvgCosz) ',a,2f17.8)" 
00266    character(*),parameter :: F02   = "('(shr_tInterp_getAvgCosz) ',a,3f17.8)" 
00267    character(*),parameter :: F03   = "('(shr_tInterp_getAvgCosz) ',2a,3(i9.8,i6))" 
00268 
00269 !-------------------------------------------------------------------------------
00270 ! Computes time avg cosz over interval [LB,UB]
00271 !-------------------------------------------------------------------------------
00272 
00273    lsize = size(lonr)
00274    allocate(cosz(lsize))
00275    if (lsize < 1 .or. size(latr) /= lsize .or. size(tavCosz) /= lsize) then
00276       call shr_sys_abort(subname//' ERROR: lon lat tavCosz sizes disagree')
00277    endif
00278 
00279    !--- get LB & UB dates ---
00280    call shr_cal_date2eday(ymd1,eday1)
00281    call shr_cal_date2eday(ymd2,eday2)
00282    reday1 = eday1 + tod1/SHR_CONST_CDAY
00283    reday2 = eday2 + tod2/SHR_CONST_CDAY
00284    if (reday1 > reday2) call shr_sys_abort(subname//'ERROR: lower-bound > upper-bound')
00285 
00286    !--- compute time average ---
00287    tavCosz = 0.0_SHR_KIND_R8 ! initialize partial sum
00288    n       = 0               ! number of samples in t-avg
00289    reday   = reday1          ! mid [LB,UB] interval t-step starts at LB 
00290    ymd     = ymd1
00291    tod     = tod1
00292    do while( reday < reday2) ! mid-interval t-steps thru interval [LB,UB]
00293 
00294       !--- advance to next time in [LB,UB] ---
00295       n = n + 1
00296       ymd0 = ymd
00297       tod0 = tod
00298       call shr_cal_advDateInt(dt,'seconds',ymd0,tod0,ymd,tod)
00299       call shr_cal_date2eday(ymd ,eday ) ! convert to elapsed days
00300       reday = eday + tod/SHR_CONST_CDAY  ! add in partial day
00301 
00302       !--- get next cosz value for t-avg ---
00303       call shr_tInterp_getCosz(cosz,lonr,latr,ymd,tod,eccen,mvelpp,lambm0,obliqr)
00304       tavCosz = tavCosz + cosz ! add to partial sum
00305 
00306    end do
00307    tavCosz = tavCosz/real(n,SHR_KIND_R8) ! form t-avg
00308 
00309    deallocate( cosz )
00310 
00311 end subroutine shr_tInterp_getAvgCosz
00312 
00313 !===============================================================================
00314 !BOP ===========================================================================
00315 !
00316 ! !IROUTINE: shr_tInterp_getCosz -- calculate cosine(solar-zenith angle).
00317 !
00318 ! !DESCRIPTION:
00319 !
00320 !    Calculate the cos(solar-zenith angle).
00321 !
00322 ! !REVISION HISTORY:
00323 !     2010-Apr - B. Kauffman - returns cosz
00324 !     2009-Oct - T. Craig - added
00325 !
00326 ! !INTERFACE: ------------------------------------------------------------------
00327 
00328 subroutine shr_tInterp_getCosz(cosz,lonr,latr,ymd,tod,eccen,mvelpp,lambm0,obliqr)
00329 
00330    implicit none
00331 
00332 ! !USES:
00333 
00334 ! !INPUT/OUTPUT PARAMETERS:
00335 
00336    real(SHR_KIND_R8)   ,intent(out)   :: cosz(:)    ! cos(zenith angle)
00337    real(SHR_KIND_R8)   ,intent(in)    :: latr(:)    ! latitude
00338    real(SHR_KIND_R8)   ,intent(in)    :: lonr(:)    ! longitude
00339    integer(SHR_KIND_IN),intent(in)    :: ymd,tod    ! date of interest
00340    real(SHR_KIND_R8)   ,intent(in)    :: eccen      ! orb param
00341    real(SHR_KIND_R8)   ,intent(in)    :: mvelpp     ! orb param
00342    real(SHR_KIND_R8)   ,intent(in)    :: lambm0     ! orb param
00343    real(SHR_KIND_R8)   ,intent(in)    :: obliqr     ! orb param
00344 
00345 !EOP
00346 
00347    !----- local  ------
00348    integer(SHR_KIND_IN) :: n
00349    integer(SHR_KIND_IN) :: lsize
00350    real(SHR_KIND_R8)    :: calday      ! julian days
00351    real(SHR_KIND_R8)    :: declin,eccf ! orb params
00352 
00353    real(SHR_KIND_R8),parameter :: solZenMin = 0.001_SHR_KIND_R8 ! min solar zenith angle
00354 
00355    !----- formats -----
00356    character(*),parameter :: subName = "(shr_tInterp_getCosz) "
00357    character(*),parameter :: F00   = "('(shr_tInterp_getCosz) ',8a)" 
00358    character(*),parameter :: F01   = "('(shr_tInterp_getCosz) ',a,2f17.8)" 
00359    character(*),parameter :: F02   = "('(shr_tInterp_getCosz) ',a,3f17.8)" 
00360    character(*),parameter :: F03   = "('(shr_tInterp_getCosz) ',2a,3(i9.8,i6))" 
00361 
00362 !-------------------------------------------------------------------------------
00363 ! Returns cos(zenith angle)
00364 !-------------------------------------------------------------------------------
00365 
00366    lsize = size(lonr)
00367    if (lsize < 1 .or. size(latr) /= lsize .or. size(cosz) /= lsize) then
00368       call shr_sys_abort(subname//' ERROR: lon lat cosz sizes disagree')
00369    endif
00370 
00371    call shr_cal_date2julian(ymd,tod,calday)
00372    call shr_orb_decl( calday, eccen, mvelpp, lambm0, obliqr, declin, eccf)
00373    do n = 1,lsize
00374       cosz(n) = max(solZenMin, shr_orb_cosz( calday, latr(n), lonr(n), declin ))
00375    end do
00376 
00377 end subroutine shr_tInterp_getCosz
00378 
00379 !===============================================================================
00380 !BOP ===========================================================================
00381 !
00382 ! !IROUTINE: shr_tInterp_setAbort -- Set local shr_tInterp abort flag
00383 !
00384 ! !DESCRIPTION:
00385 !     Set local shr_tInterp abort flag, true = abort, false = print and continue
00386 !     \newline
00387 !     call shr\_tInterp\_setAbort(.false.)
00388 !
00389 ! !REVISION HISTORY:
00390 !     2005-Apr-10  - T. Craig - first prototype
00391 !
00392 ! !INTERFACE: ------------------------------------------------------------------
00393 
00394 subroutine shr_tInterp_setAbort(flag)
00395 
00396   implicit none
00397 
00398 ! !INPUT/OUTPUT PARAMETERS:
00399 
00400   logical,intent(in) :: flag
00401 
00402 !EOP
00403 
00404   !--- formats ---
00405   character(*),parameter :: subName =   "(shr_tInterp_setAbort) "
00406   character(*),parameter :: F00     = "('(shr_tInterp_setAbort) ',a) "
00407 
00408 !-------------------------------------------------------------------------------
00409 !
00410 !-------------------------------------------------------------------------------
00411 
00412   doabort = flag
00413 
00414 end subroutine shr_tInterp_setAbort
00415 
00416 !===============================================================================
00417 !XXBOP ===========================================================================
00418 !
00419 ! !IROUTINE: shr_tInterp_abort -- local interface for abort
00420 !
00421 ! !DESCRIPTION:
00422 !     Local interface for shr\_tInterp abort calls
00423 !     \newline
00424 !     call shr\_tInterp\_abort(subName//' ERROR illegal option')
00425 !
00426 ! !REVISION HISTORY:
00427 !     2005-Apr-10  - T. Craig - first prototype
00428 !
00429 ! !INTERFACE: ------------------------------------------------------------------
00430 
00431 subroutine shr_tInterp_abort(string)
00432 
00433   implicit none
00434 
00435 ! !INPUT/OUTPUT PARAMETERS:
00436 
00437   character(*),optional,intent(in) :: string
00438 
00439 !XXEOP
00440 
00441   !--- formats ---
00442   character(SHR_KIND_CL) :: lstring
00443   character(*),parameter :: subName =   "(shr_tInterp_abort) "
00444   character(*),parameter :: F00     = "('(shr_tInterp_abort) ',a) "
00445 
00446 !-------------------------------------------------------------------------------
00447 !
00448 !-------------------------------------------------------------------------------
00449 
00450   lstring = ''
00451   if (present(string)) lstring = string
00452 
00453   if (doabort) then
00454     call shr_sys_abort(lstring)
00455   else
00456     write(s_logunit,F00) trim(lstring)
00457   endif
00458 
00459 end subroutine shr_tInterp_abort
00460 
00461 !===============================================================================
00462 !BOP ===========================================================================
00463 !
00464 ! !IROUTINE: shr_tInterp_getDebug -- Get local shr_tInterp debug level
00465 !
00466 ! !DESCRIPTION:
00467 !     Get local shr_tInterp debug level, 0 = production
00468 !     \newline
00469 !     call shr\_tInterp\_getDebug(level)
00470 !
00471 ! !REVISION HISTORY:
00472 !     2005-Jun-14  - B. Kauffman
00473 !
00474 ! !INTERFACE: ------------------------------------------------------------------
00475 
00476 subroutine shr_tInterp_getDebug(level)
00477 
00478   implicit none
00479 
00480 ! !INPUT/OUTPUT PARAMETERS:
00481 
00482   integer,intent(out) :: level
00483 
00484 !EOP
00485 
00486   !--- formats ---
00487   character(*),parameter :: subName =   "(shr_tInterp_getDebug) "
00488   character(*),parameter :: F00     = "('(shr_tInterp_getDebug) ',a) "
00489 
00490 !-------------------------------------------------------------------------------
00491 !
00492 !-------------------------------------------------------------------------------
00493 
00494   level = debug
00495 
00496 end subroutine shr_tInterp_getDebug
00497 
00498 !===============================================================================
00499 !BOP ===========================================================================
00500 !
00501 ! !IROUTINE: shr_tInterp_setDebug -- Set local shr_tInterp debug level
00502 !
00503 ! !DESCRIPTION:
00504 !     Set local shr_tInterp debug level, 0 = production
00505 !     \newline
00506 !     call shr\_tInterp\_setDebug(2)
00507 !
00508 ! !REVISION HISTORY:
00509 !     2005-Apr-10  - T. Craig - first prototype
00510 !
00511 ! !INTERFACE: ------------------------------------------------------------------
00512 
00513 subroutine shr_tInterp_setDebug(iflag)
00514 
00515   implicit none
00516 
00517 ! !INPUT/OUTPUT PARAMETERS:
00518 
00519   integer,intent(in) :: iflag
00520 
00521 !EOP
00522 
00523   !--- formats ---
00524   character(*),parameter :: subName =   "(shr_tInterp_setDebug) "
00525   character(*),parameter :: F01     = "('(shr_tInterp_setDebug) ',a,i3) "
00526 
00527 !-------------------------------------------------------------------------------
00528 !
00529 !-------------------------------------------------------------------------------
00530 
00531   debug = iflag
00532   if (debug>0 .and. s_loglev > 0) write(s_logunit,F01) "DEBUG: level changed to ",debug
00533 
00534 end subroutine shr_tInterp_setDebug
00535 
00536 !===============================================================================
00537 !===============================================================================
00538 
00539 end module shr_tInterp_mod