LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - interpolation.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 16 | 118 | 13.6 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 1 | 8 | 12.5 %

              Line data    Source code
           1             : !-------------------------------------------------------------------------------
           2             : !$Id$
           3             : !===============================================================================
           4             : module interpolation
           5             :
           6             :   use clubb_precision, only: &
           7             :     core_rknd ! Variable(s)
           8             :
           9             :   implicit none
          10             :
          11             :   private ! Default Scope
          12             :
          13             :   public :: lin_interpolate_two_points, binary_search, zlinterp_fnc, &
          14             :     lin_interpolate_on_grid, linear_interp_factor, mono_cubic_interp, plinterp_fnc, &
          15             :     pvertinterp
          16             :
          17             :   contains
          18             :
          19             : !-------------------------------------------------------------------------------
          20           0 :   function lin_interpolate_two_points( height_int, height_high, height_low, &
          21             :     var_high, var_low )
          22             :
          23             : ! Description:
          24             : ! This function computes a linear interpolation of the value of variable.
          25             : ! Given two known values of a variable at two height values, the value
          26             : ! of that variable at a height between those two height levels (rather
          27             : ! than a height outside of those two height levels) is computed.
          28             : !
          29             : ! Here is a diagram:
          30             : !
          31             : !  ################################ Height high, know variable value
          32             : !
          33             : !
          34             : !
          35             : !  -------------------------------- Height to be interpolated to; linear interpolation
          36             : !
          37             : !
          38             : !
          39             : !
          40             : !
          41             : !  ################################ Height low, know variable value
          42             : !
          43             : !
          44             : ! FORMULA:
          45             : !
          46             : ! variable(@ Height interpolation) =
          47             : !
          48             : ! [ (variable(@ Height high) - variable(@ Height low)) / (Height high - Height low) ]
          49             : ! * (Height interpolation - Height low)  +  variable(@ Height low)
          50             :
          51             : ! Comments from WRF-HOC, Brian Griffin.
          52             :
          53             : ! References:
          54             : ! None
          55             : !-------------------------------------------------------------------------------
          56             :
          57             :     use clubb_precision, only: &
          58             :         core_rknd ! Variable(s)
          59             :
          60             :     use constants_clubb, only: fstderr ! Constant
          61             :
          62             :     implicit none
          63             :
          64             :     ! Input Variables
          65             :
          66             :     real( kind = core_rknd ), intent(in) :: &
          67             :       height_int,  & ! Height to be interpolated to     [m]
          68             :       height_high, & ! Height above the interpolation   [m]
          69             :       height_low,  & ! Height below the interpolation   [m]
          70             :       var_high,    & ! Variable above the interpolation [units vary]
          71             :       var_low        ! Variable below the interpolation [units vary]
          72             :
          73             :     ! Output Variables
          74             :     real( kind = core_rknd ) :: lin_interpolate_two_points
          75             :
          76             :     ! Check for valid input
          77           0 :     if ( abs(height_low - height_high) < 1.0e-12_core_rknd ) then
          78           0 :       write(fstderr,*) "lin_interpolate_two_points: height_high and height_low cannot be equal."
          79           0 :       error stop
          80             :     end if
          81             :
          82             :     ! Compute linear interpolation
          83             :
          84             :     lin_interpolate_two_points = ( ( height_int - height_low )/( height_high - height_low ) ) &
          85           0 :       * ( var_high - var_low ) + var_low
          86             :
          87             :     return
          88             :   end function lin_interpolate_two_points
          89             :
          90             :   !-------------------------------------------------------------------------------------------------
          91           0 :   elemental real( kind = core_rknd ) function linear_interp_factor( factor, var_high, var_low )
          92             :   ! Description:
          93             :   !   Determines the coefficient for a linear interpolation
          94             :   !
          95             :   ! References:
          96             :   !   None
          97             :   !-------------------------------------------------------------------------------------------------
          98             :
          99             :     use clubb_precision, only: &
         100             :         core_rknd ! Variable(s)
         101             :
         102             :     implicit none
         103             :
         104             :     real( kind = core_rknd ), intent(in) :: &
         105             :       factor,   & ! Factor                           [units vary]
         106             :       var_high, & ! Variable above the interpolation [units vary]
         107             :       var_low     ! Variable below the interpolation [units vary]
         108             :
         109           0 :     linear_interp_factor = factor * ( var_high - var_low ) + var_low
         110             :
         111             :     return
         112             :   end function linear_interp_factor
         113             :   !-------------------------------------------------------------------------------------------------
         114           0 :   function mono_cubic_interp &
         115             :     ( z_in, km1, k00, kp1, kp2, zm1, z00, zp1, zp2, fm1, f00, fp1, fp2 ) result ( f_out )
         116             :
         117             :   ! Description:
         118             :   !   Steffen's monotone cubic interpolation method
         119             :   !   Returns monotone cubic interpolated value between x00 and xp1
         120             :
         121             :   ! Original Author:
         122             :   !   Takanobu Yamaguchi
         123             :   !   tak.yamaguchi@noaa.gov
         124             :   !
         125             :   !   This version has been modified slightly for CLUBB's coding standards and
         126             :   !   adds the 3/2 from eqn 21. -dschanen 26 Oct 2011
         127             :   !   We have also added a quintic polynomial option.
         128             :   !
         129             :   ! References:
         130             :   !   M. Steffen, Astron. Astrophys. 239, 443-450 (1990)
         131             :   !-------------------------------------------------------------------------------------------------
         132             :
         133             :     use constants_clubb, only: &
         134             :         eps ! Constant
         135             :
         136             :     use clubb_precision, only: &
         137             :         core_rknd ! Constant
         138             :
         139             :     use model_flags, only: &
         140             :         l_quintic_poly_interp ! Variable(s)
         141             :
         142             :     implicit none
         143             :
         144             :     ! External
         145             :     intrinsic :: sign, abs, min
         146             :
         147             :     ! Input Variables
         148             :     real( kind = core_rknd ), intent(in) :: &
         149             :       z_in ! The altitude to be interpolated to [m]
         150             :
         151             :     ! k-levels;  their meaning depends on whether we're extrapolating or interpolating
         152             :     integer, intent(in) :: &
         153             :       km1, k00, kp1, kp2
         154             :
         155             :     real( kind = core_rknd ), intent(in) :: &
         156             :       zm1, z00, zp1, zp2, & ! The altitudes for km1, k00, kp1, kp2      [m]
         157             :       fm1, f00, fp1, fp2    ! The field at km1, k00, kp1, and kp2       [units vary]
         158             :
         159             :     ! Output Variables
         160             :     real( kind = core_rknd ) :: f_out ! The interpolated field
         161             :
         162             :     ! Local Variables
         163             :     real( kind = core_rknd ) :: &
         164             :       hm1, h00, hp1, &
         165             :       sm1, s00, sp1, &
         166             :       p00, pp1, &
         167             :       dfdx00, dfdxp1, &
         168             :       c1, c2, c3, c4, &
         169             :       w00, wp1, &
         170             :       coef1, coef2, &
         171             :       zprime, beta, alpha, zn
         172             :
         173             :     ! ---- Begin Code ----
         174             :
         175           0 :     coef1 = 1.0_core_rknd
         176           0 :     coef2 = 1.0_core_rknd
         177             :
         178           0 :     if ( km1 <= k00 ) then
         179           0 :       hm1 = z00 - zm1
         180           0 :       h00 = zp1 - z00
         181           0 :       hp1 = zp2 - zp1
         182             :
         183           0 :       if ( km1 == k00 ) then
         184           0 :         s00 = ( fp1 - f00 ) / ( zp1 - z00 )
         185           0 :         sp1 = ( fp2 - fp1 ) / ( zp2 - zp1 )
         186           0 :         dfdx00 = s00
         187           0 :         pp1 = ( s00 * hp1 + sp1 * h00 ) / ( h00 + hp1 )
         188             :         dfdxp1 = coef1*( sign( 1.0_core_rknd, s00 ) + sign( 1.0_core_rknd, sp1 ) ) &
         189           0 :           * min( abs( s00 ), abs( sp1 ), coef2*0.5_core_rknd*abs( pp1 ) )
         190             :
         191           0 :       else if ( kp1 == kp2 ) then
         192           0 :         sm1 = ( f00 - fm1 ) / ( z00 - zm1 )
         193           0 :         s00 = ( fp1 - f00 ) / ( zp1 - z00 )
         194           0 :         p00 = ( sm1 * h00 + s00 * hm1 ) / ( hm1 + h00 )
         195             :         dfdx00 = coef1*( sign( 1.0_core_rknd, sm1 ) + sign( 1.0_core_rknd, s00 ) ) &
         196           0 :           * min( abs( sm1 ), abs( s00 ), coef2*0.5_core_rknd*abs( p00 ) )
         197           0 :         dfdxp1 = s00
         198             :
         199             :       else
         200           0 :         sm1 = ( f00 - fm1 ) / ( z00 - zm1 )
         201           0 :         s00 = ( fp1 - f00 ) / ( zp1 - z00 )
         202           0 :         sp1 = ( fp2 - fp1 ) / ( zp2 - zp1 )
         203           0 :         p00 = ( sm1 * h00 + s00 * hm1 ) / ( hm1 + h00 )
         204           0 :         pp1 = ( s00 * hp1 + sp1 * h00 ) / ( h00 + hp1 )
         205             :         dfdx00 = coef1*( sign( 1.0_core_rknd, sm1 ) + sign( 1.0_core_rknd, s00 ) ) &
         206           0 :           * min( abs( sm1 ), abs( s00 ), coef2*0.5_core_rknd*abs( p00 ) )
         207             :         dfdxp1 = coef1*( sign( 1.0_core_rknd, s00 ) + sign( 1.0_core_rknd, sp1 ) ) &
         208           0 :           * min( abs( s00 ), abs( sp1 ), coef2*0.5_core_rknd*abs( pp1 ) )
         209             :
         210             :       end if
         211             :
         212           0 :       c1 = ( dfdx00 + dfdxp1 - 2._core_rknd * s00 ) / ( h00 ** 2 )
         213           0 :       c2 = ( 3._core_rknd * s00 - 2._core_rknd * dfdx00 - dfdxp1 ) / h00
         214           0 :       c3 = dfdx00
         215           0 :       c4 = f00
         216             :
         217             :       if ( .not. l_quintic_poly_interp ) then
         218             :
         219             :         ! Old formula
         220             :         !f_out = c1 * ( (z_in - z00)**3 ) + c2 * ( (z_in - z00)**2 ) + c3 * (z_in - z00) + c4
         221             :
         222             :         ! Faster nested multiplication
         223           0 :         zprime = z_in - z00
         224           0 :         f_out =  c4 + zprime*( c3 + zprime*( c2 + ( zprime*c1 ) ) )
         225             :
         226             :       else
         227             :
         228             :        ! Use a quintic polynomial interpolation instead instead of the Steffen formula.
         229             :        ! Unlike the formula above, this formula does not guarantee monotonicity.
         230             :
         231             :         beta = 120._core_rknd * ( (fp1-f00) - 0.5_core_rknd * h00 * (dfdx00 + dfdxp1) )
         232             :
         233             :         ! Prevent an underflow by using a linear interpolation
         234             :         if ( abs( beta ) < eps ) then
         235             :           f_out = lin_interpolate_two_points( z00, zp1, zm1, &
         236             :                            fp1, fm1 )
         237             :
         238             :         else
         239             :           alpha = (6._core_rknd/beta) * h00 * (dfdxp1-dfdx00) + 0.5_core_rknd
         240             :           zn = (z_in-z00)/h00
         241             :
         242             :           f_out = ( &
         243             :                   (( (beta/20._core_rknd)*zn - (beta*(1._core_rknd+alpha) &
         244             :                     / 12._core_rknd)) * zn + (beta*alpha/6._core_rknd)) &
         245             :                     * zn**2 + dfdx00*h00 &
         246             :                   ) * zn + f00
         247             :         end if ! beta < eps
         248             :       end if ! ~quintic_polynomial
         249             :
         250             :     else
         251             :       ! Linear extrapolation
         252           0 :       wp1 = ( z_in - z00 ) / ( zp1 - z00 )
         253           0 :       w00 = 1._core_rknd - wp1
         254           0 :       f_out = wp1 * fp1 + w00 * f00
         255             :
         256             :     end if
         257             :
         258             :     return
         259             :   end function mono_cubic_interp
         260             :
         261             : !-------------------------------------------------------------------------------
         262           0 :   integer function binary_search( n, array, var ) &
         263             :     result( i )
         264             :
         265             :     ! Description:
         266             :     ! This subroutine performs a binary search to find the closest value greater
         267             :     ! than or equal to var in the array.  This function returns the index of the
         268             :     ! closest value of array that is greater than or equal to var.  It returns a
         269             :     ! value of -1 if var is outside the bounds of array.
         270             :     !
         271             :     !-----------------------------------------------------------------------
         272             :
         273             :     use clubb_precision, only: &
         274             :         core_rknd                  ! Variable(s)
         275             :
         276             :     use constants_clubb, only: &
         277             :         fstderr                    ! Variable(s)
         278             :
         279             :     use error_code, only: &
         280             :         clubb_at_least_debug_level ! Error indicator
         281             :
         282             :     implicit none
         283             :
         284             :     ! Input Variables
         285             :
         286             :     ! Size of the array
         287             :     integer, intent(in) :: n
         288             :
         289             :     ! The array being searched (must be sorted from least value to greatest
         290             :     ! value).
         291             :     real( kind = core_rknd ), dimension(n), intent(in) :: array
         292             :
         293             :     ! The value being searched for
         294             :     real( kind = core_rknd ), intent(in) :: var
         295             :
         296             :     ! Bounds of the search
         297             :     integer :: high
         298             :     integer :: low
         299             :
         300             :     ! The initial value of low has been changed from 1 to 2 due to a problem
         301             :     ! that was occuring when var was close to the lower bound.
         302             :     !
         303             :     ! The lowest value in the array (which is sorted by increasing values) is
         304             :     ! found at index 1, while the highest value in the array is found at
         305             :     ! index n.  Unless the value of var exactly corresponds with one of the
         306             :     ! values found in the array, or unless the value of var is found outside of
         307             :     ! the array, the value of var will be found between two levels of the array.
         308             :     ! In this scenario, the output of function binary_search is the index of the
         309             :     ! HIGHER level.  For example, if the value of var is found between array(1)
         310             :     ! and array(2), the output of function binary_search will be 2.
         311             :     !
         312             :     ! Therefore, the lowest index of a HIGHER level in an interpolation is 2.
         313             :     ! Thus, the initial value of low has been changed to 2.  This will prevent
         314             :     ! the value of variable "i" below from becoming 1.  If the value of "i"
         315             :     ! becomes 1, the code below tries to access array(0) (which is array(i-1)
         316             :     ! when i = 1) and produces an error.
         317             :
         318           0 :     low = 2
         319             :
         320           0 :     high = n
         321             :
         322             :     ! Ensure variable is within range of array and array has more than 1 element
         323           0 :     if ( var < array(1) .or. var > array(n) .or. n < 2 ) then
         324           0 :         i = -1
         325             :         return
         326             :     end if
         327             :
         328             :     ! Special case to cover var == array(1)
         329           0 :     if ( var >= array(1) .and. var <= array(2) ) then
         330           0 :         i = 2
         331             :         return
         332             :     end if
         333             :
         334             :
         335           0 :     do while( low <= high )
         336             :
         337           0 :       i = (low + high) / 2
         338             :
         339           0 :       if ( var > array(i-1) .and. var <= array(i) ) then
         340             :
         341             :         return
         342             :
         343           0 :       elseif ( var < array(i) ) then
         344             :
         345             :         high = i - 1
         346             :
         347           0 :       elseif ( var > array(i) ) then
         348             :
         349           0 :         low = i + 1
         350             :
         351             :       endif
         352             :
         353             :     enddo
         354             :
         355             :     ! Code should not get to this point, but return -1 to be safe
         356           0 :     if ( clubb_at_least_debug_level( 0 ) ) then
         357           0 :         write(fstderr,*) "Logic error in function binary_search."
         358             :     end if
         359             :
         360           0 :     i = -1
         361             :     return
         362             :
         363             :   end function binary_search
         364             :
         365             : !-------------------------------------------------------------------------------
         366           0 :   function plinterp_fnc( dim_out, dim_src, grid_out,  &
         367           0 :                        grid_src, var_src )  &
         368           0 :   result( var_out )
         369             : ! Description:
         370             : !   Do a linear interpolation in the vertical with pressures.  Assumes
         371             : !   values that are less than lowest source point are zero and above the
         372             : !   highest source point are zero. Also assumes altitude increases linearly.
         373             : !   This function just calls zlinterp_fnc, but negates grid_out and grid_src.
         374             :
         375             : ! References:
         376             : !   function LIN_INT from WRF-HOC
         377             : !-----------------------------------------------------------------------
         378             :
         379             :     use clubb_precision, only: &
         380             :         core_rknd ! Variable(s)
         381             :
         382             :     implicit none
         383             :
         384             :     ! Input variables
         385             :     integer, intent(in) :: dim_out, dim_src
         386             :
         387             :     real( kind = core_rknd ), dimension(dim_src), intent(in) ::  &
         388             :       grid_src,  & ! [m]
         389             :       var_src      ! [units vary]
         390             :
         391             :     real( kind = core_rknd ), dimension(dim_out), intent(in) :: &
         392             :       grid_out ! [m]
         393             :
         394             :     ! Output variable
         395             :     real( kind = core_rknd ), dimension(dim_out) :: &
         396             :       var_out ! [units vary]
         397             :
         398             :     ! ---- Begin Code ----
         399             :
         400             :     var_out = zlinterp_fnc( dim_out, dim_src, -grid_out, &
         401           0 :                             -grid_src, var_src )
         402             :
         403           0 :     return
         404           0 :   end function plinterp_fnc
         405             : !-------------------------------------------------------------------------------
         406           0 :   function zlinterp_fnc( dim_out, dim_src, grid_out,  &
         407           0 :                        grid_src, var_src )  &
         408           0 :   result( var_out )
         409             : ! Description:
         410             : !   Do a linear interpolation in the vertical.  Assumes values that
         411             : !   are less than lowest source point are zero and above the highest
         412             : !   source point are zero. Also assumes altitude increases linearly.
         413             :
         414             : ! References:
         415             : !   function LIN_INT from WRF-HOC
         416             : !-----------------------------------------------------------------------
         417             :
         418             :     use clubb_precision, only: &
         419             :         core_rknd ! Variable(s)
         420             :
         421             :     implicit none
         422             :
         423             :     ! Input variables
         424             :     integer, intent(in) :: dim_out, dim_src
         425             :
         426             :     real( kind = core_rknd ), dimension(dim_src), intent(in) ::  &
         427             :       grid_src,  & ! [m]
         428             :       var_src      ! [units vary]
         429             :
         430             :     real( kind = core_rknd ), dimension(dim_out), intent(in) :: &
         431             :       grid_out ! [m]
         432             :
         433             :     ! Output variable
         434             :     real( kind = core_rknd ), dimension(dim_out) :: &
         435             :       var_out ! [units vary]
         436             :
         437             :     ! Local variables
         438             :     integer :: k, kint, km1
         439             :
         440             : !   integer :: tst, kp1
         441             :
         442             :     ! ---- Begin Code ----
         443             :
         444           0 :     k = 1
         445             :
         446           0 :     do kint = 1, dim_out, 1
         447             :
         448             :       ! Set to 0 if we're below the input data's lowest point
         449           0 :       if ( grid_out(kint) < grid_src(1) ) then
         450           0 :         var_out(kint) = 0.0_core_rknd
         451           0 :         cycle
         452             :       end if
         453             :
         454             :       ! Increment k until the level is correct
         455             : !          do while ( grid_out(kint) > grid_src(k)
         456             : !     .                .and. k < dim_src )
         457             : !            k = k + 1
         458             : !          end do
         459             :
         460             :       ! Changed so a binary search is used instead of a sequential search
         461             : !          tst = binary_search(dim_src, grid_src, grid_out(kint))
         462           0 :       k = binary_search(dim_src, grid_src, grid_out(kint))
         463             :       ! Joshua Fasching April 2008
         464             :
         465             : !          print *, "k = ", k
         466             : !          print *, "tst = ", tst
         467             : !          print *, "dim_src = ", dim_src
         468             : !          print *,"------------------------------"
         469             :
         470             :       ! If the increment leads to a level above the data, set this
         471             :       ! point and all those above it to zero
         472             :       !if( k > dim_src ) then
         473           0 :       if ( k == -1 ) then
         474           0 :         var_out(kint:dim_out) = 0.0_core_rknd
         475             :         exit
         476             :       end if
         477             :
         478           0 :       km1 = max( 1, k-1 )
         479             :       !kp1 = min( k+1, dim_src )
         480             :
         481             :       ! Interpolate
         482           0 :       var_out(kint) = lin_interpolate_two_points( grid_out(kint), grid_src(k),  &
         483           0 :         grid_src(km1), var_src(k), var_src(km1) )
         484             :
         485             : !          ( var_src(k) - var_src(km1) ) / &
         486             : !          ( grid_src(k) - grid_src(km1) ) &
         487             : !            * ( grid_out(kint) - grid_src(km1) ) + var_src(km1) &
         488             : !            Changed to use a standard function for interpolation
         489             :
         490             :       !! Note this ends up changing the results slightly because
         491             :       !the placement of variables has been changed.
         492             :
         493             : !            Joshua Fasching April 2008
         494             :
         495             :     end do ! kint = 1..dim_out
         496             :
         497           0 :     return
         498           0 :   end function zlinterp_fnc
         499             :
         500             :   !-------------------------------------------------------------------------------
         501    17870112 :   subroutine pvertinterp( nz, ngrdcol, &
         502    17870112 :                           p_mid, p_out, input_var, &
         503    17870112 :                           interp_var )
         504             :
         505             :     implicit none
         506             :
         507             :     !------------------------ Input Variables ------------------------
         508             :     integer , intent(in)  :: &
         509             :       nz, &
         510             :       ngrdcol
         511             :
         512             :     real( kind = core_rknd ), dimension(ngrdcol,nz), intent(in)  :: &
         513             :       p_mid,      & ! input level pressure levels
         514             :       input_var     ! input  array
         515             :
         516             :     real( kind = core_rknd ), intent(in)  :: &
         517             :       p_out         ! output pressure level
         518             :
         519             :     !------------------------ Output Variables ------------------------
         520             :     real( kind = core_rknd ), dimension(ngrdcol), intent(out) :: &
         521             :       interp_var    ! output array (interpolated)
         522             :
         523             :     !------------------------ Local Variables ------------------------
         524             :     integer :: &
         525             :       i, k,   & ! Loop indices
         526             :       k_upper    ! Level indices for interpolation
         527             :
         528             :     real( kind = core_rknd ) :: &
         529             :       dpu,  & ! upper level pressure difference
         530             :       dpl     ! lower level pressure difference
         531             :
         532             :     logical :: &
         533             :       l_found,  & ! true if input levels found
         534             :       l_error     ! true if error
         535             :
         536             :     !------------------------ Begin Code ------------------------
         537             :
         538             :     ! Initialize index array and logical flags
         539    17870112 :     l_error = .false.
         540             :
         541             :     ! If we've fallen through the k=1,nz-1 loop, we cannot interpolate and
         542             :     ! must extrapolate from the bottom or top data level for at least some
         543             :     ! of the longitude points.
         544             :     !$acc parallel loop gang vector default(present)
         545   298389312 :     do i = 1, ngrdcol
         546             :
         547   298389312 :       if ( p_out >= p_mid(i,1) ) then
         548             :
         549    54547740 :         interp_var(i) = input_var(i,1)
         550             :
         551   225971460 :       elseif ( p_out <= p_mid(i,nz) ) then
         552             :
         553           0 :         interp_var(i) = input_var(i,nz)
         554             :
         555             :       else
         556             :
         557             :         l_found = .false.
         558  2966159100 :         k_upper = 1
         559             :
         560             :         ! Store level indices for interpolation.
         561             :         ! If all indices for this level have been found,
         562             :         ! do the interpolation
         563  2966159100 :         do k = 1, nz-1
         564  2966159100 :           if ( p_mid(i,k) > p_out .and. p_out >= p_mid(i,k+1) ) then
         565             :             l_found = .true.
         566             :             k_upper = k
         567             :             exit
         568             :           end if
         569             :         end do
         570             :
         571             :         if ( .not. l_found ) then
         572   225971460 :           l_error = .true.
         573             :         end if
         574             :
         575   225971460 :         dpu = p_mid(i,k_upper) - p_out
         576   225971460 :         dpl = p_out - p_mid(i,k_upper+1)
         577             :         interp_var(i) = ( input_var(i,k_upper)*dpl + input_var(i,k_upper+1)*dpu ) &
         578   225971460 :                         / ( dpl + dpu )
         579             :       end if
         580             :
         581             :     end do
         582             :     !$acc end parallel loop
         583             :
         584    17870112 :     return
         585             :
         586             :   end subroutine pvertinterp
         587             :
         588             : !-------------------------------------------------------------------------------
         589           0 :   subroutine lin_interpolate_on_grid &
         590           0 :              ( nparam, xlist, tlist, xvalue, &
         591             :                tvalue )
         592             :
         593             : ! Description:
         594             : !   Linear interpolation for 25 June 1996 altocumulus case.
         595             :
         596             : !   For example, to interpolate between two temperatures in space, put
         597             : !   your spatial coordinates in x-list and your temperature values in
         598             : !   tlist.  The point in question should have its spatial value stored
         599             : !   in xvalue, and tvalue will be the temperature at that point.
         600             :
         601             : ! Author: Michael Falk for COAMPS.
         602             : !-------------------------------------------------------------------------------
         603             :
         604             :     use constants_clubb, only: fstderr ! Constant
         605             :
         606             :     use clubb_precision, only: &
         607             :         core_rknd ! Variable(s)
         608             :
         609             :     use error_code, only: &
         610             :         clubb_at_least_debug_level ! Error indicator
         611             :
         612             :     implicit none
         613             :
         614             :     ! Input Variables
         615             :     integer, intent(in) :: nparam ! Number of parameters in xlist and tlist
         616             :
         617             :     ! Input/Output Variables
         618             :     real( kind = core_rknd ), intent(inout), dimension(nparam) ::  &
         619             :       xlist,  & ! List of x-values (independent variable)
         620             :       tlist     ! List of t-values (dependent variable)
         621             :
         622             :     real( kind = core_rknd ), intent(in) ::  &
         623             :       xvalue  ! x-value at which to interpolate
         624             :
         625             :     real( kind = core_rknd ), intent(inout) ::  &
         626             :       tvalue  ! t-value solved by interpolation
         627             :
         628             :     ! Local variables
         629             :     integer ::  &
         630             :       i     ! Loop control variable
         631             :
         632             :     integer ::  &
         633             :       bottombound, & ! Index of the smaller value in the linear interpolation
         634             :       topbound       ! Index of the larger value in the linear interpolation
         635             :
         636             : !-------------------------------------------------------------------------------
         637             : !
         638             : ! Assure that the elements are in order so that the interpolation is between
         639             : ! the two closest points to the point in question.
         640             : !
         641             : !-------------------------------------------------------------------------------
         642             :
         643           0 :      if ( clubb_at_least_debug_level( 0 ) ) then
         644           0 :        do i=2,nparam
         645           0 :          if ( xlist(i) <= xlist(i-1) ) then
         646           0 :            write(fstderr,*) "xlist must be sorted for lin_interpolate_on_grid."
         647           0 :            error stop
         648             :          end if
         649             :        end do
         650             :      end if
         651             : !-------------------------------------------------------------------------------
         652             : !
         653             : ! If the point in question is larger than the largest x-value or
         654             : ! smaller than the smallest x-value, crash.
         655             : !
         656             : !-------------------------------------------------------------------------------
         657             :
         658           0 :     if ( clubb_at_least_debug_level( 0 ) ) then
         659           0 :         if ( (xvalue < xlist(1)) .or. (xvalue > xlist(nparam)) ) then
         660           0 :           write(fstderr,*) "lin_interpolate_on_grid: Value out of range"
         661           0 :           error stop
         662             :         end if
         663             :     end if
         664             :
         665             : !-------------------------------------------------------------------------------
         666             : !
         667             : ! Find the correct top and bottom bounds, do the interpolation, return c
         668             : ! the value.
         669             : !
         670             : !-------------------------------------------------------------------------------
         671             :
         672             :     topbound = -1
         673             :     bottombound = -1
         674             :
         675           0 :     do i=2,nparam
         676           0 :       if ( (xvalue >= xlist(i-1)) .and. (xvalue <= xlist(i)) ) then
         677           0 :         bottombound = i-1
         678           0 :         topbound    = i
         679             :       end if
         680             :     end do
         681             :
         682           0 :     if ( clubb_at_least_debug_level( 1 ) .and. (topbound == -1 .or. bottombound == -1) ) then
         683           0 :       write(fstderr,*) "Sanity check failed! xlist is not properly sorted"
         684           0 :       write(fstderr,*) "in lin_interpolate_on_grid."
         685             :     end if
         686             :
         687             :     tvalue =  &
         688           0 :     lin_interpolate_two_points( xvalue, xlist(topbound), xlist(bottombound),  &
         689           0 :             tlist(topbound), tlist(bottombound) )
         690             :
         691           0 :     return
         692             :   end subroutine lin_interpolate_on_grid
         693             :
         694             : end module interpolation

Generated by: LCOV version 1.14