ecearth3/sources/tm5mp/base/tmm_mf_ncep_cdc.F90

!###############################################################################
!
!  Interface to NCEP re-analysis data.
!  
!  Database
!  --------
!  
!    http://www.cdc.noaa.gov/cdc/reanalysis/reanalysis.shtml
!  
!  NetCDF files
!  ------------
!
!  global attributes:
!    Conventions                = "CDC Non-gridded" ;
!    title                      = "4x Daily Spectral Coefficients from the NMC Reanalysis for Natural Log of Pressure at the Surface" ;
!    base_date                  = 2000s, 1s, 1s ;
!    history                    = ... ;
!    description                = "Data is from NMC initialized reanalysis",
!    platform                   = "Model" ;
!    m_fastest                  = "F" ;
!    storage_type               = "Spectral" ;
!    compression_used           = "T" ;
!  
!
!  dimensions:
!    num_values = 4032 
!    level      = 28
!    num_mean   = 2 ;
!    time       = UNLIMITED ;
!
!  level data sets:
!
!    level             float    array[level]
!    units                      = "sigma_level"
!    long_name                  = "Sigma"
!    positive                   = "down"
!
!  spectral data sets:
!
!    mean             float     array[2*level,time]
!      long_name      string    "First Spectral Coefficients for Natural Log of Pressure at the Surface"
!      units          string    "-"
!      missing_value  float
!
!    add_offset       float     array[level,time]
!      long_name                "Add offset of Spectral Coefficients for Natural Log of Pressure"
!      units                    "-"
!      missing_value 
!
!    scale_factor     float     array[level,time]
!      long_name                "Scale Factor of Spectral Coefficients for Natural Log of Pressure"
!      units                    "-"
!      missing_value
!
!    time            double     array[time]
!      units                    "hours since 1-1-1 00:00:0.0" ;
!      long_name                "Time" ;
!      delta_t                  "0000-00-00 06:00:00" ;
!
!    pres            int        array[num_values,level,time]
!      long_name                "Spectral Coefficients for Natural Log of Pressure at the Surface"
!      units                    "natural log of pressure in centibars" ;
!      missing_value     
!      precision                = 0s ;
!      least_significant_digit  = 32767s ;
!      trunc_count              = 62s ;
!      trunc_type               = "Triangular" ;
!      var_desc                 = "Pressure",              "P" ;
!      dataset                  = "NMC Reanalysis",         "L" ;
!      level_desc               = "Surface",             "F" ;
!      statistic                = "Individual Obs\n",       "I" ;
!      parent_stat              = "Other\n",              "-" ;
!
!  gaussian or lon/lat data sets:
!
!    short pres(time, lat, lon) ;
!      long_name                = "4xDaily Surface Pressure" ;
!      valid_range              = 40000.f, 115000.f ;
!      actual_range             = 48540.f, 108990.f ;
!      units                    = "Pascals" ;
!      add_offset               = 367650.f ;
!      scale_factor             = 10.f ;
!      missing_value            = 32766s ;
!      precision                = -1s ;
!      least_significant_digit  = -1s ;
!      GRIB_id                  = 1s ;
!      GRIB_name                = "PRES" ;
!      var_desc                 = "Pressure\n", "CC" ;
!      dataset                  = "NMC Reanalysis\n", "L" ;
!
!  Date and time
!  -------------
!
!   Global attribute 'base_date' specifies the first date :
!     base_date = (/2000,01,01/)
!
!   Data set 'time' contains hours since year 1 or there about.
!   Substract the first element from the array to have an offset
!   in hours from midnight at base_date.
!
!  Spectral fields
!  -----------
!
!  From the 'discription':
!
!  For each latitude (going from north pole to south pole). 
!  the associated legendre functions are defined: 
!  Pbar(m,n,theta) = 
!    sqrt((2*n+1)*factorial(n-m)/(2*factorial(n+m)))*sin(theta)**m/(2**n*factorial(n)) 
!    times the (n+m)th derivative of (x**2-1)**n with respect to x=cos(theta)
!
!    note: theta = 0.5*pi - phi, where phi is latitude and theta is colatitude.
!           Therefore, cos(theta) = sin(phi) and sin(theta) = cos(phi).
!
!  where n is degree (subscript), m is order (superscript), 
!  and theta is colatitude in radians.  
!  The functions are orthogonal polynomials on the surface of the sphere and
!  are normalized so that the integral of (PBAR(n,m,theta)**2)*sin(theta) 
!  on the interval theta=0 to theta=pi equals 1.
!  Note that Pbar(0,0,theta) = sqrt(2)/2
!
!  Truncation:
!    number of data values : 4032 = 2016 complex numbers
!    number of spectral coeff = (T+1)*(T+2)/2 = 2016, thus T=62
!
!  The 'mean' are the first spectral coeff and should be added.
!
!  Packing
!  -------
!
!  From http://www.cdc.noaa.gov/PublicData/faq.html#12 :
!  "Most of the data in our netCDF files are packed. That is to say they have been
!  transformed by a scale factor and an add offset to reduce the storage needed to
!  two bytes per value. When you extract the short integers, you must unpack the
!  data to recover the correct floating point data values. Data files that contain
!  packed data will have a non-zero add offset and/or a scale factor not equal to 1.
!  The transformation is:
!    float_value = (short_value * scale_factor) + add_offset                 "
!
!
!### macro's ###################################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tmm.inc"
!
!###############################################################################

module tmm_mf_ncep_cdc

  use GO      , only : gol, goErr, goPr
  use GO      , only : TDate

  implicit none
  
  ! --- in/out ----------------------------
  
  private
  
  public  ::  TMeteoFile_ncep_cdc
  public  ::  Init, Done
  public  ::  Get
  public  ::  ReadRecord
  
  ! --- const ------------------------------
  
  character(len=*), parameter  ::  mname = 'tmm_mf_ncep_cdc'
  
  !--- type ---------------------------------
  
  type TMeteoFile_ncep_cdc
    ! field collection
    character(len=256)     ::  paramkeys
    type(TDate)            ::  trange(2)
    ! file names
    character(len=256)     ::  dir
    integer                ::  ccyy
    character(len=1)       ::  pathsep, namesep
  end type TMeteoFile_ncep_cdc
  
  
  ! --- interfaces -------------------------
  
  interface Init
    module procedure mf_Init
  end interface

  interface Done
    module procedure mf_Done
  end interface

  interface Get
    module procedure mf_Get
  end interface

  interface ReadRecord
    module procedure mf_ReadRecord
  end interface


