PROGRAM arpsread,10
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Sample program to read history data file produced by ARPS.
!
!  Link arpsread using the following command on a IBM RISC/6000
!  system, assuming HDF, NetCDF and Savi3D libraries are not
!  available:
!
!  f90 arpsread.f read3d.o nohdfio3d.o nonetio3d.o nosviio3d.o \
!      gradsio3d.o pakio3d.o outlib3d.o ibmlib3d.o
!
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'

  INTEGER :: nx, ny, nz
  INTEGER :: nstyps
!
!-----------------------------------------------------------------------
!
!  Arrays to be read in:
!
!-----------------------------------------------------------------------
!
  REAL, allocatable :: x(:)       ! The x-coord. of the physical and
                                  ! computational grid.
  REAL, allocatable :: y(:)       ! The y-coord. of the physical and
                                  ! computational grid.
  REAL, allocatable :: z(:)       ! The z-coord. of the computational
                                  ! grid. Defined at w-point.

  REAL, allocatable :: zp(:,:,:)  ! The height of the terrain.
  REAL, allocatable :: hterain(:,:) ! Terrain height.

  REAL, allocatable :: j1(:,:,:)  ! Coordinate transformation
                                  ! Jacobian -d(zp)/d(x)
  REAL, allocatable :: j2(:,:,:)  ! Coordinate transformation
                                  ! Jacobian -d(zp)/d(y)
  REAL, allocatable :: j3(:,:,:)  ! Coordinate transformation
                                  ! Jacobian  d(zp)/d(z)

  REAL, allocatable :: uprt   (:,:,:) ! Perturbation u-velocity (m/s)
  REAL, allocatable :: vprt   (:,:,:) ! Perturbation v-velocity (m/s)
  REAL, allocatable :: wprt   (:,:,:) ! Perturbation w-velocity (m/s)
  REAL, allocatable :: ptprt  (:,:,:) ! Perturbation potential temperature (K)
  REAL, allocatable :: pprt   (:,:,:) ! Perturbation pressure (Pascal)
  REAL, allocatable :: qvprt  (:,:,:) ! Perturbation water vapor specific
                                      ! humidity (kg/kg)
  REAL, allocatable :: qc     (:,:,:) ! Cloud water mixing ratio (kg/kg)
  REAL, allocatable :: qr     (:,:,:) ! Rain water mixing ratio (kg/kg)
  REAL, allocatable :: qi     (:,:,:) ! Cloud ice mixing ratio (kg/kg)
  REAL, allocatable :: qs     (:,:,:) ! Snow mixing ratio (kg/kg)
  REAL, allocatable :: qh     (:,:,:) ! Hail mixing ratio (kg/kg)
  REAL, allocatable :: tke    (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2)
  REAL, allocatable :: kmh    (:,:,:) ! Horizontal turb. mixing coef. for
                                      ! momentum. ( m**2/s )
  REAL, allocatable :: kmv    (:,:,:) ! Vertical turb. mixing coef. for
                                      ! momentum. ( m**2/s )
  REAL, allocatable :: ubar   (:,:,:) ! Base state u-velocity (m/s)
  REAL, allocatable :: vbar   (:,:,:) ! Base state v-velocity (m/s)
  REAL, allocatable :: wbar   (:,:,:) ! Base state w-velocity (m/s)
  REAL, allocatable :: ptbar  (:,:,:) ! Base state potential temperature (K)
  REAL, allocatable :: pbar   (:,:,:) ! Base state pressure (Pascal)
  REAL, allocatable :: rhobar (:,:,:) ! Base state air density (kg/m**3)
  REAL, allocatable :: qvbar  (:,:,:) ! Base state water vapor specific
                                      ! humidity (kg/kg)

  REAL, allocatable :: u      (:,:,:) ! Total u-velocity (m/s)
  REAL, allocatable :: v      (:,:,:) ! Total v-velocity (m/s)
  REAL, allocatable :: w      (:,:,:) ! Total w-velocity (m/s)
  REAL, allocatable :: qv     (:,:,:) ! Water vapor specific humidity (kg/kg)

  INTEGER, allocatable :: soiltyp (:,:,:) ! Soil type
  REAL, allocatable ::  stypfrct(:,:,:)   ! Fraction of soil types
  INTEGER, allocatable :: vegtyp(:,:)     ! Vegetation type
  REAL, allocatable ::   lai    (:,:)     ! Leaf Area Index
  REAL, allocatable ::   roufns (:,:)     ! Surface roughness
  REAL, allocatable ::   veg    (:,:)     ! Vegetation fraction

  REAL, allocatable ::   tsfc   (:,:,:)   ! Temperature at surface (K)
  REAL, allocatable ::   tsoil  (:,:,:)   ! Deep soil temperature (K)
  REAL, allocatable ::   wetsfc (:,:,:)   ! Surface soil moisture
  REAL, allocatable ::   wetdp  (:,:,:)   ! Deep soil moisture
  REAL, allocatable ::   wetcanp(:,:,:)   ! Canopy water amount
  REAL, allocatable ::   snowdpth(:,:)    ! Snow depth (m)

  REAL, allocatable ::raing(:,:)          ! Grid supersaturation rain
  REAL, allocatable ::rainc(:,:)          ! Cumulus convective rain
  REAL, allocatable ::prcrate(:,:,:)      ! 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, allocatable ::radfrc(:,:,:)      ! Radiation forcing (K/s)
  REAL, allocatable ::radsw (:,:)        ! Solar radiation reaching the surface
  REAL, allocatable ::rnflx (:,:)        ! Net radiation flux absorbed by surface

  REAL, allocatable ::usflx (:,:)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL, allocatable ::vsflx (:,:)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL, allocatable ::ptsflx(:,:)        ! Surface heat flux (K*kg/(m*s**2))
  REAL, allocatable ::qvsflx(:,:)        ! Surface moisture flux (kg/(m**2*s))

  REAL, allocatable ::tem1(:,:,:)        ! Work arrays
  REAL, allocatable ::tem2(:,:,:)        ! Work arrays
  REAL, allocatable ::tem3(:,:,:)        ! Work arrays
