!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GETARPS                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE getarps(nx_ext,ny_ext,nz_ext,                                & 1,21
           dir_extd,extdname,extdopt,extdfmt,                           &
           extdinit,extdfcst,julfname,nstyps,                           &
           iproj_ext,scale_ext,                                         &
           trlon_ext,latnot_ext,x0_ext,y0_ext,                          &
           lat_ext,lon_ext,latu_ext,lonu_ext,latv_ext,lonv_ext,         &
           p_ext,hgt_ext,zp_ext,t_ext,qv_ext,  &
           u_ext,vatu_ext,v_ext,uatv_ext,w_ext,  &
           qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                          &
           soiltyp_ext,stypfrct_ext,vegtyp_ext,                         &
           lai_ext,roufns_ext,veg_ext,                                  &
           tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext,         &
           tem1_ext,tem2_ext,istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  ARPS version.
!
!  Reads an ARPS file for processing by ext2arps, a program
!  that converts external files to ARPS variables and format.
!  This version is useful when you want to use an ARPS file
!  with a different orientation or terrain than your final
!  ARPS product so arpsr2h does not work.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster
!  November, 1995
!
!  MODIFICATION HISTORY:
!    01/16/1996 (Yuhe Liu)
!    Added three arguments to specify the "external" ARPS data file
!    names and format.
!
!  2000/08/14 (Gene Bassett)
!  Added multiple soil types, sfcdata variables and grid staggering
!  for use with arps history format of external model data.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    dir_extd      Directory name for external file
!    extdname      Prefix string of external file name
!    extdopt       Option of external data sources
!    extdfmt       Option of external data format
!
!    extdinit      Initialized time in mm-dd-yyyy:hh:mm:ss format
!    extdfcst      Forecast hour in HHH:MM:SS format
!    julfname      File name in yyjjjhhmm format
!
!  OUTPUT:
!
!    iproj_ext     Map projection number of external data
!    scale_ext     Scale factor of external data
!    trlon_ext     True longitude of external data (degrees E)
!    latnot_ext(2) True latitude(s) of external data (degrees N)
!    x0_ext        x coordinate of origin of external data
!    y0_ext        y coordinate of origin of external data
!    lat_ext       latitude of external data points (degrees N)
!    lon_ext       longitude of external data points (degrees E)
!    latu_ext      latitude of external u points (degrees N)
!    lonu_ext      longitude of external u points (degrees E)
!    latv_ext      latitude of external v points (degrees N)
!    lonv_ext      longitude of external v points (degrees E)
!    p_ext         pressure (Pascal)
!    hgt_ext       height (m)
!    zp_ext        height (m) (on arps grid)
!    t_ext         temperature (K)
!    qv_ext        specific humidity (kg/kg)
!    u_ext         u wind component (m/s) (staggered grid, true n-s component)
!    vatu_ext      v wind component at u location (m/s)
!    v_ext         v wind component (m/s) (staggered grid, true e-w component)
!    uatv_ext      u wind component at v location (m/s)
!    qc_ext        Cloud water mixing ratio (kg/kg)
!    qr_ext        Rain water mixing ratio (kg/kg)
!    qi_ext        Ice mixing ratio (kg/kg)
!    qs_ext        Snow mixing ratio (kg/kg)
!    qh_ext        Hail mixing ratio (kg/kg)
!
!    soiltyp_ext   Soil type
!    stypfrct_ext  Soil type fraction
!    vegtyp_ext    Vegetation type
!    lai_ext       Leaf Area Index
!    roufns_ext    Surface roughness
!    veg_ext       Vegetation fraction
!
!    tsfc_ext      Surface temperature
!    tsoil_ext     Soil temperature
!    wetsfc_ext    Top layer soil moisture (fraction)
!    wetdp_ext     Deep soil moisture (fraction)
!    wetcanp_ext   Water content on canopy
!
!    istatus       status indicator
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  CHARACTER (LEN=*) :: dir_extd
  CHARACTER (LEN=*) :: extdname

  INTEGER :: extdopt
  INTEGER :: extdfmt

  INTEGER :: nstyps

  CHARACTER (LEN=19) :: extdinit
  CHARACTER (LEN=9) :: extdfcst
  CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
!  External grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj_ext
  REAL :: scale_ext,trlon_ext
  REAL :: latnot_ext(2)
  REAL :: x0_ext,y0_ext
!
!-----------------------------------------------------------------------
!
!  Output external variable arrays
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx_ext,ny_ext,nz_ext

  REAL :: lat_ext(nx_ext,ny_ext)
  REAL :: lon_ext(nx_ext,ny_ext)
  REAL :: latu_ext(nx_ext,ny_ext)
  REAL :: lonu_ext(nx_ext,ny_ext)
  REAL :: latv_ext(nx_ext,ny_ext)
  REAL :: lonv_ext(nx_ext,ny_ext)
  REAL :: p_ext(nx_ext,ny_ext,nz_ext)     ! Pressure (Pascals)
  REAL :: hgt_ext(nx_ext,ny_ext,nz_ext)   ! Height (m)
  REAL :: zp_ext(nx_ext,ny_ext,nz_ext)    ! Height (m) (on arps grid)
  REAL :: t_ext(nx_ext,ny_ext,nz_ext)     ! Temperature (K)
  REAL :: qv_ext(nx_ext,ny_ext,nz_ext)    ! Specific humidity (kg/kg)
  REAL :: u_ext(nx_ext,ny_ext,nz_ext)     ! Eastward wind component (m/s)
  REAL :: v_ext(nx_ext,ny_ext,nz_ext)     ! Northward wind component (m/s)
  REAL :: uatv_ext(nx_ext,ny_ext,nz_ext)  
  REAL :: vatu_ext(nx_ext,ny_ext,nz_ext)
  REAL :: w_ext(nx_ext,ny_ext,nz_ext)     ! Vertical velocity component (m/s)
  REAL :: qc_ext(nx_ext,ny_ext,nz_ext)    ! Cloud H2O mixing ratio (kg/kg)
  REAL :: qr_ext(nx_ext,ny_ext,nz_ext)    ! Rain  H2O mixing ratio (kg/kg)
  REAL :: qi_ext(nx_ext,ny_ext,nz_ext)    ! Ice   mixing ratio (kg/kg)
  REAL :: qs_ext(nx_ext,ny_ext,nz_ext)    ! Snow  mixing ratio (kg/kg)
  REAL :: qh_ext(nx_ext,ny_ext,nz_ext)    ! Hail  mixing ratio (kg/kg)

  INTEGER soiltyp_ext (nx_ext,ny_ext,nstyps)    ! Soil type
  REAL    stypfrct_ext(nx_ext,ny_ext,nstyps)    ! Soil type
  INTEGER vegtyp_ext  (nx_ext,ny_ext)           ! Vegetation type
  REAL    lai_ext     (nx_ext,ny_ext)           ! Leaf Area Index
  REAL    roufns_ext  (nx_ext,ny_ext)           ! Surface roughness
  REAL    veg_ext     (nx_ext,ny_ext)           ! Vegetation fraction

  REAL :: tsfc_ext   (nx_ext,ny_ext,0:nstyps)  ! Temperature at surface (K)
  REAL :: tsoil_ext  (nx_ext,ny_ext,0:nstyps)  ! Deep soil temperature (K)
  REAL :: wetsfc_ext (nx_ext,ny_ext,0:nstyps)  ! Surface soil moisture (frac)
  REAL :: wetdp_ext  (nx_ext,ny_ext,0:nstyps)  ! Deep soil moisture (fraction)
  REAL :: wetcanp_ext(nx_ext,ny_ext,0:nstyps)  ! Canopy water amount
  REAL :: snowdpth_ext(nx_ext,ny_ext)          ! Snow depth (m)

  REAL :: tem1_ext(nx_ext,ny_ext,nz_ext)   ! Temporary work array
  REAL :: tem2_ext(nx_ext,ny_ext,nz_ext)   ! Temporary work array
!
!-----------------------------------------------------------------------
!
!  Other  external variable arrays
!
!-----------------------------------------------------------------------
!
  REAL :: x_ext(nx_ext)
  REAL :: y_ext(ny_ext)
  REAL :: z_ext(nz_ext)
  REAL :: xsc(nx_ext)
  REAL :: ysc(ny_ext)
!
!  REAL :: uprt_ext(nx_ext,ny_ext,nz_ext)   ! x velocity component (m/s)
!  REAL :: vprt_ext(nx_ext,ny_ext,nz_ext)   ! y velocity component (m/s)
  REAL :: ubar_ext(nx_ext,ny_ext,nz_ext)   ! Base state x velocity
                                           ! component (m/s)
  REAL :: vbar_ext(nx_ext,ny_ext,nz_ext)   ! Base state y velocity
                                           ! component (m/s)
  REAL :: wbar_ext(nx_ext,ny_ext,nz_ext)   ! Base state z velocity
                                           ! component (m/s)
  REAL :: ptbar_ext(nx_ext,ny_ext,nz_ext)  ! Base state potential
                                           ! temperature (K)
  REAL :: pbar_ext(nx_ext,ny_ext,nz_ext)   ! Base state pressure (Pascal)
  REAL :: qvbar_ext(nx_ext,ny_ext,nz_ext)  ! Base state water vapor
                                           ! mixing ratio (kg/kg)

!  real raing_ext  (nx_ext,ny_ext)               ! Grid supersaturation rain
!  real rainc_ext  (nx_ext,ny_ext)               ! Cumulus convective rain
!  real prcrate_ext(nx_ext,ny_ext,4)
                                ! precipitation rate (kg/(m**2*s))
! prcrate(1,1,1) = total precip. rate
! prcrate(1,1,2) = grid scale precip. rate
! prcrate(1,1,3) = cumulus precip. rate
! prcrate(1,1,4) = microphysics precip. rate

!  real radfrc_ext(nx,ny,nz) ! Radiation forcing (K/s)
!  real radsw_ext (nx,ny)    ! Solar radiation reaching the surface
!  real rnflx_ext (nx,ny)    ! Net radiation flux absorbed by surface

!  real usflx_ext (nx,ny)    ! Surface flux of u-momentum (kg/(m*s**2))
!  real vsflx_ext (nx,ny)    ! Surface flux of v-momentum (kg/(m*s**2))
!  real ptsflx_ext(nx,ny)    ! Surface heat flux (K*kg/(m*s**2))
!  real qvsflx_ext(nx,ny)    ! Surface moisture flux (kg/(m**2*s))

  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=3) :: fmtn
  CHARACTER (LEN=80) :: timsnd
  INTEGER :: lenrun, tmstrln
  INTEGER :: i,j,k,ldir,ireturn
  INTEGER :: ihr,imin,isec
  REAL :: govrd
  REAL :: xctr,yctr,dx_ext,dy_ext,tvbot,tvtop,tvbar

  CHARACTER (LEN=80) :: runname_org
  CHARACTER (LEN=80) :: cmnt_org(50)
  INTEGER :: nocmnt_org,mapproj_org
  INTEGER :: month_org,day_org,year_org
  INTEGER :: hour_org,minute_org,second_org

  REAL :: latnot(2)
  REAL :: umove_org,vmove_org,xgrdorg_org,ygrdorg_org
  REAL :: trulat1_org,trulat2_org,trulon_org,sclfct_org
  REAL :: latitud_org,ctrlat_org,ctrlon_org
  REAL :: dx_org,dy_org
  REAL :: lat_org,lon_org
  REAL :: dz_org,dzmin_org,zrefsfc_org,dlayer1_org,dlayer2_org
  REAL :: zflat_org,strhtune_org
  INTEGER :: strhopt_org

  CHARACTER (LEN=80) :: grdbasfn_ext
  CHARACTER (LEN=80) :: datafn_ext
  INTEGER :: nchanl_ext,lengbf_ext
  INTEGER :: lendtf_ext

  REAL :: time_ext
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
  INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Build file names
!
!-----------------------------------------------------------------------
!
  IF ( extdfcst == '         ') extdfcst='000:00:00'

  lenrun=LEN(dir_extd)
  ldir=lenrun
  CALL strlnth( dir_extd, ldir )

  IF ( ldir == 0 .OR. dir_extd(1:ldir) == ' ' ) THEN
    dir_extd = '.'
    ldir = 1
  END IF

  IF( dir_extd(ldir:ldir) /= '/' .AND.  ldir < lenrun ) THEN
    ldir = ldir + 1
    dir_extd(ldir:ldir) = '/'
  END IF

  lenrun = LEN( extdname )
  CALL strlnth( extdname, lenrun )

  IF( extdfmt == 1 ) THEN
    fmtn = 'bin'
  ELSE IF ( extdfmt == 2 ) THEN
    fmtn = 'asc'
  ELSE IF ( extdfmt == 3 ) THEN
    fmtn = 'hdf'
  ELSE IF ( extdfmt == 4 ) THEN
    fmtn = 'pak'
  ELSE IF ( extdfmt == 6 ) THEN
    fmtn = 'bn2'
  ELSE IF ( extdfmt == 7 ) THEN
    fmtn = 'net'
  ELSE IF ( extdfmt == 8 ) THEN
    fmtn = 'npk'
  ELSE IF ( extdfmt == 9 ) THEN
    fmtn = 'gad'
  ELSE IF ( extdfmt == 10 ) THEN
    fmtn = 'grb'
  ELSE
    WRITE(6,'(a,a,a)')                                                  &
        'Unknown format, ', extdfmt, '. Program stopped in GETARPS.'
    STOP
  END IF

  READ(extdfcst,'(i3,1x,i2,1x,i2)') ihr,imin,isec

  time_ext = FLOAT( (ihr*3600)+(imin*60)+isec )
  CALL cvttsnd( time_ext, timsnd, tmstrln )

  grdbasfn_ext = dir_extd(1:ldir)//extdname(1:lenrun)                   &
                                 //'.'//fmtn//'grdbas'
  lengbf_ext = ldir + lenrun + 10

  datafn_ext = dir_extd(1:ldir)//extdname(1:lenrun)                     &
                               //'.'//fmtn//timsnd(1:tmstrln)
  lendtf_ext = ldir + lenrun + 4 + tmstrln

  WRITE(6,*) 'The external grid and base file, grdbasfn = ',            &
             grdbasfn_ext(1:lendtf_ext)

  WRITE(6,*) 'The external time dependent file,  datafn = ',            &
             datafn_ext(1:lengbf_ext)
!
!-----------------------------------------------------------------------
!
!  Since the data reader will change certain parameters stored
!  in common, they need to be saved and restored in common
!  after reading is done.
!
!-----------------------------------------------------------------------
!
  runname_org=runname
  nocmnt_org=nocmnt
  IF( nocmnt > 0 ) THEN
    DO i=1,nocmnt
      cmnt_org(i)=cmnt(i)
    END DO
  END IF
  mapproj_org=mapproj
  month_org=month
  day_org=day
  year_org=year
  hour_org=hour
  minute_org=minute
  second_org=second
!
  umove_org=umove
  vmove_org=vmove
  trulat1_org=trulat1
  trulat2_org=trulat2
  trulon_org=trulon
  sclfct_org=sclfct
  latitud_org=latitud
  ctrlat_org=ctrlat
  ctrlon_org=ctrlon
  dx_org=dx
  dy_org=dy
  xgrdorg_org=xgrdorg
  ygrdorg_org=ygrdorg
  dz_org=dz
  dzmin_org=dzmin
  zrefsfc_org=zrefsfc
  dlayer1_org=dlayer1
  dlayer2_org=dlayer2
  zflat_org=zflat
  strhtune_org=strhtune
  strhopt_org=strhopt
  CALL xytoll(1,1,0.,0.,lat_org,lon_org)
!
!-----------------------------------------------------------------------
!
!  Read ARPS data file
!
!-----------------------------------------------------------------------
!
! uatv_ext & vatu_ext used as temporary variables here
!
  CALL dtaread(nx_ext,ny_ext,nz_ext,nstyps,                             &
               extdfmt, nchanl_ext, grdbasfn_ext,lengbf_ext,            &
               datafn_ext, lendtf_ext, time_ext,                        &
               x_ext,y_ext,z_ext,zp_ext,                                &
               u_ext,v_ext,w_ext,t_ext,p_ext,                     &
               qv_ext, qc_ext, qr_ext,                                  &
               qi_ext, qs_ext, qh_ext, vatu_ext,vatu_ext,vatu_ext,      &
               ubar_ext, vbar_ext, wbar_ext,                            &
               ptbar_ext, pbar_ext, vatu_ext, qvbar_ext,                &
               soiltyp_ext,stypfrct_ext,vegtyp_ext,                     &
               lai_ext,roufns_ext,veg_ext,                              &
               tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext,     &
               snowdpth_ext,                                            &
               tem1_ext(1,1,1),tem1_ext(1,1,2),tem2_ext,                &
               tem2_ext,tem1_ext(1,1,1),tem1_ext(1,1,2),                &
               tem1_ext(1,1,1),tem1_ext(1,1,2),                         &
               tem1_ext(1,1,3),tem1_ext(1,1,4),                         &
               ireturn, tem1_ext, tem2_ext, uatv_ext)
!
!-----------------------------------------------------------------------
!
!  Set maprojection parameters
!
!-----------------------------------------------------------------------
!
  iproj_ext=mapproj
  scale_ext=sclfct
  trlon_ext=trulon
  latnot_ext(1)=trulat1
  latnot_ext(2)=trulat2
  CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
  CALL lltoxy(1,1,ctrlat,ctrlon,xctr,yctr)
  dx_ext=x_ext(2)-x_ext(1)
  x0_ext=xctr - 0.5*(nx_ext-3)*dx_ext
  dy_ext=y_ext(2)-y_ext(1)
  y0_ext=yctr - 0.5*(ny_ext-3)*dy_ext
  CALL setorig(1,x0_ext,y0_ext)
!
!-----------------------------------------------------------------------
!
!  Find lat,lon of scalar points
!
!-----------------------------------------------------------------------
!
  DO i=1,nx_ext-1
    xsc(i)=0.5*(x_ext(i)+x_ext(i+1))
  END DO
  xsc(nx_ext)=2.*xsc(nx_ext-1)-xsc(nx_ext-2)
  DO j=1,ny_ext-1
    ysc(j)=0.5*(y_ext(j)+y_ext(j+1))
  END DO
  ysc(ny_ext)=2.*ysc(ny_ext-1)-ysc(ny_ext-2)

  CALL xytoll(nx_ext,ny_ext,xsc,ysc,lat_ext,lon_ext)
  CALL xytoll(nx_ext,ny_ext,x_ext,ysc,latu_ext,lonu_ext)
  CALL xytoll(nx_ext,ny_ext,xsc,y_ext,latv_ext,lonv_ext)
!
!-----------------------------------------------------------------------
!
!  Move z field onto the scalar grid.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz_ext-1
    DO j=1,ny_ext-1
      DO i=1,nx_ext-1
        hgt_ext(i,j,k)=0.5*(zp_ext(i,j,k)+zp_ext(i,j,k+1))
      END DO
    END DO
  END DO
  DO j=1,ny_ext-1
    DO i=1,nx_ext-1
      hgt_ext(i,j,nz_ext)=(2.*hgt_ext(i,j,nz_ext-1))                    &
                             -hgt_ext(i,j,nz_ext-2)
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Combine wind perturbations and mean fields.
!
!-----------------------------------------------------------------------
!
  u_ext = u_ext + ubar_ext
  v_ext = v_ext + vbar_ext
  w_ext = w_ext + wbar_ext
!
!-----------------------------------------------------------------------
!
!  Orient u & v to true north.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz_ext

    ! get u at v grid point locations
    DO j=2,ny_ext
      DO i=1,nx_ext-1
        uatv_ext(i,j,k) = 0.25*(u_ext(i,j-1,k)+u_ext(i+1,j-1,k)  &
                               +u_ext(i,j,k)+u_ext(i+1,j,k))
      END DO
    END DO
    DO i=1,nx_ext-1
      uatv_ext(i,1,k) = 0.5*(u_ext(i,1,k)+u_ext(i+1,1,k))
    END DO
    DO j=2,ny_ext
      uatv_ext(nx_ext,j,k) = 0.5*(u_ext(nx_ext,j-1,k)+u_ext(nx_ext,j,k))
    END DO
    uatv_ext(nx_ext,1,k) = u_ext(nx_ext,1,k)

    ! get v at u grid point locations
    DO j=1,ny_ext-1
      DO i=2,nx_ext
        vatu_ext(i,j,k) = 0.25*(v_ext(i-1,j,k)+v_ext(i,j,k)  &
                               +v_ext(i-1,j+1,k)+v_ext(i,j+1,k))
      END DO
    END DO
    DO j=1,ny_ext-1
      vatu_ext(1,j,k) = 0.5*(v_ext(1,j,k)+v_ext(1,j+1,k))
    END DO
    DO i=2,nx_ext
      vatu_ext(i,ny_ext,k) = 0.5*(v_ext(i-1,ny_ext,k)+v_ext(i,ny_ext,k))
    END DO
    vatu_ext(1,ny_ext,k) = v_ext(1,ny_ext,k)

    CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),vatu_ext(1,1,k),lonu_ext,     &
       tem1_ext(1,1,k),tem2_ext(1,1,k))
    u_ext(1:nx_ext,1:ny_ext,k) = tem1_ext(1:nx_ext,1:ny_ext,k)
    vatu_ext(1:nx_ext,1:ny_ext,k) = tem2_ext(1:nx_ext,1:ny_ext,k)

    CALL uvmptoe(nx_ext,ny_ext,uatv_ext(1,1,k),v_ext(1,1,k),lonv_ext,     &
       tem1_ext(1,1,k),tem2_ext(1,1,k))
    uatv_ext(1:nx_ext,1:ny_ext,k) = tem1_ext(1:nx_ext,1:ny_ext,k)
    v_ext(1:nx_ext,1:ny_ext,k) = tem2_ext(1:nx_ext,1:ny_ext,k)

  END DO
!
!-----------------------------------------------------------------------
!
!  Combine perturbations and mean fields of scalars
!  Convert potential temperature to potential temperature
!
!-----------------------------------------------------------------------
!
  DO k=2,nz_ext-1
    DO j=2,ny_ext-1
      DO i=2,nx_ext-1
        p_ext(i,j,k)=p_ext(i,j,k)+pbar_ext(i,j,k)
        t_ext(i,j,k)=t_ext(i,j,k)+ptbar_ext(i,j,k)
        t_ext(i,j,k)=t_ext(i,j,k)*((p_ext(i,j,k)/p0)**rddcp)
        qv_ext(i,j,k)=qv_ext(i,j,k)+qvbar_ext(i,j,k)
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Supply data at the edge points (zero gradient, where missing)
!
!-----------------------------------------------------------------------
!
  CALL edgfill(hgt_ext,1,nx_ext,2,nx_ext-1,1,ny_ext,2,ny_ext-1,         &
                       1,nz_ext,1,nz_ext)
  CALL edgfill(p_ext,1,nx_ext,2,nx_ext-1,1,ny_ext,2,ny_ext-1,           &
                     1,nz_ext,1,nz_ext)
  CALL edgfill(t_ext,1,nx_ext,2,nx_ext-1,1,ny_ext,2,ny_ext-1,           &
                     1,nz_ext,1,nz_ext)
  CALL edgfill(qv_ext,1,nx_ext,2,nx_ext-1,1,ny_ext,2,ny_ext-1,          &
                     1,nz_ext,2,nz_ext-1)
! u & v now not on scalar points, so don't edgfill
!  CALL edgfill(u_ext,1,nx_ext,2,nx_ext-1,1,ny_ext,2,ny_ext-1,           &
!                     1,nz_ext,2,nz_ext-1)
!  CALL edgfill(v_ext,1,nx_ext,2,nx_ext-1,1,ny_ext,2,ny_ext-1,           &
!                     1,nz_ext,2,nz_ext-1)
!
!-----------------------------------------------------------------------
!
!  Make top and bottom mass fields via hydrostatic extrapolation.
!
!-----------------------------------------------------------------------
!
  govrd=g/rd
  DO j=1,ny_ext-1
    DO i=1,nx_ext-1