contains


  ! ==============================================================


  subroutine mf_Init( mf, dir, archivekeys, paramkey, &
                               tref, t1, t2, status )
    
    use GO, only : TDate, Get, NewDate, operator(>)
    use GO, only : goVarValue
    
    ! --- in/out ----------------------------
    
    type(TMeteoFile_ncep_cdc), intent(out)  ::  mf
    character(len=*), intent(in)            ::  dir
    character(len=*), intent(in)            ::  archivekeys
    character(len=*), intent(in)            ::  paramkey
    type(TDate), intent(in)                 ::  tref, t1, t2
    integer, intent(out)                    ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_Init'

    ! --- local -------------------------------
    
    type(TDate)         ::  tend
    
    ! --- begin --------------------------------
 
    ! store directory:
    mf%dir = dir
    
    !
    ! extract fields from archivekey :
    !   mdir=spectral
    !
    mf%pathsep  = '/'
    call goVarValue( archivekeys, ';', 'pathsep', '=', mf%pathsep, status )
    if (status>0) then; write (gol,'("in ",a)') rname; call goErr; status=1; return; end if
    !
    mf%namesep  = '-'
    call goVarValue( archivekeys, ';', 'namesep', '=', mf%namesep, status )
    if (status>0) then; write (gol,'("in ",a)') rname; call goErr; status=1; return; end if
 
    ! extract year in filename from requested time range:
    call Get( t2, year=mf%ccyy )
    tend = NewDate( time6=(/mf%ccyy,12,31,18,00,00/) )
    if ( t2 > tend ) mf%ccyy = mf%ccyy + 1
    
    ! files valid for complete year:
    mf%trange(1) = NewDate( mf%ccyy-1, 12, 31, 18, 00, 00, 001 )
    mf%trange(2) = NewDate( mf%ccyy  , 12, 31, 18, 00, 00, 000 )

    ! files contain one param only:    
    mf%paramkeys   = '-'//trim(paramkey)//'-'
    
    ! dummy filename (might be used in error messages)
    mf%filename = 'ncep cdc file'

    ! ok
    status = 0
 
  end subroutine mf_Init
  
  
  ! ***

  
  subroutine mf_Done( mf, status )
  
    ! --- in/out ------------------------------------

    type(TMeteoFile_ncep_cdc), intent(inout)   ::  mf
    integer, intent(out)                       ::  status

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_Done'
    
    ! --- begin -------------------------------------
    
    ! files have been closed in ReadRecord/WriteRecord
    
    ! ok
    status = 0
    
  end subroutine mf_Done
  


  ! ***

  
  subroutine mf_Get( mf, status, trange1, trange2, paramkeys )
    
    use GO, only : TDate
    
    ! --- in/out ----------------------------
    
    type(TMeteoFile_ncep_cdc), intent(in)     ::  mf
    integer, intent(out)                      ::  status

    type(TDate), intent(out), optional        ::  trange1, trange2
    character(len=*), intent(out), optional   ::  paramkeys
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_Get'
    
    ! --- local --------------------------------
    
    ! --- begin --------------------------------

    ! time range:
    if ( present(trange1) ) trange1 = mf%trange(1)
    if ( present(trange2) ) trange2 = mf%trange(2)
    
    ! parameter names:
    if ( present(paramkeys) ) paramkeys = mf%paramkeys

    ! ok
    status = 0
 
  end subroutine mf_Get
  

  ! ***


  ! Return a field given parameter name, time, etc.
  ! Only one of grid types is filled!

  subroutine mf_ReadRecord( mf, paramkey, t1, t2, nuv, nw, &
                                gridtype, levi, &
                                lli, ll, sp_ll, &
                                ggi, gg, sp_gg, &
                                shi, sh, lnsp_sh, &
                                tmi, status )
                                
    use parray    , only : pa_Init, pa_Done
    use GO        , only : TDate
    use Grid      , only : TLevelInfo
    use Grid      , only : TllGridInfo, TggGridInfo, TshGridInfo
    use tmm_info  , only : TMeteoInfo

    ! --- in/out -------------------------------
    
    type(TMeteoFile_ncep_cdc), intent(inout)  ::  mf
    character(len=*), intent(in)              ::  paramkey
    type(TDate), intent(in)                   ::  t1, t2
    character(len=1), intent(in)              ::  nuv
    character(len=1), intent(in)              ::  nw
    character(len=2), intent(out)             ::  gridtype
    type(TLevelInfo), intent(out)             ::  levi
    type(TllGridInfo), intent(inout)          ::  lli
    real, pointer                             ::  ll(:,:,:)
    real, pointer                             ::  sp_ll(:,:)
    type(TggGridInfo), intent(inout)          ::  ggi
    real, pointer                             ::  gg(:,:)
    real, pointer                             ::  sp_gg(:)
    type(TshGridInfo), intent(inout)          ::  shi
    complex, pointer                          ::  sh(:,:)
    complex, pointer                          ::  lnsp_sh(:)
    type(TMeteoInfo), intent(out)             ::  tmi
    integer, intent(out)                      ::  status

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_ReadRecord'
    
    ! --- local -------------------------------
    
    real, pointer         ::  ll2(:,:,:)
    real, pointer         ::  gg2(:,:)
    complex, pointer      ::  sh2(:,:)
        
    ! --- begin ---------------------------------
    
    ! combined field ?
    select case ( paramkey )
   
      !
      ! total surface stress :  sstr^2 =  ewss^2 + nsss^2
      ! 
      case ( 'sstr' )
      
        ! read first field:
        call mf_ReadRecord_1( mf, 'ewss', t1, t2, nuv, nw, gridtype, levi, &
                                lli, ll, sp_ll, ggi, gg, sp_gg, shi, sh, lnsp_sh, &
                                tmi, status )
        IF_NOTOK_RETURN(status=1)  

        ! init pointer:
        call pa_Init( ll2 ) ; call pa_Init( gg2 ) ; call pa_Init( sh2 )

        ! read second field:
        call mf_ReadRecord_1( mf, 'nsss', t1, t2, nuv, nw, gridtype, levi, &
                                lli, ll2, sp_ll, ggi, gg2, sp_gg, shi, sh2, lnsp_sh, &
                                tmi, status )
         IF_NOTOK_RETURN(status=1)      

        ! process:
        select case ( gridtype )
          case ( 'll' ) ; ll = sqrt( ll**2 + ll2**2 )
          case ( 'gg' ) ; gg = sqrt( gg**2 + gg2**2 )
          case default
            write (gol,'("unsupported gridtype for sstr :",a)') gridtype; call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
        end select   

        ! clear pointers:
        call pa_Done( ll2 ) ; call pa_Done( gg2 ) ; call pa_Done( sh2 )
 
                write (gol,'("WARNING : adhoc constant value for ncep surface stress ...")'); call goPr
                gg = 0.1  ! N/m2
 
      !
      ! sea surface temperature :  skin temperture *  sea-mask
      ! 
      case ( 'sst' )
      
        ! read skin temperature:
        call mf_ReadRecord_1( mf, 'skt', t1, t2, nuv, nw, gridtype, levi, &
                                lli, ll, sp_ll, ggi, gg, sp_gg, shi, sh, lnsp_sh, &
                                tmi, status )
        IF_NOTOK_RETURN(status=1)  

        ! init pointers:
        call pa_Init( ll2 ) ; call pa_Init( gg2 ) ; call pa_Init( sh2 )

        ! read land-sea mask:
        call mf_ReadRecord_1( mf, 'lsm', t1, t2, nuv, nw, gridtype, levi, &
                                lli, ll2, sp_ll, ggi, gg2, sp_gg, shi, sh2, lnsp_sh, &
                                tmi, status )
         IF_NOTOK_RETURN(status=1)      

        ! process:
        select case ( gridtype )
          case ( 'll' ) ; ll = ll * max( 0.0, 1.0 - ll2 )
          case ( 'gg' ) ; gg = gg * max( 0.0, 1.0 - gg2 )
          case default
            write (gol,'("unsupported gridtype for sst :",a)') gridtype; call goErr
            write (gol,'("in ",a)') rname; call goErr; status=1; return
        end select   

        ! clear pointers:
        call pa_Done( ll2 ) ; call pa_Done( gg2 ) ; call pa_Done( sh2 )

      !
      ! no specials ...
      ! 
      case default
      
        call mf_ReadRecord_1( mf, paramkey, t1, t2, nuv, nw, gridtype, levi, &
                                lli, ll, sp_ll, ggi, gg, sp_gg, shi, sh, lnsp_sh, &
                                tmi, status )
         IF_NOTOK_RETURN(status=1)
      
    end select
    
    ! ok
    status = 0

  end subroutine mf_ReadRecord
  

  ! ***

  ! Return a field given parameter name, time, etc.
  ! Only one of grid types is filled!

  subroutine mf_ReadRecord_1( mf, paramkey, t1, t2, nuv, nw, &
                                gridtype, levi, &
                                lli, ll, sp_ll, &
                                ggi, gg, sp_gg, &
                                shi, sh, lnsp_sh, &
                                tmi, status )

    use PArray    , only : pa_Init, pa_Done, pa_SetShape
    use file_hdf  , only : THdfFile, TSds, Init, Done, ReadAttribute, ReadData, GetInfo
    use binas     , only : grav
    use GO        , only : TDate, NewDate, Get, rTotal, IncrDate
    use GO        , only : operator(+), operator(-), operator(/=), operator(<)
    use Grid      , only : TLevelInfo, Init, Done
    use Grid      , only : TllGridInfo, TggGridInfo, TshGridInfo
    use Grid      , only : Interpol
    use tmm_info  , only : TMeteoInfo, Init, AddHistory

    ! --- in/out -------------------------------

    type(TMeteoFile_ncep_cdc), intent(inout)  ::  mf
    character(len=*), intent(in)              ::  paramkey
    type(TDate), intent(in)                   ::  t1, t2
    character(len=1), intent(in)              ::  nuv
    character(len=1), intent(in)              ::  nw
    character(len=2), intent(out)             ::  gridtype
    type(TLevelInfo), intent(out)             ::  levi
    type(TllGridInfo), intent(inout)          ::  lli
    real, pointer                             ::  ll(:,:,:)
    real, pointer                             ::  sp_ll(:,:)
    type(TggGridInfo), intent(inout)          ::  ggi
    real, pointer                             ::  gg(:,:)
    real, pointer                             ::  sp_gg(:)
    type(TshGridInfo), intent(inout)          ::  shi
    complex, pointer                          ::  sh(:,:)
    complex, pointer                          ::  lnsp_sh(:)
    type(TMeteoInfo), intent(out)             ::  tmi
    integer, intent(out)                      ::  status

    ! --- const --------------------------------------

    character(len=*), parameter  ::  rname = mname//'/mf_ReadRecord_1'

    ! --- local -------------------------------

    logical               ::  read_2d, read_3d

    logical               ::  constant
    type(TDate)           ::  tbase
    integer               ::  year, hours
    integer               ::  tstart
    real                  ::  tfactor
    
    integer               ::  nlev, ilev
    integer               ::  data_dims(1)
    real, allocatable     ::  levels(:)
    character(len=16)     ::  levtype

    integer               ::  shT
    real, pointer         ::  gg1(:)
    complex, pointer      ::  sh1(:)

    character(len=256)    ::  fname
    character(len=64)     ::  mdir
    character(len=16)     ::  fkey
    character(len=16)     ::  gridkey
    logical               ::  exist
    type(THdfFile)        ::  hdf
    type(TSds)            ::  sds
  
    character(len=16)     ::  sds_name

    character(len=64)     ::  key

    ! --- begin ---------------------------------

    ! no fluxes through boundaries ...
    if ( nuv /= 'n' ) then
      write (gol,'("unsupported nuv key : ",a)') nuv; call goErr
      write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if
    
    ! init info; example of history:  model==msc;sh==159;nlev==60
    call Init( tmi, paramkey, 'unknown', status )
    call AddHistory( tmi, 'model==ncep/ncar reanalysis 1', status )
    call AddHistory( tmi, 'archive==CDC netcdf archive', status )
    
    
    !
    ! ~~~ setup
    !
    
    read_2d  = .false.
    read_3d  = .false.
    constant = .false.

    ! set mf_filekey to first part of file name:
    select case ( paramkey )
      !
      case ( 'LNSP'   ) ; read_2d=.true.; mdir='spectral'     ; fkey='pres.nlog.sfc'; sds_name='pres'
      case ( 'VO'     ) ; read_3d=.true.; mdir='spectral'     ; fkey='vort'         ; sds_name='vort'
      case ( 'D'      ) ; read_3d=.true.; mdir='spectral'     ; fkey='div'          ; sds_name='div'
      case ( 'Tv'     ) ; read_3d=.true.; mdir='spectral'     ; fkey='vair'         ; sds_name='vair'
      case ( 'Q'      ) ; read_3d=.true.; mdir='spectral'     ; fkey='shum'         ; sds_name='shum'
      case ( 'W'      ) ; read_3d=.true.; mdir='pressure'     ; fkey='omega'        ; sds_name='omega'
      !
      case ( 'oro'    ) ; read_2d=.true.; mdir='surface_gauss'; fkey='hgt.sfc'      ; sds_name='hgt'   ; constant=.true.
      case ( 'lsm'    ) ; read_2d=.true.; mdir='surface_gauss'; fkey='land.sfc'     ; sds_name='land'  ; constant=.true.
      case ( 'sr'     ) ; read_2d=.true.; mdir='surface_gauss'; fkey='sfcr.sfc'     ; sds_name='sfcr' 
      case ( 'sps'    ) ; read_2d=.true.; mdir='surface_gauss'; fkey='pres.sfc'     ; sds_name='pres' 
      case ( 'ci'     ) ; read_2d=.true.; mdir='surface_gauss'; fkey='icec.sfc'     ; sds_name='icec' 
      case ( 'skt'    ) ; read_2d=.true.; mdir='surface_gauss'; fkey='skt.sfc'      ; sds_name='skt' 
      case ( 'u10m'   ) ; read_2d=.true.; mdir='surface_gauss'; fkey='uwnd.10m'     ; sds_name='uwnd' 
      case ( 'v10m'   ) ; read_2d=.true.; mdir='surface_gauss'; fkey='vwnd.10m'     ; sds_name='vwnd' 
      case ( 'slhf'   ) ; read_2d=.true.; mdir='surface_gauss'; fkey='lhtfl.sfc'    ; sds_name='lhtfl'
      case ( 'sshf'   ) ; read_2d=.true.; mdir='surface_gauss'; fkey='shtfl.sfc'    ; sds_name='shtfl'
      case ( 'ewss'   ) ; read_2d=.true.; mdir='surface_gauss'; fkey='ugwd.sfc'     ; sds_name='ugwd' 
      case ( 'nsss'   ) ; read_2d=.true.; mdir='surface_gauss'; fkey='vgwd.sfc'     ; sds_name='vgwd' 
      case ( 'lsp'    ) ; read_2d=.true.; mdir='surface_gauss'; fkey='prate.sfc'    ; sds_name='prate'
      case ( 'cp'     ) ; read_2d=.true.; mdir='surface_gauss'; fkey='cprat.sfc'    ; sds_name='cprat'
      case ( 't2m'    ) ; read_2d=.true.; mdir='surface_gauss'; fkey='air.2m'       ; sds_name='air'
      case ( 'ssr'    ) ; read_2d=.true.; mdir='surface_gauss'; fkey='dswrf.sfc'    ; sds_name='dswrf'
      case ( 'sd'     ) ; read_2d=.true.; mdir='surface_gauss'; fkey='weasd.sfc'    ; sds_name='weasd'
      case ( 'swvl1'  ) ; read_2d=.true.; mdir='surface_gauss'; fkey='soilw.0-10cm' ; sds_name='soilw'
      !
      case default
        write (gol,'("unsupported paramkey `",a,"` for ncep 2d files")') paramkey; call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return
    end select

    ! time index: number of 6 hour intervals from start of year.
    if ( constant ) then
      ! index 0 in data sets ...
      tstart = 0
    else
      ! base date:
      tbase = NewDate( year=mf%ccyy, month=1, day=1 )
      ! last 6 hours of previous year ?
      if ( t2 < tbase ) then
        ! first record covers (-6,0]
        tstart = 0
      else
        ! number of hours:
        hours = nint(rTotal( t2 - tbase, 'hour' ))
        ! index = 0, 1, 2, ...
        tstart = ceiling( real(hours)/real(6.0) )
      end if
    end if


    !
    ! ~~~ 2d field
    !
    
    if ( read_2d ) then
        
      ! example file names:
      !   pres.nlog.sfc.spec.2000.nc

      select case ( mdir )
        case ( 'surface_gauss' ) ; gridtype = 'gg' ; gridkey = '.gauss'
        case ( 'spectral'      ) ; gridtype = 'sh' ; gridkey = '.spec'
        case default
          write (gol,'("do not know gridtype for mdir : ",a)') trim(mdir); call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
      end select

      ! write filename
      if ( constant ) then
        write (fname,'(a,a,a,a,a,a,".",a)') &
                  trim(mf%dir), mf%pathsep, trim(mdir), mf%namesep, trim(fkey), trim(gridkey), 'nc'
      else
        write (fname,'(a,a,a,a,a,a,".",i4.4,".",a)') &
                  trim(mf%dir), mf%pathsep, trim(mdir), mf%namesep, trim(fkey), trim(gridkey), mf%ccyy, 'nc'
      end if

      ! file exist ?
      inquire( file=fname, exist=exist )
      if ( .not. exist ) then
        write (gol,'("file not found:")'); call goErr
        write (gol,'("  ",a)') trim(fname); call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return
      end if

      ! open file:
      call Init( hdf, trim(fname), 'read', status )
      IF_NOTOK_RETURN(status=1)
      
      ! check time
      if ( constant ) then
        tfactor = 1.0
      else
        call Check_Time( hdf, tstart, t1, t2, tfactor, status )
        IF_NOTOK_RETURN(status=1)
      end if

      ! read 2d field: fill ggi/gg, or shi/sh :
      select case ( gridtype )
        case ( 'gg' )
          call pa_Init( gg1 )
          call Read_Gaussian_2d( hdf, sds_name, tstart, ggi, gg1, status )
          IF_NOTOK_RETURN(status=1)
          call pa_SetShape( gg, ggi%np, 1 )
          gg(:,1) = gg1
          call pa_Done( gg1 )
          ! apply time factor:
          gg = gg * tfactor
        case ( 'sh' )
          call pa_Init( sh1 )
          call Read_Spectral_2d( hdf, sds_name, tstart, shi, sh1, status )
          IF_NOTOK_RETURN(status=1)
          call pa_SetShape( sh, shi%np, 1 )
          sh(:,1) = sh1
          call pa_Done( sh1 )
          ! apply time factor:
          sh = sh * tfactor
        case default
          write (gol,'("unsupported grid type : ",a)') gridtype; call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
      end select

      ! info ...
      write (key,'("sh==",i4.4)') shi%T
      call AddHistory( tmi, trim(key), status )

      ! close file:
      call Done( hdf, status )
      IF_NOTOK_RETURN(status=1)

      ! unit conversion:
      select case ( paramkey )
        !
        ! factor for conversion from cbar to Pa :
        !   [Pa] = [cbar] * 1e-2 [bar/cbar] * 1e5 [Pa/bar] = [mbar] * 1e3
        ! add to first complex coeff:
        !   sp * fac  =  exp( lnsp + nlog(fac) )
        !             =  exp( {sum_i=1,n c_i p_i} + nlog(fac) )
        !             =  exp( c_1 + {sum_i=2,n c_i p_i} + nlog(fac) )
        case ( 'LNSP' )
          sh(1,1) = sh(1,1) + cmplx(log(1.0e3),0.0)
        !
        ! ncep oro in in [m], should be [m][g] = [m m/s2]
        case ( 'oro' )
          gg(:,1) = gg(:,1) * grav
        !
        ! for some probably historical reaseon, TM expects land/sea mask in % ...
        case ( 'lsm' )
          gg(:,1) = gg(:,1) * 100.0     ! 0-1  ->  %
        !
        ! adhoc surface roughness
        case ( 'sr' )
          write (gol,'("WARNING - adhoc constant value for surface roughness ...")'); call goPr
          gg(:,1) = 0.001   ! m
        !
        ! fluxes downward (ecmwf direction) instead of upward (ncep direction)
        case ( 'slhf', 'sshf' )
          gg(:,1) = - gg(:,1)
        !
        ! kg water / m2  / s   ->   m water / s
        ! With density of 998 kg water / m3  :    kg/m2/s / (kg/m3) = m/s
        case ( 'lsp', 'cp' )
          gg(:,1) = gg(:,1) / 998.0     ! m water / s
      end select

      ! dummy levels
      call Init( levi, 1, (/0.0,0.0/), (/1.0,0.0/), status )
      IF_NOTOK_RETURN(status=1)

    end if   ! read 2d ?


    !
    ! ~~~ 3d field
    !
    
    if ( read_3d ) then
        
      ! example file names:
      !
      !           div.spec.2000.nc
      ! pres.nlog.sfc.spec.2000.nc
      !          vort.spec.2000.nc

      select case ( mdir )
        case ( 'spectral' ) ; gridtype = 'sh' ; gridkey = '.spec' ; levtype = 'sigma'
        case ( 'pressure' ) ; gridtype = 'll' ; gridkey = ''      ; levtype = 'pressure'
        case default
          write (gol,'("do not know grid- and levtype for mdir : ",a)') trim(mdir); call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
      end select

      ! write filename
      write (fname,'(a,a,a,a,a,a,".",i4.4,".",a)') &
                trim(mf%dir), mf%pathsep, trim(mdir), mf%namesep, trim(fkey), trim(gridkey), mf%ccyy, 'nc'

      ! file exist ?
      inquire( file=fname, exist=exist )
      if ( .not. exist ) then
        write (gol,'("file not found:")'); call goErr
        write (gol,'("  ",a)') trim(fname); call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return
      end if

      ! open file:
      call Init( hdf, trim(fname), 'read', status )
      IF_NOTOK_RETURN(status=1)
      
      ! check time
      if ( .not. constant ) then
        call Check_Time( hdf, tstart, t1, t2, tfactor, status )
        IF_NOTOK_RETURN(status=1)
      end if

      ! extract level stuff:
      select case ( levtype )
        !
        case ( 'sigma' )
          ! open level data set:
          call Init( sds, hdf, 'level', status )
          IF_NOTOK_RETURN(status=1)
          ! extract dimensions:
          call GetInfo( sds, status, data_dims=data_dims )
          IF_NOTOK_RETURN(status=1)
          nlev = data_dims(1)
          !call ReadData( sds, sigma, status )
          !IF_NOTOK_RETURN(status=1)
          call Done( sds, status )
          IF_NOTOK_RETURN(status=1)
          ! level defintion
          select case ( nlev )
            case ( 28 )
              call Init( levi, 'nc28', status )
              IF_NOTOK_RETURN(status=1)
            case default
              write (gol,'("level definition not supported for nlev ",i4)') nlev; call goErr
              write (gol,'("in ",a)') rname; call goErr; status=1; return
          end select
        !
        case ( 'pressure' )
          ! open level data set:
          call Init( sds, hdf, 'level', status )
          IF_NOTOK_RETURN(status=1)
          ! extract dimensions:
          call GetInfo( sds, status, data_dims=data_dims )
          IF_NOTOK_RETURN(status=1)
          nlev = data_dims(1)
          ! read pressure levels:
          allocate( levels(nlev) )
          call ReadData( sds, levels, status )
          IF_NOTOK_RETURN(status=1)
          ! close data set:
          call Done( sds, status )
          IF_NOTOK_RETURN(status=1)
          ! level defintion
          write (gol,'("WARNING - adhoc implementation of pressure levels!")'); call goPr
          call Init( levi, nlev, (/levels,0.0/), (/levels*0.0,0.0/), status )
          IF_NOTOK_RETURN(status=1)
          ! clear
          deallocate( levels )
        !
        case default
          write (gol,'("level type not supported :",a)') levtype; call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
      end select
      
      ! info ...
      write (key,'("nlev==",i3.3)') nlev
      call AddHistory( tmi, trim(key), status )

      ! read 3d field: fill lli/ll, ggi/gg, or shi/sh
      select case ( gridtype )
        case ( 'll' )
          call Read_LonLat_3d( hdf, sds_name, tstart, nlev, nw, lli, ll, status )
          IF_NOTOK_RETURN(status=1)
          ! apply time factor:
          ll = ll * tfactor
        case ( 'sh' )
          call Read_Spectral_3d( hdf, sds_name, tstart, nlev, shi, sh, status )
          IF_NOTOK_RETURN(status=1)
          ! apply time factor:
          sh = sh * tfactor
        case default
          write (gol,'("unsupported grid type : ",a)') gridtype; call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
      end select

      ! info ...
      write (key,'("sh==",i4.4)') shT
      call AddHistory( tmi, trim(key), status )

      ! close file:
      call Done( hdf, status )
      IF_NOTOK_RETURN(status=1)

      ! unit conversion:
      select case ( paramkey )
        ! For some reason, the u/v/w from VO/D needs a factor -1 ...
        ! The minus is probably caused by the upwards coordinate system of ncep
        ! instead of the downward from ecmwf.
        case ( 'VO', 'D' )
          sh = - sh
      end select

      !
      ! ~~~ surface pressure
      !

      ! name of ncep file and data set:
      mdir = 'spectral'
      fkey = 'pres.nlog.sfc' ;  sds_name = 'pres'

      ! write surface pressure filename
      write (fname,'(a,"/",a,"-",a,a,".",i4.4,".",a)') &
                trim(mf%dir), trim(mdir), trim(fkey), '.spec', mf%ccyy, 'nc'

      ! file exist ?
      inquire( file=fname, exist=exist )
      if ( .not. exist ) then
        write (gol,'("file not found:")'); call goErr
        write (gol,'("  ",a)') trim(fname); call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return
      end if

      ! open file:
      call Init( hdf, trim(fname), 'read', status )
      IF_NOTOK_RETURN(status=1)

      ! check time
      call Check_Time( hdf, tstart, t1, t2, tfactor, status )
      IF_NOTOK_RETURN(status=1)

      ! read spectral field
      call Read_Spectral_2d( hdf, sds_name, tstart, shi, lnsp_sh, status )
      IF_NOTOK_RETURN(status=1)

      ! lnsp never accumulated, thus no time factor ...

      ! close file:
      call Done( hdf, status )
      IF_NOTOK_RETURN(status=1)

      ! unit conversion:
      !   sp * fac  =  exp( lnsp + nlog(fac) )
      !             =  exp( {sum_i=1,n c_i p_i} + nlog(fac) )
      !             =  exp( c_1 + {sum_i=2,n c_i p_i} + nlog(fac) )
      ! factor for conversion from cbar to Pa :
      !   [Pa] = [cbar] * 1e-2 [bar/cbar] * 1e5 [Pa/bar] = [mbar] * 1e3
      lnsp_sh(1) = lnsp_sh(1) + cmplx(log(1.0e3),0.0)
      
    end if   ! read 3d ?
    

    !
    ! ~~~ end
    !

    ! ok
    status = 0

  end subroutine mf_ReadRecord_1


  ! ***


  subroutine Check_Time( hdf, tstart, t1, t2, tfactor, status )
  
    use file_hdf  , only : THdfFile, TSds, Init, Done, ReadAttribute, ReadData
    use GO        , only : TDate, NewDate, IncrDate, wrtgol
    use GO        , only : operator(+), operator(-), rTotal
    use GO        , only :operator(/=), operator(==), operator(>=), operator(<=)

    ! --- in/out ---------------------------------------------
    
    type(THdfFile), intent(inout)    ::  hdf
    integer, intent(in)              ::  tstart
    type(TDate), intent(in)          ::  t1, t2
    real, intent(out)                ::  tfactor
    integer, intent(out)             ::  status
    
    ! --- const --------------------------------------

    character(len=*), parameter  ::  rname = mname//'/Check_Time'

    ! --- local ----------------------------------------------
    
    type(TSds)            ::  sds

    integer               ::  base_date(3)
    real(8)               ::  time1(1), time(1), dhour
    integer               ::  idh
    type(TDate)           ::  tbase
    type(TDate)           ::  trec
    character(len=32)     ::  avg_period
    integer               ::  avg_dhour
    
    ! --- begin ----------------------------------------------
          
    ! o read start date:
    call ReadAttribute( hdf, 'base_date', base_date, status )
    IF_NOTOK_RETURN(status=1)
    tbase = NewDate( year=base_date(1), month=base_date(2), day=base_date(3) )

    ! init time dataset:
    call Init( sds, hdf, 'time', status )
    IF_NOTOK_RETURN(status=1)


    ! read hour offsets:
    call ReadData( sds, time1, status, start=(/0/) )
    if (status/=0) then
      write (gol,'("reading time value for start=0")'); call goErr
      write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if
    call ReadData( sds, time, status, start=(/tstart/) )
    if (status/=0) then
      write (gol,'("reading time value for start=",i6)') tstart; call goErr
      call wrtgol( '  base_date   : ', tbase ); call goErr
      call wrtgol( '  t1          : ', t1 ); call goErr
      call wrtgol( '  t2          : ', t2 ); call goErr
      write (gol,'("  file name   : ",a)') trim(hdf%fname); call goErr
      write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if

    ! record time :    
    dhour = time(1) - time1(1)
    trec = tbase + IncrDate( hour=floor(dhour), min=int((dhour-floor(dhour))*60.0) )

    ! by default, no time factor needs to applied:
    tfactor = 1.0
    
    ! instant field or time average ?
    if ( t1 == t2 ) then
    
      ! instant time; check:  trec  ==  t1==t2  ?
      if ( trec /= t1 ) then
        write (gol,'("time of supposed record does not match with requested time:")'); call goErr
        write (gol,'("  index        : ",i6)') tstart; call goErr
        write (gol,'("  base_time    : ",i4,2i3.2)') base_date; call goErr
        call wrtgol( '  record time  : ', trec ); call goErr
        call wrtgol( '  t1           : ', t1 ); call goErr
        call wrtgol( '  t2           : ', t2 ); call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return
      end if
      
    else if ( t1 <= t2 ) then
    
      ! interval; check :   [t1,t2]  in  [trec-6,trec]  ?
    
      ! identify averaging period:
      call ReadAttribute( sds, 'avg_period', avg_period, status )
      if (status/=0) then
        write (gol,'("reading avg_period; time inteval requested while not time average fields ?")'); call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return
      end if
      select case ( avg_period )
        case ( '0000-00-00 06:00:00' ) ; avg_dhour = 6
        case default
          write (gol,'("unsupported avg_period : ",a)') trim(avg_period); call goErr
          write (gol,'("in ",a)') rname; call goErr; status=1; return
      end select
      
      ! check:
      if ( (t1 >= trec-IncrDate(hour=avg_dhour)) .and. (t2 <= trec) ) then
         
         ! ok, compute fraction:
         idh = nint(rTotal( t2 - t1, 'hour' ))
         tfactor = real(idh)/real(avg_dhour)
         
      else
      
        write (gol,'("time of supposed record does not match with requested time:")'); call goErr
        write (gol,'("  index        : ",i6)') tstart; call goErr
        write (gol,'("  base_time    : ",i4,2i3.2)') base_date; call goErr
        call wrtgol( '  record time  : ', trec ); call goErr
        write (gol,'("  avg_period   : ",a)') trim(avg_period); call goErr
        call wrtgol( '  t1           : ', t1 ); call goErr
        call wrtgol( '  t2           : ', t2 ); call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return

      end if
      
    else
      
      write (gol,'("times should be the same or interval [t1,t2] : ")'); call goErr
      call wrtgol( '  t1  : ', t1 ); call goErr
      call wrtgol( '  t2  : ', t2 ); call goErr
      write (gol,'("in ",a)') rname; call goErr; status=1; return

    end if  ! instant or interval

    ! ok
    status = 0
    
  end subroutine Check_Time
  

  ! ***
  

  subroutine Read_Spectral_2d( hdf, sds_name, tstart, shi, sh, status )
  
    use parray  , only : pa_SetShape 
    use file_hdf, only : THdfFile, TSds, Init, Done, ReadAttribute, ReadData
    use Grid    , only : TshGridInfo, Init

    ! --- in/out -------------------------------------
    
    type(THdfFile), intent(inout)     ::  hdf
    character(len=*), intent(in)      ::  sds_name
    integer, intent(in)               ::  tstart
    type(TShGridInfo),  intent(out)   ::  shi
    complex, pointer                  ::  sh(:)
    integer, intent(out)              ::  status
    
    ! --- const --------------------------------------

    character(len=*), parameter  ::  rname = mname//'/Read_Spectral_2d'

    ! --- local -------------------------------
    
    type(TSds)            ::  sds
    character(len=1)      ::  compression_used
    real                  ::  meanr(2,1)
    real                  ::  add_offset(1)
    real                  ::  scale_factor(1)
    integer               ::  trunc_count
    
    ! some cdc nc files contain a float for the packed data ...
    real, allocatable     ::  idata(:,:)

    ! --- begin ---------------------------------
    
    ! mean = first complex coeff
    call Init( sds, hdf, 'mean', status )
    IF_NOTOK_RETURN(status=1)
    call ReadData( sds, meanr, status, start=(/0,tstart/) )
    IF_NOTOK_RETURN(status=1)
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)

    ! packed ?
    call ReadAttribute( hdf, 'compression_used', compression_used, status )
    IF_NOTOK_RETURN(status=1)

    ! only packed yet ...
    if ( scan(compression_used,'Tt') == 0 ) then
      write (gol,'("only packed ncep data supported yet ...")'); call goErr
      write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if

    ! read packing offset:
    call Init( sds, hdf, 'add_offset', status )
    IF_NOTOK_RETURN(status=1)
    call ReadData( sds, add_offset, status, start=(/tstart/) )
    IF_NOTOK_RETURN(status=1)
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)

    ! read packing factor:
    call Init( sds, hdf, 'scale_factor', status )
    IF_NOTOK_RETURN(status=1)
    call ReadData( sds, scale_factor, status, start=(/tstart/) )
    IF_NOTOK_RETURN(status=1)
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)

    ! open actual data set:
    call Init( sds, hdf, sds_name, status )
    IF_NOTOK_RETURN(status=1)
    
    ! read spectral truncation:
    call ReadAttribute( sds, 'trunc_count', trunc_count, status )
    IF_NOTOK_RETURN(status=1)

    ! setup output spectral definition:
    call Init( shi, trunc_count, status )    
    IF_NOTOK_RETURN(status=1)
    
    ! data array:
    allocate( idata(shi%np*2,1) )
    
    ! read data:
    call ReadData( sds, idata, status, start=(/0,tstart/) )
    IF_NOTOK_RETURN(status=1)

    ! close data set:
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)

    ! setup output grid:
    call pa_SetShape( sh, shi%np )
    
    ! unpack and transform from realreal to complex:
    sh = transfer( ( idata(:,1) * scale_factor(1) ) + add_offset(1) , sh )

    ! add mean to first coeff:
    select case ( sds_name )
      case ( 'orog' )
        ! for some reason, adding the mean messes gives a bias in orography ..
      case default
        sh(1) = sh(1) + cmplx(meanr(1,1),meanr(2,1))
    end select

    ! convert from NCEP spectral coeff to ECMWF spectral coeff:
    sh = sh / sqrt(2.0)
    
    ! clear
    deallocate(  idata )

    ! ok
    status = 0

  end subroutine Read_Spectral_2d


  ! ***
  

  subroutine Read_Spectral_3d( hdf, sds_name, tstart, nlev, shi, sh, status )

    use parray  , only : pa_SetShape  
    use file_hdf, only : THdfFile, TSds, Init, Done, ReadAttribute, ReadData
    use grid    , only : TshGridInfo, Init

    ! --- in/out -------------------------------------
    
    type(THdfFile), intent(inout)     ::  hdf
    character(len=*), intent(in)      ::  sds_name
    integer, intent(in)               ::  tstart
    integer, intent(in)               ::  nlev
    type(TShGridInfo),  intent(out)   ::  shi
    complex, pointer                  ::  sh(:,:)
    integer, intent(out)              ::  status
    
    ! --- const --------------------------------------

    character(len=*), parameter  ::  rname = mname//'/Read_Spectral_3d'

    ! --- local -------------------------------
    
    integer               ::  ilev
    type(TSds)            ::  sds
    character(len=1)      ::  compression_used
    real                  ::  meanr(nlev*2,1)
    real                  ::  add_offset(nlev,1)
    real                  ::  scale_factor(nlev,1)
    integer, allocatable  ::  idata(:,:,:)
    integer               ::  trunc_count

    ! --- begin ---------------------------------
    
    ! mean = first complex coeff
    call Init( sds, hdf, 'mean', status )
    IF_NOTOK_RETURN(status=1)
    call ReadData( sds, meanr, status, start=(/0,tstart/) )
    IF_NOTOK_RETURN(status=1)
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)

    ! packed ?
    call ReadAttribute( hdf, 'compression_used', compression_used, status )
    IF_NOTOK_RETURN(status=1)

    ! only packed yet ...
    if ( scan(compression_used,'Tt') == 0 ) then
      write (gol,'("only packed ncep data supported yet ...")'); call goErr
      write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if

    ! read packing offset:
    call Init( sds, hdf, 'add_offset', status )
    IF_NOTOK_RETURN(status=1)
    call ReadData( sds, add_offset, status, start=(/0,tstart/) )
    IF_NOTOK_RETURN(status=1)
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)

    ! read packing factor:
    call Init( sds, hdf, 'scale_factor', status )
    IF_NOTOK_RETURN(status=1)
    call ReadData( sds, scale_factor, status, start=(/0,tstart/) )
    IF_NOTOK_RETURN(status=1)
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)

    ! open actual data set:
    call Init( sds, hdf, sds_name, status )
    IF_NOTOK_RETURN(status=1)
    
    ! read spectral truncation:
    call ReadAttribute( sds, 'trunc_count', trunc_count, status )
    IF_NOTOK_RETURN(status=1)

    ! setup output spectral definition:
    call Init( shi, trunc_count, status )    
    IF_NOTOK_RETURN(status=1)
    
    ! data array:
    allocate( idata(shi%np*2,nlev,1) )

    ! read data:
    call ReadData( sds, idata, status, start=(/0,0,tstart/) )
    IF_NOTOK_RETURN(status=1)

    ! close data set:
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)

    ! setup output grid:
    call pa_SetShape( sh, shi%np, nlev )
    
    ! unpack and transform from realreal to complex:
    do ilev = 1, nlev
      sh(:,ilev) = transfer( ( idata(:,ilev,1) * scale_factor(ilev,1) ) + add_offset(ilev,1) , sh(:,ilev) )
    end do

    ! add mean to first coeff:
    do ilev = 1, nlev
      sh(1,ilev) = sh(1,ilev) + cmplx(meanr((ilev-1)*2+1,1),meanr((ilev-1)*2+2,1))
    end do

    ! convert from NCEP spectral coeff to ECMWF spectral coeff:
    sh = sh / sqrt(2.0)

    ! clear
    deallocate( idata        )
    
    ! ok
    status = 0

  end subroutine Read_Spectral_3d
  
  ! ***
  

  subroutine Read_Gaussian_2d( hdf, sds_name, tstart, ggi, gg, status )
  
    use parray  , only : pa_SetShape 
    use file_hdf, only : THdfFile, TSds, Init, Done, ReadAttribute, ReadData, GetInfo
    use Grid    , only : TggGridInfo, Init

    ! --- in/out -------------------------------------
    
    type(THdfFile), intent(inout)     ::  hdf
    character(len=*), intent(in)      ::  sds_name
    integer, intent(in)               ::  tstart
    type(TggGridInfo),  intent(out)   ::  ggi
    real, pointer                     ::  gg(:)
    integer, intent(out)              ::  status
    
    ! --- const --------------------------------------

    character(len=*), parameter  ::  rname = mname//'/Read_Gaussian_2d'

    ! --- local -------------------------------
    
    type(TSds)            ::  sds
    integer               ::  data_dims(3)
    real                  ::  add_offset
    real                  ::  scale_factor
    integer, allocatable  ::  idata(:,:,:)

    ! --- begin ---------------------------------
    
    ! open actual data set:
    call Init( sds, hdf, sds_name, status )
    IF_NOTOK_RETURN(status=1)

    ! extract grid size:    
    call GetInfo( sds, status, data_dims=data_dims )
    IF_NOTOK_RETURN(status=1)
    
    ! setup grid definition:
    call Init( ggi, data_dims(2)/2, .false., status )
    IF_NOTOK_RETURN(status=1)
    
    ! data array:
    allocate( idata(data_dims(1),data_dims(2),1) )
    
    ! read packing:
    call ReadAttribute( sds, 'add_offset', add_offset, status )
    IF_NOTOK_RETURN(status=1)
    call ReadAttribute( sds, 'scale_factor', scale_factor, status )
    IF_NOTOK_RETURN(status=1)
    
    ! read data:
    call ReadData( sds, idata, status, start=(/0,0,tstart/) )
    IF_NOTOK_RETURN(status=1)

    ! close data set:
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)

    ! setup output grid:
    call pa_SetShape( gg, ggi%np )
    
    ! unpack and transform from rank2 to rank1 :
    gg = transfer( ( idata(:,:,1) * scale_factor ) + add_offset , gg )

    ! clear
    deallocate(  idata )

    ! ok
    status = 0

  end subroutine Read_Gaussian_2d


  ! ***
  

  subroutine Read_LonLat_3d( hdf, sds_name, tstart, nlev, nw, lli, ll, status )

    use parray  , only : pa_SetShape  
    use file_hdf, only : THdfFile, TSds, Init, Done, ReadAttribute, ReadData, GetInfo
    use grid    , only : TllGridInfo, Init

    ! --- in/out -------------------------------------
    
    type(THdfFile), intent(inout)     ::  hdf
    character(len=*), intent(in)      ::  sds_name
    integer, intent(in)               ::  tstart
    integer, intent(in)               ::  nlev
    character(len=*), intent(in)      ::  nw
    type(TllGridInfo),  intent(out)   ::  lli
    real, pointer                     ::  ll(:,:,:)
    integer, intent(out)              ::  status
    
    ! --- const --------------------------------------

    character(len=*), parameter  ::  rname = mname//'/Read_LonLat_3d'

    ! --- local -------------------------------
    
    integer               ::  j
    type(TSds)            ::  sds
    integer               ::  data_dims(4)
    real                  ::  add_offset
    real                  ::  scale_factor
    integer, allocatable  ::  idata(:,:,:,:)

    ! --- begin ---------------------------------
    
    ! open data set:
    call Init( sds, hdf, sds_name, status )
    IF_NOTOK_RETURN(status=1)
    
    ! extract grid size:    
    call GetInfo( sds, status, data_dims=data_dims )
    IF_NOTOK_RETURN(status=1)
    
    ! setup grid definition:
    call Init( lli,   0.0, 360.0/ data_dims(1)   , data_dims(1), &
                    -90.0, 180.0/(data_dims(2)-1), data_dims(2), status )
    IF_NOTOK_RETURN(status=1)
    
    ! data array:
    allocate( idata(data_dims(1),data_dims(2),nlev,1) )
    
    ! read packing:
    call ReadAttribute( sds, 'add_offset', add_offset, status )
    IF_NOTOK_RETURN(status=1)
    call ReadAttribute( sds, 'scale_factor', scale_factor, status )
    IF_NOTOK_RETURN(status=1)
    
    ! read data:
    call ReadData( sds, idata, status, start=(/0,0,0,tstart/) )
    IF_NOTOK_RETURN(status=1)

    ! close data set:
    call Done( sds, status )
    IF_NOTOK_RETURN(status=1)

    ! setup output grid:
    select case ( nw )
      case ( 'n' )
        call pa_SetShape( ll, lli%nlon, lli%nlat, nlev )
      case ( 'w' )
        call pa_SetShape( ll, lli%nlon, lli%nlat, nlev+1 )
      case default
        write (gol,'("unsupported nw : ",a)') nw; call goErr
        write (gol,'("in ",a)') rname; call goErr; status=1; return
    end select
    
    ! initial zero:
    ll = 0.0
    
    ! unpack, and transform from north->south to south->north:
    do j = 1, lli%nlat
      ll(:,j,1:nlev) = real( ( idata(:,lli%nlat+1-j,:,1) * scale_factor ) + add_offset )
    end do

    ! clear
    deallocate(  idata )

    ! ok
    status = 0

  end subroutine Read_LonLat_3d
  


