LCOV - code coverage report
Current view: | top level - utils - interpolate_data.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 221 | 494 | 44.7 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 10 | 15 | 66.7 %

              Line data    Source code
           1             : module interpolate_data
           2             : ! Description:
           3             : !   Routines for interpolation of data in latitude, longitude and time.
           4             : !
           5             : ! Author: Gathered from a number of places and put into the current format by Jim Edwards
           6             : !
           7             : ! Modules Used:
           8             : !
           9             :   use shr_kind_mod,   only: r8 => shr_kind_r8
          10             :   use cam_abortutils, only: endrun
          11             :   use cam_logfile,    only: iulog
          12             :   implicit none
          13             :   private
          14             : !
          15             : ! Public Methods:
          16             : !
          17             :
          18             :   public :: interp_type, lininterp, vertinterp, bilin, get_timeinterp_factors
          19             :   public :: lininterp_init, lininterp_finish
          20             :   type interp_type
          21             :      real(r8), pointer :: wgts(:)
          22             :      real(r8), pointer :: wgtn(:)
          23             :      integer, pointer  :: jjm(:)
          24             :      integer, pointer  :: jjp(:)
          25             :   end type interp_type
          26             :   interface lininterp
          27             :      module procedure lininterp_original
          28             :      module procedure lininterp_full1d
          29             :      module procedure lininterp1d
          30             :      module procedure lininterp2d2d
          31             :      module procedure lininterp2d1d
          32             :      module procedure lininterp3d2d
          33             :   end interface
          34             :
          35             : integer, parameter, public :: extrap_method_zero = 0
          36             : integer, parameter, public :: extrap_method_bndry = 1
          37             : integer, parameter, public :: extrap_method_cycle = 2
          38             :
          39             : contains
          40        1536 :   subroutine lininterp_full1d (arrin, yin, nin, arrout, yout, nout)
          41             :     integer, intent(in) :: nin, nout
          42             :     real(r8), intent(in) :: arrin(nin), yin(nin), yout(nout)
          43             :     real(r8), intent(out) :: arrout(nout)
          44             :     type (interp_type) :: interp_wgts
          45             :
          46        1536 :     call lininterp_init(yin, nin, yout, nout, extrap_method_bndry, interp_wgts)
          47        1536 :     call lininterp1d(arrin, nin, arrout, nout, interp_wgts)
          48        1536 :     call lininterp_finish(interp_wgts)
          49             :
          50        1536 :   end subroutine lininterp_full1d
          51             :
          52   692513712 :   subroutine lininterp_init(yin, nin, yout, nout, extrap_method, interp_wgts, &
          53             :        cyclicmin, cyclicmax)
          54             : !
          55             : ! Description:
          56             : !   Initialize a variable of type(interp_type) with weights for linear interpolation.
          57             : !       this variable can then be used in calls to lininterp1d and lininterp2d.
          58             : !   yin is a 1d array of length nin of locations to interpolate from - this array must
          59             : !       be monotonic but can be increasing or decreasing
          60             : !   yout is a 1d array of length nout of locations to interpolate to, this array need
          61             : !       not be ordered
          62             : !   extrap_method determines how to handle yout points beyond the bounds of yin
          63             : !       if 0 set values outside output grid to 0
          64             : !       if 1 set to boundary value
          65             : !       if 2 set to cyclic boundaries
          66             : !         optional values cyclicmin and cyclicmax can be used to set the bounds of the
          67             : !         cyclic mapping - these default to 0 and 360.
          68             : !
          69             :
          70             :     integer, intent(in) :: nin
          71             :     integer, intent(in) :: nout
          72             :     real(r8), intent(in) :: yin(:)           ! input mesh
          73             :     real(r8), intent(in) :: yout(:)         ! output mesh
          74             :     integer, intent(in) :: extrap_method       ! if 0 set values outside output grid to 0
          75             :                                                ! if 1 set to boundary value
          76             :                                                ! if 2 set to cyclic boundaries
          77             :     real(r8), intent(in), optional :: cyclicmin, cyclicmax
          78             :
          79             :     type (interp_type), intent(out) :: interp_wgts
          80             :     real(r8) :: cmin, cmax
          81             :     real(r8) :: extrap
          82             :     real(r8) :: dyinwrap
          83             :     real(r8) :: ratio
          84             :     real(r8) :: avgdyin
          85             :     integer :: i, j, icount
          86             :     integer :: jj
          87   692513712 :     real(r8), pointer :: wgts(:)
          88   692513712 :     real(r8), pointer :: wgtn(:)
          89   692513712 :     integer, pointer :: jjm(:)
          90   692513712 :     integer, pointer :: jjp(:)
          91             :     logical :: increasing
          92             :     !
          93             :     ! Check validity of input coordinate arrays: must be monotonically increasing,
          94             :     ! and have a total of at least 2 elements
          95             :     !
          96   692513712 :     if (nin.lt.2) then
          97           0 :        call endrun('LININTERP: Must have at least 2 input points for interpolation')
          98             :     end if
          99   692513712 :     if(present(cyclicmin)) then
         100      173376 :        cmin=cyclicmin
         101             :     else
         102             :        cmin=0_r8
         103             :     end if
         104   692513712 :     if(present(cyclicmax)) then
         105      173376 :        cmax=cyclicmax
         106             :     else
         107             :        cmax=360_r8
         108             :     end if
         109   692513712 :     if(cmax<=cmin) then
         110           0 :        call endrun('LININTERP: cyclic min value must be < max value')
         111             :     end if
         112   692513712 :     increasing=.true.
         113   692513712 :     icount = 0
         114 >12860*10^7 :     do j=1,nin-1
         115 >12860*10^7 :        if (yin(j).gt.yin(j+1)) icount = icount + 1
         116             :     end do
         117   692513712 :     if(icount.eq.nin-1) then
         118   149254768 :        increasing = .false.
         119             :        icount=0
         120             :     endif
         121   543258944 :     if (icount.gt.0) then
         122           0 :        call endrun('LININTERP: Non-monotonic input coordinate array found')
         123             :     end if
         124             :     allocate(interp_wgts%jjm(nout), &
         125             :          interp_wgts%jjp(nout), &
         126             :          interp_wgts%wgts(nout), &
         127  4847595984 :          interp_wgts%wgtn(nout))
         128             :
         129   692513712 :     jjm => interp_wgts%jjm
         130   692513712 :     jjp => interp_wgts%jjp
         131   692513712 :     wgts =>  interp_wgts%wgts
         132   692513712 :     wgtn =>  interp_wgts%wgtn
         133             :
         134             :     !
         135             :     ! Initialize index arrays for later checking
         136             :     !
         137  1395005020 :     jjm = 0
         138  1395005020 :     jjp = 0
         139             :
         140   692513712 :     extrap = 0._r8
         141   692516784 :     select case (extrap_method)
         142             :     case (extrap_method_zero)
         143             :        !
         144             :        ! For values which extend beyond boundaries, set weights
         145             :        ! such that values will be 0.
         146             :        !
         147      288768 :        do j=1,nout
         148      288768 :           if(increasing) then
         149           0 :              if (yout(j).lt.yin(1)) then
         150           0 :                 jjm(j) = 1
         151           0 :                 jjp(j) = 1
         152           0 :                 wgts(j) = 0._r8
         153           0 :                 wgtn(j) = 0._r8
         154           0 :                 extrap = extrap + 1._r8
         155           0 :              else if (yout(j).gt.yin(nin)) then
         156           0 :                 jjm(j) = nin
         157           0 :                 jjp(j) = nin
         158           0 :                 wgts(j) = 0._r8
         159           0 :                 wgtn(j) = 0._r8
         160           0 :                 extrap = extrap + 1._r8
         161             :              end if
         162             :           else
         163      285696 :              if (yout(j).gt.yin(1)) then
         164           0 :                 jjm(j) = 1
         165           0 :                 jjp(j) = 1
         166           0 :                 wgts(j) = 0._r8
         167           0 :                 wgtn(j) = 0._r8
         168           0 :                 extrap = extrap + 1._r8
         169      285696 :              else if (yout(j).lt.yin(nin)) then
         170       52224 :                 jjm(j) = nin
         171       52224 :                 jjp(j) = nin
         172       52224 :                 wgts(j) = 0._r8
         173       52224 :                 wgtn(j) = 0._r8
         174       52224 :                 extrap = extrap + 1._r8
         175             :              end if
         176             :           end if
         177             :        end do
         178             :     case (extrap_method_bndry)
         179             :        !
         180             :        ! For values which extend beyond boundaries, set weights
         181             :        ! such that values will just be copied.
         182             :        !
         183  1391821276 :        do j=1,nout
         184  1391821276 :           if(increasing) then
         185   550141308 :              if (yout(j).le.yin(1)) then
         186    19878900 :                 jjm(j) = 1
         187    19878900 :                 jjp(j) = 1
         188    19878900 :                 wgts(j) = 1._r8
         189    19878900 :                 wgtn(j) = 0._r8
         190    19878900 :                 extrap = extrap + 1._r8
         191   530262408 :              else if (yout(j).gt.yin(nin)) then
         192    47049174 :                 jjm(j) = nin
         193    47049174 :                 jjp(j) = nin
         194    47049174 :                 wgts(j) = 1._r8
         195    47049174 :                 wgtn(j) = 0._r8
         196    47049174 :                 extrap = extrap + 1._r8
         197             :              end if
         198             :           else
         199   149342704 :              if (yout(j).gt.yin(1)) then
         200         428 :                 jjm(j) = 1
         201         428 :                 jjp(j) = 1
         202         428 :                 wgts(j) = 1._r8
         203         428 :                 wgtn(j) = 0._r8
         204         428 :                 extrap = extrap + 1._r8
         205   149342276 :              else if (yout(j).le.yin(nin)) then
         206           8 :                 jjm(j) = nin
         207           8 :                 jjp(j) = nin
         208           8 :                 wgts(j) = 1._r8
         209           8 :                 wgtn(j) = 0._r8
         210           8 :                 extrap = extrap + 1._r8
         211             :              end if
         212             :           end if
         213             :        end do
         214             :     case (extrap_method_cycle)
         215             :        !
         216             :        ! For values which extend beyond boundaries, set weights
         217             :        ! for circular boundaries
         218             :        !
         219      173376 :        dyinwrap = yin(1) + (cmax-cmin) - yin(nin)
         220      173376 :        avgdyin = abs(yin(nin)-yin(1))/(nin-1._r8)
         221      173376 :        ratio = dyinwrap/avgdyin
         222      173376 :        if (ratio < 0.9_r8 .or. ratio > 1.1_r8) then
         223           0 :           write(iulog,*) 'Lininterp: Bad dyinwrap value =',dyinwrap,&
         224           0 :                ' avg=', avgdyin, yin(1),yin(nin)
         225           0 :           call endrun('interpolate_data')
         226             :        end if
         227             :
         228   695408688 :        do j=1,nout
         229     2894976 :           if(increasing) then
         230     2721600 :              if (yout(j) <= yin(1)) then
         231           0 :                 jjm(j) = nin
         232           0 :                 jjp(j) = 1
         233           0 :                 wgts(j) = (yin(1)-yout(j))/dyinwrap
         234           0 :                 wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap
         235     2721600 :              else if (yout(j) > yin(nin)) then
         236        7876 :                 jjm(j) = nin
         237        7876 :                 jjp(j) = 1
         238        7876 :                 wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap
         239        7876 :                 wgtn(j) = (yout(j)-yin(nin))/dyinwrap
         240             :              end if
         241             :           else
         242           0 :              if (yout(j) > yin(1)) then
         243           0 :                 jjm(j) = nin
         244           0 :                 jjp(j) = 1
         245           0 :                 wgts(j) = (yin(1)-yout(j))/dyinwrap
         246           0 :                 wgtn(j) = (yout(j)+(cmax-cmin) - yin(nin))/dyinwrap
         247           0 :              else if (yout(j) <= yin(nin)) then
         248           0 :                 jjm(j) = nin
         249           0 :                 jjp(j) = 1
         250           0 :                 wgts(j) = (yin(1)+(cmax-cmin)-yout(j))/dyinwrap
         251           0 :                 wgtn(j) = (yout(j)+(cmax-cmin)-yin(nin))/dyinwrap
         252             :              end if
         253             :
         254             :           endif
         255             :        end do
         256             :     end select
         257             :
         258             :     !
         259             :     ! Loop though output indices finding input indices and weights
         260             :     !
         261   692513712 :     if(increasing) then
         262  1096121852 :        do j=1,nout
         263 47076308796 :           do jj=1,nin-1
         264 46533049852 :              if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then
         265   485926958 :                 jjm(j) = jj
         266   485926958 :                 jjp(j) = jj + 1
         267   485926958 :                 wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj))
         268   485926958 :                 wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj))
         269   485926958 :                 exit
         270             :              end if
         271             :           end do
         272             :        end do
         273             :     else
         274   298883168 :        do j=1,nout
         275  4775225084 :           do jj=1,nin-1
         276  4625970316 :              if (yout(j).le.yin(jj) .and. yout(j).gt.yin(jj+1)) then
         277   149575740 :                 jjm(j) = jj
         278   149575740 :                 jjp(j) = jj + 1
         279   149575740 :                 wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj))
         280   149575740 :                 wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj))
         281   149575740 :                 exit
         282             :              end if
         283             :           end do
         284             :        end do
         285             :     end if
         286             :
         287             :     !
         288             :     ! Check that interp/extrap points have been found for all outputs
         289             :     !
         290   692513712 :     icount = 0
         291  1395005020 :     do j=1,nout
         292   702491308 :        if (jjm(j).eq.0 .or. jjp(j).eq.0) icount = icount + 1
         293   702491308 :        ratio=wgts(j)+wgtn(j)
         294  1395005020 :        if((ratio<0.9_r8.or.ratio>1.1_r8).and.extrap_method.ne.0) then
         295           0 :           write(iulog,*) j, wgts(j),wgtn(j),jjm(j),jjp(j), increasing,extrap_method
         296           0 :           call endrun('Bad weight computed in LININTERP_init')
         297             :        end if
         298             :     end do
         299   692513712 :     if (icount.gt.0) then
         300           0 :        call endrun('LININTERP: Point found without interp indices')
         301             :     end if
         302             :
         303   692513712 :   end subroutine lininterp_init
         304             :
         305 18888844864 :   subroutine lininterp1d (arrin, n1, arrout, m1, interp_wgts)
         306             :     !-----------------------------------------------------------------------
         307             :     !
         308             :     ! Purpose: Do a linear interpolation from input mesh to output
         309             :     !          mesh with weights as set in lininterp_init.
         310             :     !
         311             :     !
         312             :     ! Author: Jim Edwards
         313             :     !
         314             :     !-----------------------------------------------------------------------
         315             :     !-----------------------------------------------------------------------
         316             :     implicit none
         317             :     !-----------------------------------------------------------------------
         318             :     !
         319             :     ! Arguments
         320             :     !
         321             :     integer, intent(in) :: n1                 ! number of input latitudes
         322             :     integer, intent(in) :: m1                ! number of output latitudes
         323             :
         324             :     real(r8), intent(in) :: arrin(n1)    ! input array of values to interpolate
         325             :     type(interp_type), intent(in) :: interp_wgts
         326             :     real(r8), intent(out) :: arrout(m1) ! interpolated array
         327             :
         328             :     !
         329             :     ! Local workspace
         330             :     !
         331             :     integer j                ! latitude indices
         332 18888844864 :     integer, pointer :: jjm(:)
         333 18888844864 :     integer, pointer :: jjp(:)
         334             :
         335 18888844864 :     real(r8), pointer :: wgts(:)
         336 18888844864 :     real(r8), pointer :: wgtn(:)
         337             :
         338             :
         339 18888844864 :     jjm => interp_wgts%jjm
         340 18888844864 :     jjp => interp_wgts%jjp
         341 18888844864 :     wgts =>  interp_wgts%wgts
         342 18888844864 :     wgtn =>  interp_wgts%wgtn
         343             :
         344             :     !
         345             :     ! Do the interpolation
         346             :     !
         347 38055588992 :     do j=1,m1
         348 38055588992 :       arrout(j) = arrin(jjm(j))*wgts(j) + arrin(jjp(j))*wgtn(j)
         349             :     end do
         350             :
         351 18888844864 :     return
         352 18888844864 :   end subroutine lininterp1d
         353             :
         354       61920 :   subroutine lininterp2d2d(arrin, n1, n2, arrout, m1, m2, wgt1, wgt2)
         355             :     implicit none
         356             :     !-----------------------------------------------------------------------
         357             :     !
         358             :     ! Arguments
         359             :     !
         360             :     integer, intent(in) :: n1, n2, m1, m2
         361             :     real(r8), intent(in) :: arrin(n1,n2)    ! input array of values to interpolate
         362             :     type(interp_type), intent(in) :: wgt1, wgt2
         363             :     real(r8), intent(out) :: arrout(m1,m2) ! interpolated array
         364             :     !
         365             :     ! locals
         366             :     !
         367             :     integer i,j                ! indices
         368       61920 :     integer, pointer :: iim(:), jjm(:)
         369       61920 :     integer, pointer :: iip(:), jjp(:)
         370             :
         371       61920 :     real(r8), pointer :: wgts1(:), wgts2(:)
         372       61920 :     real(r8), pointer :: wgtn1(:), wgtn2(:)
         373             :
         374      123840 :     real(r8) :: arrtmp(n1,m2)
         375             :
         376             :
         377       61920 :     jjm => wgt2%jjm
         378       61920 :     jjp => wgt2%jjp
         379       61920 :     wgts2 => wgt2%wgts
         380       61920 :     wgtn2 => wgt2%wgtn
         381             :
         382       61920 :     iim => wgt1%jjm
         383       61920 :     iip => wgt1%jjp
         384       61920 :     wgts1 => wgt1%wgts
         385       61920 :     wgtn1 => wgt1%wgtn
         386             :
         387     5820480 :     do j=1,m2
         388   213128640 :       do i=1,n1
         389   213066720 :         arrtmp(i,j) = arrin(i,jjm(j))*wgts2(j) + arrin(i,jjp(j))*wgtn2(j)
         390             :       end do
         391             :     end do
         392             :
         393     5820480 :     do j=1,m2
         394    94488540 :       do i=1,m1
         395    94426620 :         arrout(i,j) = arrtmp(iim(i),j)*wgts1(i) + arrtmp(iip(i),j)*wgtn1(i)
         396             :       end do
         397             :     end do
         398             :
         399       61920 :   end subroutine lininterp2d2d
         400  4328237264 :   subroutine lininterp2d1d(arrin, n1, n2, arrout, m1, wgt1, wgt2, fldname)
         401             :     implicit none
         402             :     !-----------------------------------------------------------------------
         403             :     !
         404             :     ! Arguments
         405             :     !
         406             :     integer, intent(in) :: n1, n2, m1
         407             :     real(r8), intent(in) :: arrin(n1,n2)    ! input array of values to interpolate
         408             :     type(interp_type), intent(in) :: wgt1, wgt2
         409             :     real(r8), intent(out) :: arrout(m1) ! interpolated array
         410             :     character(len=*), intent(in), optional :: fldname(:)
         411             :     !
         412             :     ! locals
         413             :     !
         414             :     integer i                ! indices
         415  4328237264 :     integer, pointer :: iim(:), jjm(:)
         416  4328237264 :     integer, pointer :: iip(:), jjp(:)
         417             :
         418  4328237264 :     real(r8), pointer :: wgts(:), wgte(:)
         419  4328237264 :     real(r8), pointer :: wgtn(:), wgtw(:)
         420             :
         421  4328237264 :     jjm => wgt2%jjm
         422  4328237264 :     jjp => wgt2%jjp
         423  4328237264 :     wgts => wgt2%wgts
         424  4328237264 :     wgtn => wgt2%wgtn
         425             :
         426  4328237264 :     iim => wgt1%jjm
         427  4328237264 :     iip => wgt1%jjp
         428  4328237264 :     wgtw => wgt1%wgts
         429  4328237264 :     wgte => wgt1%wgtn
         430             :
         431  8658203680 :     do i=1,m1
         432 17319865664 :        arrout(i) = arrin(iim(i),jjm(i))*wgtw(i)*wgts(i)+arrin(iip(i),jjm(i))*wgte(i)*wgts(i) + &
         433 25978069344 :                    arrin(iim(i),jjp(i))*wgtw(i)*wgtn(i)+arrin(iip(i),jjp(i))*wgte(i)*wgtn(i)
         434             :     end do
         435             :
         436             :
         437  4328237264 :   end subroutine lininterp2d1d
         438      123840 :   subroutine lininterp3d2d(arrin, n1, n2, n3, arrout, m1, len1, wgt1, wgt2)
         439             :     implicit none
         440             :     !-----------------------------------------------------------------------
         441             :     !
         442             :     ! Arguments
         443             :     !
         444             :     integer, intent(in) :: n1, n2, n3, m1, len1   ! m1 is to len1 as ncols is to pcols
         445             :     real(r8), intent(in) :: arrin(n1,n2,n3)    ! input array of values to interpolate
         446             :     type(interp_type), intent(in) :: wgt1, wgt2
         447             :     real(r8), intent(out) :: arrout(len1, n3) ! interpolated array
         448             :
         449             :     !
         450             :     ! locals
         451             :     !
         452             :     integer i, k               ! indices
         453      123840 :     integer, pointer :: iim(:), jjm(:)
         454      123840 :     integer, pointer :: iip(:), jjp(:)
         455             :
         456      123840 :     real(r8), pointer :: wgts(:), wgte(:)
         457      123840 :     real(r8), pointer :: wgtn(:), wgtw(:)
         458             :
         459      123840 :     jjm => wgt2%jjm
         460      123840 :     jjp => wgt2%jjp
         461      123840 :     wgts => wgt2%wgts
         462      123840 :     wgtn => wgt2%wgtn
         463             :
         464      123840 :     iim => wgt1%jjm
         465      123840 :     iip => wgt1%jjp
         466      123840 :     wgtw => wgt1%wgts
         467      123840 :     wgte => wgt1%wgtn
         468             :
         469     8792640 :     do k=1,n3
         470   144872640 :        do i=1,m1
         471   680400000 :           arrout(i,k) = arrin(iim(i),jjm(i),k)*wgtw(i)*wgts(i)+arrin(iip(i),jjm(i),k)*wgte(i)*wgts(i) + &
         472   825148800 :                arrin(iim(i),jjp(i),k)*wgtw(i)*wgtn(i)+arrin(iip(i),jjp(i),k)*wgte(i)*wgtn(i)
         473             :        end do
         474             :     end do
         475             :
         476      123840 :   end subroutine lininterp3d2d
         477             :
         478             :
         479             :
         480             :
         481   692513712 :   subroutine lininterp_finish(interp_wgts)
         482             :     type(interp_type) :: interp_wgts
         483             :
         484           0 :     deallocate(interp_wgts%jjm, &
         485           0 :          interp_wgts%jjp, &
         486           0 :          interp_wgts%wgts, &
         487   692513712 :          interp_wgts%wgtn)
         488             :
         489             :     nullify(interp_wgts%jjm, &
         490   692513712 :          interp_wgts%jjp, &
         491   692513712 :          interp_wgts%wgts, &
         492   692513712 :          interp_wgts%wgtn)
         493   692513712 :   end subroutine lininterp_finish
         494             :
         495           0 :   subroutine lininterp_original (arrin, yin, nlev, nlatin, arrout, &
         496           0 :        yout, nlatout)
         497             :     !-----------------------------------------------------------------------
         498             :     !
         499             :     ! Purpose: Do a linear interpolation from input mesh defined by yin to output
         500             :     !          mesh defined by yout.  Where extrapolation is necessary, values will
         501             :     !          be copied from the extreme edge of the input grid.  Vectorization is over
         502             :     !          the vertical (nlev) dimension.
         503             :     !
         504             :     ! Method: Check validity of input, then determine weights, then do the N-S interpolation.
         505             :     !
         506             :     ! Author: Jim Rosinski
         507             :     ! Modified: Jim Edwards so that there is no requirement of monotonacity on the yout array
         508             :     !
         509             :     !-----------------------------------------------------------------------
         510             :     implicit none
         511             :     !-----------------------------------------------------------------------
         512             :     !
         513             :     ! Arguments
         514             :     !
         515             :     integer, intent(in) :: nlev                   ! number of vertical levels
         516             :     integer, intent(in) :: nlatin                 ! number of input latitudes
         517             :     integer, intent(in) :: nlatout                ! number of output latitudes
         518             :
         519             :     real(r8), intent(in) :: arrin(nlev,nlatin)    ! input array of values to interpolate
         520             :     real(r8), intent(in) :: yin(nlatin)           ! input mesh
         521             :     real(r8), intent(in) :: yout(nlatout)         ! output mesh
         522             :
         523             :     real(r8), intent(out) :: arrout(nlev,nlatout) ! interpolated array
         524             :     !
         525             :     ! Local workspace
         526             :     !
         527             :     integer j, jj              ! latitude indices
         528             :     integer jjprev             ! latitude indices
         529             :     integer k                  ! level index
         530             :     integer icount             ! number of values
         531             :
         532             :     real(r8) extrap            ! percent grid non-overlap
         533             :     !
         534             :     ! Dynamic
         535             :     !
         536           0 :     integer :: jjm(nlatout)
         537           0 :     integer :: jjp(nlatout)
         538             :
         539           0 :     real(r8) :: wgts(nlatout)
         540           0 :     real(r8) :: wgtn(nlatout)
         541             :     !
         542             :     ! Check validity of input coordinate arrays: must be monotonically increasing,
         543             :     ! and have a total of at least 2 elements
         544             :     !
         545           0 :     if (nlatin.lt.2) then
         546           0 :        call endrun('LININTERP: Must have at least 2 input points for interpolation')
         547             :     end if
         548             :
         549           0 :     icount = 0
         550           0 :     do j=1,nlatin-1
         551           0 :        if (yin(j).gt.yin(j+1)) icount = icount + 1
         552             :     end do
         553             :
         554             :
         555           0 :     if (icount.gt.0) then
         556           0 :        call endrun('LININTERP: Non-monotonic coordinate array(s) found')
         557             :     end if
         558             :     !
         559             :     ! Initialize index arrays for later checking
         560             :     !
         561           0 :     do j=1,nlatout
         562           0 :        jjm(j) = 0
         563           0 :        jjp(j) = 0
         564             :     end do
         565             :     !
         566             :     ! For values which extend beyond N and S boundaries, set weights
         567             :     ! such that values will just be copied.
         568             :     !
         569             :     extrap = 0._r8
         570             :
         571           0 :     do j=1,nlatout
         572           0 :        if (yout(j).le.yin(1)) then
         573           0 :           jjm(j) = 1
         574           0 :           jjp(j) = 1
         575           0 :           wgts(j) = 1._r8
         576           0 :           wgtn(j) = 0._r8
         577             :           extrap=extrap+1._r8
         578           0 :        else if (yout(j).gt.yin(nlatin)) then
         579           0 :           jjm(j) = nlatin
         580           0 :           jjp(j) = nlatin
         581           0 :           wgts(j) = 1._r8
         582           0 :           wgtn(j) = 0._r8
         583             :           extrap=extrap+1._r8
         584             :        endif
         585             :     end do
         586             :
         587             :     !
         588             :     ! Loop though output indices finding input indices and weights
         589             :     !
         590           0 :     do j=1,nlatout
         591           0 :        do jj=1,nlatin-1
         592           0 :           if (yout(j).gt.yin(jj) .and. yout(j).le.yin(jj+1)) then
         593           0 :              jjm(j) = jj
         594           0 :              jjp(j) = jj + 1
         595           0 :              wgts(j) = (yin(jj+1)-yout(j))/(yin(jj+1)-yin(jj))
         596           0 :              wgtn(j) = (yout(j)-yin(jj))/(yin(jj+1)-yin(jj))
         597           0 :              exit
         598             :           end if
         599             :        end do
         600             :     end do
         601             :     !
         602             :     ! Check that interp/extrap points have been found for all outputs
         603             :     !
         604             :     icount = 0
         605           0 :     do j=1,nlatout
         606           0 :        if (jjm(j).eq.0 .or. jjp(j).eq.0) then
         607           0 :           icount = icount + 1
         608             :        end if
         609             :     end do
         610           0 :     if (icount.gt.0) then
         611           0 :        call endrun('LININTERP: Point found without interp indices')
         612             :     end if
         613             :     !
         614             :     ! Do the interpolation
         615             :     !
         616           0 :     do j=1,nlatout
         617           0 :        do k=1,nlev
         618           0 :           arrout(k,j) = arrin(k,jjm(j))*wgts(j) + arrin(k,jjp(j))*wgtn(j)
         619             :        end do
         620             :     end do
         621             :
         622           0 :     return
         623             :   end subroutine lininterp_original
         624             :
         625             :
         626           0 :   subroutine bilin (arrin, xin, yin, nlondin, nlonin, &
         627           0 :        nlevdin, nlev, nlatin, arrout, xout, &
         628           0 :        yout, nlondout, nlonout, nlevdout, nlatout)
         629             :     !-----------------------------------------------------------------------
         630             :     !
         631             :     ! Purpose:
         632             :     !
         633             :     ! Do a bilinear interpolation from input mesh defined by xin, yin to output
         634             :     ! mesh defined by xout, yout.  Circularity is assumed in the x-direction so
         635             :     ! input x-grid must be in degrees east and must start from Greenwich.  When
         636             :     ! extrapolation is necessary in the N-S direction, values will be copied
         637             :     ! from the extreme edge of the input grid.  Vectorization is over the
         638             :     ! longitude dimension.  Input grid is assumed rectangular. Output grid
         639             :     ! is assumed ragged, i.e. xout is a 2-d mesh.
         640             :     !
         641             :     ! Method: Interpolate first in longitude, then in latitude.
         642             :     !
         643             :     ! Author: Jim Rosinski
         644             :     !
         645             :     !-----------------------------------------------------------------------
         646             :     use shr_kind_mod,     only: r8 => shr_kind_r8
         647             :     use cam_abortutils,   only: endrun
         648             :     !-----------------------------------------------------------------------
         649             :     implicit none
         650             :     !-----------------------------------------------------------------------
         651             :     !
         652             :     ! Input arguments
         653             :     !
         654             :     integer, intent(in) :: nlondin                        ! longitude dimension of input grid
         655             :     integer, intent(in) :: nlonin                         ! number of real longitudes (input)
         656             :     integer, intent(in) :: nlevdin                        ! vertical dimension of input grid
         657             :     integer, intent(in) :: nlev                           ! number of vertical levels
         658             :     integer, intent(in) :: nlatin                         ! number of input latitudes
         659             :     integer, intent(in) :: nlatout                        ! number of output latitudes
         660             :     integer, intent(in) :: nlondout                       ! longitude dimension of output grid
         661             :     integer, intent(in) :: nlonout(nlatout)               ! number of output longitudes per lat
         662             :     integer, intent(in) :: nlevdout                       ! vertical dimension of output grid
         663             :
         664             :     real(r8), intent(in) :: arrin(nlondin,nlevdin,nlatin) ! input array of values to interpolate
         665             :     real(r8), intent(in) :: xin(nlondin)                  ! input x mesh
         666             :     real(r8), intent(in) :: yin(nlatin)                   ! input y mesh
         667             :     real(r8), intent(in) :: xout(nlondout,nlatout)        ! output x mesh
         668             :     real(r8), intent(in) :: yout(nlatout)                 ! output y mesh
         669             :     !
         670             :     ! Output arguments
         671             :     !
         672             :     real(r8), intent(out) :: arrout(nlondout,nlevdout,nlatout) ! interpolated array
         673             :     !
         674             :     ! Local workspace
         675             :     !
         676             :     integer :: i, ii, iw, ie, iiprev ! longitude indices
         677             :     integer :: j, jj, js, jn, jjprev ! latitude indices
         678             :     integer :: k                     ! level index
         679             :     integer :: icount                ! number of bad values
         680             :
         681             :     real(r8) :: extrap               ! percent grid non-overlap
         682             :     real(r8) :: dxinwrap             ! delta-x on input grid for 2-pi
         683             :     real(r8) :: avgdxin              ! avg input delta-x
         684             :     real(r8) :: ratio                ! compare dxinwrap to avgdxin
         685             :     real(r8) :: sum                  ! sum of weights (used for testing)
         686             :     !
         687             :     ! Dynamic
         688             :     !
         689           0 :     integer :: iim(nlondout)         ! interpolation index to the left
         690           0 :     integer :: iip(nlondout)         ! interpolation index to the right
         691           0 :     integer :: jjm(nlatout)          ! interpolation index to the south
         692           0 :     integer :: jjp(nlatout)          ! interpolation index to the north
         693             :
         694           0 :     real(r8) :: wgts(nlatout)        ! interpolation weight to the north
         695           0 :     real(r8) :: wgtn(nlatout)        ! interpolation weight to the north
         696           0 :     real(r8) :: wgte(nlondout)       ! interpolation weight to the north
         697           0 :     real(r8) :: wgtw(nlondout)       ! interpolation weight to the north
         698           0 :     real(r8) :: igrid(nlonin)        ! interpolation weight to the north
         699             :     !
         700             :     ! Check validity of input coordinate arrays: must be monotonically increasing,
         701             :     ! and have a total of at least 2 elements
         702             :     !
         703           0 :     if (nlonin < 2 .or. nlatin < 2) then
         704           0 :        call endrun ('BILIN: Must have at least 2 input points for interpolation')
         705             :     end if
         706             :
         707           0 :     if (xin(1) < 0._r8 .or. xin(nlonin) > 360._r8) then
         708           0 :        call endrun ('BILIN: Input x-grid must be between 0 and 360')
         709             :     end if
         710             :
         711           0 :     icount = 0
         712           0 :     do i=1,nlonin-1
         713           0 :        if (xin(i) >= xin(i+1)) icount = icount + 1
         714             :     end do
         715             :
         716           0 :     do j=1,nlatin-1
         717           0 :        if (yin(j) >= yin(j+1)) icount = icount + 1
         718             :     end do
         719             :
         720           0 :     do j=1,nlatout-1
         721           0 :        if (yout(j) >= yout(j+1)) icount = icount + 1
         722             :     end do
         723             :
         724           0 :     do j=1,nlatout
         725           0 :        do i=1,nlonout(j)-1
         726           0 :           if (xout(i,j) >= xout(i+1,j)) icount = icount + 1
         727             :        end do
         728             :     end do
         729             :
         730           0 :     if (icount > 0) then
         731           0 :        call endrun ('BILIN: Non-monotonic coordinate array(s) found')
         732             :     end if
         733             :
         734           0 :     if (yout(nlatout) <= yin(1) .or. yout(1) >= yin(nlatin)) then
         735           0 :        call endrun ('BILIN: No overlap in y-coordinates')
         736             :     end if
         737             :
         738           0 :     do j=1,nlatout
         739           0 :        if (xout(1,j) < 0._r8 .or. xout(nlonout(j),j) > 360._r8) then
         740           0 :           call endrun ('BILIN: Output x-grid must be between 0 and 360')
         741             :        end if
         742             :
         743           0 :        if (xout(nlonout(j),j) <= xin(1) .or.  &
         744           0 :             xout(1,j)          >= xin(nlonin)) then
         745           0 :           call endrun ('BILIN: No overlap in x-coordinates')
         746             :        end if
         747             :     end do
         748             :     !
         749             :     ! Initialize index arrays for later checking
         750             :     !
         751           0 :     do j=1,nlatout
         752           0 :        jjm(j) = 0
         753           0 :        jjp(j) = 0
         754             :     end do
         755             :     !
         756             :     ! For values which extend beyond N and S boundaries, set weights
         757             :     ! such that values will just be copied.
         758             :     !
         759           0 :     do js=1,nlatout
         760           0 :        if (yout(js) > yin(1)) exit
         761           0 :        jjm(js) = 1
         762           0 :        jjp(js) = 1
         763           0 :        wgts(js) = 1._r8
         764           0 :        wgtn(js) = 0._r8
         765             :     end do
         766             :
         767           0 :     do jn=nlatout,1,-1
         768           0 :        if (yout(jn) <= yin(nlatin)) exit
         769           0 :        jjm(jn) = nlatin
         770           0 :        jjp(jn) = nlatin
         771           0 :        wgts(jn) = 1._r8
         772           0 :        wgtn(jn) = 0._r8
         773             :     end do
         774             :     !
         775             :     ! Loop though output indices finding input indices and weights
         776             :     !
         777           0 :     jjprev = 1
         778           0 :     do j=js,jn
         779           0 :        do jj=jjprev,nlatin-1
         780           0 :           if (yout(j) > yin(jj) .and. yout(j) <= yin(jj+1)) then
         781           0 :              jjm(j) = jj
         782           0 :              jjp(j) = jj + 1
         783           0 :              wgts(j) = (yin(jj+1) - yout(j)) / (yin(jj+1) - yin(jj))
         784           0 :              wgtn(j) = (yout(j)   - yin(jj)) / (yin(jj+1) - yin(jj))
         785           0 :              goto 30
         786             :           end if
         787             :        end do
         788           0 :        call endrun ('BILIN: Failed to find interp values')
         789           0 : 30     jjprev = jj
         790             :     end do
         791             :
         792           0 :     dxinwrap = xin(1) + 360._r8 - xin(nlonin)
         793             :     !
         794             :     ! Check for sane dxinwrap values.  Allow to differ no more than 10% from avg
         795             :     !
         796           0 :     avgdxin = (xin(nlonin)-xin(1))/(nlonin-1._r8)
         797           0 :     ratio = dxinwrap/avgdxin
         798           0 :     if (ratio < 0.9_r8 .or. ratio > 1.1_r8) then
         799           0 :        write(iulog,*)'BILIN: Insane dxinwrap value =',dxinwrap,' avg=', avgdxin
         800           0 :        call endrun
         801             :     end if
         802             :     !
         803             :     ! Check that interp/extrap points have been found for all outputs, and that
         804             :     ! they are reasonable (i.e. within 32-bit roundoff)
         805             :     !
         806           0 :     icount = 0
         807           0 :     do j=1,nlatout
         808           0 :        if (jjm(j) == 0 .or. jjp(j) == 0) icount = icount + 1
         809           0 :        sum = wgts(j) + wgtn(j)
         810           0 :        if (sum < 0.99999_r8 .or. sum > 1.00001_r8) icount = icount + 1
         811           0 :        if (wgts(j) < 0._r8 .or. wgts(j) > 1._r8) icount = icount + 1
         812           0 :        if (wgtn(j) < 0._r8 .or. wgtn(j) > 1._r8) icount = icount + 1
         813             :     end do
         814             :
         815           0 :     if (icount > 0) then
         816           0 :        call endrun ('BILIN: something bad in latitude indices or weights')
         817             :     end if
         818             :     !
         819             :     ! Do the bilinear interpolation
         820             :     !
         821           0 :     do j=1,nlatout
         822             :        !
         823             :        ! Initialize index arrays for later checking
         824             :        !
         825           0 :        do i=1,nlondout
         826           0 :           iim(i) = 0
         827           0 :           iip(i) = 0
         828             :        end do
         829             :        !
         830             :        ! For values which extend beyond E and W boundaries, set weights such that
         831             :        ! values will be interpolated between E and W edges of input grid.
         832             :        !
         833           0 :        do iw=1,nlonout(j)
         834           0 :           if (xout(iw,j) > xin(1)) exit
         835           0 :           iim(iw) = nlonin
         836           0 :           iip(iw) = 1
         837           0 :           wgtw(iw) = (xin(1)        - xout(iw,j))   /dxinwrap
         838           0 :           wgte(iw) = (xout(iw,j)+360._r8 - xin(nlonin))/dxinwrap
         839             :        end do
         840             :
         841           0 :        do ie=nlonout(j),1,-1
         842           0 :           if (xout(ie,j) <= xin(nlonin)) exit
         843           0 :           iim(ie) = nlonin
         844           0 :           iip(ie) = 1
         845           0 :           wgtw(ie) = (xin(1)+360._r8 - xout(ie,j))   /dxinwrap
         846           0 :           wgte(ie) = (xout(ie,j)    - xin(nlonin))/dxinwrap
         847             :        end do
         848             :        !
         849             :        ! Loop though output indices finding input indices and weights
         850             :        !
         851             :        iiprev = 1
         852           0 :        do i=iw,ie
         853           0 :           do ii=iiprev,nlonin-1
         854           0 :              if (xout(i,j) > xin(ii) .and. xout(i,j) <= xin(ii+1)) then
         855           0 :                 iim(i) = ii
         856           0 :                 iip(i) = ii + 1
         857           0 :                 wgtw(i) = (xin(ii+1) - xout(i,j)) / (xin(ii+1) - xin(ii))
         858           0 :                 wgte(i) = (xout(i,j) - xin(ii))   / (xin(ii+1) - xin(ii))
         859           0 :                 goto 60
         860             :              end if
         861             :           end do
         862           0 :           call endrun ('BILIN: Failed to find interp values')
         863           0 : 60        iiprev = ii
         864             :        end do
         865             :
         866           0 :        icount = 0
         867           0 :        do i=1,nlonout(j)
         868           0 :           if (iim(i) == 0 .or. iip(i) == 0) icount = icount + 1
         869           0 :           sum = wgtw(i) + wgte(i)
         870           0 :           if (sum < 0.99999_r8 .or. sum > 1.00001_r8) icount = icount + 1
         871           0 :           if (wgtw(i) < 0._r8 .or. wgtw(i) > 1._r8) icount = icount + 1
         872           0 :           if (wgte(i) < 0._r8 .or. wgte(i) > 1._r8) icount = icount + 1
         873             :        end do
         874             :
         875           0 :        if (icount > 0) then
         876           0 :           write(iulog,*)'BILIN: j=',j,' Something bad in longitude indices or weights'
         877           0 :           call endrun
         878             :        end if
         879             :        !
         880             :        ! Do the interpolation, 1st in longitude then latitude
         881             :        !
         882           0 :        do k=1,nlev
         883           0 :           do i=1,nlonin
         884           0 :              igrid(i) = arrin(i,k,jjm(j))*wgts(j) + arrin(i,k,jjp(j))*wgtn(j)
         885             :           end do
         886             :
         887           0 :           do i=1,nlonout(j)
         888           0 :              arrout(i,k,j) = igrid(iim(i))*wgtw(i) + igrid(iip(i))*wgte(i)
         889             :           end do
         890             :        end do
         891             :     end do
         892             :
         893             :
         894           0 :     return
         895             :   end subroutine bilin
         896             :
         897             : !=========================================================================================
         898             :
         899     5969088 : subroutine vertinterp(ncol, ncold, nlev, pin, pout, arrin, arrout, &
         900           0 :                       extrapolate, ln_interp, ps, phis, tbot)
         901             :
         902             :    ! Vertically interpolate input array to output pressure level.  The
         903             :    ! interpolation is linear in either pressure (the default) or ln pressure.
         904             :    !
         905             :    ! If above the top boundary then return top boundary value.
         906             :    !
         907             :    ! If below bottom boundary then use bottom boundary value, or optionally
         908             :    ! for T or Z use the extrapolation algorithm from ECMWF (which was taken
         909             :    ! from the NCL code base).
         910             :
         911             :    use physconst, only: rair, rga
         912             :
         913             :    !------------------------------Arguments--------------------------------
         914             :    integer,  intent(in)  :: ncol              ! number active columns
         915             :    integer,  intent(in)  :: ncold             ! column dimension
         916             :    integer,  intent(in)  :: nlev              ! vertical dimension
         917             :    real(r8), intent(in)  :: pin(ncold,nlev)   ! input pressure levels
         918             :    real(r8), intent(in)  :: pout              ! output pressure level
         919             :    real(r8), intent(in)  :: arrin(ncold,nlev) ! input  array
         920             :    real(r8), intent(out) :: arrout(ncold)     ! output array (interpolated)
         921             :
         922             :    character(len=1), optional, intent(in) :: extrapolate ! either 'T' or 'Z'
         923             :    logical,          optional, intent(in) :: ln_interp   ! set true to interolate
         924             :                                                          ! in ln(pressure)
         925             :    real(r8),         optional, intent(in) :: ps(ncold)   ! surface pressure
         926             :    real(r8),         optional, intent(in) :: phis(ncold) ! surface geopotential
         927             :    real(r8),         optional, intent(in) :: tbot(ncold) ! temperature at bottom level
         928             :
         929             :    !---------------------------Local variables-----------------------------
         930             :    real(r8) :: alpha
         931             :    logical  :: linear_interp
         932             :    logical  :: do_extrapolate_T, do_extrapolate_Z
         933             :
         934             :    integer  :: i,k              ! indices
         935     5969088 :    integer  :: kupper(ncold)    ! Level indices for interpolation
         936             :    real(r8) :: dpu              ! upper level pressure difference
         937             :    real(r8) :: dpl              ! lower level pressure difference
         938     5969088 :    logical  :: found(ncold)     ! true if input levels found
         939             :    logical  :: error            ! error flag
         940             :    !----------------------------------------------------------------------------
         941             :
         942     2984544 :    alpha = 0.0065_r8*rair*rga
         943             :
         944     2984544 :    do_extrapolate_T = .false.
         945     2984544 :    do_extrapolate_Z = .false.
         946     2984544 :    if (present(extrapolate)) then
         947             :
         948           0 :       if (extrapolate == 'T') do_extrapolate_T = .true.
         949           0 :       if (extrapolate == 'Z') do_extrapolate_Z = .true.
         950             :
         951           0 :       if (.not. do_extrapolate_T .and. .not. do_extrapolate_Z) then
         952           0 :          call endrun ('VERTINTERP: extrapolate must be T or Z')
         953             :       end if
         954           0 :       if (.not. present(ps)) then
         955           0 :          call endrun ('VERTINTERP: ps required for extrapolation')
         956             :       end if
         957           0 :       if (.not. present(phis)) then
         958           0 :          call endrun ('VERTINTERP: phis required for extrapolation')
         959             :       end if
         960           0 :       if (do_extrapolate_Z) then
         961           0 :          if (.not. present(tbot)) then
         962           0 :             call endrun ('VERTINTERP: extrapolate must be T or Z')
         963             :          end if
         964             :       end if
         965             :    end if
         966             :
         967     2984544 :    linear_interp = .true.
         968     2984544 :    if (present(ln_interp)) then
         969           0 :       if (ln_interp) linear_interp = .false.
         970             :    end if
         971             :
         972    49834944 :    do i = 1, ncol
         973    46850400 :       found(i)  = .false.
         974    49834944 :       kupper(i) = 1
         975             :    end do
         976     2984544 :    error = .false.
         977             :
         978             :    ! Find indices of upper pressure bound
         979   280547136 :    do k = 1, nlev - 1
         980  4637634336 :       do i = 1, ncol
         981  4634649792 :          if ((.not. found(i)) .and. pin(i,k)<pout .and. pout<=pin(i,k+1)) then
         982    46850400 :             found(i)  = .true.
         983    46850400 :             kupper(i) = k
         984             :          end if
         985             :       end do
         986             :    end do
         987             :
         988    49834944 :    do i = 1, ncol
         989             :
         990    49834944 :       if (pout <= pin(i,1)) then
         991             :
         992             :          ! if above top pressure level then use value at top (no interpolation)
         993           0 :          arrout(i) = arrin(i,1)
         994             :
         995    46850400 :       else if (pout >= pin(i,nlev)) then
         996             :
         997           0 :          if (do_extrapolate_T) then
         998             :             ! use ECMWF algorithm to extrapolate T
         999           0 :             arrout(i) = extrapolate_T()
        1000             :
        1001           0 :          else if (do_extrapolate_Z) then
        1002             :             ! use ECMWF algorithm to extrapolate Z
        1003           0 :             arrout(i) = extrapolate_Z()
        1004             :
        1005             :          else
        1006             :             ! use bottom value
        1007           0 :             arrout(i) = arrin(i,nlev)
        1008             :          end if
        1009             :
        1010    46850400 :       else if (found(i)) then
        1011    46850400 :          if (linear_interp) then
        1012             :             ! linear interpolation in p
        1013    46850400 :             dpu = pout - pin(i,kupper(i))
        1014    46850400 :             dpl = pin(i,kupper(i)+1) - pout
        1015    46850400 :             arrout(i) = (arrin(i,kupper(i))*dpl + arrin(i,kupper(i)+1)*dpu)/(dpl + dpu)
        1016             :          else
        1017             :             ! linear interpolation in ln(p)
        1018           0 :             arrout(i) = arrin(i,kupper(i)) + (arrin(i,kupper(i)+1) - arrin(i,kupper(i)))* &
        1019           0 :                                              log(pout/pin(i,kupper(i))) /                 &
        1020           0 :                                              log(pin(i,kupper(i)+1)/pin(i,kupper(i)))
        1021             :          end if
        1022             :       else
        1023             :          error = .true.
        1024             :       end if
        1025             :    end do
        1026             :
        1027             :    ! Error check
        1028     5969088 :    if (error) then
        1029           0 :       call endrun ('VERTINTERP: ERROR FLAG')
        1030             :    end if
        1031             :
        1032             : contains
        1033             :
        1034             :    !======================================================================================
        1035             :
        1036           0 :    real(r8) function extrapolate_T()
        1037             :
        1038             :       ! local variables
        1039             :       real(r8) :: sixth
        1040             :       real(r8) :: tstar
        1041             :       real(r8) :: hgt
        1042             :       real(r8) :: alnp
        1043             :       real(r8) :: t0
        1044             :       real(r8) :: tplat
        1045             :       real(r8) :: tprime0
        1046             :       !-------------------------------------------------------------------------
        1047             :
        1048           0 :       sixth  = 1._r8 / 6._r8
        1049           0 :       tstar = arrin(i,nlev)*(1._r8 + alpha*(ps(i)/pin(i,nlev) - 1._r8))
        1050           0 :       hgt   = phis(i)*rga
        1051             :
        1052           0 :       if (hgt < 2000._r8) then
        1053           0 :          alnp = alpha*log(pout/ps(i))
        1054             :       else
        1055           0 :          t0 = tstar + 0.0065_r8*hgt
        1056           0 :          tplat = min(t0, 298._r8)
        1057             :
        1058           0 :          if (hgt <= 2500._r8) then
        1059             :             tprime0 = 0.002_r8*((2500._r8 - hgt)*t0 + &
        1060           0 :                                 (hgt - 2000._r8)*tplat)
        1061             :          else
        1062             :             tprime0 = tplat
        1063             :          end if
        1064             :
        1065           0 :          if (tprime0 < tstar) then
        1066             :             alnp = 0._r8
        1067             :          else
        1068           0 :             alnp = rair*(tprime0 - tstar)/phis(i) * log(pout/ps(i))
        1069             :          end if
        1070             :
        1071             :       end if
        1072             :
        1073           0 :       extrapolate_T = tstar*(1._r8 + alnp + 0.5_r8*alnp**2 + sixth*alnp**3)
        1074             :
        1075     2984544 :    end function extrapolate_T
        1076             :
        1077             :    !======================================================================================
        1078             :
        1079           0 :    real(r8) function extrapolate_Z()
        1080             :
        1081             :       ! local variables
        1082             :       real(r8) :: sixth
        1083             :       real(r8) :: hgt
        1084             :       real(r8) :: tstar
        1085             :       real(r8) :: t0
        1086             :       real(r8) :: alph
        1087             :       real(r8) :: alnp
        1088             :       !-------------------------------------------------------------------------
        1089             :
        1090           0 :       sixth  = 1._r8 / 6._r8
        1091           0 :       hgt   = phis(i)*rga
        1092           0 :       tstar = tbot(i)*(1._r8 + alpha*(ps(i)/pin(i,nlev) - 1._r8))
        1093           0 :       t0 = tstar + 0.0065_r8*hgt
        1094             :
        1095           0 :       if (tstar <= 290.5_r8 .and. t0 > 290.5_r8) then
        1096           0 :          alph = rair/phis(i) * (290.5_r8 - tstar)
        1097             :
        1098           0 :       else if (tstar > 290.5_r8 .and. t0 > 290.5_r8) then
        1099           0 :          alph  = 0._r8
        1100           0 :          tstar = 0.5_r8*(290.5_r8 + tstar)
        1101             :
        1102             :       else
        1103             :          alph = alpha
        1104             :
        1105             :       end if
        1106             :
        1107           0 :       if (tstar < 255._r8) then
        1108           0 :          tstar = 0.5_r8*(tstar + 255._r8)
        1109             :       end if
        1110           0 :       alnp = alph*log(pout/ps(i))
        1111             :
        1112             :       extrapolate_Z = hgt - rair*tstar*rga*log(pout/ps(i))* &
        1113           0 :                       (1._r8 + 0.5_r8*alnp + sixth*alnp**2)
        1114             :
        1115           0 :    end function extrapolate_Z
        1116             :
        1117             : end subroutine vertinterp
        1118             :
        1119             : !=========================================================================================
        1120             :
        1121      741888 :   subroutine get_timeinterp_factors (cycflag, np1, cdayminus, cdayplus, cday, &
        1122             :        fact1, fact2, str)
        1123             :     !---------------------------------------------------------------------------
        1124             :     !
        1125             :     ! Purpose: Determine time interpolation factors (normally for a boundary dataset)
        1126             :     !          for linear interpolation.
        1127             :     !
        1128             :     ! Method:  Assume 365 days per year.  Output variable fact1 will be the weight to
        1129             :     !          apply to data at calendar time "cdayminus", and fact2 the weight to apply
        1130             :     !          to data at time "cdayplus".  Combining these values will produce a result
        1131             :     !          valid at time "cday".  Output arguments fact1 and fact2 will be between
        1132             :     !          0 and 1, and fact1 + fact2 = 1 to roundoff.
        1133             :     !
        1134             :     ! Author:  Jim Rosinski
        1135             :     !
        1136             :     !---------------------------------------------------------------------------
        1137             :     implicit none
        1138             :     !
        1139             :     ! Arguments
        1140             :     !
        1141             :     logical, intent(in) :: cycflag             ! flag indicates whether dataset is being cycled yearly
        1142             :
        1143             :     integer, intent(in) :: np1                 ! index points to forward time slice matching cdayplus
        1144             :
        1145             :     real(r8), intent(in) :: cdayminus          ! calendar day of rearward time slice
        1146             :     real(r8), intent(in) :: cdayplus           ! calendar day of forward time slice
        1147             :     real(r8), intent(in) :: cday               ! calenar day to be interpolated to
        1148             :     real(r8), intent(out) :: fact1             ! time interpolation factor to apply to rearward time slice
        1149             :     real(r8), intent(out) :: fact2             ! time interpolation factor to apply to forward time slice
        1150             :
        1151             :     character(len=*), intent(in) :: str        ! string to be added to print in case of error (normally the callers name)
        1152             :     !
        1153             :     ! Local workspace
        1154             :     !
        1155             :     real(r8) :: deltat                         ! time difference (days) between cdayminus and cdayplus
        1156             :     real(r8), parameter :: daysperyear = 365._r8  ! number of days in a year
        1157             :     !
        1158             :     ! Initial sanity checks
        1159             :     !
        1160      741888 :     if (np1 == 1 .and. .not. cycflag) then
        1161           0 :        call endrun ('GETFACTORS:'//str//' cycflag false and forward month index = Jan. not allowed')
        1162             :     end if
        1163             :
        1164      741888 :     if (np1 < 1) then
        1165           0 :        call endrun ('GETFACTORS:'//str//' input arg np1 must be > 0')
        1166             :     end if
        1167             :
        1168      741888 :     if (cycflag) then
        1169      741888 :        if ((cday < 1._r8) .or. (cday > (daysperyear+1._r8))) then
        1170           0 :           write(iulog,*) 'GETFACTORS:', str, ' bad cday=',cday
        1171           0 :           call endrun ()
        1172             :        end if
        1173             :     else
        1174           0 :        if (cday < 1._r8) then
        1175           0 :           write(iulog,*) 'GETFACTORS:', str, ' bad cday=',cday
        1176           0 :           call endrun ()
        1177             :        end if
        1178             :     end if
        1179             :     !
        1180             :     ! Determine time interpolation factors.  Account for December-January
        1181             :     ! interpolation if dataset is being cycled yearly.
        1182             :     !
        1183      741888 :     if (cycflag .and. np1 == 1) then                     ! Dec-Jan interpolation
        1184      741888 :        deltat = cdayplus + daysperyear - cdayminus
        1185      741888 :        if (cday > cdayplus) then                         ! We are in December
        1186           0 :           fact1 = (cdayplus + daysperyear - cday)/deltat
        1187           0 :           fact2 = (cday - cdayminus)/deltat
        1188             :        else                                              ! We are in January
        1189      741888 :           fact1 = (cdayplus - cday)/deltat
        1190      741888 :           fact2 = (cday + daysperyear - cdayminus)/deltat
        1191             :        end if
        1192             :     else
        1193           0 :        deltat = cdayplus - cdayminus
        1194           0 :        fact1 = (cdayplus - cday)/deltat
        1195           0 :        fact2 = (cday - cdayminus)/deltat
        1196             :     end if
        1197             :
        1198      741888 :     if (.not. valid_timeinterp_factors (fact1, fact2)) then
        1199           0 :        write(iulog,*) 'GETFACTORS: ', str, ' bad fact1 and/or fact2=', fact1, fact2
        1200           0 :        call endrun ()
        1201             :     end if
        1202             :
        1203      741888 :     return
        1204             :   end subroutine get_timeinterp_factors
        1205             :
        1206      741888 :   logical function valid_timeinterp_factors (fact1, fact2)
        1207             :     !---------------------------------------------------------------------------
        1208             :     !
        1209             :     ! Purpose: check sanity of time interpolation factors to within 32-bit roundoff
        1210             :     !
        1211             :     !---------------------------------------------------------------------------
        1212             :     implicit none
        1213             :
        1214             :     real(r8), intent(in) :: fact1, fact2           ! time interpolation factors
        1215             :
        1216      741888 :     valid_timeinterp_factors = .true.
        1217             :
        1218             :     ! The fact1 .ne. fact1 and fact2 .ne. fact2 comparisons are to detect NaNs.
        1219             :     if (abs(fact1+fact2-1._r8) > 1.e-6_r8 .or. &
        1220             :          fact1 > 1.000001_r8 .or. fact1 < -1.e-6_r8 .or. &
        1221             :          fact2 > 1.000001_r8 .or. fact2 < -1.e-6_r8 .or. &
        1222      741888 :          fact1 .ne. fact1 .or. fact2 .ne. fact2) then
        1223             :
        1224           0 :        valid_timeinterp_factors = .false.
        1225             :     end if
        1226             :
        1227             :     return
        1228             :   end function valid_timeinterp_factors
        1229             :
        1230           0 : end module interpolate_data

Generated by: LCOV version 1.14