!
      t_ext(i,j,1)=2.*t_ext(i,j,2)-t_ext(i,j,3)
      tvbot=t_ext(i,j,1) * ( 1.0 + 0.608*qv_ext(i,j,1) )
      tvtop=t_ext(i,j,2) * ( 1.0 + 0.608*qv_ext(i,j,2) )
      tvbar=0.5*(tvtop+tvbot)
      p_ext(i,j,1)=p_ext(i,j,2)                                         &
                   *EXP(govrd*(hgt_ext(i,j,2)-hgt_ext(i,j,1))/tvbar)
!
      t_ext(i,j,nz_ext)=2.*t_ext(i,j,nz_ext-1)-t_ext(i,j,nz_ext-2)
      tvbot=t_ext(i,j,nz_ext-1)*(1.0 + 0.608*qv_ext(i,j,nz_ext-1))
      tvtop=t_ext(i,j,nz_ext)  *(1.0 + 0.608*qv_ext(i,j,nz_ext))
      tvbar=0.5*(tvtop+tvbot)
      p_ext(i,j,nz_ext)=p_ext(i,j,nz_ext-1)                             &
                   *EXP(govrd*(hgt_ext(i,j,nz_ext-1)-                   &
                        hgt_ext(i,j,nz_ext))/tvbar)
    END DO
  END DO

  CALL edgfill(p_ext,1,nx_ext,1,nx_ext-1,1,ny_ext,1,ny_ext-1,           &
                     1,nz_ext,1,nz_ext)
  CALL edgfill(t_ext,1,nx_ext,1,nx_ext-1,1,ny_ext,1,ny_ext-1,1          &
                      ,nz_ext,1,nz_ext)
!
!-----------------------------------------------------------------------
!
!  Reset info in common to original values
!
!-----------------------------------------------------------------------
!
  runname=runname_org
  nocmnt=nocmnt_org
  IF( nocmnt > 0 ) THEN
    DO i=1,nocmnt
      cmnt(i)=cmnt_org(i)
    END DO
  END IF
  mapproj=mapproj_org
  month=month_org
  day=day_org
  year=year_org
  hour=hour_org
  minute=minute_org
  second=second_org
!
  umove=umove_org
  vmove=vmove_org
  xgrdorg=xgrdorg_org
  ygrdorg=ygrdorg_org
  trulat1=trulat1_org
  trulat2=trulat2_org
  trulon=trulon_org
  sclfct=sclfct_org
  latitud=latitud_org
  ctrlat=ctrlat_org
  ctrlon=ctrlon_org
  dx=dx_org
  dy=dy_org
  dz=dz_org
  dzmin=dzmin_org
  zrefsfc=zrefsfc_org
  dlayer1=dlayer1_org
  dlayer2=dlayer2_org
  zflat=zflat_org
  strhtune=strhtune_org
  strhopt=strhopt_org
  xgrdorg=xgrdorg_org
  ygrdorg=ygrdorg_org
!
!-----------------------------------------------------------------------
!
!  Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
  latnot(1)=trulat1
  latnot(2)=trulat2
  CALL setmapr(mapproj,sclfct,latnot,trulon)
  CALL setorig(2,lat_org,lon_org)
  istatus=1
  RETURN
!
!-----------------------------------------------------------------------
!
!  Error destination
!
!-----------------------------------------------------------------------
!
  WRITE(6,'(a)')  ' Error reading last field, returning'
  RETURN
END SUBROUTINE getarps
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GETNMCRUC87                ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE getnmcruc87(nx_ext,ny_ext,nz_ext,                            & 1,12
           dir_extd,extdname,extdopt,extdfmt,                           &
           extdinit,extdfcst,julfname,                                  &
           iproj_ext,scale_ext,                                         &
           trlon_ext,latnot_ext,x0_ext,y0_ext,                          &
           lat_ext,lon_ext,                                             &
           p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,                      &
           qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                          &
           tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext,         &
           trn_ext,psfc_ext,                                            &
           istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  ARPS version.
!
!  Reads an ARPS file for processing by ext2arps, a program
!  that converts external files to ARPS variables and format.
!  This version is useful when you want to use an ARPS file
!  with a different orientation or terrain than your final
!  ARPS product so arpsr2h does not work.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu and Jinxing Zong
!  01/18/1996
!
!  MODIFICATION HISTORY:
!    6/5/1997 Jinxing Zong and Yuhe Liu
!    Virtualization of temperature is accounted for in variable
!    conversion. Values of constants cp, rd/cp and g used in RUC are
!    adopted in making the variable conversion. Subroutine TV2TQ and
!    function ESW are added to facilitate the conversion.
!
!    01/26/1998 (Yuhe Liu)
!    Removed function esw, instead, used unified function f_es in
!    file thermolib3d.f
!
!    1999/11/30 (Gene Bassett)
!    Changed deep soil temperature and moisture to be an average from
!    0-100cm.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    dir_extd      Directory name for external file
!    extdname      Prefix string of external file name
!    extdopt       Option of external data sources
!    extdfmt       Option of external data format

!    extdinit      Initialized time in mm-dd-yyyy:hh:mm:ss format
!    extdfcst      Forecast hour in HHH:MM:SS format
!    julfname      File name in yyjjjhhmm format
!
!  OUTPUT:
!
!    iproj_ext     Map projection number of external data
!    scale_ext     Scale factor of external data
!    trlon_ext     True longitude of external data (degrees E)
!    latnot_ext(2) True latitude(s) of external data (degrees N)
!    x0_ext        x coordinate of origin of external data
!    y0_ext        y coordinate of origin of external data
!    lat_ext       latitude of external data points (degrees N)
!    lon_ext       longitude of external data points (degrees E)
!    p_ext         pressure (Pascal)
!    hgt_ext       height (m)
!    t_ext         temperature (K)
!    qv_ext        specific humidity (kg/kg)
!    u_ext         u wind component (m/s)
!    v_ext         v wind component (m/s)
!    qc_ext        Cloud water mixing ratio (kg/kg)
!    qr_ext        Rain water mixing ratio (kg/kg)
!    qi_ext        Ice mixing ratio (kg/kg)
!    qs_ext        Snow mixing ratio (kg/kg)
!    qh_ext        Hail mixing ratio (kg/kg)
!
!    tsfc_ext      Surface temperature
!    tsoil_ext     Soil temperature
!    wetsfc_ext    Top layer soil moisture (fraction)
!    wetdp_ext     Deep soil moisture (fraction)
!    wetcanp_ext   Water content on canopy
!
!    trn_ext       External terrain (m)
!    psfc_ext      Surface pressure (Pa)
!
!    istatus       status indicator
!
!  WORK ARRAYS:
!
!    var_grb3d     Arrays to store the GRIB 3-D variables:
!                  var_grb3d(nxgrb,nygrb,nzgrb,1,1) - pressure (Pa)
!                  var_grb3d(nxgrb,nygrb,nzgrb,2,1) - Montgomery stream
!                                                     function (m**2/s**2)
!                  var_grb3d(nxgrb,nygrb,nzgrb,3,1) - Potential
!                                                     temperature (K)
!                  var_grb3d(nxgrb,nygrb,nzgrb,4,1) - Condensation pressure
!                                                     of lifted parcel (Pa)
!                  var_grb3d(nxgrb,nygrb,nzgrb,5,1) - u wind (m/s)
!                  var_grb3d(nxgrb,nygrb,nzgrb,6,1) - v wind (m/s)
!                  var_grb3d(nxgrb,nygrb,nzgrb,1,2) - soil temp. (K)
!                  var_grb3d(nxgrb,nygrb,nzgrb,2,2) - soil moist.
!                                                     (fraction)
!
!    var_grb2d     Arrays to store the GRIB 2-D variables:
!                  var_grb2d(nxgrb,nygrb,1) - Ground surface temperature (K)
!                  var_grb2d(nxgrb,nygrb,2) - Geometric height (m)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'gribcst.inc'

  CHARACTER (LEN=*) :: dir_extd
  CHARACTER (LEN=*) :: extdname

  INTEGER :: extdopt
  INTEGER :: extdfmt

  CHARACTER (LEN=19) :: extdinit
  CHARACTER (LEN=9) :: extdfcst
  CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
!  External grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj_ext
  REAL :: scale_ext,trlon_ext
  REAL :: latnot_ext(2)
  REAL :: x0_ext,y0_ext
  REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
!  Output external variable arrays
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx_ext, ny_ext, nz_ext

  REAL :: lat_ext(nx_ext,ny_ext)
  REAL :: lon_ext(nx_ext,ny_ext)
  REAL :: p_ext  (nx_ext,ny_ext,nz_ext)   ! Pressure (Pascals)
  REAL :: hgt_ext(nx_ext,ny_ext,nz_ext)   ! Height (m)
  REAL :: t_ext  (nx_ext,ny_ext,nz_ext)   ! Temperature (K)
  REAL :: qv_ext (nx_ext,ny_ext,nz_ext)   ! Specific humidity (kg/kg)
  REAL :: u_ext  (nx_ext,ny_ext,nz_ext)   ! Eastward wind component
  REAL :: v_ext  (nx_ext,ny_ext,nz_ext)   ! Northward wind component
  REAL :: qc_ext (nx_ext,ny_ext,nz_ext)   ! Cloud H2O mixing ratio (kg/kg)
  REAL :: qr_ext (nx_ext,ny_ext,nz_ext)   ! Rain  H2O mixing ratio (kg/kg)
  REAL :: qi_ext (nx_ext,ny_ext,nz_ext)   ! Ice   mixing ratio (kg/kg)
  REAL :: qs_ext (nx_ext,ny_ext,nz_ext)   ! Snow  mixing ratio (kg/kg)
  REAL :: qh_ext (nx_ext,ny_ext,nz_ext)   ! Hail  mixing ratio (kg/kg)

  REAL :: tsfc_ext   (nx_ext,ny_ext)      ! Temperature at surface (K)
  REAL :: tsoil_ext  (nx_ext,ny_ext)      ! Deep soil temperature (K)
  REAL :: wetsfc_ext (nx_ext,ny_ext)      ! Surface soil moisture (fraction)
  REAL :: wetdp_ext  (nx_ext,ny_ext)      ! Deep soil moisture (fraction)
  REAL :: wetcanp_ext(nx_ext,ny_ext)      ! Canopy water amount

  REAL :: trn_ext    (nx_ext,ny_ext)          ! Geometrical heights
  REAL :: psfc_ext   (nx_ext,ny_ext)      ! Surface pressure (Pa)

!
!-----------------------------------------------------------------------
!
!  Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: var_grb2d(:,:,:,:)   ! GRIB variables
  REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
  INTEGER, ALLOCATABLE :: var_lev3d(:,:,:)  ! Levels (hybrid) for
                                            ! each 3-D variable
  REAL, ALLOCATABLE :: rcdata(:)            ! temporary data array
!
!-----------------------------------------------------------------------
!
!  Other  external variable arrays
!
!-----------------------------------------------------------------------
!
  REAL :: x_ext(nx_ext)
  REAL :: y_ext(ny_ext)

  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Original grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj
  REAL :: scale,trlon,x0,y0
  REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: gribfile
  CHARACTER (LEN=13) :: gribtime
  INTEGER :: i,j,k
  INTEGER :: iyr,imo,iday,myr, jldy
  INTEGER :: ihr,imin,isec
  INTEGER :: ifhr,ifmin,ifsec
  INTEGER :: grbflen, len1, lenrun

  INTEGER :: m,n,nz1,max_nr2d,max_nr3d

  REAL :: pilev

  REAL :: tv_ext, tvc_ext
  REAL :: rovcp_p, cpd_p, g0_p, rd_p

  INTEGER :: chklev, lvscan, kk, jj

  REAL :: tema, temb

  INTEGER :: iret             ! Return flag

  REAL, ALLOCATABLE :: utmp(:,:), vtmp(:,:)
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  GRIB grid information
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=42) :: gridesc ! Grid description

  INTEGER :: iproj_grb    ! Map projection indicator
  INTEGER :: gthin        ! Indicator of whether the grid is "thinned"

  INTEGER :: ni_grb       ! Number of points along x-axis
  INTEGER :: nj_grb       ! Number of points along y-axis
  INTEGER :: np_grb       ! Total number of horizontal grid points

  INTEGER :: nk_grb       ! Number of vertical parameters
  REAL :: zk_grb(nz_ext)   ! Vertical coordinate parameters

  INTEGER :: npeq         ! Number of lat circles from pole to equator
  INTEGER :: nit(nz_ext)   ! Number of x-points for thinned grid

  REAL :: pi_grb          ! x-coordinate of pole point
  REAL :: pj_grb          ! y-coordinate of pole point
  INTEGER :: ipole        ! Projection center flag

  REAL :: di_grb          ! x-direction increment or grid length
  REAL :: dj_grb          ! y-direction increment or grid length

  REAL :: latsw           ! Latitude  of South West corner point
  REAL :: lonsw           ! Longitude of South West corner point
  REAL :: latne           ! Latitude  of North East corner point
  REAL :: lonne           ! Longitude of North East corner point

  REAL :: lattru1         ! Latitude (1st) at which projection is true
  REAL :: lattru2         ! Latitude (2nd) at which projection is true
  REAL :: lontrue         ! Longitude      at which projection is true

  REAL :: latrot          ! Latitude  of southern pole of rotation
  REAL :: lonrot          ! Longitude of southern pole of rotation
  REAL :: angrot          ! Angle of rotation

  REAL :: latstr          ! Latitude  of the pole of stretching
  REAL :: lonstr          ! Longitude of the pole of stretching
  REAL :: facstr          ! Stretching factor

  INTEGER :: scanmode     ! Scanning indicator
  INTEGER :: iscan        ! x-direction   scanning indicator
  INTEGER :: jscan        ! y-direction   scanning indicator
  INTEGER :: kscan        ! FORTRAN index scanning indicator

  INTEGER :: ires         ! Resolution direction increments indicator
  INTEGER :: iearth       ! Earth shape indicator: spherical or oblate?
  INTEGER :: icomp        ! (u,v) components decomposition indicator

  INTEGER :: jpenta       ! J-Pentagonal resolution parameter
  INTEGER :: kpenta       ! K-Pentagonal resolution parameter
  INTEGER :: mpenta       ! M-Pentagonal resolution parameter
  INTEGER :: ispect       ! Spectral representation type
  INTEGER :: icoeff       ! Spectral coefficient storage mode

  REAL :: xp_grb          ! X coordinate of sub-satellite point
  REAL :: yp_grb          ! Y coordinate of sub-satellite point
  REAL :: xo_grb          ! X coordinate of image sector origin
  REAL :: yo_grb          ! Y coordinate of image sector origin
  REAL :: zo_grb          ! Camera altitude from center of Earth
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
  ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
  ALLOCATE(rcdata(nx_ext*ny_ext))
  ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))
  ALLOCATE(utmp(nx_ext,ny_ext))
  ALLOCATE(vtmp(nx_ext,ny_ext))

  IF ( extdfcst == '         ') extdfcst='000:00:00'

  READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)')                  &
      iyr,imo,iday,ihr,imin,isec

  CALL julday(iyr,imo,iday,jldy)

  myr=MOD(iyr,100)
  ifhr=0
  ifmin=0
  ifsec=0

  READ (extdfcst,'(i3)',ERR=4,END=4) ifhr

  4     CONTINUE

  WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)')                      &
         iyr,imo,iday,ihr,'f',ifhr

  len1=LEN(dir_extd)
  grbflen=len1

  CALL strlnth( dir_extd, grbflen )

  IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
    dir_extd = '.'
    grbflen=1
  END IF

  IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
    grbflen=grbflen+1
    dir_extd(grbflen:grbflen)='/'
  END IF

  lenrun = LEN( extdname )
  CALL strlnth( extdname, lenrun )

  gribfile = dir_extd(1:grbflen)//extdname(1:lenrun)                    &
                                //'.'//gribtime(1:13)
  grbflen = grbflen + lenrun + 14
!
!-----------------------------------------------------------------------
!
!  RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
  gridtyp   = ruc87grid
  mproj_grb = ruc87proj

  n2dvars  = ruc87nvs2d
  n2dlvtps = ruc87nlvt2d

  DO m=1,n2dlvtps
    DO n=1,n2dvars
      var_id2d(n,m) = ruc87var_id2d(n,m)
    END DO
    levtyp2d(m) = ruc87levs2d(m)
  END DO

  n3dvars  = ruc87nvs3d
  n3dlvtps = ruc87nlvt3d

  DO m=1,n3dlvtps
    DO n=1,n3dvars
      var_id3d(n,m) = ruc87var_id3d(n,m)
    END DO
    levtyp3d(m) = ruc87levs3d(m)
  END DO

  CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime,        &
                gridesc, iproj_grb, gthin,                              &
                ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit,          &
                pi_grb,pj_grb,ipole, di_grb,dj_grb,                     &
                latsw,lonsw, latne,lonne,                               &
                latrot,lonrot,angrot,                                   &
                latstr,lonstr,facstr,                                   &
                lattru1,lattru2,lontrue,                                &
                scanmode, iscan,jscan,kscan,                            &
                ires,iearth,icomp,                                      &
                jpenta,kpenta,mpenta,ispect,icoeff,                     &
                xp_grb,yp_grb, xo_grb,yo_grb,zo_grb,                    &
                rcdata,var_grb2d,var_grb3d,var_lev3d,iret)

  max_nr2d = 0
  DO n=1,n2dvars
    DO m=1,n2dlvtps
      max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
    END DO
  END DO

  max_nr3d = 0
  DO n=1,n3dvars
    max_nr3d = MAX( max_nr3d, var_nr3d(n,1) )
  END DO

  IF ( max_nr3d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 3-D variables was found in the GRIB file',                  &
        'Program stopped in GETNMCRUC.'
    STOP
  END IF

  IF ( max_nr2d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 2-D variables was found in the GRIB file'
  END IF

!  write (6,'(/a7,2x,6(i7))')
!    :      'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)

!  DO 60 k=1,max_nr3d
!     write (6,'(i5,4x,6(i7))')
!    :      k,(var_lev3d(k,n,1),n=1,n3dvars)

  DO k=1,max_nr3d
    DO n=2,n3dvars
      IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
        WRITE (6,'(a)')                                                 &
            'Variables were not at the same level.',                    &
            'Program stopped in GETNMCRUC.'
        STOP
      END IF
    END DO
  END DO

  IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN      ! Center in NP
    iproj_ext = 1
  ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -1
  ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN  ! Center in NP
    iproj_ext = 2
  ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -2
  ELSE IF ( iproj_grb == 1 ) THEN
    iproj_ext = 3
  ELSE IF ( iproj_grb == 0 ) THEN
    iproj_ext = 4
  ELSE
    WRITE (6,'(a)')                                                     &
        'Unknown map projection. Set to non-projection.'
    iproj_ext = 0
  END IF

  scale_ext = 1.0
  latnot_ext(1) = lattru1
  latnot_ext(2) = lattru2
  trlon_ext = lontrue

  dx_ext = di_grb
  dy_ext = dj_grb

  CALL getmapr(iproj,scale,latnot,trlon,x0,y0)
  CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
  CALL lltoxy(1,1,latsw,lonsw,x0_ext,y0_ext)

  DO i=1,nx_ext
    x_ext(i)=x0_ext+(i-1)*dx_ext
  END DO

  DO j=1,ny_ext
    y_ext(j)=y0_ext+(j-1)*dy_ext
  END DO

  CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!
!-----------------------------------------------------------------------
!
!  Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
  DO j=1,ny_ext
    DO i=1,nx_ext
      IF ( var_nr2d(1,1) == 0 ) THEN
        tsfc_ext   (i,j) = -999.0
      ELSE
        tsfc_ext   (i,j) = var_grb2d(i,j,1,1)
      END IF

      IF ( var_nr2d(2,1) == 0 ) THEN
        trn_ext    (i,j) = -999.0
      ELSE
        trn_ext    (i,j) = var_grb2d(i,j,2,1)
      END IF

      IF ( var_nr3d(1,2) == 0 ) THEN
        tsfc_ext   (i,j) = var_grb3d(i,j,1,1,2)
             ! note that this is the 5cm value and not the surface value
        tsoil_ext  (i,j) = 0.1 * var_grb3d(i,j,1,1,2)  & !   5cm
                     + 0.2 * var_grb3d(i,j,2,1,2)  & !  20cm
                     + 0.4 * var_grb3d(i,j,3,1,2)  & !  40cm
                     + 0.3 * var_grb3d(i,j,4,1,2)  ! 160cm
        wetsfc_ext (i,j) = var_grb3d(i,j,1,2,2)
             ! note that this is the 5cm value and not the surface value
        wetdp_ext  (i,j) = 0.1 * var_grb3d(i,j,1,2,2)  & !   5cm
                     + 0.2 * var_grb3d(i,j,2,2,2)  & !  20cm
                     + 0.4 * var_grb3d(i,j,3,2,2)  & !  40cm
                     + 0.3 * var_grb3d(i,j,4,2,2)  ! 160cm
        wetcanp_ext(i,j) = var_grb2d(i,j,3,1)*1.e-3  ! in meters
      ELSE
        tsfc_ext   (i,j) = -999.0
        tsoil_ext  (i,j) = -999.0
        wetsfc_ext (i,j) = -999.0
        wetdp_ext  (i,j) = -999.0
        wetcanp_ext(i,j) = -999.0
      END IF

      psfc_ext   (i,j) = -999.0
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!

  cpd_p   = 1004.686    ! cp in RUC
  rovcp_p = 0.285714    ! rd/cp used in RUC
  g0_p    = 9.80665     ! gravity in RUC

  nz1 = MIN(var_nr3d(1,1),nz_ext)

  IF ( var_lev3d(1,1,1) < var_lev3d(nz1,1,1) ) THEN  ! 1st level at sfc
    chklev = 1
    lvscan = 0
  ELSE
    chklev = -1
    lvscan = nz1+1
  END IF

  DO k=1,nz1
    kk = chklev * k + lvscan
    DO j=1,ny_ext
      DO i=1,nx_ext

        p_ext(i,j,kk) = var_grb3d(i,j,k,1,1)       ! Pressure (Pa)

        pilev = (p_ext(i,j,kk)/100000.)**rovcp_p
        tv_ext = var_grb3d(i,j,k,3,1)*pilev        ! Virtual Temperature (K)
        hgt_ext(i,j,kk) = (var_grb3d(i,j,k,2,1)-cpd_p*tv_ext)/g0_p
                                ! Height (m) from Mongmery function with
! M = CpTv + gz

        u_ext(i,j,kk) = var_grb3d(i,j,k,5,1)       ! u wind (m/s)
        v_ext(i,j,kk) = var_grb3d(i,j,k,6,1)       ! v wind (m/s)

        tvc_ext = var_grb3d(i,j,k,3,1)                                  &
                     * (var_grb3d(i,j,k,4,1)/100000.)**rovcp_p
                              ! Virtual Temperature (K) at LCL
        CALL tv2tq( tvc_ext, var_grb3d(i,j,k,4,1), tema, temb )
        t_ext(i,j,kk)  = tema                                           &
            * (p_ext(i,j,kk)/var_grb3d(i,j,k,4,1))**rovcp_p  ! Temperature (K)
        qv_ext(i,j,kk) = temb                      ! Specific humidity
        qc_ext(i,j,kk) = -999.
        qr_ext(i,j,kk) = -999.
        qi_ext(i,j,kk) = -999.
        qs_ext(i,j,kk) = -999.
        qh_ext(i,j,kk) = -999.
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Rotate winds to be relative to true north.
!  The RUC data are sent as grid-relative.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz1
!2001-05-16 GMB: Having umap & uear (or vmap & vear) point to
!the same array causes numerical errors when optimizing.
    CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k),               &
                 lon_ext,utmp,vtmp)
    u_ext(:,:,k) = utmp(:,:)
    v_ext(:,:,k) = vtmp(:,:)
  END DO

!
!-----------------------------------------------------------------------
!
!  Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
  CALL setmapr(iproj,scale,latnot,trlon)
  CALL setorig(1,x0,y0)

  istatus = 1

  DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)
  DEALLOCATE(utmp)
  DEALLOCATE(vtmp)

  RETURN
