LCOV - code coverage report
Current view: | top level - utils - table_interp_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 27 | 55 | 49.1 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 3 | 6 | 50.0 %

              Line data    Source code
           1             : !----------------------------------------------------------------------------
           2             : ! Utility module used for interpolation of aerosol optics table
           3             : !  NOTE: Results will be set to table edges for interpolations beyond
           4             : !        the edges -- no extropolations
           5             : !----------------------------------------------------------------------------
           6             : module table_interp_mod
           7             :   use shr_kind_mod, only: r8=>shr_kind_r8
           8             :
           9             :   implicit none
          10             :
          11             :   private
          12             :   public :: table_interp
          13             :   public :: table_interp_wghts
          14             :   public :: table_interp_calcwghts
          15             :
          16             :   ! overload the interpolation routines
          17             :   interface table_interp
          18             :      module procedure interp1d
          19             :      module procedure interp2d
          20             :      module procedure interp4d
          21             :   end interface table_interp
          22             :
          23             :   ! interpolation weights and indices
          24             :   type :: table_interp_wghts
          25             :      real(r8) :: wt1
          26             :      real(r8) :: wt2
          27             :      integer  :: ix1
          28             :      integer  :: ix2
          29             :   end type table_interp_wghts
          30             :
          31             : contains
          32             :
          33             :   !--------------------------------------------------------------------------
          34             :   ! 1-D interpolation
          35             :   !--------------------------------------------------------------------------
          36           0 :   pure function interp1d( ncol, nxs, xwghts, tbl ) result(res)
          37             :
          38             :     integer, intent(in)  :: ncol                         ! number of model columns
          39             :     integer, intent(in)  :: nxs                          ! table size
          40             :     real(r8), intent(in) :: tbl(nxs)                     ! table values to be interpolated
          41             :     type(table_interp_wghts), intent(in) :: xwghts(ncol) ! interpolation weights and indices
          42             :
          43             :     real(r8) :: res(ncol)
          44             :
          45             :     integer :: i
          46             :
          47           0 :     do i = 1,ncol
          48             :
          49           0 :        res(i) = xwghts(i)%wt1*tbl(xwghts(i)%ix1) &
          50           0 :               + xwghts(i)%wt2*tbl(xwghts(i)%ix2)
          51             :
          52             :     end do
          53             :
          54           0 :   end function interp1d
          55             :
          56             :   !--------------------------------------------------------------------------
          57             :   ! 2-D interpolation
          58             :   !--------------------------------------------------------------------------
          59 64394521344 :   pure function interp2d( ncoef, ncol, nxs, nys, xwghts, ywghts, tbl ) result(res)
          60             :
          61             :     integer, intent(in)  :: ncoef                        ! number chebyshev coefficients
          62             :     integer, intent(in)  :: ncol                         ! number of model columns
          63             :     integer, intent(in)  :: nxs                          ! table x-dimension size
          64             :     integer, intent(in)  :: nys                          ! table y-dimension size
          65             :     real(r8), intent(in) :: tbl(ncoef,nxs,nys)           ! table values to be interpolated
          66             :     type(table_interp_wghts), intent(in) :: xwghts(ncol) ! x interpolation weights and indices
          67             :     type(table_interp_wghts), intent(in) :: ywghts(ncol) ! y interpolation weights and indices
          68             :
          69             :     real(r8) :: res(ncoef,ncol)
          70             :
          71 16098630336 :     real(r8) :: fx(ncoef,2)
          72             :
          73             :     integer :: i
          74             :
          75 >26880*10^7 :     do i = 1,ncol
          76             :
          77             :        ! interp x dir
          78 >75813*10^7 :        fx(:,1) = xwghts(i)%wt1*tbl(:,xwghts(i)%ix1,ywghts(i)%ix1) & ! @ y1
          79 >22743*10^8 :                + xwghts(i)%wt2*tbl(:,xwghts(i)%ix2,ywghts(i)%ix1)
          80 >25271*10^7 :        fx(:,2) = xwghts(i)%wt1*tbl(:,xwghts(i)%ix1,ywghts(i)%ix2) & ! @ y2
          81 >17689*10^8 :                + xwghts(i)%wt2*tbl(:,xwghts(i)%ix2,ywghts(i)%ix2)
          82             :
          83             :        ! interp y dir
          84 >15323*10^8 :        res(:,i) = ywghts(i)%wt1*fx(:,1) + ywghts(i)%wt2*fx(:,2)
          85             :
          86             :     end do
          87             :
          88 16098630336 :   end function interp2d
          89             :
          90             :   !--------------------------------------------------------------------------
          91             :   ! 4-D interpolation
          92             :   !--------------------------------------------------------------------------
          93           0 :   pure function interp4d( ncol, nxs, nys, nzs, nts, xwghts, ywghts, zwghts, twghts, tbl ) result(res)
          94             :
          95             :     integer, intent(in)  :: ncol                         ! number of model columns
          96             :     integer, intent(in)  :: nxs                          ! table x-dimension size
          97             :     integer, intent(in)  :: nys                          ! table y-dimension size
          98             :     integer, intent(in)  :: nzs                          ! table z-dimension size
          99             :     integer, intent(in)  :: nts                          ! table t-dimension size
         100             :     real(r8), intent(in) :: tbl(nxs,nys,nzs,nts)         ! table values to be interpolated
         101             :     type(table_interp_wghts), intent(in) :: xwghts(ncol) ! x interpolation weights and indices
         102             :     type(table_interp_wghts), intent(in) :: ywghts(ncol) ! y interpolation weights and indices
         103             :     type(table_interp_wghts), intent(in) :: zwghts(ncol) ! z interpolation weights and indices
         104             :     type(table_interp_wghts), intent(in) :: twghts(ncol) ! t interpolation weights and indices
         105             :
         106             :     real(r8) :: res(ncol)
         107             :
         108             :     real(r8) :: fx(8)
         109             :     real(r8) :: fy(4)
         110             :     real(r8) :: fz(2)
         111             :
         112             :     integer :: i
         113             :
         114           0 :     do i = 1,ncol
         115             :
         116             :        ! interp x dir
         117           0 :        fx(1) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix1,zwghts(i)%ix1,twghts(i)%ix1) & ! @ y1, z1, t1
         118           0 :              + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix1,zwghts(i)%ix1,twghts(i)%ix1)
         119           0 :        fx(2) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix2,zwghts(i)%ix1,twghts(i)%ix1) & ! @ y2, z1, t1
         120           0 :              + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix2,zwghts(i)%ix1,twghts(i)%ix1)
         121             :
         122           0 :        fx(3) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix1,zwghts(i)%ix2,twghts(i)%ix1) & ! @ y1, z2, t1
         123           0 :              + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix1,zwghts(i)%ix2,twghts(i)%ix1)
         124             :        fx(4) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix2,zwghts(i)%ix2,twghts(i)%ix1) & ! @ y2, z2, t1
         125           0 :              + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix2,zwghts(i)%ix2,twghts(i)%ix1)
         126             :
         127           0 :        fx(5) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix1,zwghts(i)%ix1,twghts(i)%ix2) & ! @ y1, z1, t2
         128           0 :              + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix1,zwghts(i)%ix1,twghts(i)%ix2)
         129             :        fx(6) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix2,zwghts(i)%ix1,twghts(i)%ix2) & ! @ y2, z1, t2
         130           0 :              + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix2,zwghts(i)%ix1,twghts(i)%ix2)
         131             :
         132             :        fx(7) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix1,zwghts(i)%ix2,twghts(i)%ix2) & ! @ y1, z2, t2
         133           0 :              + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix1,zwghts(i)%ix2,twghts(i)%ix2)
         134             :        fx(8) = xwghts(i)%wt1*tbl(xwghts(i)%ix1,ywghts(i)%ix2,zwghts(i)%ix2,twghts(i)%ix2) & ! @ y2, z2, t2
         135           0 :              + xwghts(i)%wt2*tbl(xwghts(i)%ix2,ywghts(i)%ix2,zwghts(i)%ix2,twghts(i)%ix2)
         136             :
         137             :        ! interp y dir
         138           0 :        fy(1) = ywghts(i)%wt1*fx(1) + ywghts(i)%wt2*fx(2) ! @ z1, t1
         139           0 :        fy(2) = ywghts(i)%wt1*fx(3) + ywghts(i)%wt2*fx(4) ! @ z2, t1
         140           0 :        fy(3) = ywghts(i)%wt1*fx(5) + ywghts(i)%wt2*fx(6) ! @ z1, t2
         141           0 :        fy(4) = ywghts(i)%wt1*fx(7) + ywghts(i)%wt2*fx(8) ! @ z2, t2
         142             :
         143             :        ! interp z dir
         144           0 :        fz(1) = zwghts(i)%wt1*fy(1) + zwghts(i)%wt2*fy(2) ! @ t1
         145           0 :        fz(2) = zwghts(i)%wt1*fy(3) + zwghts(i)%wt2*fy(4) ! @ t2
         146             :
         147             :        ! interp t dir
         148           0 :        res(i) = twghts(i)%wt1*fz(1) + twghts(i)%wt2*fz(2)
         149             :
         150             :     end do
         151             :
         152           0 :   end function interp4d
         153             :
         154             :   !--------------------------------------------------------------------------
         155             :   ! determines interpolation weights and indices for given values at the model columns
         156             :   !--------------------------------------------------------------------------
         157 16653755520 :   pure function table_interp_calcwghts( ngrid, xgrid, ncols, xcols ) result(wghts)
         158             :
         159             :     integer,  intent(in) :: ngrid        ! number of grid point values
         160             :     real(r8), intent(in) :: xgrid(ngrid) ! grid point values
         161             :     integer,  intent(in) :: ncols        ! number of model columns
         162             :     real(r8), intent(in) :: xcols(ncols) ! values at the model columns
         163             :
         164             :     type(table_interp_wghts) :: wghts(ncols) ! interpolations weights at the model columns
         165             :
         166             :     integer :: i
         167 16653755520 :     real(r8) :: xs(ncols)
         168             :
         169 >27807*10^7 :     xs(:) = xcols(:)
         170             :
         171             :     ! do not extrapolate beyond the edges of the table
         172 >27807*10^7 :     where(xs < xgrid(1))
         173 16653755520 :        xs = xgrid(1)
         174             :     end where
         175 >27807*10^7 :     where(xs > xgrid(ngrid))
         176             :        xs = xgrid(ngrid)
         177             :     end where
         178             :
         179 >27807*10^7 :     do i = 1,ncols
         180 >26142*10^7 :        wghts(i)%ix2 = find_index(ngrid,xgrid,xs(i))
         181 >26142*10^7 :        wghts(i)%ix1 = wghts(i)%ix2 - 1
         182 >26142*10^7 :        wghts(i)%wt1 = (xgrid(wghts(i)%ix2)-xs(i)) &
         183 >52285*10^7 :                      /(xgrid(wghts(i)%ix2)-xgrid(wghts(i)%ix1))
         184 >27807*10^7 :        wghts(i)%wt2 = 1._r8 - wghts(i)%wt1
         185             :     end do
         186             :
         187 16653755520 :   end function table_interp_calcwghts
         188             :
         189             :   ! private methods
         190             :   !--------------------------------------------------------------------------
         191             :   !--------------------------------------------------------------------------
         192             :   ! determines last index of grid vals of which is greater then or equal to
         193             :   ! value vx
         194             :   !--------------------------------------------------------------------------
         195 >26142*10^7 :   pure function find_index( nvals, vals, vx ) result(res)
         196             :     integer,  intent(in) :: nvals
         197             :     real(r8), intent(in) :: vals(nvals)
         198             :     real(r8), intent(in) :: vx
         199             :     integer :: res
         200             :
         201             :     integer :: ndx
         202             :
         203 >26142*10^7 :     res = -1
         204             :
         205 >11667*10^8 :     find_ndx: do ndx = 2, nvals
         206 >11667*10^8 :        if (vals(ndx)>=vx) then
         207             :           res = ndx
         208             :           exit find_ndx
         209             :        end if
         210             :     end do find_ndx
         211             :
         212 >26142*10^7 :   end function find_index
         213             :
         214           0 : end module table_interp_mod

Generated by: LCOV version 1.14