PROGRAM arpsread,7
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Sample program to read history data files produced by ARPS.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!    12/12/2003: Rewrote based on arpscvt. Input file names will be 
!    specified in arpsread.input.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz,nzsoil   ! Grid dimensions.
  INTEGER :: nstyps            ! Maximum number of soil types.

  INTEGER :: hinfmt,nhisfile_max,nhisfile,lengbf,nf,lenfil
  PARAMETER (nhisfile_max=200)
  CHARACTER (LEN=132) :: grdbasfn,hisfile(nhisfile_max)
!
!-----------------------------------------------------------------------
!  Include files:
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
! INCLUDE 'grid.inc'
! INCLUDE 'indtflg.inc'
! INCLUDE 'exbc.inc'
!
!-----------------------------------------------------------------------
!  Arrays to be read in:
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: x(:)   ! The x-coord. of the physical and
                              ! computational grid. Defined at u-point.
  REAL, ALLOCATABLE :: y(:)   ! The y-coord. of the physical and
                              ! computational grid. Defined at v-point.
  REAL, ALLOCATABLE :: z(:)   ! The z-coord. of the computational
                              ! grid. Defined at w-point.

  REAL, ALLOCATABLE :: zp(:,:,:)  ! The height of the terrain.
  REAL, ALLOCATABLE :: zpsoil(:,:,:)  ! The height of the terrain.

  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(:,:,:)    ! Soil type
  INTEGER, ALLOCATABLE :: vegtyp(:,:)     ! Vegetation type
  REAL, ALLOCATABLE :: lai    (:,:)   ! Leaf Area Index
  REAL, ALLOCATABLE :: roufns (:,:)   ! Surface roughness
  REAL, ALLOCATABLE :: veg    (:,:)   ! Vegetation fraction

  REAL, ALLOCATABLE :: tsoil (:,:,:,:) ! Soil temperature (K)
  REAL, ALLOCATABLE :: qsoil (:,:,:,:) ! Soil moisture (m**3/m**3) 
  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 :: radswnet (:,:) ! Net solar radiation, SWin - SWout
  REAL, ALLOCATABLE :: radlwin  (:,:) ! Incoming longwave radiation

  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 :: reflectivity(:,:,:)  ! simulated reflectivity (dBZ)
!
!-----------------------------------------------------------------------
! Temporary working arrays:
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: tem1(:,:,:)
  REAL, ALLOCATABLE :: tem2(:,:,:)
  REAL, ALLOCATABLE :: tem3(:,:,:)