END SUBROUTINE getnmcruc87
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GETNMCRUC211               ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE getnmcruc211(nx_ext,ny_ext,nz_ext,                           & 1,12
           dir_extd,extdname,extdopt,extdfmt,                           &
           extdinit,extdfcst,julfname,                                  &
           iproj_ext,scale_ext,                                         &
           trlon_ext,latnot_ext,x0_ext,y0_ext,                          &
           lat_ext,lon_ext,                                             &
           p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,                      &
           qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                          &
           tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext,         &
           trn_ext,psfc_ext, t_2m_ext, rh_2m_ext, u_10m_ext,            &
           v_10m_ext, istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  ARPS version.
!
!  Reads an ARPS file for processing by ext2arps, a program
!  that converts external files to ARPS variables and format.
!  This version is useful when you want to use an ARPS file
!  with a different orientation or terrain than your final
!  ARPS product so arpsr2h does not work.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Dennis Moon, SSESCO
!  06/19/1996
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    dir_extd      Directory name for external file
!    extdname      Prefix string of external file name
!    extdopt       Option of external data sources
!    extdfmt       Option of external data format

!    extdinit      Initialized time in mm-dd-yyyy:hh:mm:ss format
!    extdfcst      Forecast hour in HHH:MM:SS format
!    julfname      File name in yyjjjhhmm format
!
!  OUTPUT:
!
!    iproj_ext     Map projection number of external data
!    scale_ext     Scale factor of external data
!    trlon_ext     True longitude of external data (degrees E)
!    latnot_ext(2) True latitude(s) of external data (degrees N)
!    x0_ext        x coordinate of origin of external data
!    y0_ext        y coordinate of origin of external data
!    lat_ext       latitude of external data points (degrees N)
!    lon_ext       longitude of external data points (degrees E)
!    p_ext         pressure (Pascal)
!    hgt_ext       height (m)
!    t_ext         temperature (K)
!    qv_ext        specific humidity (kg/kg)
!    u_ext         u wind component (m/s)
!    v_ext         v wind component (m/s)
!    qc_ext        Cloud water mixing ratio (kg/kg)
!    qr_ext        Rain water mixing ratio (kg/kg)
!    qi_ext        Ice mixing ratio (kg/kg)
!    qs_ext        Snow mixing ratio (kg/kg)
!    qh_ext        Hail mixing ratio (kg/kg)
!
!    tsfc_ext      Surface temperature
!    tsoil_ext     Soil temperature
!    wetsfc_ext    Top layer soil moisture (fraction)
!    wetdp_ext     Deep soil moisture (fraction)
!    wetcanp_ext   Water content on canopy
!
!    trn_ext       External terrain (m)
!    psfc_ext      Surface pressure (Pa)
!
!    T_2m_ext      Temperature at 2m AGL
!    rh_2m_ext     Specific Humidity at 2m AGL
!    U_10m_ext     U at 10m AGL
!    V_10m_ext     V at 10m AGL
!
!    istatus       status indicator
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'gribcst.inc'

  CHARACTER (LEN=*) :: dir_extd
  CHARACTER (LEN=*) :: extdname

  INTEGER :: extdopt
  INTEGER :: extdfmt

  CHARACTER (LEN=19) :: extdinit
  CHARACTER (LEN=9) :: extdfcst
  CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
!  External grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj_ext
  REAL :: scale_ext,trlon_ext
  REAL :: latnot_ext(2)
  REAL :: x0_ext,y0_ext
  REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
!  Output external variable arrays
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx_ext,ny_ext,nz_ext

  REAL :: lat_ext(nx_ext,ny_ext)
  REAL :: lon_ext(nx_ext,ny_ext)
  REAL :: p_ext  (nx_ext,ny_ext,nz_ext)   ! Pressure (Pascals)
  REAL :: hgt_ext(nx_ext,ny_ext,nz_ext)   ! Height (m)
  REAL :: t_ext  (nx_ext,ny_ext,nz_ext)   ! Temperature (K)
  REAL :: qv_ext (nx_ext,ny_ext,nz_ext)   ! Specific humidity (kg/kg)
  REAL :: u_ext  (nx_ext,ny_ext,nz_ext)   ! Eastward wind component
  REAL :: v_ext  (nx_ext,ny_ext,nz_ext)   ! Northward wind component
  REAL :: qc_ext (nx_ext,ny_ext,nz_ext)   ! Cloud H2O mixing ratio (kg/kg)
  REAL :: qr_ext (nx_ext,ny_ext,nz_ext)   ! Rain  H2O mixing ratio (kg/kg)
  REAL :: qi_ext (nx_ext,ny_ext,nz_ext)   ! Ice   mixing ratio (kg/kg)
  REAL :: qs_ext (nx_ext,ny_ext,nz_ext)   ! Snow  mixing ratio (kg/kg)
  REAL :: qh_ext (nx_ext,ny_ext,nz_ext)   ! Hail  mixing ratio (kg/kg)

  REAL :: tsfc_ext   (nx_ext,ny_ext)      ! Temperature at surface (K)
  REAL :: tsoil_ext  (nx_ext,ny_ext)      ! Deep soil temperature (K)
  REAL :: wetsfc_ext (nx_ext,ny_ext)      ! Surface soil moisture (fraction)
  REAL :: wetdp_ext  (nx_ext,ny_ext)      ! Deep soil moisture (fraction)
  REAL :: wetcanp_ext(nx_ext,ny_ext)      ! Canopy water amount

  REAL :: trn_ext    (nx_ext,ny_ext)          ! Geometrical heights
  REAL :: psfc_ext   (nx_ext,ny_ext)      ! Surface pressure (Pa)

  REAL :: t_2m_ext (nx_ext,ny_ext)
  REAL :: rh_2m_ext(nx_ext,ny_ext)
  REAL :: u_10m_ext(nx_ext,ny_ext)
  REAL :: v_10m_ext(nx_ext,ny_ext)
!
!----------------------------------------------------------------------
!
!  Other  external variable arrays
!
!----------------------------------------------------------------------
!
  REAL :: x_ext(nx_ext)
  REAL :: y_ext(ny_ext)
  REAL :: z_ext(nz_ext)

  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: var_grb2d(:,:,:,:)   ! GRIB variables
  REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
  INTEGER, ALLOCATABLE :: var_lev3d(:,:,:)  ! Levels (hybrid) for
                                            ! each 3-D variable
  REAL, ALLOCATABLE :: rcdata(:)            ! temporary data array
!
!-----------------------------------------------------------------------
!
!  Original grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj
  REAL :: scale,trlon,x0,y0
  REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: gribfile
  CHARACTER (LEN=13) :: gribtime
  CHARACTER (LEN=10) :: runstr
  CHARACTER (LEN=3) :: fmtn
  INTEGER :: ichr,bchar,echar
  INTEGER :: i,j,k,ldir,ireturn
  INTEGER :: iyr,imo,iday,myr, jldy
  INTEGER :: ihr,imin,isec
  INTEGER :: ifhr,ifmin,ifsec
  INTEGER :: grbflen, len1, lenrun

  INTEGER :: m,n,nz1,max_nr2d,max_nr3d

  REAL :: qvsat, pilev, qvsatice

  REAL :: tema, temb

  INTEGER :: chklev, lvscan, kk, jj

  INTEGER :: iret             ! Return flag

  REAL, ALLOCATABLE :: utmp(:,:), vtmp(:,:)
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  GRIB grid information
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=42) :: gridesc ! Grid description

  INTEGER :: iproj_grb    ! Map projection indicator
  INTEGER :: gthin        ! Indicator of whether the grid is "thinned"

  INTEGER :: ni_grb       ! Number of points along x-axis
  INTEGER :: nj_grb       ! Number of points along y-axis
  INTEGER :: np_grb       ! Total number of horizontal grid points

  INTEGER :: nk_grb       ! Number of vertical parameters
  REAL :: zk_grb(nz_ext)  ! Vertical coordinate parameters

  INTEGER :: npeq         ! Number of lat circles from pole to equator
  INTEGER :: nit(nz_ext)  ! Number of x-points for thinned grid

  REAL :: pi_grb          ! x-coordinate of pole point
  REAL :: pj_grb          ! y-coordinate of pole point
  INTEGER :: ipole        ! Projection center flag

  REAL :: di_grb          ! x-direction increment or grid length
  REAL :: dj_grb          ! y-direction increment or grid length

  REAL :: lrb             ! y-direction increment or grid length

  REAL :: latsw           ! Latitude  of South West corner point
  REAL :: lonsw           ! Longitude of South West corner point
  REAL :: latne           ! Latitude  of North East corner point
  REAL :: lonne           ! Longitude of North East corner point

  REAL :: lattru1         ! Latitude (1st) at which projection is true
  REAL :: lattru2         ! Latitude (2nd) at which projection is true
  REAL :: lontrue         ! Longitude      at which projection is true

  REAL :: latrot          ! Latitude  of southern pole of rotation
  REAL :: lonrot          ! Longitude of southern pole of rotation
  REAL :: angrot          ! Angle of rotation

  REAL :: latstr          ! Latitude  of the pole of stretching
  REAL :: lonstr          ! Longitude of the pole of stretching
  REAL :: facstr          ! Stretching factor

  INTEGER :: scanmode     ! Scanning indicator
  INTEGER :: iscan        ! x-direction   scanning indicator
  INTEGER :: jscan        ! y-direction   scanning indicator
  INTEGER :: kscan        ! FORTRAN index scanning indicator

  INTEGER :: ires         ! Resolution direction increments indicator
  INTEGER :: iearth       ! Earth shape indicator: spherical or oblate?
  INTEGER :: icomp        ! (u,v) components decomposition indicator

  INTEGER :: jpenta       ! J-Pentagonal resolution parameter
  INTEGER :: kpenta       ! K-Pentagonal resolution parameter
  INTEGER :: mpenta       ! M-Pentagonal resolution parameter
  INTEGER :: ispect       ! Spectral representation type
  INTEGER :: icoeff       ! Spectral coefficient storage mode

  REAL :: xp_grb          ! X coordinate of sub-satellite point
  REAL :: yp_grb          ! Y coordinate of sub-satellite point
  REAL :: xo_grb          ! X coordinate of image sector origin
  REAL :: yo_grb          ! Y coordinate of image sector origin
  REAL :: zo_grb          ! Camera altitude from center of Earth
!
!-----------------------------------------------------------------------
!
!  Function f_qvsatl and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_qvsatl
!fpp$ expand (f_desdt)
!fpp$ expand (f_qvsatl)
!dir$ inline always f_desdt, f_qvsatl
!*$*  inline routine (f_desdt, f_qvsatl)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
  ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
  ALLOCATE(rcdata(nx_ext*ny_ext))
  ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))
  ALLOCATE(utmp(nx_ext,ny_ext))
  ALLOCATE(vtmp(nx_ext,ny_ext))

  IF ( extdfcst == '         ') extdfcst='000:00:00'

  READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)')                  &
      iyr,imo,iday,ihr,imin,isec

  CALL julday(iyr,imo,iday,jldy)

  myr=MOD(iyr,100)
  ifhr=0
  ifmin=0
  ifsec=0

  READ (extdfcst,'(i3)',ERR=4,END=4) ifhr

  4     CONTINUE

  WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)')                      &
         iyr,imo,iday,ihr,'f',ifhr

  len1=LEN(dir_extd)
  grbflen=len1

  CALL strlnth( dir_extd, grbflen )

  IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
    dir_extd = '.'
    grbflen=1
  END IF
  IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
    grbflen=grbflen+1
    dir_extd(grbflen:grbflen)='/'
  END IF

  lenrun = LEN( extdname )
  CALL strlnth( extdname, lenrun )

  gribfile =                                                            &
      dir_extd(1:grbflen)//extdname(1:lenrun)//'.'//gribtime(1:13)
  grbflen = grbflen + lenrun + 14

!
!-----------------------------------------------------------------------
!
!  RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
  gridtyp   = ruc211grid
  mproj_grb = ruc211proj

  n2dvars  = ruc211nvs2d
  n2dlvtps = ruc211nlvt2d

  DO k=1,n2dlvtps
    DO n=1,n2dvars
      var_id2d(n,k) = ruc211var_id2d(n,k)
    END DO
    levtyp2d(k) = ruc211levs2d(k)
  END DO

  n3dvars  = ruc211nvs3d
  n3dlvtps = ruc211nlvt3d

  DO m=1,n3dlvtps
    DO n=1,n3dvars
      var_id3d(n,m) = ruc211var_id3d(n,m)
    END DO
    levtyp3d(m) = ruc211levs3d(m)
  END DO

  CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime,        &
                gridesc, iproj_grb, gthin,                              &
                ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit,          &
                pi_grb,pj_grb,ipole, di_grb,dj_grb,                     &
                latsw,lonsw, latne,lonne,                               &
                latrot,lonrot,angrot,                                   &
                latstr,lonstr,facstr,                                   &
                lattru1,lattru2,lontrue,                                &
                scanmode, iscan,jscan,kscan,                            &
                ires,iearth,icomp,                                      &
                jpenta,kpenta,mpenta,ispect,icoeff,                     &
                xp_grb,yp_grb, xo_grb,yo_grb,zo_grb,                    &
                rcdata,var_grb2d,var_grb3d,var_lev3d,iret)

  max_nr2d = 0
  DO n=1,n2dvars
    DO m=1,n2dlvtps
      max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
    END DO
  END DO

  max_nr3d = 0
  DO n=1,n3dvars
    max_nr3d = MAX( max_nr3d, var_nr3d(n,1))
  END DO

  IF ( max_nr3d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 3-D variables was found in the GRIB file',                  &
        'Program stopped in GETNMCRUC211.'
    STOP
  END IF

  IF ( max_nr2d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 2-D variables was found in the GRIB file'
  END IF

!  write (6,'(/a7,2x,6(i7))')
!    :  'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)

!  DO 60 k=1,max_nr3d
!    write (6,'(i5,4x,6(i7))')
!    :  k,(var_lev3d(k,n,1),n=1,n3dvars)

  DO k=1,max_nr3d
    DO n=2,n3dvars
      IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
        WRITE (6,'(a)')                                                 &
            'Variables were not at the same level.',                    &
            'Program stopped in GETNMCRUC211.'
        STOP
      END IF
    END DO
  END DO

  IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN      ! Center in NP
    iproj_ext = 1
  ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -1
  ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN  ! Center in NP
    iproj_ext = 2
  ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -2
  ELSE IF ( iproj_grb == 1 ) THEN
    iproj_ext = 3
  ELSE IF ( iproj_grb == 0 ) THEN
    iproj_ext = 4
  ELSE
    WRITE (6,'(a)')                                                     &
        'Unknown map projection. Set to non-projection.'
    iproj_ext = 0
  END IF

  scale_ext = 1.0
  latnot_ext(1) = lattru1
  latnot_ext(2) = lattru2
  trlon_ext = lontrue

  dx_ext = di_grb
  dy_ext = dj_grb

  CALL getmapr(iproj,scale,latnot,trlon,x0,y0)
  CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
  CALL lltoxy(1,1,latsw,lonsw,x0_ext,y0_ext)

  DO i=1,nx_ext
    x_ext(i)=x0_ext+(i-1)*dx_ext
  END DO

  DO j=1,ny_ext
    y_ext(j)=y0_ext+(j-1)*dy_ext
  END DO

  CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!
!-----------------------------------------------------------------------
!
!  Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
  DO j=1,ny_ext
    DO i=1,nx_ext
      IF ( var_nr2d(1,1) == 0 ) THEN
        psfc_ext   (i,j) = -999.0
      ELSE
        psfc_ext   (i,j) = var_grb2d(i,j,1,1)
      END IF

      IF ( var_nr2d(2,1) == 0 ) THEN
        trn_ext    (i,j) = -999.0
      ELSE
        trn_ext    (i,j) = var_grb2d(i,j,2,1)
      END IF

      IF( var_nr2d(1,2) == 0.) THEN
        t_2m_ext(i,j)= -999.0
      ELSE
        t_2m_ext(i,j)= var_grb2d(i,j,1,2)
      END IF

!   at this point we are reading in the relative humidity
!   later we'll convert to specific humidity

      IF( var_nr2d(2,2) == 0.) THEN
        rh_2m_ext(i,j)= -999.0
      ELSE
        rh_2m_ext(i,j)= var_grb2d(i,j,2,2)
      END IF

      IF( var_nr2d(3,2) == 0.) THEN
        u_10m_ext(i,j)= -999.0
      ELSE
        u_10m_ext(i,j)= var_grb2d(i,j,3,2)
      END IF

      IF( var_nr2d(4,2) == 0.) THEN
        v_10m_ext(i,j)= -999.0
      ELSE
        v_10m_ext(i,j)= var_grb2d(i,j,4,2)
      END IF

      tsfc_ext   (i,j) = -999.0
      tsoil_ext  (i,j) = -999.0
      wetsfc_ext (i,j) = -999.0
      wetdp_ext  (i,j) = -999.0
      wetcanp_ext(i,j) = -999.0

    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!
  nz1 = MIN(var_nr3d(1,1),nz_ext)

  IF ( var_lev3d(1,1,1) > var_lev3d(nz1,1,1) ) THEN  ! 1st level at sfc
    chklev = 1
    lvscan = 0
  ELSE
    chklev = -1
    lvscan = nz1+1
  END IF

  DO k=1,nz1
    kk = chklev * k + lvscan
    DO j=1,ny_ext
      DO i=1,nx_ext
        p_ext(i,j,kk) = 100.0 * FLOAT(var_lev3d(k,1,1)) ! Pressure
        t_ext(i,j,kk) = var_grb3d(i,j,k,2,1)    ! Temperature (K)
        u_ext(i,j,kk) = var_grb3d(i,j,k,4,1)    ! u wind (m/s)
        v_ext(i,j,kk) = var_grb3d(i,j,k,5,1)    ! v wind (m/s)
        hgt_ext(i,j,kk) = var_grb3d(i,j,k,1,1)    ! height (m)

!   check for portions of constant pressure grids that are below the surface

!       IF(((kk == 1) .OR. (p_ext(i,j,kk) > psfc_ext(i,j))) .AND.       &
!             (psfc_ext(i,j) > 0.) ) THEN
!         p_ext(i,j,kk)= psfc_ext(i,j) - 1. * (kk - 1)
!         t_ext(i,j,kk)= t_2m_ext(i,j)
!         u_ext(i,j,kk)= u_10m_ext(i,j)
!         v_ext(i,j,kk)= v_10m_ext(i,j)
!         hgt_ext(i,j,kk)= trn_ext(i,j) + 0.1 * (kk - 1)
!       END IF

        IF( (p_ext(i,j,kk) > 0.0) .AND. (t_ext(i,j,kk) > 0.0) ) THEN
          qvsat = f_qvsatl( p_ext(i,j,kk), t_ext(i,j,kk) )
          qv_ext(i,j,kk)= var_grb3d(i,j,k,3,1) * qvsat * 0.01

          IF((kk == 1) .OR. (p_ext(i,j,kk) > psfc_ext(i,j)))THEN
            qv_ext(i,j,kk)= rh_2m_ext(i,j) * qvsat * 0.01
          END IF
          qc_ext(i,j,kk)= 0.0
          qi_ext(i,j,kk)= 0.0
        ELSE
          qv_ext(i,j,kk)= 0.0
          qi_ext(i,j,kk)= 0.0
          qc_ext(i,j,kk)= 0.0
        END IF

        qr_ext (i,j,kk) = -999.
        qs_ext (i,j,kk) = -999.
        qh_ext (i,j,kk) = -999.
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Rotate winds to be relative to true north.
!  The RUCawips data are sent as grid-relative.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz1
!2001-05-16 GMB: Having umap & uear (or vmap & vear) point to
!the same array causes numerical errors when optimizing.
    CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k),               &
                 lon_ext,utmp,vtmp)
    u_ext(:,:,k) = utmp(:,:)
    v_ext(:,:,k) = vtmp(:,:)
  END DO

!
!-----------------------------------------------------------------------
!
!  Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
  CALL setmapr(iproj,scale,latnot,trlon)
  CALL setorig(1,x0,y0)

  istatus = 1
  DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)
  DEALLOCATE(utmp)
  DEALLOCATE(vtmp)

  RETURN