!
!-----------------------------------------------------------------------
!
!  Misc. internal variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: hinfmt, nchin, istatus
  INTEGER :: lengbf, lenfil, ireturn
  CHARACTER (LEN=80) :: filename
  CHARACTER (LEN=80) :: grdbasfn
  REAL :: time

  hinfmt = 1  ! Data format set to 1 for unformatted binary

  grdbasfn = 'may20.grbgrdbas'  ! File containing base state and grid arrays
  lengbf = len_trim(grdbasfn)

  filename = 'may20.grb003600'  ! History data file at 3600 s.
  lenfil = len_trim(filename )
!
!-----------------------------------------------------------------------
!
!  Obtain the grid dimensions from input data.
!
!-----------------------------------------------------------------------
!

  CALL get_dims_from_data(hinfmt,grdbasfn(1:lengbf),                    &
       nx,ny,nz,nstyps, ireturn)

  IF( ireturn /= 0 ) THEN
    PRINT*,'Problem occured when trying to get dimensions from data.'
    PRINT*,'Program stopped.'
    STOP
  END IF

  WRITE(6,'(3(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz

!
!
!  Allocate arrays
!
!-----------------------------------------------------------------------
!
  allocate(x      (nx),stat=istatus)
  allocate(y      (ny),stat=istatus)
  allocate(z      (nz),stat=istatus)
  allocate(zp     (nx,ny,nz),stat=istatus)
  allocate(hterain(nx,ny),stat=istatus)
  allocate(j1     (nx,ny,nz),stat=istatus)
  allocate(j2     (nx,ny,nz),stat=istatus)
  allocate(j3     (nx,ny,nz),stat=istatus)
  allocate(uprt   (nx,ny,nz),stat=istatus)
  allocate(vprt   (nx,ny,nz),stat=istatus)
  allocate(wprt   (nx,ny,nz),stat=istatus)
  allocate(ptprt  (nx,ny,nz),stat=istatus)
  allocate(pprt   (nx,ny,nz),stat=istatus)
  allocate(qvprt  (nx,ny,nz),stat=istatus)
  allocate(qc     (nx,ny,nz),stat=istatus)
  allocate(qr     (nx,ny,nz),stat=istatus)
  allocate(qi     (nx,ny,nz),stat=istatus)
  allocate(qs     (nx,ny,nz),stat=istatus)
  allocate(qh     (nx,ny,nz),stat=istatus)
  allocate(tke    (nx,ny,nz),stat=istatus)
  allocate(kmh    (nx,ny,nz),stat=istatus)
  allocate(kmv    (nx,ny,nz),stat=istatus)
  allocate(ubar   (nx,ny,nz),stat=istatus)
  allocate(vbar   (nx,ny,nz),stat=istatus)
  allocate(wbar   (nx,ny,nz),stat=istatus)
  allocate(ptbar  (nx,ny,nz),stat=istatus)
  allocate(pbar   (nx,ny,nz),stat=istatus)
  allocate(rhobar (nx,ny,nz),stat=istatus)
  allocate(qvbar  (nx,ny,nz),stat=istatus)
  allocate(u      (nx,ny,nz),stat=istatus)
  allocate(v      (nx,ny,nz),stat=istatus)
  allocate(w      (nx,ny,nz),stat=istatus)
  allocate(qv     (nx,ny,nz),stat=istatus)
  allocate(soiltyp (nx,ny,nstyps),stat=istatus)
  allocate(stypfrct(nx,ny,nstyps),stat=istatus)
  allocate(vegtyp(nx,ny),stat=istatus)
  allocate(lai    (nx,ny),stat=istatus)
  allocate(roufns (nx,ny),stat=istatus)
  allocate(veg    (nx,ny),stat=istatus)
  allocate(tsfc   (nx,ny,0:nstyps),stat=istatus)
  allocate(tsoil  (nx,ny,0:nstyps),stat=istatus)
  allocate(wetsfc (nx,ny,0:nstyps),stat=istatus)
  allocate(wetdp  (nx,ny,0:nstyps),stat=istatus)
  allocate(wetcanp(nx,ny,0:nstyps),stat=istatus)
  allocate(snowdpth(nx,ny),stat=istatus)
  allocate(raing(nx,ny),stat=istatus)
  allocate(rainc(nx,ny),stat=istatus)
  allocate(prcrate(nx,ny,4),stat=istatus)
  allocate(radfrc(nx,ny,nz),stat=istatus)
  allocate(radsw (nx,ny),stat=istatus)
  allocate(rnflx (nx,ny),stat=istatus)
  allocate(usflx (nx,ny),stat=istatus)
  allocate(vsflx (nx,ny),stat=istatus)
  allocate(ptsflx(nx,ny),stat=istatus)
  allocate(qvsflx(nx,ny),stat=istatus)
  allocate(tem1(nx,ny,nz),stat=istatus)
  allocate(tem2(nx,ny,nz),stat=istatus)
  allocate(tem3(nx,ny,nz),stat=istatus)

  CALL dtaread(nx,ny,nz,nstyps,                                         &
       hinfmt, nchin,grdbasfn(1:16),lengbf,                             &
       filename(1:16),lenfil,time,                                      &
       x,y,z,zp, uprt ,vprt ,wprt ,ptprt, pprt ,                        &
       qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                          &
       ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,                    &
       soiltyp,stypfrct,vegtyp,lai,roufns,veg,                          &
       tsfc, tsoil, wetsfc, wetdp, wetcanp,snowdpth,                    &
       raing,rainc,prcrate,                                             &
       radfrc,radsw,rnflx,                                              &
       usflx,vsflx,ptsflx,qvsflx,                                       &
       ireturn, tem1,tem2,tem3)

  curtim = time


  STOP

END PROGRAM arpsread