!
!-----------------------------------------------------------------------
! Misc. internal variables
!-----------------------------------------------------------------------
!
  INTEGER :: nchin

  INTEGER :: grdbas
  INTEGER :: i,j,k,ireturn

  REAL :: time

  INTEGER :: nfile

  INTEGER :: length, ierr, istatus

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!  Get the names of the input data files.
!-----------------------------------------------------------------------
!
  CALL get_input_file_names(hinfmt,grdbasfn,hisfile,nhisfile)

  lengbf = len_trim(grdbasfn)

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

  IF (nstyps <= 0) nstyps = 1
  nstyp = nstyps

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

  WRITE(6,'(4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz ,'nzsoil=',nzsoil  
 
  ALLOCATE(x      (nx))       
  ALLOCATE(y      (ny))       
  ALLOCATE(z      (nz))       
  ALLOCATE(zp     (nx,ny,nz)) 
  ALLOCATE(zpsoil (nx,ny,nzsoil)) 

  ALLOCATE(uprt (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "uprt")

  ALLOCATE(vprt (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "vprt")

  ALLOCATE(wprt (nx,ny,nz),STAT=istatus)
  CALL check_alloc_status(istatus, "wprt")

  ALLOCATE(ptprt  (nx,ny,nz)) 
  ALLOCATE(pprt   (nx,ny,nz)) 
  ALLOCATE(qvprt  (nx,ny,nz)) 
  ALLOCATE(qc     (nx,ny,nz)) 
  ALLOCATE(qr     (nx,ny,nz)) 
  ALLOCATE(qi     (nx,ny,nz)) 
  ALLOCATE(qs     (nx,ny,nz)) 
  ALLOCATE(qh     (nx,ny,nz)) 


  ALLOCATE(tke    (nx,ny,nz)) 
  ALLOCATE(kmh    (nx,ny,nz)) 
  ALLOCATE(kmv    (nx,ny,nz)) 
  ALLOCATE(ubar   (nx,ny,nz)) 
  ALLOCATE(vbar   (nx,ny,nz)) 
  ALLOCATE(wbar   (nx,ny,nz)) 
  ALLOCATE(ptbar  (nx,ny,nz)) 
  ALLOCATE(pbar   (nx,ny,nz)) 
  ALLOCATE(rhobar (nx,ny,nz)) 

  ALLOCATE(qvbar  (nx,ny,nz)) 
  ALLOCATE(u      (nx,ny,nz)) 
  ALLOCATE(v      (nx,ny,nz)) 
  ALLOCATE(w      (nx,ny,nz)) 
  ALLOCATE(qv     (nx,ny,nz)) 

  ALLOCATE(soiltyp (nx,ny,nstyps))    
  ALLOCATE(stypfrct(nx,ny,nstyps))    
  ALLOCATE(vegtyp (nx,ny))     
  ALLOCATE(lai    (nx,ny))    
  ALLOCATE(roufns (nx,ny))    
  ALLOCATE(veg    (nx,ny))    

  ALLOCATE(tsoil (nx,ny,nzsoil,0:nstyps))    
  ALLOCATE(qsoil (nx,ny,nzsoil,0:nstyps))    
  ALLOCATE(wetcanp(nx,ny,0:nstyps))    
  ALLOCATE(snowdpth(nx,ny))            

  ALLOCATE(raing(nx,ny))         
  ALLOCATE(rainc(nx,ny))         
  ALLOCATE(prcrate(nx,ny,4))     

  ALLOCATE(radsw (nx,ny))        
  ALLOCATE(rnflx (nx,ny))        
  ALLOCATE(radswnet (nx,ny))        
  ALLOCATE(radlwin (nx,ny))        
  
  ALLOCATE(radfrc(nx,ny,nz))     

  ALLOCATE(usflx (nx,ny))        
  ALLOCATE(vsflx (nx,ny))        
  ALLOCATE(ptsflx(nx,ny))        
  ALLOCATE(qvsflx(nx,ny))        
  ALLOCATE(tem1(nx,ny,nz))
  ALLOCATE(tem2(nx,ny,nz))
  ALLOCATE(tem3(nx,ny,nz))

  ALLOCATE(reflectivity(nx,ny,nz))

  x      =0.0
  y      =0.0
  z      =0.0
  zp     =0.0
  zpsoil =0.0

  uprt   =0.0
  vprt   =0.0
  wprt   =0.0
  ptprt  =0.0
  pprt   =0.0
  qvprt  =0.0
  qc     =0.0
  qr     =0.0
  qi     =0.0
  qs     =0.0
  qh     =0.0
  tke    =0.0
  kmh    =0.0
  kmv    =0.0
  ubar   =0.0
  vbar   =0.0
  wbar   =0.0
  ptbar  =0.0
  pbar   =0.0
  rhobar =0.0
  qvbar  =0.0
  u      =0.0
  v      =0.0
  w      =0.0
  qv     =0.0

  soiltyp =0.0
  stypfrct=0.0
  vegtyp =0.0
  lai    =0.0
  roufns =0.0
  veg    =0.0

  tsoil  =0.0
  qsoil  =0.0
  wetcanp=0.0
  snowdpth=0.0

  raing=0.0
  rainc=0.0
  prcrate=0.0

  radfrc=0.0
  radsw =0.0
  rnflx =0.0
  radswnet = 0.0
  radlwin = 0.0

  usflx =0.0
  vsflx =0.0
  ptsflx=0.0
  qvsflx=0.0
  tem1=0.0
  tem2=0.0
  tem3=0.0
!
!-----------------------------------------------------------------------
!  Get the name of the input data set.
!-----------------------------------------------------------------------
!
  ldirnam=LEN(dirname)
  CALL strlnth( dirname , ldirnam)

  lengbf=len_trim(grdbasfn)
  WRITE(6,'(/a,a)')' The grid/base name is ', grdbasfn(1:lengbf)

  DO nfile = 1,nhisfile

    lenfil=len_trim(hisfile(nfile))
    WRITE(6,'(/a,a,a)')                                                 &
        ' Data set ', trim(hisfile(nfile)),' to be converted.'
!
!-----------------------------------------------------------------------
!  Read all input data arrays
!-----------------------------------------------------------------------

    CALL dtaread(nx,ny,nz,nzsoil,nstyps,                                &
                 hinfmt, nchin,grdbasfn(1:lengbf),lengbf,               &
                 hisfile(nfile)(1:lenfil),lenfil,time,                  &
                 x,y,z,zp,zpsoil, 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,                &
                 tsoil, qsoil, wetcanp,snowdpth,                        &
                 raing,rainc,prcrate,                                   &
                 radfrc,radsw,rnflx,radswnet,radlwin,                   &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 ireturn, tem1, tem2, tem3)

    curtim = time

    if( ireturn /= 0 ) then 
       WRITE(6,'(1x,a,i2,/1x,a)')                                       &
      'Data read was unsuccessful. ireturn =', ireturn,' Job stopped.'
       STOP
    endif

!
!   Do whatever you want to do with the data here
!

  ENDDO

  STOP
END PROGRAM ARPSREAD