END SUBROUTINE getnmcruc211
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GETREANALT62               ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE getreanalt62(nx_ext,ny_ext,nz_ext,                           & 1,4
           dir_extd,extdname,extdopt,extdfmt,                           &
           extdinit,extdfcst,julfname,                                  &
           iproj_ext,scale_ext,                                         &
           trlon_ext,latnot_ext,x0_ext,y0_ext,                          &
           lat_ext,lon_ext,                                             &
           p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,                      &
           qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                          &
           tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext,         &
           trn_ext,psfc_ext, istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!
!  Reads a Global ReAnalaysis file for processing by ext2arps,
!  a program that converts external files to ARPS variables and format.
!  This version is useful when you want to use an ARPS file
!  with a different orientation or terrain than your final
!  ARPS product so arpsr2h does not work.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Dennis Moon, SSESCO
!  06/19/1998
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    dir_extd      Directory name for external file
!    extdname      Prefix string of external file name
!    extdopt       Option of external data sources
!    extdfmt       Option of external data format
!
!    extdinit      Initialized time in mm-dd-yyyy:hh:mm:ss format
!    extdfcst      Forecast hour in HHH:MM:SS format
!    julfname      File name in yyjjjhhmm format
!
!  OUTPUT:
!
!    iproj_ext     Map projection number of external data
!    scale_ext     Scale factor of external data
!    trlon_ext     True longitude of external data (degrees E)
!    latnot_ext(2) True latitude(s) of external data (degrees N)
!    x0_ext        x coordinate of origin of external data
!    y0_ext        y coordinate of origin of external data
!    lat_ext       latitude of external data points (degrees N)
!    lon_ext       longitude of external data points (degrees E)
!    p_ext         pressure (Pascal)
!    hgt_ext       height (m)
!    t_ext         temperature (K)
!    qv_ext        specific humidity (kg/kg)
!    u_ext         u wind component (m/s)
!    v_ext         v wind component (m/s)
!    qc_ext        Cloud water mixing ratio (kg/kg)
!    qr_ext        Rain water mixing ratio (kg/kg)
!    qi_ext        Ice mixing ratio (kg/kg)
!    qs_ext        Snow mixing ratio (kg/kg)
!    qh_ext        Hail mixing ratio (kg/kg)
!
!    tsfc_ext      Surface temperature
!    tsoil_ext     Soil temperature
!    wetsfc_ext    Top layer soil moisture (fraction)
!    wetdp_ext     Deep soil moisture (fraction)
!    wetcanp_ext   Water content on canopy
!
!    trn_ext       External terrain (m)
!    psfc_ext      Surface pressure (Pa)
!
!    istatus       status indicator
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'gribcst.inc'

  CHARACTER (LEN=*) :: dir_extd
  CHARACTER (LEN=*) :: extdname

  INTEGER :: extdopt
  INTEGER :: extdfmt

  CHARACTER (LEN=19) :: extdinit
  CHARACTER (LEN=9) :: extdfcst
  CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
!  External grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj_ext
  REAL :: scale_ext,trlon_ext
  REAL :: latnot_ext(2)
  REAL :: x0_ext,y0_ext
  REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
!  Output external variable arrays
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx_ext,ny_ext,nz_ext

  REAL :: lat_ext(nx_ext,ny_ext)
  REAL :: lon_ext(nx_ext,ny_ext)
  REAL :: p_ext  (nx_ext,ny_ext,nz_ext)   ! Pressure (Pascals)
  REAL :: hgt_ext(nx_ext,ny_ext,nz_ext)   ! Height (m)
  REAL :: t_ext  (nx_ext,ny_ext,nz_ext)   ! Temperature (K)
  REAL :: qv_ext (nx_ext,ny_ext,nz_ext)   ! Specific humidity (kg/kg)
  REAL :: u_ext  (nx_ext,ny_ext,nz_ext)   ! Eastward wind component
  REAL :: v_ext  (nx_ext,ny_ext,nz_ext)   ! Northward wind component
  REAL :: qc_ext (nx_ext,ny_ext,nz_ext)   ! Cloud H2O mixing ratio (kg/kg)
  REAL :: qr_ext (nx_ext,ny_ext,nz_ext)   ! Rain  H2O mixing ratio (kg/kg)
  REAL :: qi_ext (nx_ext,ny_ext,nz_ext)   ! Ice   mixing ratio (kg/kg)
  REAL :: qs_ext (nx_ext,ny_ext,nz_ext)   ! Snow  mixing ratio (kg/kg)
  REAL :: qh_ext (nx_ext,ny_ext,nz_ext)   ! Hail  mixing ratio (kg/kg)

  REAL :: tsfc_ext   (nx_ext,ny_ext)      ! Temperature at surface (K)
  REAL :: tsoil_ext  (nx_ext,ny_ext)      ! Deep soil temperature (K)
  REAL :: wetsfc_ext (nx_ext,ny_ext)      ! Surface soil moisture (fraction)
  REAL :: wetdp_ext  (nx_ext,ny_ext)      ! Deep soil moisture (fraction)
  REAL :: wetcanp_ext(nx_ext,ny_ext)      ! Canopy water amount

  REAL :: trn_ext    (nx_ext,ny_ext)          ! Geometrical heights
  REAL :: psfc_ext   (nx_ext,ny_ext)      ! Surface pressure (Pa)

!
!----------------------------------------------------------------------
!
!  Other  external variable arrays
!
!----------------------------------------------------------------------
!
  REAL :: x_ext(nx_ext)
  REAL :: y_ext(ny_ext)
  REAL :: z_ext(nz_ext)

  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: var_grb2d(:,:,:,:)   ! GRIB variables
  REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
  INTEGER, ALLOCATABLE :: var_lev3d(:,:,:)  ! Levels (hybrid) for
                                            ! each 3-D variable
  REAL, ALLOCATABLE :: rcdata(:)            ! temporary data array
!
!-----------------------------------------------------------------------
!
!  Original grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj
  REAL :: scale,trlon,x0,y0
  REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: gribfile
  CHARACTER (LEN=13) :: gribtime
  CHARACTER (LEN=10) :: runstr
  CHARACTER (LEN=3) :: fmtn
  INTEGER :: ichr,bchar,echar
  INTEGER :: i,j,k,ldir,ireturn
  INTEGER :: iyr,imo,iday,myr, jldy
  INTEGER :: ihr,imin,isec
  INTEGER :: ifhr,ifmin,ifsec
  INTEGER :: grbflen, len1, lenrun

  INTEGER :: m,n,nz1,max_nr2d,max_nr3d

  REAL :: qvsat, pilev, qvsatice

  REAL :: tema, temb, qvavg, mixrat, tvirtbar, tmpval
  REAL :: gausslat(200)

  INTEGER :: chklev, lvscan, kk, jj

  INTEGER :: iret             ! Return flag
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  GRIB grid information
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=42) :: gridesc ! Grid description

  INTEGER :: iproj_grb    ! Map projection indicator
  INTEGER :: gthin        ! Indicator of whether the grid is "thinned"

  INTEGER :: ni_grb       ! Number of points along x-axis
  INTEGER :: nj_grb       ! Number of points along y-axis
  INTEGER :: np_grb       ! Total number of horizontal grid points

  INTEGER :: nk_grb       ! Number of vertical parameters
  REAL :: zk_grb(nz_ext)  ! Vertical coordinate parameters

  INTEGER :: npeq         ! Number of lat circles from pole to equator
  INTEGER :: nit(nz_ext)  ! Number of x-points for thinned grid

  REAL :: pi_grb          ! x-coordinate of pole point
  REAL :: pj_grb          ! y-coordinate of pole point
  INTEGER :: ipole        ! Projection center flag

  REAL :: di_grb          ! x-direction increment or grid length
  REAL :: dj_grb          ! y-direction increment or grid length

  REAL :: lrb             ! y-direction increment or grid length

  REAL :: latsw           ! Latitude  of South West corner point
  REAL :: lonsw           ! Longitude of South West corner point
  REAL :: latne           ! Latitude  of North East corner point
  REAL :: lonne           ! Longitude of North East corner point

  REAL :: lattru1         ! Latitude (1st) at which projection is true
  REAL :: lattru2         ! Latitude (2nd) at which projection is true
  REAL :: lontrue         ! Longitude      at which projection is true

  REAL :: latrot          ! Latitude  of southern pole of rotation
  REAL :: lonrot          ! Longitude of southern pole of rotation
  REAL :: angrot          ! Angle of rotation

  REAL :: latstr          ! Latitude  of the pole of stretching
  REAL :: lonstr          ! Longitude of the pole of stretching
  REAL :: facstr          ! Stretching factor

  INTEGER :: scanmode     ! Scanning indicator
  INTEGER :: iscan        ! x-direction   scanning indicator
  INTEGER :: jscan        ! y-direction   scanning indicator
  INTEGER :: kscan        ! FORTRAN index scanning indicator

  INTEGER :: ires         ! Resolution direction increments indicator
  INTEGER :: iearth       ! Earth shape indicator: spherical or oblate?
  INTEGER :: icomp        ! (u,v) components decomposition indicator

  INTEGER :: jpenta       ! J-Pentagonal resolution parameter
  INTEGER :: kpenta       ! K-Pentagonal resolution parameter
  INTEGER :: mpenta       ! M-Pentagonal resolution parameter
  INTEGER :: ispect       ! Spectral representation type
  INTEGER :: icoeff       ! Spectral coefficient storage mode

  REAL :: xp_grb          ! X coordinate of sub-satellite point
  REAL :: yp_grb          ! Y coordinate of sub-satellite point
  REAL :: xo_grb          ! X coordinate of image sector origin
  REAL :: yo_grb          ! Y coordinate of image sector origin
  REAL :: zo_grb          ! Camera altitude from center of Earth
!
!-----------------------------------------------------------------------
!
!  T62 Gaussian grid latitudes
!
!-----------------------------------------------------------------------
!
  DATA gausslat /  88.5420,  86.6532,  84.7532,  82.8508,  80.9474,     &
                   79.0435,  77.1393,  75.2351,  73.3307,  71.4262,     &
                   69.5217,  67.6171,  65.7125,  63.8079,  61.9033,     &
                   59.9986,  58.0940,  56.1893,  54.2846,  52.3799,     &
                   50.4752,  48.5705,  46.6658,  44.7611,  42.8564,     &
                   40.9517,  39.0470,  37.1422,  35.2375,  33.3328,     &
                   31.4281,  29.5234,  27.6186,  25.7139,  23.8092,     &
                   21.9044,  19.9997,  18.0950,  16.1902,  14.2855,     &
                   12.3808,  10.4760,   8.5713,   6.6666,   4.7618,     &
                    2.8571,   0.9524,  -0.9524,  -2.8571,  -4.7618,     &
                   -6.6666,  -8.5713, -10.4760, -12.3808, -14.2855,     &
                  -16.1902, -18.0950, -19.9997, -21.9044, -23.8092,     &
                  -25.7139, -27.6186, -29.5234, -31.4281, -33.3328,     &
                  -35.2375, -37.1422, -39.0470, -40.9517, -42.8564,     &
                  -44.7611, -46.6658, -48.5705, -50.4752, -52.3799,     &
                  -54.2846, -56.1893, -58.0940, -59.9986, -61.9033,     &
                  -63.8079, -65.7125, -67.6171, -69.5217, -71.4262,     &
                  -73.3307, -75.2351, -77.1393, -79.0435, -80.9474,     &
                  -82.8508, -84.7532, -86.6532, -88.5420, 106*0.0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
  ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
  ALLOCATE(rcdata(nx_ext*ny_ext))
  ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))

  IF ( extdfcst == '         ') extdfcst='000:00:00'

  READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)')                  &
      iyr,imo,iday,ihr,imin,isec

  CALL julday(iyr,imo,iday,jldy)

  myr=MOD(iyr,100)
  ifhr=0
  ifmin=0
  ifsec=0

  READ (extdfcst,'(i3)',ERR=4,END=4) ifhr

  4     CONTINUE

  WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)')                      &
         iyr,imo,iday,ihr,'f',ifhr

  len1=LEN(dir_extd)
  grbflen=len1

  CALL strlnth( dir_extd, grbflen )

  IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
    dir_extd = '.'
    grbflen=1
  END IF
  IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
    grbflen=grbflen+1
    dir_extd(grbflen:grbflen)='/'
  END IF

  lenrun = LEN( extdname )
  CALL strlnth( extdname, lenrun )

  gribfile =                                                            &
      dir_extd(1:grbflen)//extdname(1:lenrun)//'.'//gribtime(1:13)
  grbflen = grbflen + lenrun + 14

!
!-----------------------------------------------------------------------
!
!  RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
  gridtyp   = reanalt62grid
  mproj_grb = reanalt62proj

  n2dvars  = reanalt62nvs2d
  n2dlvtps = reanalt62nlvt2d

  DO k=1,n2dlvtps
    DO n=1,n2dvars
      var_id2d(n,k) = reanalt62var_id2d(n,k)
    END DO
    levtyp2d(k) = reanalt62levs2d(k)
  END DO

  n3dvars  = reanalt62nvs3d
  n3dlvtps = reanalt62nlvt3d

  DO m=1,n3dlvtps
    DO n=1,n3dvars
      var_id3d(n,m) = reanalt62var_id3d(n,m)
    END DO
    levtyp3d(m) = reanalt62levs3d(m)
  END DO

  CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime,        &
                gridesc, iproj_grb, gthin,                              &
                ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit,          &
                pi_grb,pj_grb,ipole, di_grb,dj_grb,                     &
                latsw,lonsw, latne,lonne,                               &
                latrot,lonrot,angrot,                                   &
                latstr,lonstr,facstr,                                   &
                lattru1,lattru2,lontrue,                                &
                scanmode, iscan,jscan,kscan,                            &
                ires,iearth,icomp,                                      &
                jpenta,kpenta,mpenta,ispect,icoeff,                     &
                xp_grb,yp_grb, xo_grb,yo_grb,zo_grb,                    &
                rcdata,var_grb2d,var_grb3d,var_lev3d,iret)

  max_nr2d = 0
  DO n=1,n2dvars
    DO m=1,n2dlvtps
      max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
    END DO
  END DO

  max_nr3d = 0
  DO n=1,n3dvars
    max_nr3d = MAX( max_nr3d, var_nr3d(n,1))
  END DO

  IF ( max_nr3d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 3-D variables was found in the GRIB file',                  &
        'Program stopped in GETREANALT62.'
    STOP
  END IF

  IF ( max_nr2d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 2-D variables was found in the GRIB file'
  END IF

!  write (6,'(/a7,2x,6(i7))')
!    :  'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)

!  DO 60 k=1,max_nr3d
!    write (6,'(i5,4x,6(i7))')
!    :  k,(var_lev3d(k,n,1),n=1,n3dvars)

  DO k=1,max_nr3d
    DO n=2,n3dvars
      IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
        WRITE (6,'(a)')                                                 &
            'Variables were not at the same level.',                    &
            'Program stopped in GETREANALT62.'
        STOP
      END IF
    END DO
  END DO

  IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN      ! Center in NP
    iproj_ext = 1
  ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -1
  ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN  ! Center in NP
    iproj_ext = 2
  ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -2
  ELSE IF ( iproj_grb == 1 ) THEN
    iproj_ext = 3
  ELSE IF ( (iproj_grb == 0) .OR. (iproj_grb == 4) ) THEN
    iproj_ext = 4
  ELSE
    WRITE (6,'(a)')                                                     &
        'Unknown map projection. Set to non-projection.'
    iproj_ext = 0
  END IF

  scale_ext = 1.0
  latnot_ext(1) = lattru1
  latnot_ext(2) = lattru2
  trlon_ext = lontrue

  latnot_ext(1)= 0.
  trlon_ext= 0.

  dx_ext = di_grb
  dy_ext = dj_grb

  DO i=1, nx_ext
    DO j=1, ny_ext
      lat_ext(i,j)= gausslat(j)
      lon_ext(i,j)= lonsw + (i-1) * dx_ext
    END DO
  END DO

!   swap the latitude values so that +j moves to the north

  DO i=1, nx_ext
    DO j=1, ny_ext/2
      tmpval= lat_ext(i,ny_ext-j+1)
      lat_ext(i,ny_ext-j+1)= lat_ext(i,j)
      lat_ext(i,j)= tmpval
    END DO
  END DO

!   call getmapr(iproj,scale,latnot,trlon,x0,y0)
!   call setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
!   call lltoxy(1,1,LatSW,LonSW,x0_ext,y0_ext)

!   DO 100 i=1,nx_ext
!     x_ext(i)=x0_ext+(i-1)*dx_ext
!100   CONTINUE

!   DO 110 j=1,ny_ext
!     y_ext(j)=y0_ext+(j-1)*dy_ext
!110   CONTINUE

!   CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!
!-----------------------------------------------------------------------
!
!  Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
  DO j=1,ny_ext
    DO i=1,nx_ext
      IF ( var_nr2d(1,1) == 0 ) THEN
        psfc_ext   (i,j) = -999.0
      ELSE
        psfc_ext   (i,j) = var_grb2d(i,j,1,1)
      END IF

      IF ( var_nr2d(2,1) == 0 ) THEN
        trn_ext    (i,j) = -999.0
      ELSE
        trn_ext    (i,j) = var_grb2d(i,j,2,1)
      END IF

      tsfc_ext   (i,j) = -999.0
      tsoil_ext  (i,j) = -999.0
      wetsfc_ext (i,j) = -999.0
      wetdp_ext  (i,j) = -999.0
      wetcanp_ext(i,j) = -999.0

    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!
  nz1 = MIN(var_nr3d(1,1),nz_ext)

  IF ( var_lev3d(1,1,1) > var_lev3d(nz1,1,1) ) THEN  ! 1st level at sfc
    chklev = 1
    lvscan = 0
  ELSE
    chklev = -1
    lvscan = nz1+1
  END IF

  DO k=1,nz1
    kk = chklev * k + lvscan
    DO j=1,ny_ext
      DO i=1,nx_ext
        p_ext(i,j,kk) = 1.e-4 * FLOAT(var_lev3d(k,1,1))                 &
                         * psfc_ext(i,j)        ! Pressure
        t_ext(i,j,kk) = var_grb3d(i,j,k,1,1)    ! Temperature (K)
        qv_ext(i,j,kk)= var_grb3d(i,j,k,2,1)    ! Specific Humidity
        u_ext(i,j,kk) = var_grb3d(i,j,k,3,1)    ! u wind (m/s)
        v_ext(i,j,kk) = var_grb3d(i,j,k,4,1)    ! v wind (m/s)

!   calculate height from sigma-P values

        IF( k == 1) THEN
          mixrat= 1.0 / ( 1./qv_ext(i,j,kk) - 1.)
          tvirtbar= t_ext(i,j,kk) * ( 1.0 + .61* mixrat)
          hgt_ext(i,j,kk)= trn_ext(i,j) + 29.3 * tvirtbar *             &
                          LOG(psfc_ext(i,j)/p_ext(i,j,kk))
        ELSE
          qvavg= 0.5 * (qv_ext(i,j,kk) + qv_ext(i,j,kk-1))
          mixrat= 1.0 / (1./qvavg - .1)
          tvirtbar= 0.5 * (t_ext(i,j,kk) + t_ext(i,j,kk-1)) *           &
                          ( 1.0 + .61 * mixrat)
          hgt_ext(i,j,kk)= hgt_ext(i,j,kk-1) + 29.3 * tvirtbar *        &
                          LOG(p_ext(i,j,kk-1)/p_ext(i,j,kk))
        END IF

        qc_ext(i,j,kk)= 0.0
        qi_ext(i,j,kk)= 0.0

        qr_ext (i,j,kk) = -999.
        qs_ext (i,j,kk) = -999.
        qh_ext (i,j,kk) = -999.
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  No need to rotate winds to be relative to true north.
!  The Re-analysis grid winds are true N-S, E-W
!
!-----------------------------------------------------------------------
!
!   DO 300 k=1,nz1
!     CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k),
!  :               lon_ext,u_ext(1,1,k),v_ext(1,1,k))
!300   CONTINUE

!
!-----------------------------------------------------------------------
!
!  Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
!   call setmapr(iproj,scale,latnot,trlon)
!   call setorig(1,x0,y0)

  istatus = 1

  DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)

  RETURN
END SUBROUTINE getreanalt62

!
!##################################################################
!##################################################################
!######                                                      ######
!######                   SUBROUTINE TV2TQ                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE tv2tq(tv,p,t,q) 1,4
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  REAL :: tv       ! virtual temperature (K)
  REAL :: p        ! pressure (pa)
  REAL :: t        ! temperature (K)
  REAL :: q        ! specific humidity (kg/kg)

  REAL :: tg,qg,fac,tvg,dtvdt,tgnew
  REAL :: es, es1

  INTEGER :: it
!
!-----------------------------------------------------------------------
!
!  Function and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_es

!fpp$ expand (f_es)
!dir$ inline always f_es
!*$*  inline routine (f_es)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  it = 0
!
! -- TG - guess at non-virtual temperature
!
  es = f_es(p,tv)
  tg = tv/( 1.0 + 0.6078 * 0.62197*es                                   &
                / ( p - es ) )
  10    CONTINUE
  it = it + 1
!
! -- QG - guess at specific humidity
!
  es = f_es(p,tg)
  qg = 0.62197*es/( p - es )
  fac = 1.+0.6078*qg
!
! -- TVG - guess at virtual temperature
!
  tvg = tg*fac
!
! -- DTVDT - d(Tv)/dT estimated over 0.1 K interval
!    from 1st term Taylor series expansion
!
  es1 = f_es(p,tg+0.1)
  dtvdt = fac + tg * 0.6078 * 0.62197*p/(p-es)**2                       &
       * ( es1 - es ) / 0.1
  tgnew = tg + (tv-tvg)/dtvdt
!
!    write(12,*) ' '
!    write(12,*) 'IT=',IT,'  TV=',TV,'  TVG=',TVG,
!    :              '  TV-TVG=',TV-TVG,'  DTVDT=',DTVDT
!
  IF (ABS(tv-tvg) < 1.0E-4) GO TO 20
  tg = tgnew
  GO TO 10

  20    CONTINUE

  t = tgnew
  es = f_es(p,t)
  q = 0.62197*es / ( p - es )

  RETURN
END SUBROUTINE tv2tq
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GETNMCETA212               ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. soiltyp_ext