end module tmm_mf_ncep_cdc


!program test
!
!  use file_hdf
!  use grid
!
!  ! --- const ---------------------------
!
!  integer, parameter ::  T = 62
!
!  character(len=*), parameter  ::  fdir = '/misc/p71/co2/ncep.reanalysis/spectral/'
!  
!  ! --- local ----------------------------------
!
!  type(THdfFile)  ::  nc, hdf
!  type(TSds)      ::  sds
!  integer         ::  status
!
!  real, allocatable     ::  rr(:)
!  real                  ::  mean(2,1)
!  real                  ::  add_offset(1)
!  real                  ::  scale_factor(1)
!  real, allocatable     ::  idata(:,:)
!
!  type(TshGridInfo)     ::  shi
!  complex, allocatable  ::  cc(:)
!  type(TshGrid)         ::  shc
!
!  integer  ::  tstart
!
!  type(TllGridInfo)  ::  lli
!  real   ::  ll(120,90)
!
!  ! --- begin ------------------------------------
!
!  print *, ''
!  print *, 'test: start'
!
!  call Init( shi, T, status )
!  if (status/=0) stop 'error'
!
!  allocate( cc(shi%np) )
!
!  allocate( idata(shi%np*2,1) )
!  allocate( rr(shi%np*2) )
!
!  !
!  ! time
!  !
!
!  tstart = 0
!
!  !
!  ! read
!  !
!  
!  call Init( nc, fdir//'pres.nlog.sfc.spec.2000.nc', 'read', status )
!  if (status/=0) stop 'error'
!
!  call Init( sds, nc, 'mean', status )
!  if (status/=0) stop 'error'
!  call ReadData( sds, mean, status, start=(/0,tstart/) )
!  if (status/=0) stop 'error'
!  call Done( sds, status )
!  if (status/=0) stop 'error'
!
!  call Init( sds, nc, 'add_offset', status )
!  if (status/=0) stop 'error'
!  call ReadData( sds, add_offset, status, start=(/tstart/) )
!  if (status/=0) stop 'error'
!  call Done( sds, status )
!  if (status/=0) stop 'error'
!
!  call Init( sds, nc, 'scale_factor', status )
!  if (status/=0) stop 'error'
!  call ReadData( sds, scale_factor, status, start=(/tstart/) )
!  if (status/=0) stop 'error'
!  call Done( sds, status )
!  if (status/=0) stop 'error'
!
!  call Init( sds, nc, 'pres', status )
!  if (status/=0) stop 'error'
!  call ReadData( sds, idata, status, start=(/0,tstart/) )
!  if (status/=0) stop 'error'
!  call Done( sds, status )
!  if (status/=0) stop 'error'
!
!  call Done( nc, status )
!  if (status/=0) stop 'error'
!
!  !
!  ! sh grid
!  !
!
!  ! From http://www.cdc.noaa.gov/PublicData/faq.html#12 :
!  !
!  ! Most of the data in our netCDF files are packed. That is to say they have been
!  ! transformed by a scale factor and an add offset to reduce the storage needed to
!  ! two bytes per value. When you extract the short integers, you must unpack the
!  ! data to recover the correct floating point data values. Data files that contain
!  ! packed data will have a non-zero add offset and/or a scale factor not equal to 1.
!  !
!  ! The transformation is:
!  !  float_value = (short_value * scale_factor) + add_offset 
!  !
!
!  print *, mean(1,1)
!  print *, add_offset
!  print *, scale_factor
!  print *, sum(idata)/size(idata)
!
!  rr = ( idata(:,1) * scale_factor(1) ) + add_offset(1)
!
!  print *, minval(rr), maxval(rr)
!  print *, rr(1), exp(rr(1))
!
!  ! convert to complex coeff:
!  cc = transfer( rr, cc )
!
!  ! convert from NOAA spectral coeff to ECMWF spectral coeff:
!  cc = cc / sqrt(2.0)
!
!  !
!  ! convert units from  nlog(bar)  to  nlog(Pa)
!  !  1 bar = 1e5 Pa
!  !  sp_Pa = exp( sp_nlog_cbar ) * 1e5
!  !        = exp( sp_nlog_cbar + nlog(1e5) )
!  !        = exp( sum cnm pnm  + nlog(1e5) )
!  !
!  ! add conversion offset to real part first complex coeff,
!  ! which represent the first global constant mode of lnsp :
!  cc(1) = cc(1) + cmplx(log(1.0e5),0.0)
!
!  call Init( shc )
!  call Set( shc, T, cc )
!
!  !
!  ! interpol
!  !
!
!  call Init( lli, -180.0+3.0/2, 3.0, 120, -90.0+2.0/2, 2.0, 90, status )
!  if (status/=0) stop 'error'
!
!  call Interpol( shc, lli, ll )
!
!  call Done( lli, status )
!  if (status/=0) stop 'error'
!
!  !
!  ! dump
!  !
!
!  call Init( hdf, 'lnsp.hdf', 'create', status )
!  if (status/=0) stop 'error'
!  call Init( sds, hdf, 'lnsp', shape(ll), 'real(4)', status ) 
!  if (status/=0) stop 'error'
!  call WriteData( sds, ll, status )
!  if (status/=0) stop 'error'
!  call Done( sds, status )
!  if (status/=0) stop 'error'
!  call Done( hdf, status )
!  if (status/=0) stop 'error'
!
!  !
!  ! done
!  !
!
!  deallocate( rr )
!  deallocate( idata )
!
!  deallocate( cc )
!  call Done( shc )
!  call Done( shi )
!
!  print *, 'test: end'
!  print *, ''
!
!end program test