LCOV - code coverage report
Current view: | top level - control - getinterpnetcdfdata.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 107 | 0.0 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 0 | 2 | 0.0 %

              Line data    Source code
           1             : module getinterpnetcdfdata
           2             :
           3             : ! Description:
           4             : !   Routines for extracting a column from a netcdf file
           5             : !
           6             : ! Author:
           7             : !
           8             : ! Modules Used:
           9             : !
          10             :   use cam_abortutils, only: endrun
          11             :   use pmgrid,         only: plev
          12             :   use cam_logfile,    only: iulog
          13             :
          14             :   implicit none
          15             :   private
          16             : !
          17             : ! Public Methods:
          18             : !
          19             :   public getinterpncdata
          20             :
          21             : contains
          22             :
          23           0 : subroutine getinterpncdata( NCID, camlat, camlon, TimeIdx, &
          24             :    varName, have_surfdat, surfdat, fill_ends, scm_crm_mode, &
          25           0 :    press, npress, ps, hyam, hybm, outData, STATUS )
          26             :
          27             : !     getinterpncdata: extracts the entire level dimension for a
          28             : !     particular lat,lon,time from a netCDF file
          29             : !     and interpolates it onto the input pressure levels, placing
          30             : !     result in outData, and the error status inx STATUS
          31             :
          32             :    use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8
          33             :    use shr_scam_mod, only: shr_scam_GetCloseLatLon
          34             :    use netcdf
          35             :    implicit none
          36             : !-----------------------------------------------------------------------
          37             :
          38             :
          39             : !     ---------- inputs ------------
          40             :
          41             :    integer, intent(in)  :: NCID          ! NetCDF ID
          42             :    integer, intent(in)  :: TimeIdx       ! time index
          43             :    real(r8), intent(in) :: camlat,camlon ! target lat and lon to be extracted
          44             :    logical, intent(in)  :: have_surfdat  ! is surfdat provided
          45             :    logical, intent(in)  :: fill_ends     ! extrapolate the end values
          46             :    logical, intent(in)  :: scm_crm_mode  ! scam column radiation mode
          47             :    integer, intent(in)  :: npress        ! number of dataset pressure levels
          48             :    real(r8), intent(in) :: press(npress) ! dataset pressure levels
          49             :    real(r8), intent(in) :: ps            ! surface pressure
          50             :    real(r8), intent(in) :: hyam(:)       ! dataset hybrid midpoint pressure levels
          51             :    real(r8), intent(in) :: hybm(:)       ! dataset hybrid midpoint pressure levels
          52             :
          53             : !     ---------- outputs ----------
          54             :
          55             :    real(r8), intent(inout)  :: outData(:)    ! interpolated output data
          56             :    integer, intent(out)   :: STATUS        ! return status of netcdf calls
          57             :
          58             : !     -------  locals ---------
          59             :
          60             :    real(r8)  surfdat       ! surface value to be added before interpolation
          61             :    integer nlev          ! number of levels in dataset
          62             :    integer     latIdx        ! latitude index
          63             :    integer     lonIdx        ! longitude index
          64           0 :    real(r8),allocatable :: tmp(:)
          65             :    real(r8)        closelat,closelon
          66             :    real(r8)        dx, dy, m             ! slope for interpolation of surfdat
          67             :    integer     varID
          68             :    integer     var_ndims
          69             :    integer     dims_set
          70             :    integer     i
          71             :    integer     var_dimIDs( NF90_MAX_VAR_DIMS )
          72             :    integer     start( NF90_MAX_VAR_DIMS )
          73             :    integer     count( NF90_MAX_VAR_DIMS )
          74             :
          75             :    character   varName*(*)
          76             :    character   dim_name*( 256 )
          77             :    real(r8)        missing_val
          78             :    logical     usable_var
          79             :
          80             : !     -------  code ---------
          81             :
          82           0 :    call shr_scam_GetCloseLatLon(ncid,camlat,camlon,closelat,closelon,latidx,lonidx)
          83             :
          84             : !
          85             : ! Check mode: double or single precision
          86             : !
          87             :
          88             : !
          89             : ! Get var ID.  Check to make sure it's there.
          90             : !
          91           0 :    STATUS = NF90_INQ_VARID( NCID, varName, varID )
          92             :
          93           0 :    if ( STATUS .NE. NF90_NOERR )  return
          94             :
          95             : !
          96             : ! Check the var variable's information with what we are expecting
          97             : ! it to be.
          98             : !
          99             :
         100           0 :    STATUS = NF90_INQUIRE_VARIABLE( NCID, varID, ndims=var_ndims )
         101           0 :    if ( var_ndims .GT. 4 ) then
         102           0 :       write(iulog,* ) 'ERROR - extractdata.F: The input var',varName, &
         103           0 :          'has', var_ndims, 'dimensions'
         104           0 :       STATUS = -1
         105             :    endif
         106             :
         107             : !
         108             : !     surface variables
         109             : !
         110           0 :    if ( var_ndims .EQ. 0 ) then
         111           0 :       STATUS = NF90_GET_VAR( NCID, varID, outData )
         112           0 :       return
         113             :    endif
         114             :
         115           0 :    STATUS = NF90_INQUIRE_VARIABLE( NCID, varID, dimids=var_dimIDs )
         116           0 :    if ( STATUS .NE. NF90_NOERR ) then
         117           0 :       write(iulog,* ) 'ERROR - extractdata.F:Cant get dimension IDs for', varName
         118           0 :       return
         119             :    endif
         120             : !
         121             : !     Initialize the start and count arrays
         122             : !
         123           0 :    dims_set = 0
         124           0 :    nlev = 1
         125           0 :    do i =  var_ndims, 1, -1
         126             :
         127           0 :       usable_var = .false.
         128           0 :       STATUS = NF90_INQUIRE_DIMENSION( NCID, var_dimIDs( i ), dim_name )
         129             :
         130           0 :       if ( dim_name .EQ. 'lat' ) then
         131           0 :          start( i ) =  latIdx
         132           0 :          count( i ) = 1           ! Extract a single value
         133           0 :          dims_set = dims_set + 1
         134           0 :          usable_var = .true.
         135             :       endif
         136             :
         137           0 :       if ( dim_name .EQ. 'lon' .or. dim_name .EQ. 'ncol' .or. dim_name .EQ. 'ncol_d' ) then
         138           0 :          start( i ) = lonIdx
         139           0 :          count( i ) = 1           ! Extract a single value
         140           0 :          dims_set = dims_set + 1
         141           0 :          usable_var = .true.
         142             :       endif
         143             :
         144           0 :       if ( dim_name .EQ. 'lev' ) then
         145           0 :          STATUS = NF90_INQUIRE_DIMENSION( NCID, var_dimIDs( i ), len=nlev )
         146           0 :          start( i ) = 1
         147           0 :          count( i ) = nlev       ! Extract all levels
         148           0 :          dims_set = dims_set + 1
         149           0 :          usable_var = .true.
         150             :       endif
         151             :
         152           0 :       if ( dim_name .EQ. 'ilev' ) then
         153           0 :          STATUS = NF90_INQUIRE_DIMENSION( NCID, var_dimIDs( i ), len=nlev )
         154           0 :          start( i ) = 1
         155           0 :          count( i ) = nlev        ! Extract all levels
         156           0 :          dims_set = dims_set + 1
         157           0 :          usable_var = .true.
         158             :       endif
         159             :
         160           0 :       if ( dim_name .EQ. 'time' .OR. dim_name .EQ. 'tsec' ) then
         161           0 :          start( i ) = TimeIdx
         162           0 :          count( i ) = 1           ! Extract a single value
         163           0 :          dims_set = dims_set + 1
         164             :          usable_var = .true.
         165             :       endif
         166             :
         167           0 :       if ( usable_var .EQV. .false. ) then
         168           0 :          write(iulog,* )'ERROR - extractdata.F: The input var ',varName, &
         169           0 :             ' has an unusable dimension ', dim_name
         170           0 :          STATUS = 1
         171             :       endif
         172             :    end do
         173             :
         174           0 :    if ( dims_set .NE. var_ndims ) then
         175           0 :       write(iulog,* )'ERROR - extractdata.F: Could not find all the', &
         176           0 :          ' dimensions for input var ', varName
         177           0 :       write(iulog,* )'Found ',dims_set, ' of ',var_ndims
         178           0 :       STATUS = 1
         179             :    endif
         180             :
         181           0 :    allocate(tmp(nlev+1))
         182             :
         183           0 :    STATUS = NF90_GET_VAR( NCID, varID, tmp, start, count )
         184             :
         185           0 :    if ( STATUS .NE. NF90_NOERR ) then
         186           0 :       write(iulog,* )'ERROR - extractdata.F: Could not get data for input var ', varName
         187           0 :       return
         188             :    endif
         189             :
         190           0 :    if ( nlev .eq. 1 ) then
         191           0 :       outdata(1) = tmp(1)
         192           0 :       return                 ! no need to do interpolation
         193             :    endif
         194             : !   if ( use_camiop .and. nlev.eq.plev) then
         195           0 :    if ( nlev.eq.plev .or. nlev.eq.plev+1) then
         196           0 :       outData(:nlev)= tmp(:nlev)! no need to do interpolation
         197             :    else
         198             : !
         199             : !     add the surface data if available, else
         200             : !     fill in missing surface data by extrapolation
         201             : !
         202           0 :       if(.not.scm_crm_mode) then
         203           0 :          if ( have_surfdat ) then
         204           0 :             tmp(npress) = surfdat
         205             :          else
         206           0 :             dy = press(npress-1) - press(npress-2)
         207           0 :             dx = tmp(npress-1) - tmp(npress-2)
         208           0 :             if ( dx .ne. 0.0_r8 ) then
         209           0 :                m = dy/dx
         210           0 :                tmp(npress) = ((press(npress) - press(npress-1)) / m ) + tmp(npress-1)
         211             :             else
         212           0 :                tmp(npress) = tmp(npress-1)
         213             :             endif
         214           0 :             surfdat = tmp(npress)
         215             :          endif
         216             :       endif
         217             :
         218             : #if DEBUG > 1
         219             : !
         220             : !     check data for missing values
         221             : !
         222             :
         223             :       STATUS = NF90_GET_ATT( NCID, varID, 'missing_value', missing_val )
         224             :       if ( STATUS .NE. NF90_NOERR ) then
         225             :          missing_val = -9999999.0_r8
         226             :       endif
         227             : !
         228             : ! reset status to zero
         229             : !
         230             :       STATUS = 0
         231             : !
         232             :       do i=1, npress
         233             :          if ( tmp(i) .eq. missing_val ) then
         234             :             write(iulog,*) 'ERROR - missing value found in ', varname
         235             :             write(iulog,*) 'time,lat,lon,lev = ' ,timeidx, latidx, lonidx, i
         236             :             stop
         237             :          endif
         238             :       enddo
         239             : #endif
         240             : !
         241           0 :       call interplevs( tmp(:npress), press, npress, ps, fill_ends, hyam, hybm, outdata )
         242             :
         243             :    endif
         244             :
         245           0 :    deallocate(tmp)
         246           0 :    return
         247           0 :  end subroutine getinterpncdata
         248             :
         249           0 : subroutine interplevs( inputdata,   dplevs,   nlev, &
         250           0 :                        ps, fill_ends, hyam, hybm, outdata)
         251             :
         252             :    use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8
         253             :    use interpolate_data, only: lininterp
         254             :    implicit none
         255             :
         256             : !
         257             : !     WARNING: ps, siga and sigb must be initialized before calling this routine
         258             : !
         259             :
         260             : !------------------------------Commons----------------------------------
         261             :
         262             : !-----------------------------------------------------------------------
         263             :
         264             :
         265             : !     ------- inputs -----------
         266             :    integer, intent(in) :: nlev                 ! num press levels in dataset
         267             :
         268             :    real(r8), intent(in) :: ps                  ! surface pressure
         269             :    real(r8), intent(in) :: hyam(:)             ! a midpoint pressure
         270             :    real(r8), intent(in) :: hybm(:)             ! b midpoint pressure
         271             :    real(r8), intent(in) :: inputdata(nlev)     ! data from netcdf dataset
         272             :    real(r8), intent(in) :: dplevs(nlev)        ! input data pressure levels
         273             :
         274             :    logical, intent(in) :: fill_ends            ! fill in missing end values(used for
         275             :                                                ! global model datasets)
         276             :
         277             :
         278             : ! ------- outputs ----------
         279             :    real(r8), intent(inout) :: outdata(:)      ! interpolated column data
         280             :
         281             : ! ------- locals -----------
         282             :
         283             :    real(r8) mplevs( PLEV )
         284             :    real(r8) interpdata( PLEV )
         285             :
         286             :
         287             :    integer dstart_lev, dend_lev
         288             :    integer mstart_lev, mend_lev
         289             :    integer data_nlevs, model_nlevs, i
         290             :    integer STATUS
         291             :
         292             : !
         293             : !     Initialize  model_pressure_levels.  ps should be set in the calling
         294             : !     routine to the value in the dataset
         295             : !
         296           0 :    do i = 1, plev
         297           0 :       mplevs( i ) = 1000.0_r8 * hyam( i ) + ps * hybm( i ) / 100.0_r8
         298             :    end do
         299             : !
         300             : !     the following algorithm assumes that pressures are increasing in the
         301             : !     arrays
         302             : !
         303             : !
         304             : !     Find the data pressure levels that are just outside the range
         305             : !     of the model pressure levels, and that contain valid values
         306             : !
         307             :    dstart_lev = 1
         308           0 :    do i= 1, nlev
         309           0 :       if ( dplevs(i) .LE. mplevs(1) ) dstart_lev  = i
         310             :    end do
         311             :
         312             :    dend_lev = nlev
         313           0 :    do i= nlev, 1, -1
         314           0 :       if ( dplevs(i) .GE. mplevs(plev) ) then
         315           0 :          dend_lev  = i
         316             :       endif
         317             :    end do
         318             : !
         319             : !     Find the model pressure levels that are just inside the range
         320             : !     of the data pressure levels
         321             : !
         322             :    mstart_lev = 1
         323           0 :    do i=plev, 1, -1
         324           0 :       if ( mplevs( i ) .GE. dplevs( dstart_lev ) )  mstart_lev = i
         325             :    end do
         326             :
         327             :    mend_lev = plev
         328           0 :    do i=1,plev
         329           0 :       if ( mplevs( i ) .LE. dplevs( dend_lev ) ) mend_lev = i
         330             :    end do
         331             :
         332           0 :    data_nlevs = dend_lev - dstart_lev +1
         333           0 :    model_nlevs = mend_lev - mstart_lev +1
         334             :
         335           0 :    call lininterp (inputdata(dstart_lev:dend_lev),dplevs(dstart_lev:dend_lev),data_nlevs, &
         336           0 :                    interpdata,mplevs(mstart_lev:mend_lev),model_nlevs)
         337             : !
         338             : !     interpolate data onto the model pressure levels
         339             : !
         340             : !!$   call lininterp (inputdata,dplevs,nlev, &
         341             : !!$                outdata(:plev),mplevs,plev)
         342           0 :    do i=1 , model_nlevs
         343           0 :       outdata( i+mstart_lev-1 ) = interpdata( i )
         344             :    end do
         345             : !
         346             : !     fill in the missing end values
         347             : !           (usually  done if this is global model dataset)
         348             : !
         349           0 :    if ( fill_ends ) then
         350           0 :       do i=1, mstart_lev
         351           0 :          outdata(i) = inputdata(1)
         352             :       end do
         353           0 :       do i= mend_lev, plev
         354           0 :          outdata(i) = inputdata(nlev)
         355             :       end do
         356             :    end if
         357             :
         358           0 :    return
         359             : end subroutine interplevs
         360             : end module getinterpnetcdfdata

Generated by: LCOV version 1.14