SUBROUTINE getnmceta212(nx_ext,ny_ext,nz_ext,                           & 1,11
           dir_extd,extdname,extdopt,extdfmt,                           &
           extdinit,extdfcst,julfname,                                  &
           iproj_ext,scale_ext,                                         &
           trlon_ext,latnot_ext,x0_ext,y0_ext,                          &
           lat_ext,lon_ext,                                             &
           p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,                      &
           qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                          &
           tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext,         &
           snowdpth_ext,trn_ext,psfc_ext,soiltyp_ext,                   &
           istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Reads a NMC GRIB (Grid #212, 40km) ETA file for processing by
!  ext2arps, a program that converts external files to ARPS variables
!  and format.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!  02/22/1996
!
!  MODIFICATION HISTORY:
!
!  1999/08/03 (Gene Bassett)
!  Modified the deep soil moisture and temperature to be an average of
!  the three soil layers covering 0 to 100 cm in depth.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    dir_extd      Directory name for external file
!    extdname      Prefix string of external file name
!    extdopt       Option of external data sources
!    extdfmt       Option of external data format

!    extdinit      Initialized time in mm-dd-yyyy:hh:mm:ss format
!    extdfcst      Forecast hour in HHH:MM:SS format
!    julfname      File name in yyjjjhhmm format
!
!  OUTPUT:
!
!    iproj_ext     Map projection number of external data
!    scale_ext     Scale factor of external data
!    trlon_ext     True longitude of external data (degrees E)
!    latnot_ext(2) True latitude(s) of external data (degrees N)
!    x0_ext        x coordinate of origin of external data
!    y0_ext        y coordinate of origin of external data
!    lat_ext       latitude of external data points (degrees N)
!    lon_ext       longitude of external data points (degrees E)
!    p_ext         pressure (Pascal)
!    hgt_ext       height (m)
!    t_ext         temperature (K)
!    qv_ext        specific humidity (kg/kg)
!    u_ext         u wind component (m/s)
!    v_ext         v wind component (m/s)
!    qc_ext        Cloud water mixing ratio (kg/kg)
!    qr_ext        Rain water mixing ratio (kg/kg)
!    qi_ext        Ice mixing ratio (kg/kg)
!    qs_ext        Snow mixing ratio (kg/kg)
!    qh_ext        Hail mixing ratio (kg/kg)
!
!    tsfc_ext      Surface temperature
!    tsoil_ext     Soil temperature
!    wetsfc_ext    Top layer soil moisture (fraction)
!    wetdp_ext     Deep soil moisture (fraction)
!    wetcanp_ext   Water content on canopy
!
!    trn_ext       External terrain (m)
!    psfc_ext      Surface pressure (Pa)
!
!    istatus       status indicator
!
!  WORK ARRAYS:
!
!    var_grb3d     Arrays to store the GRIB 3-D variables:
!                  var_grb3d(nxgrb,nygrb,nzgrb,1,1) - Temperature (K)
!                  var_grb3d(nxgrb,nygrb,nzgrb,2,1) - Specific humidity
!                                                     (kg/kg)
!                  var_grb3d(nxgrb,nygrb,nzgrb,3,1) - u wind (m/s)
!                  var_grb3d(nxgrb,nygrb,nzgrb,4,1) - v wind (m/s)
!                  var_grb3d(nxgrb,nygrb,nzgrb,5,1) - Geopotential
!                                                     height (gpm)
!                  var_grb3d(nxgrb,nygrb,nzgrb,6,1) - Pressure vertical
!                                                     velocity (Pa/s)
!                                                     (if applied)
!                  var_grb3d(nxgrb,nygrb,nzgrb,1,2) - soil temp. (K)
!                  var_grb3d(nxgrb,nygrb,nzgrb,2,2) - vol. soil moist.
!                                                     (m**3/m**3)
!
!    var_grb2d     Arrays to store the GRIB 2-D variables:
!                  var_grb2d(nxgrb,nygrb,1) - Surface pressure (Pa)
!                  var_grb2d(nxgrb,nygrb,2) - Geopotential height (gpm)
!                  var_grb2d(nxgrb,nygrb,3) - Surface temperature (K)
!                  var_grb2d(nxgrb,nygrb,4) - Plant canopy surface
!                                             water (kg/m**2)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'gribcst.inc'

  CHARACTER (LEN=*) :: dir_extd
  CHARACTER (LEN=*) :: extdname

  INTEGER :: extdopt
  INTEGER :: extdfmt

  CHARACTER (LEN=19) :: extdinit
  CHARACTER (LEN=9) :: extdfcst
  CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
!  External grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj_ext
  REAL :: scale_ext,trlon_ext
  REAL :: latnot_ext(2)
  REAL :: x0_ext,y0_ext
  REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
!  Output external variable arrays
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx_ext,ny_ext,nz_ext

  REAL :: lat_ext(nx_ext,ny_ext)
  REAL :: lon_ext(nx_ext,ny_ext)
  REAL :: p_ext  (nx_ext,ny_ext,nz_ext)   ! Pressure (Pascals)
  REAL :: hgt_ext(nx_ext,ny_ext,nz_ext)   ! Height (m)
  REAL :: t_ext  (nx_ext,ny_ext,nz_ext)   ! Temperature (K)
  REAL :: qv_ext (nx_ext,ny_ext,nz_ext)   ! Specific humidity (kg/kg)
  REAL :: u_ext  (nx_ext,ny_ext,nz_ext)   ! Eastward wind component
  REAL :: v_ext  (nx_ext,ny_ext,nz_ext)   ! Northward wind component
  REAL :: qc_ext (nx_ext,ny_ext,nz_ext)   ! Cloud H2O mixing ratio (kg/kg)
  REAL :: qr_ext (nx_ext,ny_ext,nz_ext)   ! Rain  H2O mixing ratio (kg/kg)
  REAL :: qi_ext (nx_ext,ny_ext,nz_ext)   ! Ice   mixing ratio (kg/kg)
  REAL :: qs_ext (nx_ext,ny_ext,nz_ext)   ! Snow  mixing ratio (kg/kg)
  REAL :: qh_ext (nx_ext,ny_ext,nz_ext)   ! Hail  mixing ratio (kg/kg)

  REAL :: tsfc_ext   (nx_ext,ny_ext)      ! Temperature at surface (K)
  REAL :: tsoil_ext  (nx_ext,ny_ext)      ! Deep soil temperature (K)
  REAL :: wetsfc_ext (nx_ext,ny_ext)      ! Surface soil moisture (fraction)
  REAL :: wetdp_ext  (nx_ext,ny_ext)      ! Deep soil moisture (fraction)
  REAL :: wetcanp_ext(nx_ext,ny_ext)      ! Canopy water amount
  REAL :: snowdpth_ext(nx_ext,ny_ext)     ! Snow depth (m)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2001-03-14 GMB
  INTEGER :: soiltyp_ext(nx_ext,ny_ext)   ! Soil type
  !MARK FIXME continue ...

  REAL :: trn_ext    (nx_ext,ny_ext)      ! External terrain (m)
  REAL :: psfc_ext   (nx_ext,ny_ext)      ! Surface pressure (Pa)
!
!-----------------------------------------------------------------------
!
!  Other  external variable arrays
!
!-----------------------------------------------------------------------
!
  REAL :: x_ext(nx_ext)
  REAL :: y_ext(ny_ext)

  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: var_grb2d(:,:,:,:)   ! GRIB variables
  REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
  INTEGER, ALLOCATABLE :: var_lev3d(:,:,:)  ! Levels (hybrid) for
                                            ! each 3-D variable
  REAL, ALLOCATABLE :: rcdata(:)            ! temporary data array
!
!-----------------------------------------------------------------------
!
!  Original grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj
  REAL :: scale,trlon,x0,y0
  REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: gribfile
  CHARACTER (LEN=13) :: gribtime
  INTEGER :: i,j,k,kk
  INTEGER :: iyr,imo,iday,myr, jldy
  INTEGER :: ihr,imin,isec
  INTEGER :: ifhr,ifmin,ifsec
  INTEGER :: grbflen, len1, lenrun

  INTEGER :: m,n,nz1,max_nr2d,max_nr3d

  REAL :: govrd

  INTEGER :: chklev, lvscan

  INTEGER :: iret             ! Return flag

  REAL, ALLOCATABLE :: utmp(:,:), vtmp(:,:)
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  GRIB grid information
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=42) :: gridesc ! Grid description

  INTEGER :: iproj_grb    ! Map projection indicator
  INTEGER :: gthin        ! Indicator of whether the grid is "thinned"

  INTEGER :: ni_grb       ! Number of points along x-axis
  INTEGER :: nj_grb       ! Number of points along y-axis
  INTEGER :: np_grb       ! Total number of horizontal grid points

  INTEGER :: nk_grb       ! Number of vertical parameters
  REAL :: zk_grb(nz_ext)  ! Vertical coordinate parameters

  INTEGER :: npeq         ! Number of lat circles from pole to equator
  INTEGER :: nit(nz_ext)  ! Number of x-points for thinned grid

  REAL :: pi_grb          ! x-coordinate of pole point
  REAL :: pj_grb          ! y-coordinate of pole point
  INTEGER :: ipole        ! Projection center flag

  REAL :: di_grb          ! x-direction increment or grid length
  REAL :: dj_grb          ! y-direction increment or grid length

  REAL :: latsw           ! Latitude  of South West corner point
  REAL :: lonsw           ! Longitude of South West corner point
  REAL :: latne           ! Latitude  of North East corner point
  REAL :: lonne           ! Longitude of North East corner point

  REAL :: lattru1         ! Latitude (1st) at which projection is true
  REAL :: lattru2         ! Latitude (2nd) at which projection is true
  REAL :: lontrue         ! Longitude      at which projection is true

  REAL :: latrot          ! Latitude  of southern pole of rotation
  REAL :: lonrot          ! Longitude of southern pole of rotation
  REAL :: angrot          ! Angle of rotation

  REAL :: latstr          ! Latitude  of the pole of stretching
  REAL :: lonstr          ! Longitude of the pole of stretching
  REAL :: facstr          ! Stretching factor

  INTEGER :: scanmode     ! Scanning indicator
  INTEGER :: iscan        ! x-direction   scanning indicator
  INTEGER :: jscan        ! y-direction   scanning indicator
  INTEGER :: kscan        ! FORTRAN index scanning indicator

  INTEGER :: ires         ! Resolution direction increments indicator
  INTEGER :: iearth       ! Earth shape indicator: spherical or oblate?
  INTEGER :: icomp        ! (u,v) components decomposition indicator

  INTEGER :: jpenta       ! J-Pentagonal resolution parameter
  INTEGER :: kpenta       ! K-Pentagonal resolution parameter
  INTEGER :: mpenta       ! M-Pentagonal resolution parameter
  INTEGER :: ispect       ! Spectral representation type
  INTEGER :: icoeff       ! Spectral coefficient storage mode

  REAL :: xp_grb          ! X coordinate of sub-satellite point
  REAL :: yp_grb          ! Y coordinate of sub-satellite point
  REAL :: xo_grb          ! X coordinate of image sector origin
  REAL :: yo_grb          ! Y coordinate of image sector origin
  REAL :: zo_grb          ! Camera altitude from center of Earth

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
  ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
  ALLOCATE(rcdata(nx_ext*ny_ext))
  ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))
  ALLOCATE(utmp(nx_ext,ny_ext))
  ALLOCATE(vtmp(nx_ext,ny_ext))

  IF(extdfcst == '         ') extdfcst='000:00:00'

  READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)')                  &
       iyr,imo,iday,ihr,imin,isec

  CALL julday(iyr,imo,iday,jldy)

  myr=MOD(iyr,100)
  ifhr=0
  ifmin=0
  ifsec=0

  READ(extdfcst,'(i3)',ERR=4,END=4) ifhr

4 CONTINUE

  WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)')                      &
         iyr,imo,iday,ihr,'f',ifhr

  len1=LEN(dir_extd)
  grbflen=len1

  CALL strlnth( dir_extd, grbflen )

  IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
    dir_extd = '.'
    grbflen=1
  END IF

  IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
    grbflen=grbflen+1
    dir_extd(grbflen:grbflen)='/'
  END IF

  lenrun = LEN( extdname )
  CALL strlnth( extdname, lenrun )

  gribfile = dir_extd(1:grbflen)//extdname(1:lenrun)                    &
                                //'.'//gribtime(1:13)
  grbflen = grbflen + lenrun + 14
!
!-----------------------------------------------------------------------
!
!  RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
  govrd = g/rd

  gridtyp  = eta212grid
  mproj_grb = eta212proj

  n2dvars  = eta212nvs2d
  n2dlvtps = eta212nlvt2d

  DO k=1,n2dlvtps
    DO n=1,n2dvars
      var_id2d(n,k) = eta212var_id2d(n,k)
    END DO
    levtyp2d(k) = eta212levs2d(k)
  END DO

  n3dvars  = eta212nvs3d
  n3dlvtps = eta212nlvt3d

  DO m=1,n3dlvtps
    DO n=1,n3dvars
      var_id3d(n,m) = eta212var_id3d(n,m)
    END DO
    levtyp3d(m) = eta212levs3d(m)
  END DO

  CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime,        &
                gridesc, iproj_grb, gthin,                              &
                ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit,          &
                pi_grb,pj_grb,ipole, di_grb,dj_grb,                     &
                latsw,lonsw, latne,lonne,                               &
                latrot,lonrot,angrot,                                   &
                latstr,lonstr,facstr,                                   &
                lattru1,lattru2,lontrue,                                &
                scanmode, iscan,jscan,kscan,                            &
                ires,iearth,icomp,                                      &
                jpenta,kpenta,mpenta,ispect,icoeff,                     &
                xp_grb,yp_grb, xo_grb,yo_grb,zo_grb,                    &
                rcdata,var_grb2d,var_grb3d,var_lev3d,iret)

  max_nr2d = 0
  DO n=1,n2dvars
    DO m=1,n2dlvtps
      max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
    END DO
  END DO

  max_nr3d = 0
  DO n=1,n3dvars
    max_nr3d = MAX( max_nr3d, var_nr3d(n,1) )
  END DO

  IF ( max_nr3d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 3-D variable was found in the GRIB file',                   &
        'Program stopped in GETNMCETA212.'
    STOP
  END IF

  IF ( max_nr2d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 2-D variables was found in the GRIB file'
  END IF

!  write (6,'(/a7,2x,6(i7))')
!    :    'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)

!  DO 60 k=1,max_nr3d
!    write (6,'(/i5,4x,6(i7))')
!    :    k,(var_lev3d(k,n,1),n=1,n3dvars)

  DO k=1,max_nr3d
    DO n=2,n3dvars
      IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
        WRITE (6,'(a)')                                                 &
            'Variables were not at the same level.',                    &
            'Program stopped in GETNMCETA212.'
        STOP
      END IF
    END DO
  END DO

  IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN      ! Center in NP
    iproj_ext = 1
  ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -1
  ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN  ! Center in NP
    iproj_ext = 2
  ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -2
  ELSE IF ( iproj_grb == 1 ) THEN
    iproj_ext = 3
  ELSE IF ( iproj_grb == 0 ) THEN
    iproj_ext = 4
  ELSE
    WRITE (6,'(a)')                                                     &
        'Unknown map projection. Set to non-projection.'
    iproj_ext = 0
  END IF

  scale_ext = 1.0
  latnot_ext(1) = lattru1
  latnot_ext(2) = lattru2
  trlon_ext = lontrue

  dx_ext = di_grb
  dy_ext = dj_grb

  CALL getmapr(iproj,scale,latnot,trlon,x0,y0)
  CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
  CALL lltoxy(1,1,latsw,lonsw,x0_ext,y0_ext)

  DO i=1,nx_ext
    x_ext(i)=x0_ext+(i-1)*dx_ext
  END DO

  DO j=1,ny_ext
    y_ext(j)=y0_ext+(j-1)*dy_ext
  END DO

  CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!
!-----------------------------------------------------------------------
!
!  Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
  DO j=1,ny_ext
    DO i=1,nx_ext
      IF ( var_nr2d(1,1) == 0 ) THEN
        psfc_ext   (i,j) = -999.0
      ELSE
        psfc_ext   (i,j) = var_grb2d(i,j,1,1) * 100.0
      END IF
      IF ( var_nr2d(2,1) == 0 ) THEN
        trn_ext    (i,j) = -999.0
      ELSE
        trn_ext    (i,j) = var_grb2d(i,j,2,1)
      END IF

      IF ( var_nr3d(1,2) == 0 ) THEN
        tsfc_ext  (i,j) = -999.0
        tsoil_ext (i,j) = -999.0
        wetsfc_ext(i,j) = -999.0
        wetdp_ext (i,j) = -999.0
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
        soiltyp_ext (i,j) = 0
      ELSE
        tsfc_ext  (i,j) = var_grb2d(i,j,3,1)       ! sfc temp.

        IF ( nint(var_grb2d(i,j,5,1)) == 1 ) THEN  ! soil temp over land
          tsoil_ext (i,j) = 0.1 * var_grb3d(i,j,1,1,2) & !   0-10cm
                      + 0.3 * var_grb3d(i,j,2,1,2) & !  10-40cm
                      + 0.6 * var_grb3d(i,j,3,1,2) ! 40-100cm
        ELSE                                       ! sfc temp over sea
          tsoil_ext (i,j) = var_grb2d(i,j,3,1)
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
          soiltyp_ext (i,j) = 13
        END IF

        wetsfc_ext(i,j) = var_grb3d(i,j,1,2,2)
        wetdp_ext(i,j)  = 0.1 * var_grb3d(i,j,1,2,2) & !   0-10cm
                    + 0.3 * var_grb3d(i,j,2,2,2) & !  10-40cm
                    + 0.6 * var_grb3d(i,j,3,2,2) ! 40-100cm
      END IF

      IF ( var_nr2d(4,1) == 0 ) THEN
        wetcanp_ext(i,j) = -999.0
      ELSE
        wetcanp_ext(i,j) = var_grb2d(i,j,4,1)*1.e-3     ! in meter
      END IF

      IF ( var_nr2d(6,1) == 0 ) THEN
        snowdpth_ext(i,j) = -999
      ELSE

!      Convert water equiv. of accum. snow depth (kg/m**2) to meters
!      (where 1 meter liquid water is set equivqlent to 10 meters snow).
!          0.01 = 10. (m snow/m liquid) / (1000 kg/m**3)

        snowdpth_ext(i,j) = 0.01 * var_grb2d(i,j,6,1)  ! in meters

      END IF

    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!
  nz1 = MIN(var_nr3d(1,1),nz_ext)

  IF ( var_lev3d(1,1,1) > var_lev3d(nz1,1,1) ) THEN  ! 1st level at sfc
    chklev = 1
    lvscan = 0
  ELSE
    chklev = -1
    lvscan = nz1+1
  END IF

  DO k=1,nz1
    kk = chklev * k + lvscan
    DO j=1,ny_ext
      DO i=1,nx_ext
!      p_ext  (i,j,kk) = psfc_ext(i,j)
!    :                         * float(var_lev3d(k,1))/10000.
              ! Pressure derived from sigma ccordinates

        p_ext  (i,j,kk) = 100.0 * FLOAT(var_lev3d(k,1,1)) ! Pressure
        t_ext  (i,j,kk) = var_grb3d(i,j,k,1,1)    ! Temperature (K)
        qv_ext (i,j,kk) = var_grb3d(i,j,k,2,1)    ! Spec. humidity
                                                  ! (kg/kg)
        u_ext  (i,j,kk) = var_grb3d(i,j,k,3,1)    ! u wind (m/s)
        v_ext  (i,j,kk) = var_grb3d(i,j,k,4,1)    ! v wind (m/s)
        hgt_ext(i,j,kk) = var_grb3d(i,j,k,5,1)    ! height (m)

        qc_ext (i,j,kk) = -999.
        qr_ext (i,j,kk) = -999.
        qi_ext (i,j,kk) = -999.
        qs_ext (i,j,kk) = -999.
        qh_ext (i,j,kk) = -999.
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate the height by using hytrostatical equation
!
!-----------------------------------------------------------------------
!
!   DO 220 j=1,ny_ext
!   DO 220 i=1,nx_ext
!     tema = 0.5 * ( tsfc_ext (i,j) + t_ext (i,j,1) )
!     temb = 0.5 * ( qvsfc_ext(i,j) + qv_ext(i,j,1) )
!     temc = tema * (1.0 + 0.608*temb)         ! Virtual temperature
!     tema = 0.5 * ( psfc_ext (i,j) + p_ext (i,j,1) )
!     temb = temc / tema
!
!     hgt_ext(i,j,1) = trn_ext(i,j) + govrd * temb      ! Height (m)
!  :                 * ( psfc_ext(i,j) - p_ext(i,j,1) )
!220   CONTINUE
!
!   DO 230 k=2,nz1
!   DO 230 j=1,ny_ext
!   DO 230 i=1,nx_ext
!     tema = 0.5 * ( t_ext (i,j,k-1) + t_ext (i,j,k) )
!     temb = 0.5 * ( qv_ext(i,j,k-1) + qv_ext(i,j,k) )
!     temc = tema * (1.0 + 0.608*temb)         ! Virtual temperature
!     tema = 0.5 * ( p_ext (i,j,k-1) + p_ext (i,j,k) )
!     temb = temc / tema
!
!     hgt_ext(i,j,k) = hgt_ext(i,j,k-1) + govrd * temb   ! Height (m)
!  :                 * ( p_ext (i,j,k-1) - p_ext (i,j,k) )
!230   CONTINUE
!
!
!-----------------------------------------------------------------------
!
!  Rotate winds to be relative to true north.
!  The RUC data are sent as grid-relative.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz1
!2001-05-16 GMB: Having umap & uear (or vmap & vear) point to
!the same array causes numerical errors when optimizing.
    CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k),               &
                 lon_ext,utmp,vtmp)
    u_ext(:,:,k) = utmp(:,:)
    v_ext(:,:,k) = vtmp(:,:)
  END DO

!
!-----------------------------------------------------------------------
!
!  Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
  CALL setmapr(iproj,scale,latnot,trlon)
  CALL setorig(1,x0,y0)

  istatus = 1

  DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)
  DEALLOCATE(utmp)
  DEALLOCATE(vtmp)

  RETURN
END SUBROUTINE getnmceta212
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  SUBROUTINE SWAPJ                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE swapj(nx,ny,nz,varval)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  This routine swaps values in the y direction so that the
!  reanalysis grid goes S to N in the +j direction
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Dennis Moon, SSESCO
!  06/19/1998
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!

  IMPLICIT NONE

  INTEGER*4 nx,ny,nz,i,j,k
  REAL*4 varval(nx,ny,nz), tmpval
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=1,nz
    DO i=1,nx
      DO j=1,ny/2
        tmpval = varval(i,ny-j+1,k)
        varval(i,ny-j+1,k) = varval(i,j,k)
        varval(i,j,k) = tmpval
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE swapj
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE GETNMCRUCN236               ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE getnmcrucn236(nx_ext,ny_ext,nz_ext,                          & 1,11
           dir_extd,extdname,extdopt,extdfmt,                           &
           extdinit,extdfcst,julfname,                                  &
           iproj_ext,scale_ext,                                         &
           trlon_ext,latnot_ext,x0_ext,y0_ext,                          &
           lat_ext,lon_ext,                                             &
           p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,                      &
           qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                          &
           tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext,         &
           trn_ext,psfc_ext,snowdpth_ext,                               &
           istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Reads in a GRIB file containing native coordinate RUC2 data
!  (Grid 236) and extracts/converts selected variables for use by
!  EXT2ARPS.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Eric Kemp
!  09/17/1999
!  Based on subroutine GETNMCRUC87.
!
!  MODIFICATION HISTORY:
!
!  11/01/1999 Eric Kemp
!  Corrected several bugs in variable retrievals, and added
!  retrieval of RUC2 hydrometeor fields.
!
!  1999/11/30 Gene Bassett
!  Corrected RUC2 soil moisture and changed deep soil moisture &
!  temperature to be an average from 0 to 100cm.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    dir_extd      Directory name for external file
!    extdname      Prefix string of external file name
!    extdopt       Option of external data sources
!    extdfmt       Option of external data format

!    extdinit      Initialized time in mm-dd-yyyy:hh:mm:ss format
!    extdfcst      Forecast hour in HHH:MM:SS format
!    julfname      File name in yyjjjhhmm format
!
!  OUTPUT:
!
!    iproj_ext     Map projection number of external data
!    scale_ext     Scale factor of external data
!    trlon_ext     True longitude of external data (degrees E)
!    latnot_ext(2) True latitude(s) of external data (degrees N)
!    x0_ext        x coordinate of origin of external data
!    y0_ext        y coordinate of origin of external data
!    lat_ext       latitude of external data points (degrees N)
!    lon_ext       longitude of external data points (degrees E)
!    p_ext         pressure (Pascal)
!    hgt_ext       height (m)
!    t_ext         temperature (K)
!    qv_ext        specific humidity (kg/kg)
!    u_ext         u wind component (m/s)
!    v_ext         v wind component (m/s)
!    qc_ext        Cloud water mixing ratio (kg/kg)
!    qr_ext        Rain water mixing ratio (kg/kg)
!    qi_ext        Ice mixing ratio (kg/kg)
!    qs_ext        Snow mixing ratio (kg/kg)
!    qh_ext        Hail mixing ratio (kg/kg)
!
!    tsfc_ext      Surface temperature
!    tsoil_ext     Soil temperature
!    wetsfc_ext    Top layer soil moisture (fraction)
!    wetdp_ext     Deep soil moisture (fraction)
!    wetcanp_ext   Water content on canopy
!
!    trn_ext       External terrain (m)
!    psfc_ext      Surface pressure (Pa)
!    snowdpth_ext  Snow depth (m)
!
!    istatus       status indicator
!
!  WORK ARRAYS:
!
!    var_grb3d     Arrays to store the GRIB 3-D variables:
!                  var_grb3d(nxgrb,nygrb,nzgrb,1,1) - pressure (Pa)
!                  var_grb3d(nxgrb,nygrb,nzgrb,2,1) - height (m)
!                  var_grb3d(nxgrb,nygrb,nzgrb,3,1) - Virtual Potential
!                                                     temperature (K)
!                  var_grb3d(nxgrb,nygrb,nzgrb,4,1) - Water vapor
!                                                     mixing ratio
!                                                     (kg/kg)
!                  var_grb3d(nxgrb,nygrb,nzgrb,5,1) - u wind (m/s)
!                  var_grb3d(nxgrb,nygrb,nzgrb,6,1) - v wind (m/s)
!                  var_grb3d(nxgrb,nygrb,nzgrb,7,1)
!                        - cloud water mixing ratio (kg/kg)
!                  var_grb3d(nxgrb,nygrb,nzgrb,8,1)
!                        - rain water mixing ratio (kg/kg)
!                  var_grb3d(nxgrb,nygrb,nzgrb,9,1)
!                        - ice mixing ratio (kg/kg)
!                  var_grb3d(nxgrb,nygrb,nzgrb,10,1)
!                        - snow mixing ratio (kg/kg)
!                  var_grb3d(nxgrb,nygrb,nzgrb,11,1)
!                        - graupel mixing ratio (kg/kg)
!                  var_grb3d(nxgrb,nygrb,nzgrb,1,2) - soil temp. (K)
!                  var_grb3d(nxgrb,nygrb,nzgrb,2,2) - soil moist.
!                                                    (fraction)
!
!    var_grb2d     Arrays to store the GRIB 2-D variables:
!                  var_grb2d(nxgrb,nygrb,1) - Canopy Water
!                                             (kg/m**2)
!                  var_grb2d(nxgrb,nygrb,2) - Snow depth (m)
!                  var_grb2d(nxgrb,nygrb,3) - Surface soil temperature (K)
!                  var_grb2d(nxgrb,nygrb,4) - Surface soil moisture
!                                             (fraction)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'gribcst.inc'

  CHARACTER (LEN=*) :: dir_extd
  CHARACTER (LEN=*) :: extdname

  INTEGER :: extdopt
  INTEGER :: extdfmt

  CHARACTER (LEN=19) :: extdinit
  CHARACTER (LEN=9) :: extdfcst
  CHARACTER (LEN=9) :: julfname
!
!
!-----------------------------------------------------------------------
!
!  External grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj_ext
  REAL :: scale_ext,trlon_ext
  REAL :: latnot_ext(2)
  REAL :: x0_ext,y0_ext
  REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
!  Output external variable arrays
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx_ext, ny_ext, nz_ext

  REAL :: lat_ext(nx_ext,ny_ext)
  REAL :: lon_ext(nx_ext,ny_ext)
  REAL :: p_ext  (nx_ext,ny_ext,nz_ext)   ! Pressure (Pascals)
  REAL :: hgt_ext(nx_ext,ny_ext,nz_ext)   ! Height (m)
  REAL :: t_ext  (nx_ext,ny_ext,nz_ext)   ! Temperature (K)
  REAL :: qv_ext (nx_ext,ny_ext,nz_ext)   ! Specific humidity (kg/kg)
  REAL :: u_ext  (nx_ext,ny_ext,nz_ext)   ! Eastward wind component
  REAL :: v_ext  (nx_ext,ny_ext,nz_ext)   ! Northward wind component
  REAL :: qc_ext (nx_ext,ny_ext,nz_ext)   ! Cloud H2O mixing ratio (kg/kg)
  REAL :: qr_ext (nx_ext,ny_ext,nz_ext)   ! Rain  H2O mixing ratio (kg/kg)
  REAL :: qi_ext (nx_ext,ny_ext,nz_ext)   ! Ice   mixing ratio (kg/kg)
  REAL :: qs_ext (nx_ext,ny_ext,nz_ext)   ! Snow  mixing ratio (kg/kg)
  REAL :: qh_ext (nx_ext,ny_ext,nz_ext)   ! Hail  mixing ratio (kg/kg)

  REAL :: tsfc_ext   (nx_ext,ny_ext)      ! Temperature at surface (K)
  REAL :: tsoil_ext  (nx_ext,ny_ext)      ! Deep soil temperature (K)
  REAL :: wetsfc_ext (nx_ext,ny_ext)      ! Surface soil moisture (fraction)
  REAL :: wetdp_ext  (nx_ext,ny_ext)      ! Deep soil moisture (fraction)
  REAL :: wetcanp_ext(nx_ext,ny_ext)      ! Canopy water amount

  REAL :: trn_ext    (nx_ext,ny_ext)          ! Geometrical heights
  REAL :: psfc_ext   (nx_ext,ny_ext)      ! Surface pressure (Pa)
  REAL :: snowdpth_ext (nx_ext,ny_ext)    ! Snow depth (m)
!
!-----------------------------------------------------------------------
!
!  Other  external variable arrays
!
!-----------------------------------------------------------------------
!
  REAL :: x_ext(nx_ext)
  REAL :: y_ext(ny_ext)

  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: var_grb2d(:,:,:,:)   ! GRIB variables
  REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
  INTEGER, ALLOCATABLE :: var_lev3d(:,:,:)  ! Levels (hybrid) for
                                            ! each 3-D variable
  REAL, ALLOCATABLE :: rcdata(:)            ! temporary data array
!
!-----------------------------------------------------------------------
!
!  Original grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj
  REAL :: scale,trlon,x0,y0
  REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: gribfile
  CHARACTER (LEN=13) :: gribtime
  INTEGER :: i,j,k
  INTEGER :: iyr,imo,iday,myr, jldy
  INTEGER :: ihr,imin,isec
  INTEGER :: ifhr,ifmin,ifsec
  INTEGER :: grbflen, len1, lenrun

  INTEGER :: m,n,nz1,max_nr2d,max_nr3d

  REAL :: pilev

  REAL :: tv_ext, tvc_ext
  REAL :: rovcp_p, cpd_p, g0_p, rd_p

  INTEGER :: chklev, lvscan, kk, jj

  REAL :: tema, temb

  REAL :: a,b

  INTEGER :: iret             ! Return flag

  REAL, ALLOCATABLE :: utmp(:,:), vtmp(:,:)
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  GRIB grid information
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=42) :: gridesc ! Grid description

  INTEGER :: iproj_grb    ! Map projection indicator
  INTEGER :: gthin        ! Indicator of whether the grid is "thinned"

  INTEGER :: ni_grb       ! Number of points along x-axis
  INTEGER :: nj_grb       ! Number of points along y-axis
  INTEGER :: np_grb       ! Total number of horizontal grid points

  INTEGER :: nk_grb       ! Number of vertical parameters
  REAL :: zk_grb(nz_ext)  ! Vertical coordinate parameters

  INTEGER :: npeq         ! Number of lat circles from pole to equator
  INTEGER :: nit(nz_ext)  ! Number of x-points for thinned grid

  REAL :: pi_grb          ! x-coordinate of pole point
  REAL :: pj_grb          ! y-coordinate of pole point
  INTEGER :: ipole        ! Projection center flag

  REAL :: di_grb          ! x-direction increment or grid length
  REAL :: dj_grb          ! y-direction increment or grid length

  REAL :: latsw           ! Latitude  of South West corner point
  REAL :: lonsw           ! Longitude of South West corner point
  REAL :: latne           ! Latitude  of North East corner point
  REAL :: lonne           ! Longitude of North East corner point

  REAL :: lattru1         ! Latitude (1st) at which projection is true
  REAL :: lattru2         ! Latitude (2nd) at which projection is true
  REAL :: lontrue         ! Longitude      at which projection is true

  REAL :: latrot          ! Latitude  of southern pole of rotation
  REAL :: lonrot          ! Longitude of southern pole of rotation
  REAL :: angrot          ! Angle of rotation

  REAL :: latstr          ! Latitude  of the pole of stretching
  REAL :: lonstr          ! Longitude of the pole of stretching
  REAL :: facstr          ! Stretching factor

  INTEGER :: scanmode     ! Scanning indicator
  INTEGER :: iscan        ! x-direction   scanning indicator
  INTEGER :: jscan        ! y-direction   scanning indicator
  INTEGER :: kscan        ! FORTRAN index scanning indicator

  INTEGER :: ires         ! Resolution direction increments indicator
  INTEGER :: iearth       ! Earth shape indicator: spherical or oblate?
  INTEGER :: icomp        ! (u,v) components decomposition indicator

  INTEGER :: jpenta       ! J-Pentagonal resolution parameter
  INTEGER :: kpenta       ! K-Pentagonal resolution parameter
  INTEGER :: mpenta       ! M-Pentagonal resolution parameter
  INTEGER :: ispect       ! Spectral representation type
  INTEGER :: icoeff       ! Spectral coefficient storage mode

  REAL :: xp_grb          ! X coordinate of sub-satellite point
  REAL :: yp_grb          ! Y coordinate of sub-satellite point
  REAL :: xo_grb          ! X coordinate of image sector origin
  REAL :: yo_grb          ! Y coordinate of image sector origin
  REAL :: zo_grb          ! Camera altitude from center of Earth

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
  ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
  ALLOCATE(rcdata(nx_ext*ny_ext))
  ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))
  ALLOCATE(utmp(nx_ext,ny_ext))
  ALLOCATE(vtmp(nx_ext,ny_ext))

  IF ( extdfcst == '         ') extdfcst='000:00:00'

  READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)')                  &
      iyr,imo,iday,ihr,imin,isec

  CALL julday(iyr,imo,iday,jldy)

  myr=MOD(iyr,100)
  ifhr=0
  ifmin=0
  ifsec=0

  READ (extdfcst,'(i3)',ERR=4,END=4) ifhr

  4     CONTINUE

  WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)')                      &
         iyr,imo,iday,ihr,'f',ifhr

  len1=LEN(dir_extd)
  grbflen=len1

  CALL strlnth( dir_extd, grbflen )

  IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
    dir_extd = '.'
    grbflen=1
  END IF

  IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
    grbflen=grbflen+1
    dir_extd(grbflen:grbflen)='/'
  END IF

  lenrun = LEN( extdname )
  CALL strlnth( extdname, lenrun )

  gribfile = dir_extd(1:grbflen)//extdname(1:lenrun)                    &
                                //'.'//gribtime(1:13)
  grbflen = grbflen + lenrun + 14
!
!-----------------------------------------------------------------------
!
!  RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
  gridtyp   = rucn236grid
  mproj_grb = rucn236proj

  n2dvars  = rucn236nvs2d
  n2dlvtps = rucn236nlvt2d

  DO m=1,n2dlvtps
    DO n=1,n2dvars
      var_id2d(n,m) = rucn236var_id2d(n,m)
    END DO
    levtyp2d(m) = rucn236levs2d(m)
  END DO

  n3dvars  = rucn236nvs3d
  n3dlvtps = rucn236nlvt3d

  DO m=1,n3dlvtps
    DO n=1,n3dvars
      var_id3d(n,m) = rucn236var_id3d(n,m)
    END DO
    levtyp3d(m) = rucn236levs3d(m)
  END DO

  CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime,        &
                gridesc, iproj_grb, gthin,                              &
                ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit,          &
                pi_grb,pj_grb,ipole, di_grb,dj_grb,                     &
                latsw,lonsw, latne,lonne,                               &
                latrot,lonrot,angrot,                                   &
                latstr,lonstr,facstr,                                   &
                lattru1,lattru2,lontrue,                                &
                scanmode, iscan,jscan,kscan,                            &
                ires,iearth,icomp,                                      &
                jpenta,kpenta,mpenta,ispect,icoeff,                     &
                xp_grb,yp_grb, xo_grb,yo_grb,zo_grb,                    &
                rcdata,var_grb2d,var_grb3d,var_lev3d,iret)

  max_nr2d = 0
  DO n=1,n2dvars
    DO m=1,n2dlvtps
      max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
    END DO
  END DO

  max_nr3d = 0
  DO n=1,n3dvars
    max_nr3d = MAX( max_nr3d, var_nr3d(n,1) )
  END DO

  IF ( max_nr3d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 3-D variables was found in the GRIB file',                  &
        'Program stopped in GETNMCRUC.'
    STOP
  END IF

  IF ( max_nr2d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 2-D variables was found in the GRIB file'
  END IF

!  write (6,'(/a7,2x,6(i7))')
!    :      'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)

!  DO 60 k=1,max_nr3d
!     write (6,'(i5,4x,6(i7))')
!    :      k,(var_lev3d(k,n,1),n=1,n3dvars)

  DO k=1,max_nr3d
    DO n=2,n3dvars
      IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
        WRITE (6,'(a)')                                                 &
            'Variables were not at the same level.',                    &
            'Program stopped in GETNMCRUC.'
        WRITE(6,*)                                                      &
            'var_lev3d(k,1,1) = ',var_lev3d(k,1,1),                     &
            'var_lev3d(k,n,1) = ',var_lev3d(k,n,1)
        STOP
      END IF
    END DO
  END DO

  IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN      ! Center in NP
    iproj_ext = 1
  ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -1
  ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN  ! Center in NP
    iproj_ext = 2
  ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -2
  ELSE IF ( iproj_grb == 1 ) THEN
    iproj_ext = 3
  ELSE IF ( iproj_grb == 0 ) THEN
    iproj_ext = 4
  ELSE
    WRITE (6,'(a)')                                                     &
        'Unknown map projection. Set to non-projection.'
    iproj_ext = 0
  END IF

  scale_ext = 1.0
  latnot_ext(1) = lattru1
  latnot_ext(2) = lattru2
  trlon_ext = lontrue

  dx_ext = di_grb
  dy_ext = dj_grb

  CALL getmapr(iproj,scale,latnot,trlon,x0,y0)
  CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
  CALL lltoxy(1,1,latsw,lonsw,x0_ext,y0_ext)

  DO i=1,nx_ext
    x_ext(i)=x0_ext+(i-1)*dx_ext
  END DO

  DO j=1,ny_ext
    y_ext(j)=y0_ext+(j-1)*dy_ext
  END DO

  CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!
!-----------------------------------------------------------------------
!
!  Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
  DO j=1,ny_ext
    DO i=1,nx_ext

      tsfc_ext(i,j) = -999.0

      trn_ext(i,j) = -999.0

      IF ( var_nr3d(1,2) == 0 ) THEN
        tsfc_ext   (i,j) = -999.0
        tsoil_ext  (i,j) = -999.0
        wetsfc_ext (i,j) = -999.0
        wetdp_ext  (i,j) = -999.0
        wetcanp_ext(i,j) = -999.0
      ELSE
!      tsfc_ext   (i,j) = var_grb3d(i,j,1,1,2)
        tsfc_ext   (i,j) = var_grb2d(i,j,3,1)
        tsoil_ext  (i,j) = 0.1 * var_grb3d(i,j,1,1,2)  & !   5cm
                     + 0.2 * var_grb3d(i,j,2,1,2)  & !  20cm
                     + 0.4 * var_grb3d(i,j,3,1,2)  & !  40cm
                     + 0.3 * var_grb3d(i,j,4,1,2)  ! 160cm
!      wetsfc_ext (i,j) = var_grb3d(i,j,1,2,2)
        wetsfc_ext (i,j) = var_grb2d(i,j,4,1)
        wetdp_ext  (i,j) = 0.1 * var_grb3d(i,j,1,2,2)  & !   5cm
                     + 0.2 * var_grb3d(i,j,2,2,2)  & !  20cm
                     + 0.4 * var_grb3d(i,j,3,2,2)  & !  40cm
                     + 0.3 * var_grb3d(i,j,4,2,2)  ! 160cm
        wetcanp_ext(i,j) = var_grb2d(i,j,1,1)*1.e-3  ! in meters
      END IF

      psfc_ext   (i,j) = -999.0

      IF ( var_nr2d(2,1) == 0 ) THEN
        snowdpth_ext(i,j) = -999
      ELSE
        snowdpth_ext(i,j) = var_grb2d(i,j,2,1)  ! in meters
      END IF

    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!

  cpd_p   = 1004.686    ! cp in RUC
  rovcp_p = 0.285714    ! rd/cp used in RUC
  g0_p    = 9.80665     ! gravity in RUC

  nz1 = MIN(var_nr3d(1,1),nz_ext)

  IF ( var_lev3d(1,1,1) < var_lev3d(nz1,1,1) ) THEN  ! 1st level at sfc
    chklev = 1
    lvscan = 0
  ELSE
    chklev = -1
    lvscan = nz1+1
  END IF

  DO k=1,nz1
    kk = chklev * k + lvscan
    DO j=1,ny_ext
      DO i=1,nx_ext

        p_ext(i,j,kk) = var_grb3d(i,j,k,1,1)       ! Pressure (Pa)

        hgt_ext(i,j,kk) = var_grb3d(i,j,k,2,1)     ! Height (m)

        u_ext(i,j,kk) = var_grb3d(i,j,k,5,1)       ! u wind (m/s)
        v_ext(i,j,kk) = var_grb3d(i,j,k,6,1)       ! v wind (m/s)

        a = REAL(100000)/var_grb3d(i,j,k,1,1)
        a = a**rovcp_p
        tvc_ext = var_grb3d(i,j,k,3,1)/a ! Virtual Temperature
        b = 0.61*var_grb3d(i,j,k,4,1)
        b = REAL(1) + b
        t_ext(i,j,kk) = tvc_ext/b ! Temperature (K)

        a = var_grb3d(i,j,k,4,1)*var_grb3d(i,j,k,1,1)
        a = a/(0.622 - var_grb3d(i,j,k,4,1))
        qv_ext(i,j,kk) = a*0.622/var_grb3d(i,j,k,1,1) ! SpecificHumidity
!
!-----------------------------------------------------------------------
!
!      Retrieve hydrometeor data.
!
!-----------------------------------------------------------------------
!
        IF (var_nr3d(7,1) == 0) THEN
          qc_ext(i,j,kk) = -999.
        ELSE
          qc_ext(i,j,kk) = var_grb3d(i,j,k,7,1)
        END IF

        IF (var_nr3d(8,1) == 0) THEN
          qr_ext(i,j,kk) = -999.
        ELSE
          qr_ext(i,j,kk) = var_grb3d(i,j,k,8,1)
        END IF

        IF (var_nr3d(9,1) == 0) THEN
          qi_ext(i,j,kk) = -999.
        ELSE
          qi_ext(i,j,kk) = var_grb3d(i,j,k,9,1)
        END IF

        IF (var_nr3d(10,1) == 0) THEN
          qs_ext(i,j,kk) = -999.
        ELSE
          qs_ext(i,j,kk) = var_grb3d(i,j,k,10,1)
        END IF

        IF (var_nr3d(11,1) == 0) THEN
          qh_ext(i,j,kk) = -999.
        ELSE
          qh_ext(i,j,kk) = var_grb3d(i,j,k,11,1)
        END IF

      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Rotate winds to be relative to true north.
!  The RUC data are sent as grid-relative.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz1
!2001-05-16 GMB: Having umap & uear (or vmap & vear) point to
!the same array causes numerical errors when optimizing.
    CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k),               &
                 lon_ext,utmp,vtmp)
    u_ext(:,:,k) = utmp(:,:)
    v_ext(:,:,k) = vtmp(:,:)
  END DO

!
!-----------------------------------------------------------------------
!
!  Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
  CALL setmapr(iproj,scale,latnot,trlon)
  CALL setorig(1,x0,y0)

  istatus = 1

  DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)
  DEALLOCATE(utmp)
  DEALLOCATE(vtmp)

  RETURN
END SUBROUTINE getnmcrucn236
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE GETNMCRUCP236               ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE getnmcrucp236(nx_ext,ny_ext,nz_ext,                          & 1,12
           dir_extd,extdname,extdopt,extdfmt,                           &
           extdinit,extdfcst,julfname,                                  &
           iproj_ext,scale_ext,                                         &
           trlon_ext,latnot_ext,x0_ext,y0_ext,                          &
           lat_ext,lon_ext,                                             &
           p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,                      &
           qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                          &
           tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext,         &
           trn_ext,psfc_ext, t_2m_ext, rh_2m_ext,                       &
           u_10m_ext, v_10m_ext, snowdpth_ext,                          &
           istatus)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Reads in a GRIB file containing isobaric RUC2 data (Grid 236)
!  and extracts/converts selected variables for use by EXT2ARPS.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Eric Kemp
!  09/17/1999
!  Based on subroutine GETNMCRUC211.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    dir_extd      Directory name for external file
!    extdname      Prefix string of external file name
!    extdopt       Option of external data sources
!    extdfmt       Option of external data format

!    extdinit      Initialized time in mm-dd-yyyy:hh:mm:ss format
!    extdfcst      Forecast hour in HHH:MM:SS format
!    julfname      File name in yyjjjhhmm format
!
!  OUTPUT:
!
!    iproj_ext     Map projection number of external data
!    scale_ext     Scale factor of external data
!    trlon_ext     True longitude of external data (degrees E)
!    latnot_ext(2) True latitude(s) of external data (degrees N)
!    x0_ext        x coordinate of origin of external data
!    y0_ext        y coordinate of origin of external data
!    lat_ext       latitude of external data points (degrees N)
!    lon_ext       longitude of external data points (degrees E)
!    p_ext         pressure (Pascal)
!    hgt_ext       height (m)
!    t_ext         temperature (K)
!    qv_ext        specific humidity (kg/kg)
!    u_ext         u wind component (m/s)
!    v_ext         v wind component (m/s)
!    qc_ext        Cloud water mixing ratio (kg/kg)
!    qr_ext        Rain water mixing ratio (kg/kg)
!    qi_ext        Ice mixing ratio (kg/kg)
!    qs_ext        Snow mixing ratio (kg/kg)
!    qh_ext        Hail mixing ratio (kg/kg)
!
!    tsfc_ext      Surface temperature
!    tsoil_ext     Soil temperature
!    wetsfc_ext    Top layer soil moisture (fraction)
!    wetdp_ext     Deep soil moisture (fraction)
!    wetcanp_ext   Water content on canopy
!
!    trn_ext       External terrain (m)
!    psfc_ext      Surface pressure (Pa)
!
!    T_2m_ext      Temperature at 2m AGL
!    rh_2m_ext     Specific Humidity at 2m AGL
!    U_10m_ext     U at 10m AGL
!    V_10m_ext     V at 10m AGL
!
!    snowdpth_ext  Snow depth (m)
!
!    istatus       status indicator
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'gribcst.inc'

  CHARACTER (LEN=*) :: dir_extd
  CHARACTER (LEN=*) :: extdname

  INTEGER :: extdopt
  INTEGER :: extdfmt

  CHARACTER (LEN=19) :: extdinit
  CHARACTER (LEN=9) :: extdfcst
  CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
!  External grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj_ext
  REAL :: scale_ext,trlon_ext
  REAL :: latnot_ext(2)
  REAL :: x0_ext,y0_ext
  REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
!  Output external variable arrays
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx_ext,ny_ext,nz_ext

  REAL :: lat_ext(nx_ext,ny_ext)
  REAL :: lon_ext(nx_ext,ny_ext)
  REAL :: p_ext  (nx_ext,ny_ext,nz_ext)   ! Pressure (Pascals)
  REAL :: hgt_ext(nx_ext,ny_ext,nz_ext)   ! Height (m)
  REAL :: t_ext  (nx_ext,ny_ext,nz_ext)   ! Temperature (K)
  REAL :: qv_ext (nx_ext,ny_ext,nz_ext)   ! Specific humidity (kg/kg)
  REAL :: u_ext  (nx_ext,ny_ext,nz_ext)   ! Eastward wind component
  REAL :: v_ext  (nx_ext,ny_ext,nz_ext)   ! Northward wind component
  REAL :: qc_ext (nx_ext,ny_ext,nz_ext)   ! Cloud H2O mixing ratio
                                          ! (kg/kg)
  REAL :: qr_ext (nx_ext,ny_ext,nz_ext)   ! Rain  H2O mixing ratio
                                          ! (kg/kg)
  REAL :: qi_ext (nx_ext,ny_ext,nz_ext)   ! Ice   mixing ratio (kg/kg)
  REAL :: qs_ext (nx_ext,ny_ext,nz_ext)   ! Snow  mixing ratio (kg/kg)
  REAL :: qh_ext (nx_ext,ny_ext,nz_ext)   ! Hail  mixing ratio (kg/kg)

  REAL :: tsfc_ext   (nx_ext,ny_ext)      ! Temperature at surface (K)
  REAL :: tsoil_ext  (nx_ext,ny_ext)      ! Deep soil temperature (K)
  REAL :: wetsfc_ext (nx_ext,ny_ext)      ! Surface soil moisture (fraction)
  REAL :: wetdp_ext  (nx_ext,ny_ext)      ! Deep soil moisture (fraction)
  REAL :: wetcanp_ext(nx_ext,ny_ext)      ! Canopy water amount

  REAL :: trn_ext    (nx_ext,ny_ext)          ! Geometrical heights
  REAL :: psfc_ext   (nx_ext,ny_ext)      ! Surface pressure (Pa)

  REAL :: t_2m_ext (nx_ext,ny_ext)
  REAL :: rh_2m_ext(nx_ext,ny_ext)
  REAL :: u_10m_ext(nx_ext,ny_ext)
  REAL :: v_10m_ext(nx_ext,ny_ext)

  REAL :: snowdpth_ext(nx_ext,ny_ext)     ! Snow depth (m)
!
!----------------------------------------------------------------------
!
!  Other  external variable arrays
!
!----------------------------------------------------------------------
!
  REAL :: x_ext(nx_ext)
  REAL :: y_ext(ny_ext)
  REAL :: z_ext(nz_ext)

  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: var_grb2d(:,:,:,:)   ! GRIB variables
  REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
  INTEGER, ALLOCATABLE :: var_lev3d(:,:,:)  ! Levels (hybrid) for
                                            ! each 3-D variable
  REAL, ALLOCATABLE :: rcdata(:)            ! temporary data array
!
!-----------------------------------------------------------------------
!
!  Original grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj
  REAL :: scale,trlon,x0,y0
  REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: gribfile
  CHARACTER (LEN=13) :: gribtime
  CHARACTER (LEN=10) :: runstr
  CHARACTER (LEN=3) :: fmtn
  INTEGER :: ichr,bchar,echar
  INTEGER :: i,j,k,ldir,ireturn
  INTEGER :: iyr,imo,iday,myr, jldy
  INTEGER :: ihr,imin,isec
  INTEGER :: ifhr,ifmin,ifsec
  INTEGER :: grbflen, len1, lenrun

  INTEGER :: m,n,nz1,max_nr2d,max_nr3d

  REAL :: qvsat, pilev, qvsatice

  REAL :: tema, temb

  INTEGER :: chklev, lvscan, kk, jj

  INTEGER :: iret             ! Return flag

  REAL, ALLOCATABLE :: utmp(:,:), vtmp(:,:)
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!
!-----------------------------------------------------------------------
!
!  GRIB grid information
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=42) :: gridesc ! Grid description

  INTEGER :: iproj_grb    ! Map projection indicator
  INTEGER :: gthin        ! Indicator of whether the grid is "thinned"

  INTEGER :: ni_grb       ! Number of points along x-axis
  INTEGER :: nj_grb       ! Number of points along y-axis
  INTEGER :: np_grb       ! Total number of horizontal grid points

  INTEGER :: nk_grb       ! Number of vertical parameters
  REAL :: zk_grb(nz_ext)  ! Vertical coordinate parameters

  INTEGER :: npeq         ! Number of lat circles from pole to equator
  INTEGER :: nit(nz_ext)  ! Number of x-points for thinned grid

  REAL :: pi_grb          ! x-coordinate of pole point
  REAL :: pj_grb          ! y-coordinate of pole point
  INTEGER :: ipole        ! Projection center flag

  REAL :: di_grb          ! x-direction increment or grid length
  REAL :: dj_grb          ! y-direction increment or grid length

  REAL :: lrb             ! y-direction increment or grid length

  REAL :: latsw           ! Latitude  of South West corner point
  REAL :: lonsw           ! Longitude of South West corner point
  REAL :: latne           ! Latitude  of North East corner point
  REAL :: lonne           ! Longitude of North East corner point

  REAL :: lattru1         ! Latitude (1st) at which projection is true
  REAL :: lattru2         ! Latitude (2nd) at which projection is true
  REAL :: lontrue         ! Longitude      at which projection is true

  REAL :: latrot          ! Latitude  of southern pole of rotation
  REAL :: lonrot          ! Longitude of southern pole of rotation
  REAL :: angrot          ! Angle of rotation

  REAL :: latstr          ! Latitude  of the pole of stretching
  REAL :: lonstr          ! Longitude of the pole of stretching
  REAL :: facstr          ! Stretching factor

  INTEGER :: scanmode     ! Scanning indicator
  INTEGER :: iscan        ! x-direction   scanning indicator
  INTEGER :: jscan        ! y-direction   scanning indicator
  INTEGER :: kscan        ! FORTRAN index scanning indicator

  INTEGER :: ires         ! Resolution direction increments indicator
  INTEGER :: iearth       ! Earth shape indicator: spherical or oblate?
  INTEGER :: icomp        ! (u,v) components decomposition indicator

  INTEGER :: jpenta       ! J-Pentagonal resolution parameter
  INTEGER :: kpenta       ! K-Pentagonal resolution parameter
  INTEGER :: mpenta       ! M-Pentagonal resolution parameter
  INTEGER :: ispect       ! Spectral representation type
  INTEGER :: icoeff       ! Spectral coefficient storage mode

  REAL :: xp_grb          ! X coordinate of sub-satellite point
  REAL :: yp_grb          ! Y coordinate of sub-satellite point
  REAL :: xo_grb          ! X coordinate of image sector origin
  REAL :: yo_grb          ! Y coordinate of image sector origin
  REAL :: zo_grb          ! Camera altitude from center of Earth
!
!
!-----------------------------------------------------------------------
!
!  Function f_qvsatl and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_qvsatl
!fpp$ expand (f_desdt)
!fpp$ expand (f_qvsatl)
!dir$ inline always f_desdt, f_qvsatl
!*$*  inline routine (f_desdt, f_qvsatl)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
  ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
  ALLOCATE(rcdata(nx_ext*ny_ext))
  ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))
  ALLOCATE(utmp(nx_ext,ny_ext))
  ALLOCATE(vtmp(nx_ext,ny_ext))

  IF ( extdfcst == '         ') extdfcst='000:00:00'

  READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)')                  &
      iyr,imo,iday,ihr,imin,isec

  CALL julday(iyr,imo,iday,jldy)

  myr=MOD(iyr,100)
  ifhr=0
  ifmin=0
  ifsec=0

  READ (extdfcst,'(i3)',ERR=4,END=4) ifhr

  4     CONTINUE

  WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)')                      &
         iyr,imo,iday,ihr,'f',ifhr

  len1=LEN(dir_extd)
  grbflen=len1

  CALL strlnth( dir_extd, grbflen )

  IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
    dir_extd = '.'
    grbflen=1
  END IF
  IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
    grbflen=grbflen+1
    dir_extd(grbflen:grbflen)='/'
  END IF

  lenrun = LEN( extdname )
  CALL strlnth( extdname, lenrun )

  gribfile =                                                            &
      dir_extd(1:grbflen)//extdname(1:lenrun)//'.'//gribtime(1:13)
  grbflen = grbflen + lenrun + 14

!
!-----------------------------------------------------------------------
!
!  RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
  gridtyp   = rucp236grid
  mproj_grb = rucp236proj

  n2dvars  = rucp236nvs2d
  n2dlvtps = rucp236nlvt2d

  DO k=1,n2dlvtps
    DO n=1,n2dvars
      var_id2d(n,k) = rucp236var_id2d(n,k)
    END DO
    levtyp2d(k) = rucp236levs2d(k)
  END DO

  n3dvars  = rucp236nvs3d
  n3dlvtps = rucp236nlvt3d

  DO m=1,n3dlvtps
    DO n=1,n3dvars
      var_id3d(n,m) = rucp236var_id3d(n,m)
    END DO
    levtyp3d(m) = rucp236levs3d(m)
  END DO

  CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime,        &
                gridesc, iproj_grb, gthin,                              &
                ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit,          &
                pi_grb,pj_grb,ipole, di_grb,dj_grb,                     &
                latsw,lonsw, latne,lonne,                               &
                latrot,lonrot,angrot,                                   &
                latstr,lonstr,facstr,                                   &
                lattru1,lattru2,lontrue,                                &
                scanmode, iscan,jscan,kscan,                            &
                ires,iearth,icomp,                                      &
                jpenta,kpenta,mpenta,ispect,icoeff,                     &
                xp_grb,yp_grb, xo_grb,yo_grb,zo_grb,                    &
                rcdata,var_grb2d,var_grb3d,var_lev3d,iret)

  max_nr2d = 0
  DO n=1,n2dvars
    DO m=1,n2dlvtps
      max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
    END DO
  END DO

  max_nr3d = 0
  DO n=1,n3dvars
    max_nr3d = MAX( max_nr3d, var_nr3d(n,1))
  END DO

  IF ( max_nr3d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 3-D variables was found in the GRIB file',                  &
        'Program stopped in GETNMCRUCP236.'
    STOP
  END IF

  IF ( max_nr2d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 2-D variables was found in the GRIB file'
  END IF

!  write (6,'(/a7,2x,6(i7))')
!    :  'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)

!  DO 60 k=1,max_nr3d
!    write (6,'(i5,4x,6(i7))')
!    :  k,(var_lev3d(k,n,1),n=1,n3dvars)

  DO k=1,max_nr3d
    DO n=2,n3dvars
      IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
        WRITE (6,'(a)')                                                 &
            'Variables were not at the same level.',                    &
            'Program stopped in GETNMCRUCP236.'
        STOP
      END IF
    END DO
  END DO

  IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN      ! Center in NP
    iproj_ext = 1
  ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -1
  ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN  ! Center in NP
    iproj_ext = 2
  ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -2
  ELSE IF ( iproj_grb == 1 ) THEN
    iproj_ext = 3
  ELSE IF ( iproj_grb == 0 ) THEN
    iproj_ext = 4
  ELSE
    WRITE (6,'(a)')                                                     &
        'Unknown map projection. Set to non-projection.'
    iproj_ext = 0
  END IF

  scale_ext = 1.0
  latnot_ext(1) = lattru1
  latnot_ext(2) = lattru2
  trlon_ext = lontrue

  dx_ext = di_grb
  dy_ext = dj_grb

  CALL getmapr(iproj,scale,latnot,trlon,x0,y0)
  CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
  CALL lltoxy(1,1,latsw,lonsw,x0_ext,y0_ext)

  DO i=1,nx_ext
    x_ext(i)=x0_ext+(i-1)*dx_ext
  END DO

  DO j=1,ny_ext
    y_ext(j)=y0_ext+(j-1)*dy_ext
  END DO

  CALL xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ext,lon_ext)
!
!-----------------------------------------------------------------------
!
!  Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
  DO j=1,ny_ext
    DO i=1,nx_ext
      IF ( var_nr2d(1,1) == 0 ) THEN
        psfc_ext   (i,j) = -999.0
      ELSE
        psfc_ext   (i,j) = var_grb2d(i,j,1,1)
      END IF

      IF ( var_nr2d(2,1) == 0 ) THEN
        trn_ext    (i,j) = -999.0
      ELSE
        trn_ext    (i,j) = var_grb2d(i,j,2,1)
      END IF

      IF( var_nr2d(1,2) == 0.) THEN
        t_2m_ext(i,j)= -999.0
      ELSE
        t_2m_ext(i,j)= var_grb2d(i,j,1,2)
      END IF

!   at this point we are reading in the relative humidity
!   later we'll convert to specific humidity

      IF( var_nr2d(2,2) == 0.) THEN
        rh_2m_ext(i,j)= -999.0
      ELSE
        rh_2m_ext(i,j)= var_grb2d(i,j,2,2)
      END IF

      IF( var_nr2d(3,2) == 0.) THEN
        u_10m_ext(i,j)= -999.0
      ELSE
        u_10m_ext(i,j)= var_grb2d(i,j,3,2)
      END IF

      IF( var_nr2d(4,2) == 0.) THEN
        v_10m_ext(i,j)= -999.0
      ELSE
        v_10m_ext(i,j)= var_grb2d(i,j,4,2)
      END IF


      tsfc_ext   (i,j) = -999.0
      tsoil_ext  (i,j) = -999.0
      wetsfc_ext (i,j) = -999.0
      wetdp_ext  (i,j) = -999.0
      wetcanp_ext(i,j) = -999.0


      IF ( var_nr2d(3,1) == 0 ) THEN
        snowdpth_ext(i,j) = -999
      ELSE
        snowdpth_ext(i,j) = var_grb2d(i,j,3,1)   ! in meters
      END IF

    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!
  nz1 = MIN(var_nr3d(1,1),nz_ext)

  IF ( var_lev3d(1,1,1) > var_lev3d(nz1,1,1) ) THEN  ! 1st level at
                                                     !sfc
    chklev = 1
    lvscan = 0
  ELSE
    chklev = -1
    lvscan = nz1+1
  END IF

  DO k=1,nz1
    kk = chklev * k + lvscan
    DO j=1,ny_ext
      DO i=1,nx_ext
        p_ext(i,j,kk) = 100.0 * FLOAT(var_lev3d(k,1,1)) ! Pressure
        t_ext(i,j,kk) = var_grb3d(i,j,k,2,1)    ! Temperature (K)
        u_ext(i,j,kk) = var_grb3d(i,j,k,4,1)    ! u wind (m/s)
        v_ext(i,j,kk) = var_grb3d(i,j,k,5,1)    ! v wind (m/s)
        hgt_ext(i,j,kk) = var_grb3d(i,j,k,1,1)    ! height (m)

!   check for portions of constant pressure grids that are below the
!   surface
! 
!       IF(((kk == 1) .OR. (p_ext(i,j,kk) > psfc_ext(i,j))) .AND.       &
!             (psfc_ext(i,j) > 0.) ) THEN
!         p_ext(i,j,kk)= psfc_ext(i,j) - 1. * (kk - 1)
!         t_ext(i,j,kk)= t_2m_ext(i,j)
!         u_ext(i,j,kk)= u_10m_ext(i,j)
!         v_ext(i,j,kk)= v_10m_ext(i,j)
!         hgt_ext(i,j,kk)= trn_ext(i,j) + 0.1 * (kk - 1)
!       END IF

        IF( (p_ext(i,j,kk) > 0.0) .AND. (t_ext(i,j,kk) > 0.0) ) THEN
          qvsat = f_qvsatl( p_ext(i,j,kk), t_ext(i,j,kk) )
          qv_ext(i,j,kk)= var_grb3d(i,j,k,3,1) * qvsat * 0.01

          IF((kk == 1) .OR. (p_ext(i,j,kk) > psfc_ext(i,j)))THEN
            qv_ext(i,j,kk)= rh_2m_ext(i,j) * qvsat * 0.01
          END IF
          qc_ext(i,j,kk)= 0.0
          qi_ext(i,j,kk)= 0.0
        ELSE
          qv_ext(i,j,kk)= 0.0
          qi_ext(i,j,kk)= 0.0
          qc_ext(i,j,kk)= 0.0
        END IF

        qr_ext (i,j,kk) = -999.
        qs_ext (i,j,kk) = -999.
        qh_ext (i,j,kk) = -999.
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Rotate winds to be relative to true north.
!  The RUCawips data are sent as grid-relative.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz1
!2001-05-16 GMB: Having umap & uear (or vmap & vear) point to
!the same array causes numerical errors when optimizing.
    CALL uvmptoe(nx_ext,ny_ext,u_ext(1,1,k),v_ext(1,1,k),               &
                 lon_ext,utmp,vtmp)
    u_ext(:,:,k) = utmp(:,:)
    v_ext(:,:,k) = vtmp(:,:)
  END DO

!
!-----------------------------------------------------------------------
!
!  Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
  CALL setmapr(iproj,scale,latnot,trlon)
  CALL setorig(1,x0,y0)

  istatus = 1

  DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)
  DEALLOCATE(utmp)
  DEALLOCATE(vtmp)

  RETURN
END SUBROUTINE getnmcrucp236
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GETNCEPAVN3                ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE getncepavn3(nx_ext,ny_ext,nz_ext,                            & 1,5
           dir_extd,extdname,extdopt,extdfmt,                           &
           extdinit,extdfcst,julfname,                                  &
           iproj_ext,scale_ext,                                         &
           trlon_ext,latnot_ext,x0_ext,y0_ext,                          &
           lat_ext,lon_ext,                                             &
           p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,                      &
           qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                          &
           tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext,         &
           snowdpth_ext,trn_ext,psfc_ext,                               &
           istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Reads a NCEP AVN GRIB (Grid #3, 1x1 degree) data file for
!  processing by ext2arps.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Donghai Wang and Yuhe Liu
!  05/19/1999
!
!  MODIFICATION HISTORY:
!    03/23/2000 (Donghai Wang)
!    Fixed some bugs and modified the code for ARPS4.5.1 version.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    dir_extd      Directory name for external file
!    extdname      Prefix string of external file name
!    extdopt       Option of external data sources
!    extdfmt       Option of external data format

!    extdinit      Initialized time in mm-dd-yyyy:hh:mm:ss format
!    extdfcst      Forecast hour in HHH:MM:SS format
!    julfname      File name in yyjjjhhmm format
!
!  OUTPUT:
!
!    iproj_ext     Map projection number of external data
!    scale_ext     Scale factor of external data
!    trlon_ext     True longitude of external data (degrees E)
!    latnot_ext(2) True latitude(s) of external data (degrees N)
!    x0_ext        x coordinate of origin of external data
!    y0_ext        y coordinate of origin of external data
!    lat_ext       latitude of external data points (degrees N)
!    lon_ext       longitude of external data points (degrees E)
!    p_ext         pressure (Pascal)
!    hgt_ext       height (m)
!    t_ext         temperature (K)
!    qv_ext        specific humidity (kg/kg)
!    u_ext         u wind component (m/s)
!    v_ext         v wind component (m/s)
!    qc_ext        Cloud water mixing ratio (kg/kg)
!    qr_ext        Rain water mixing ratio (kg/kg)
!    qi_ext        Ice mixing ratio (kg/kg)
!    qs_ext        Snow mixing ratio (kg/kg)
!    qh_ext        Hail mixing ratio (kg/kg)
!
!    tsfc_ext      Surface temperature
!    tsoil_ext     Soil temperature
!    wetsfc_ext    Top layer soil moisture
!    wetdp_ext     Deep soil moisture
!    wetcanp_ext   Water content on canopy
!
!    trn_ext       External terrain (m)
!    psfc_ext      Surface pressure (Pa)
!
!    istatus       status indicator
!
!  WORK ARRAYS:
!
!    var_grb3d     Arrays to store the GRIB 3-D variables:
!                  var_grb3d(nxgrb,nygrb,nzgrb,1,1) - Temperature (K)
!                  var_grb3d(nxgrb,nygrb,nzgrb,2,1) - Specific humidity
!                                                     (kg/kg)
!                  var_grb3d(nxgrb,nygrb,nzgrb,3,1) - u wind (m/s)
!                  var_grb3d(nxgrb,nygrb,nzgrb,4,1) - v wind (m/s)
!                  var_grb3d(nxgrb,nygrb,nzgrb,5,1) - Geopotential
!                                                     height (gpm)
!                  var_grb3d(nxgrb,nygrb,nzgrb,6,1) - Pressure vertical
!                                                     velocity (Pa/s)
!                                                     (if applied)
!                  var_grb3d(nxgrb,nygrb,nzgrb,1,2) - soil temp. (K)
!                  var_grb3d(nxgrb,nygrb,nzgrb,2,2) - vol. soil moist.
!                                                     (m**3/m**3)
!
!    var_grb2d     Arrays to store the GRIB 2-D variables:
!                  var_grb2d(nxgrb,nygrb,1) - Surface pressure (Pa)
!                  var_grb2d(nxgrb,nygrb,2) - Geopotential height (gpm)
!                  var_grb2d(nxgrb,nygrb,3) - Surface temperature (K)
!                  var_grb2d(nxgrb,nygrb,4) - Plant canopy surface
!                                             water (kg/m**2)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'gribcst.inc'

  CHARACTER (LEN=*) :: dir_extd
  CHARACTER (LEN=*) :: extdname

  INTEGER :: extdopt
  INTEGER :: extdfmt

  CHARACTER (LEN=19) :: extdinit
  CHARACTER (LEN=9) :: extdfcst
  CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
!  External grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj_ext
  REAL :: scale_ext,trlon_ext
  REAL :: latnot_ext(2)
  REAL :: x0_ext,y0_ext
  REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
!  Output external variable arrays
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx_ext,ny_ext,nz_ext

  REAL :: lat_ext(nx_ext,ny_ext)
  REAL :: lon_ext(nx_ext,ny_ext)
  REAL :: p_ext  (nx_ext,ny_ext,nz_ext)   ! Pressure (Pascals)
  REAL :: hgt_ext(nx_ext,ny_ext,nz_ext)   ! Height (m)
  REAL :: t_ext  (nx_ext,ny_ext,nz_ext)   ! Temperature (K)
  REAL :: qv_ext (nx_ext,ny_ext,nz_ext)   ! Specific humidity (kg/kg)
  REAL :: u_ext  (nx_ext,ny_ext,nz_ext)   ! Eastward wind component
  REAL :: v_ext  (nx_ext,ny_ext,nz_ext)   ! Northward wind component
  REAL :: qc_ext (nx_ext,ny_ext,nz_ext)   ! Cloud H2O mixing ratio (kg/kg)
  REAL :: qr_ext (nx_ext,ny_ext,nz_ext)   ! Rain  H2O mixing ratio (kg/kg)
  REAL :: qi_ext (nx_ext,ny_ext,nz_ext)   ! Ice   mixing ratio (kg/kg)
  REAL :: qs_ext (nx_ext,ny_ext,nz_ext)   ! Snow  mixing ratio (kg/kg)
  REAL :: qh_ext (nx_ext,ny_ext,nz_ext)   ! Hail  mixing ratio (kg/kg)

  REAL :: tsfc_ext   (nx_ext,ny_ext)      ! Temperature at surface (K)
  REAL :: tsoil_ext  (nx_ext,ny_ext)      ! Deep soil temperature (K)
  REAL :: wetsfc_ext (nx_ext,ny_ext)      ! Surface soil moisture
  REAL :: wetdp_ext  (nx_ext,ny_ext)      ! Deep soil moisture
  REAL :: wetcanp_ext(nx_ext,ny_ext)      ! Canopy water amount
  REAL :: snowdpth_ext(nx_ext,ny_ext)     ! Snow depth (m)

  REAL :: trn_ext    (nx_ext,ny_ext)      ! External terrain (m)
  REAL :: psfc_ext   (nx_ext,ny_ext)      ! Surface pressure (Pa)
!
!-----------------------------------------------------------------------
!
!  Other  external variable arrays
!
!-----------------------------------------------------------------------
!
  REAL :: x_ext(nx_ext)
  REAL :: y_ext(ny_ext)

  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  Work arrays for storing grib data
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: var_grb2d(:,:,:,:)   ! GRIB variables
  REAL, ALLOCATABLE :: var_grb3d(:,:,:,:,:) ! GRIB 3-D variables
  INTEGER, ALLOCATABLE :: var_lev3d(:,:,:)  ! Levels (hybrid) for
                                            ! each 3-D variable
  REAL, ALLOCATABLE :: rcdata(:)            ! temporary data array
!
!-----------------------------------------------------------------------
!
!  Original grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj
  REAL :: scale,trlon,x0,y0
  REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: gribfile
  CHARACTER (LEN=13) :: gribtime
  INTEGER :: i,j,k,kk
  INTEGER :: iyr,imo,iday,myr, jldy
  INTEGER :: ihr,imin,isec
  INTEGER :: ifhr,ifmin,ifsec
  INTEGER :: grbflen, len1, lenrun

  INTEGER :: m,n,nz1,max_nr2d,max_nr3d,min_nr3d,nz2

  REAL :: govrd

  INTEGER :: chklev, lvscan

  INTEGER :: iret             ! Return flag
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  GRIB grid information
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=42) :: gridesc ! Grid description

  INTEGER :: iproj_grb    ! Map projection indicator
  INTEGER :: gthin        ! Indicator of whether the grid is "thinned"

  INTEGER :: ni_grb       ! Number of points along x-axis
  INTEGER :: nj_grb       ! Number of points along y-axis
  INTEGER :: np_grb       ! Total number of horizontal grid points

  INTEGER :: nk_grb       ! Number of vertical parameters
  REAL :: zk_grb(nz_ext)  ! Vertical coordinate parameters

  INTEGER :: npeq         ! Number of lat circles from pole to equator
  INTEGER :: nit(nz_ext)  ! Number of x-points for thinned grid

  REAL :: pi_grb          ! x-coordinate of pole point
  REAL :: pj_grb          ! y-coordinate of pole point
  INTEGER :: ipole        ! Projection center flag

  REAL :: di_grb          ! x-direction increment or grid length
  REAL :: dj_grb          ! y-direction increment or grid length

  REAL :: latsw           ! Latitude  of South West corner point
  REAL :: lonsw           ! Longitude of South West corner point
  REAL :: latne           ! Latitude  of North East corner point
  REAL :: lonne           ! Longitude of North East corner point

  REAL :: lattru1         ! Latitude (1st) at which projection is true
  REAL :: lattru2         ! Latitude (2nd) at which projection is true
  REAL :: lontrue         ! Longitude      at which projection is true

  REAL :: latrot          ! Latitude  of southern pole of rotation
  REAL :: lonrot          ! Longitude of southern pole of rotation
  REAL :: angrot          ! Angle of rotation

  REAL :: latstr          ! Latitude  of the pole of stretching
  REAL :: lonstr          ! Longitude of the pole of stretching
  REAL :: facstr          ! Stretching factor

  INTEGER :: scanmode     ! Scanning indicator
  INTEGER :: iscan        ! x-direction   scanning indicator
  INTEGER :: jscan        ! y-direction   scanning indicator
  INTEGER :: kscan        ! FORTRAN index scanning indicator

  INTEGER :: ires         ! Resolution direction increments indicator
  INTEGER :: iearth       ! Earth shape indicator: spherical or oblate?
  INTEGER :: icomp        ! (u,v) components decomposition indicator

  INTEGER :: jpenta       ! J-Pentagonal resolution parameter
  INTEGER :: kpenta       ! K-Pentagonal resolution parameter
  INTEGER :: mpenta       ! M-Pentagonal resolution parameter
  INTEGER :: ispect       ! Spectral representation type
  INTEGER :: icoeff       ! Spectral coefficient storage mode

  REAL :: xp_grb          ! X coordinate of sub-satellite point
  REAL :: yp_grb          ! Y coordinate of sub-satellite point
  REAL :: xo_grb          ! X coordinate of image sector origin
  REAL :: yo_grb          ! Y coordinate of image sector origin
  REAL :: zo_grb          ! Camera altitude from center of Earth
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ALLOCATE(var_grb2d(nx_ext,ny_ext,n2dvs,n2dlvt))
  ALLOCATE(var_grb3d(nx_ext,ny_ext,nz_ext,n3dvs,n3dlvt))
  ALLOCATE(rcdata(nx_ext*ny_ext))
  ALLOCATE(var_lev3d(nz_ext,n3dvs,n3dlvt))

  IF(extdfcst == '         ') extdfcst='000:00:00'

  READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)')                  &
       iyr,imo,iday,ihr,imin,isec

  CALL julday(iyr,imo,iday,jldy)

  myr=MOD(iyr,100)
  ifhr=0
  ifmin=0
  ifsec=0

  READ(extdfcst,'(i3)',ERR=4,END=4) ifhr

  4     CONTINUE

  WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)')                      &
         iyr,imo,iday,ihr,'f',ifhr

  len1=LEN(dir_extd)
  grbflen=len1

  CALL strlnth( dir_extd, grbflen )

  IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
    dir_extd = '.'
    grbflen=1
  END IF

  IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
    grbflen=grbflen+1
    dir_extd(grbflen:grbflen)='/'
  END IF

  lenrun = LEN( extdname )
  CALL strlnth( extdname, lenrun )

  gribfile = dir_extd(1:grbflen)//extdname(1:lenrun)                    &
                                //'.'//gribtime(1:13)
  grbflen = grbflen + lenrun + 14
!
!-----------------------------------------------------------------------
!
!  RDNMCGRB reads NMC GRIB data
!
!-----------------------------------------------------------------------
!
  gridtyp  = avn3grid
  mproj_grb = avn3proj

  n2dvars  = avn3nvs2d
  n2dlvtps = avn3nlvt2d

  DO k=1,n2dlvtps
    DO n=1,n2dvars
      var_id2d(n,k) = avn3var_id2d(n,k)
    END DO
    levtyp2d(k) = avn3levs2d(k)
  END DO

  n3dvars  = avn3nvs3d
  n3dlvtps = avn3nlvt3d

  DO m=1,n3dlvtps
    DO n=1,n3dvars
      var_id3d(n,m) = avn3var_id3d(n,m)
    END DO
    levtyp3d(m) = avn3levs3d(m)
  END DO

  CALL rdnmcgrb(nx_ext,ny_ext,nz_ext,gribfile,grbflen, gribtime,        &
                gridesc, iproj_grb, gthin,                              &
                ni_grb,nj_grb,np_grb, nk_grb,zk_grb, npeq,nit,          &
                pi_grb,pj_grb,ipole, di_grb,dj_grb,                     &
                latsw,lonsw, latne,lonne,                               &
                latrot,lonrot,angrot,                                   &
                latstr,lonstr,facstr,                                   &
                lattru1,lattru2,lontrue,                                &
                scanmode, iscan,jscan,kscan,                            &
                ires,iearth,icomp,                                      &
                jpenta,kpenta,mpenta,ispect,icoeff,                     &
                xp_grb,yp_grb, xo_grb,yo_grb,zo_grb,                    &
                rcdata,var_grb2d,var_grb3d,var_lev3d,iret)

  max_nr2d = 0
  DO n=1,n2dvars
    DO m=1,n2dlvtps
      max_nr2d = MAX( max_nr2d, var_nr2d(n,m) )
    END DO
  END DO

  max_nr3d = 0
  min_nr3d = nz_ext
  DO n=1,n3dvars
    max_nr3d = MAX( max_nr3d, var_nr3d(n,1) )
    min_nr3d = MIN( min_nr3d, var_nr3d(n,1) )
  END DO

  IF ( max_nr3d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 3-D variable was found in the GRIB file',                   &
         'Program stopped in GETNCEPAVN3.'
    STOP
  END IF

  IF ( max_nr2d == 0 ) THEN
    WRITE (6,'(a)')                                                     &
        'No 2-D variables was found in the GRIB file'
  END IF

!  write (6,'(/a7,2x,6(i7))')
!    :    'Lev\\VID',(var_id3d(n,1),n=1,n3dvars)

!  DO 60 k=1,max_nr3d
!    var_lev3d(k,5,1) = var_lev3d(k,1,1)
!    var_lev3d(k,6,1) = var_lev3d(k,1,1)
!    write (6,'(/i5,4x,6(i7))')
!    :    k,(var_lev3d(k,n,1),n=1,n3dvars)

  DO k=1,min_nr3d
    DO n=2,n3dvars
      IF ( var_lev3d(k,1,1) /= var_lev3d(k,n,1) ) THEN
        WRITE (6,'(a)')                                                 &
            'Variables were not at the same level.',                    &
            'Program stopped in GETNCEPAVN3.'
        STOP
      END IF
    END DO
  END DO

  IF ( iproj_grb == 5 .AND. ipole == 0 ) THEN      ! Center in NP
    iproj_ext = 1
  ELSE IF ( iproj_grb == 5 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -1
  ELSE IF ( iproj_grb == 3 .AND. ipole == 0 ) THEN  ! Center in NP
    iproj_ext = 2
  ELSE IF ( iproj_grb == 3 .AND. ipole == 1 ) THEN  ! Center in SP
    iproj_ext = -2
  ELSE IF ( iproj_grb == 1 ) THEN
    iproj_ext = 3
  ELSE IF ( iproj_grb == 0 ) THEN
    iproj_ext = 4
  ELSE
    WRITE (6,'(a)')                                                     &
	'Unknown map projection. Set to non-projection.'
    iproj_ext = 0
  END IF

  scale_ext = 1.0
  latnot_ext(1) = lattru1
  latnot_ext(2) = lattru2
  trlon_ext = lontrue

  dx_ext = di_grb
  dy_ext = dj_grb

  DO i=1, nx_ext
    DO j=1, ny_ext
      lon_ext(i,j)= lonsw + (i-1) * dx_ext
      lat_ext(i,j)= latsw + (j-1) * dy_ext
    END DO
  END DO

  PRINT *,'LatSW = ',latsw,' LonSW = ',lonsw
!
!-----------------------------------------------------------------------
!
!  Retrieve 2-D variables
!
!-----------------------------------------------------------------------
!
  DO j=1,ny_ext
    DO i=1,nx_ext
      IF ( var_nr2d(1,1) == 0 ) THEN
        psfc_ext   (i,j) = -999.0
      ELSE
        psfc_ext   (i,j) = var_grb2d(i,j,1,1) * 100.0
      END IF
      IF ( var_nr2d(2,1) == 0 ) THEN
        trn_ext    (i,j) = -999.0
      ELSE
        trn_ext    (i,j) = var_grb2d(i,j,2,1)/g
      END IF

      IF ( var_nr3d(1,2) == 0 ) THEN
        tsfc_ext  (i,j) = -999.0
        tsoil_ext (i,j) = -999.0
        wetsfc_ext(i,j) = -999.0
        wetdp_ext (i,j) = -999.0
      ELSE
        tsfc_ext  (i,j) = var_grb2d(i,j,3,1)       ! sfc temp.

        IF ( nint(var_grb2d(i,j,5,1)) == 1 ) THEN  ! soil temp over land
          tsoil_ext (i,j) = var_grb3d(i,j,1,1,2)

          IF ( tsoil_ext (i,j) <= 200. ) THEN
            tsoil_ext (i,j) = tsfc_ext(i,j)
          END IF

          wetsfc_ext(i,j) = var_grb3d(i,j,2,2,2)
          wetdp_ext(i,j)  = var_grb3d(i,j,1,2,2)
        ELSE                                       ! sfc temp over sea

          tsoil_ext (i,j) = tsfc_ext(i,j)

          wetsfc_ext(i,j) = 1.0
          wetdp_ext(i,j)  = 1.0
        END IF

      END IF

      IF ( var_nr2d(4,1) == 0 ) THEN
        wetcanp_ext(i,j) = -999.0
      ELSE
        wetcanp_ext(i,j) = var_grb2d(i,j,4,1)*1.e-3     ! in meter
      END IF

      IF ( var_nr2d(6,1) == 0 ) THEN
        snowdpth_ext(i,j) = -999.
      ELSE

!      Convert water equiv. of accum. snow depth (kg/m**2) to meters
!      (where 1 meter liquid water is set equivqlent to 10 meters snow).
!          0.01 = 10. (m snow/m liquid) / (1000 kg/m**3)

        snowdpth_ext(i,j) = 0.01 * var_grb2d(i,j,6,1)  ! in meters

      END IF

    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Retrieve 3-D variables
!
!-----------------------------------------------------------------------
!
  nz1 = MIN(var_nr3d(1,1),nz_ext)

  IF ( var_lev3d(1,1,1) > var_lev3d(nz1,1,1) ) THEN  ! 1st level at sfc
    chklev = 1
    lvscan = 0
  ELSE
    chklev = -1
    lvscan = nz1+1
  END IF

  DO k=1,nz1
    kk = chklev * k + lvscan
    DO j=1,ny_ext
      DO i=1,nx_ext
        p_ext  (i,j,kk) = 100.0 * FLOAT(var_lev3d(k,1,1)) ! Pressure
        hgt_ext(i,j,kk) = var_grb3d(i,j,k,1,1)
        u_ext  (i,j,kk) = var_grb3d(i,j,k,2,1)    ! u wind (m/s)
        v_ext  (i,j,kk) = var_grb3d(i,j,k,3,1)    ! v wind (m/s)

        t_ext  (i,j,kk) = var_grb3d(i,j,k,4,1)    ! Temperature (K)

        qc_ext (i,j,kk) = -999.
        qr_ext (i,j,kk) = -999.
        qi_ext (i,j,kk) = -999.
        qs_ext (i,j,kk) = -999.
        qh_ext (i,j,kk) = -999.
      END DO
    END DO
  END DO

  CALL getqvs(nx_ext,ny_ext,nz1, 1,nx_ext,1,ny_ext,1,nz1,               &
              p_ext, t_ext, qv_ext )

  nz2 = MIN( nz1, min_nr3d )

  DO k=1,nz2
    kk = chklev * k + lvscan
    DO j=1,ny_ext
      DO i=1,nx_ext
        qv_ext(i,j,kk) = 0.01*var_grb3d(i,j,k,5,1)*qv_ext(i,j,kk)
      END DO
    END DO
  END DO

  IF ( nz2 < nz1 ) THEN
    DO k=nz2+1,nz1
      kk = chklev * k + lvscan
      DO j=1,ny_ext
        DO i=1,nx_ext
          qv_ext(i,j,kk) = 0.0
          var_lev3d(k,5,1) = var_lev3d(k,1,1)
        END DO
      END DO
    END DO
  END IF

  istatus = 1

  DEALLOCATE(var_grb2d,var_grb3d,rcdata,var_lev3d)

  RETURN
END SUBROUTINE getncepavn3
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GET_AVN_BIN                ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE get_avn_bin(nx_ext,ny_ext,nz_ext,                            & 1,4
           dir_extd,extdname,extdinit,extdfcst,julfname,                &
           iproj_ext,scale_ext,trlon_ext,latnot_ext,x0_ext,y0_ext,      &
           lat_ext,lon_ext,p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,      &
           qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                          &
           tsfc_ext,tsoil_ext,wetsfc_ext,wetdp_ext,wetcanp_ext,         &
           snowdpth_ext,trn_ext,psfc_ext,                               &
           istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Reads and pass out a section of NCEP AVN GRIB
!  (Grid #3, 1x1 degree) data file (created by program EXTRACT_AVN).
!
!-----------------------------------------------------------------------
!
!  AUTHOR: M. Xue
!  07/25/2000
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    dir_extd      Directory name for external file
!    extdname      Prefix string of external file name
!    extdinit      Initialized time in mm-dd-yyyy:hh:mm:ss format
!    extdfcst      Forecast hour in HHH:MM:SS format
!    julfname      File name in yyjjjhhmm format
!
!  OUTPUT:
!
!    iproj_ext     Map projection number of external data
!    scale_ext     Scale factor of external data
!    trlon_ext     True longitude of external data (degrees E)
!    latnot_ext(2) True latitude(s) of external data (degrees N)
!    x0_ext        x coordinate of origin of external data
!    y0_ext        y coordinate of origin of external data
!    lat_ext       latitude of external data points (degrees N)
!    lon_ext       longitude of external data points (degrees E)
!    p_ext         pressure (Pascal)
!    hgt_ext       height (m)
!    t_ext         temperature (K)
!    qv_ext        specific humidity (kg/kg)
!    u_ext         u wind component (m/s)
!    v_ext         v wind component (m/s)
!    qc_ext        Cloud water mixing ratio (kg/kg)
!    qr_ext        Rain water mixing ratio (kg/kg)
!    qi_ext        Ice mixing ratio (kg/kg)
!    qs_ext        Snow mixing ratio (kg/kg)
!    qh_ext        Hail mixing ratio (kg/kg)
!
!    tsfc_ext      Surface temperature
!    tsoil_ext     Soil temperature
!    wetsfc_ext    Top layer soil moisture
!    wetdp_ext     Deep soil moisture
!    wetcanp_ext   Water content on canopy
!
!    trn_ext       External terrain (m)
!    psfc_ext      Surface pressure (Pa)
!
!    istatus       status indicator
!
!  WORK ARRAYS:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx_ext,ny_ext,nz_ext

  CHARACTER (LEN=*) :: dir_extd
  CHARACTER (LEN=*) :: extdname
  CHARACTER (LEN=19) :: extdinit
  CHARACTER (LEN=9) :: extdfcst
  CHARACTER (LEN=9) :: julfname
!
!-----------------------------------------------------------------------
!
!  External grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj_ext
  REAL :: scale_ext,trlon_ext
  REAL :: latnot_ext(2)
  REAL :: x0_ext,y0_ext
  REAL :: dx_ext,dy_ext
!
!-----------------------------------------------------------------------
!
!  Output external variable arrays
!
!-----------------------------------------------------------------------
!
  REAL :: lat_ext(nx_ext,ny_ext)
  REAL :: lon_ext(nx_ext,ny_ext)
  REAL :: p_ext  (nx_ext,ny_ext,nz_ext)   ! Pressure (Pascals)
  REAL :: hgt_ext(nx_ext,ny_ext,nz_ext)   ! Height (m)
  REAL :: t_ext  (nx_ext,ny_ext,nz_ext)   ! Temperature (K)
  REAL :: qv_ext (nx_ext,ny_ext,nz_ext)   ! Specific humidity (kg/kg)
  REAL :: u_ext  (nx_ext,ny_ext,nz_ext)   ! Eastward wind component
  REAL :: v_ext  (nx_ext,ny_ext,nz_ext)   ! Northward wind component
  REAL :: qc_ext (nx_ext,ny_ext,nz_ext)   ! Cloud H2O mixing ratio (kg/kg)
  REAL :: qr_ext (nx_ext,ny_ext,nz_ext)   ! Rain  H2O mixing ratio (kg/kg)
  REAL :: qi_ext (nx_ext,ny_ext,nz_ext)   ! Ice   mixing ratio (kg/kg)
  REAL :: qs_ext (nx_ext,ny_ext,nz_ext)   ! Snow  mixing ratio (kg/kg)
  REAL :: qh_ext (nx_ext,ny_ext,nz_ext)   ! Hail  mixing ratio (kg/kg)

  REAL :: tsfc_ext   (nx_ext,ny_ext)      ! Temperature at surface (K)
  REAL :: tsoil_ext  (nx_ext,ny_ext)      ! Deep soil temperature (K)
  REAL :: wetsfc_ext (nx_ext,ny_ext)      ! Surface soil moisture
  REAL :: wetdp_ext  (nx_ext,ny_ext)      ! Deep soil moisture
  REAL :: wetcanp_ext(nx_ext,ny_ext)      ! Canopy water amount
  REAL :: snowdpth_ext(nx_ext,ny_ext)     ! Snow depth (m)

  REAL :: trn_ext    (nx_ext,ny_ext)      ! External terrain (m)
  REAL :: psfc_ext   (nx_ext,ny_ext)      ! Surface pressure (Pa)

  REAL :: temp_ext   (nx_ext,ny_ext)      ! Automatic work array

!-----------------------------------------------------------------------
!
!  Local variables
!
!-----------------------------------------------------------------------

  INTEGER :: istatus,ierr
  INTEGER :: nunit, idummy
  REAL :: rdummy
  CHARACTER (LEN=10) :: var_label

  CHARACTER (LEN=80) :: gribfile
  CHARACTER (LEN=13) :: gribtime
  INTEGER :: i,j,k
  INTEGER :: iyr,imo,iday,myr, jldy
  INTEGER :: ihr,imin,isec
  INTEGER :: ifhr,ifmin,ifsec
  INTEGER :: grbflen, len1, lenrun

  INTEGER :: iuout, ivout, ipout, ihout,itout,                          &
          iqvout, itsfcout,itsoilout,iwsfcout,iwdpout,                  &
          iwcnpout,isnowout,itrnout,ipsfcout,                                &
          iugrdout,ivgrdout,itgrdout,iptgrdout,irhgrdout,ipmslout

  INTEGER :: nx_ext_in, ny_ext_in, nz_ext_in
  REAL :: latbgn,latend,lonbgn,lonend, del_lat, del_lon
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF(extdfcst == '         ') extdfcst='000:00:00'

  READ (extdinit,'(i4,1x,i2,1x,i2,1x,i2,1x,i2,1x,i2)')                  &
       iyr,imo,iday,ihr,imin,isec

  CALL julday(iyr,imo,iday,jldy)

  myr=MOD(iyr,100)
  ifhr=0
  ifmin=0
  ifsec=0

  READ(extdfcst,'(i3)',ERR=4,END=4) ifhr

  4     CONTINUE

  WRITE (gribtime,'(i4.4,i2.2,i2.2,i2.2,a1,i2.2)')                      &
         iyr,imo,iday,ihr,'f',ifhr

  len1=len_trim(dir_extd)
  grbflen=len1

  IF( grbflen == 0 .OR. dir_extd(1:grbflen) == ' ' ) THEN
    dir_extd = '.'
    grbflen=1
  END IF

  IF( dir_extd(grbflen:grbflen) /= '/'.AND.grbflen < len1 ) THEN
    grbflen=grbflen+1
    dir_extd(grbflen:grbflen)='/'
  END IF

  lenrun = len_trim( extdname )

  gribfile = dir_extd(1:grbflen)//extdname(1:lenrun)                    &
                                //'.'//gribtime(1:13)
  grbflen = grbflen + lenrun + 14

  psfc_ext   = -999.0
  trn_ext    = -999.0
  tsfc_ext   = -999.0
  tsoil_ext  = -999.0
  wetsfc_ext = -999.0
  wetdp_ext  = -999.0
  wetcanp_ext= -999.0
  snowdpth_ext=-999.0

  u_ext   = -999.0
  v_ext   = -999.0
  p_ext   = -999.0
  hgt_ext = -999.0
  t_ext   = -999.0
  qv_ext  = -999.0

  CALL getunit( nunit)
  CALL asnctl ('NEWLOCAL', 1, ierr)
  CALL asnfile(trim(gribfile), '-F f77 -N ieee', ierr)

  OPEN(UNIT=nunit,FILE=trim(gribfile),                                  &
       STATUS='old',FORM='unformatted',IOSTAT=istatus)

  IF( istatus /= 0 ) THEN
    WRITE(6,'(1x,a,a,/1x,i3,a)')                                        &
        'Error occured when opening file ',trim(gribfile),              &
        'using FORTRAN unit ',nunit,' Program stopped.'
    STOP
  END IF

  PRINT*,'To read file ',trim(gribfile)

  READ (nunit,ERR=910,END=920) nx_ext_in,ny_ext_in, nz_ext_in

  IF( nx_ext_in /= nx_ext.OR.ny_ext_in /= ny_ext .OR. nz_ext_in /= nz_ext ) THEN
    WRITE(6,'(1x,a)')                                                   &
         ' Dimensions in GET_AVN_BIN inconsistent with data.'
    WRITE(6,'(1x,a,3I15)') ' Read were: ',                              &
         nx_ext_in, ny_ext_in, nz_ext_in
    WRITE(6,'(1x,a,3I15)') ' Expected:  ',                              &
         nx_ext, ny_ext, nz_ext
    WRITE(6,'(1x,a)')                                                   &
         ' Program aborted in GET_AVN_BIN.'
    STOP
  END IF

  READ (nunit,ERR=910,END=920) lonbgn,lonend,latbgn,latend
  READ (nunit,ERR=910,END=920) del_lon, del_lat

  READ (nunit,ERR=910,END=920)                                          &
        iuout, ivout, ipout, ihout,itout,                               &
        iqvout, itsfcout,itsoilout,iwsfcout,iwdpout,                    &
        iwcnpout,isnowout,itrnout,ipsfcout,iugrdout,                        &
        ivgrdout,itgrdout,iptgrdout,irhgrdout,ipmslout,                        &
        iproj_ext,idummy,idummy, idummy, idummy,                        &
        idummy, idummy,  idummy, idummy, idummy

  READ (nunit,ERR=910,END=920)                                          &
        scale_ext,trlon_ext,latnot_ext(1),latnot_ext(2),                &
        x0_ext, y0_ext, rdummy, rdummy, rdummy, rdummy,                 &
        rdummy, rdummy, rdummy, rdummy, rdummy,                         &
        rdummy, rdummy, rdummy, rdummy, rdummy,                         &
        rdummy, rdummy, rdummy, rdummy, rdummy,                         &
        rdummy, rdummy, rdummy, rdummy, rdummy

  IF( iuout == 1 ) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) u_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array u_ext.'
  END IF

  IF( ivout == 1  ) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) v_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array v_ext.'
  END IF

  IF( ipout == 1  ) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) p_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array p_ext.'
  END IF

  IF( ihout == 1  ) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) hgt_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array hgt_ext.'
  END IF

  IF( itout == 1  ) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) t_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array t_ext.'
  END IF

  IF( iqvout== 1  ) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) qv_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array qv_ext.'
  END IF

  IF( itsfcout==1 ) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) tsfc_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array tsfc_ext.'
  END IF

  IF( itsoilout==1) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) tsoil_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array tsoil_ext.'
  END IF

  IF( iwsfcout ==1) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) wetsfc_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array wetsfc_ext.'
  END IF

  IF( iwdpout ==1 ) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) wetdp_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array wetdp_ext.'
  END IF


  IF( iwcnpout ==1) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) wetcanp_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array wetcanp_ext.'
  END IF

  IF( isnowout ==1) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) snowdpth_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array snowdpth_ext.'
  END IF

  IF( itrnout ==1 ) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) trn_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array trn_ext.'
  END IF

  IF( ipsfcout ==1) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) psfc_ext
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array psfc_ext.'
  END IF

  IF( iugrdout ==1) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) temp_ext ! discarded
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array temp_ext.'
  END IF

  IF( ivgrdout ==1) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) temp_ext ! discarded
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array temp_ext.'
  END IF

  IF( itgrdout ==1) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) temp_ext ! discarded
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array temp_ext.'
  END IF

  IF( irhgrdout ==1) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) temp_ext ! discarded
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array temp_ext.'
  END IF

  IF( iptgrdout ==1) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) temp_ext ! discarded
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array temp_ext.'
  END IF

  IF( ipmslout ==1) THEN
    READ (nunit,ERR=910,END=920) var_label
    READ (nunit,ERR=910,END=920) temp_ext ! discarded
    WRITE(6,'(a,a,a)')'Read ',var_label,' into array temp_ext.'
  END IF

  CLOSE (UNIT=nunit)
  CALL retunit( nunit)

  PRINT*,'Finished reading file ',trim(gribfile)
  PRINT*,' '

  dx_ext = del_lon
  dy_ext = del_lat

  IF( lonend > 180.0 ) THEN
    lonend = lonend - 360.0
  END IF

  IF( lonbgn > 180.0 ) THEN
    lonbgn = lonbgn - 360.0
  END IF

  IF( lonbgn > lonend ) lonbgn = lonbgn - 360.0

  DO i=1,nx_ext
    DO j=1,ny_ext
      lon_ext(i,j)= lonbgn + (i-1) * dx_ext
      lat_ext(i,j)= latbgn + (j-1) * dy_ext - 90.0
    END DO
  END DO

  istatus = 1

  RETURN
!
!-----------------------------------------------------------------------
!
!  Error during read
!
!----------------------------------------------------------------------
!

  910   CONTINUE
  WRITE(6,'(/a/)') ' Error reading data in GET_AVN_BIN.'
  istatus =2
  RETURN
!
!-----------------------------------------------------------------------
!
!  End-of-file during read
!
!----------------------------------------------------------------------
!

  920   CONTINUE
  WRITE(6,'(/a/)') ' End of file reached in GET_AVN_BIN.'
  istatus =3
  RETURN
END SUBROUTINE get_avn_bin