PROGRAM arpsensic,25
!
!##################################################################
!##################################################################
!######                                                      ######
!######                  PROGRAM ARPSENSIC                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Generate one perturbation initial condition files from two
!  sets of ARPS output files and write it out for the use in
!  ENSEMBLE forecast.
!  The idea is based on the SLAF (Scaled Lagged Average Forecast).
!   The procedure is as follows:
!    1, read file a;
!    2, read in file b;
!    3, find the difference bwtween a and b, store it in a. a=a-b
!    4, read in the base file c or use b as the control file (iread=0).
!    5, generate the perturbation  c=c+a/n (b is used as c in code)
!    6, n is input as the variable iorder; it can be ... -2,-1,+1,+2 ...
!
!  It shares with the model the include file dims.inc for
!  definition of dimensions and domain size. a,b,c have the same
!  dimensions and the same grid structure.
!
!  Parameters grdout,varout,mstout,iceout and trbout should be input
!  with the same values as in the data dump subroutines in the model.
!
!  AUTHOR: Dingchen Hou
! History: Apr. 30, 1998: developed from the framework of ARPSDIFF.
!       Sep. 15, 1999: modified to include soil variables in
!              perturbation change input to namelist format.
!       Feb-Apr, 2002: (F.KONG) Major modifications to:
!              - accommodate BGM (Breeding Fast Growing Mode) IC
!                generation, including the initial random perturbation
!                procedure (with inibred = 1 and specified iseed).
!                During the regular breeding cycles, certain scale
!                factors are calculated to control the amplitude of
!                the growing perturbations (with iensopt = 1). When
!                generating BGM IC, the two 12h forecasts from the
!                paired breeding members are specified as a and b,
!                and the analysis data (c) must be read in (iread=1).
!
!                When generating initial BGM IC, the only one data
!                needed is the analysis (both a and b)
!
!              - read in domain config directly from data
!
!-----------------------------------------------------------------------
!
!  DATA ARRAYS READ IN:
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       z coordinate of grid points in physical space (m)
!
!    u        x component of velocity (m/s)
!    v        y component of velocity (m/s)
!    w        vertical component of velocity in Cartesian
!             coordinates (m/s).
!
!    ptprt    perturbation potential temperature (K)
!    pprt     perturbation pressure (Pascal)
!
!    qv       water vapor mixing ratio (kg/kg)
!    qc       Cloud water mixing ratio (kg/kg)
!    qr       Rainwater mixing ratio (kg/kg)
!    qi       Cloud ice mixing ratio (kg/kg)
!    qs       Snow mixing ratio (kg/kg)
!    qh       Hail mixing ratio (kg/kg)
!
!    ubar     Base state x velocity component (m/s)
!    vbar     Base state y velocity component (m/s)
!    wbar     Base state z velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhobar   Base state density (kg/m**3)
!    qvbar    Base state water vapor mixing ratio (kg/kg)
!
!    soiltyp  Soil type
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Temperature at surface (K)
!    tsoil    Deep soil temperature (K)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!
!  CALCULATED DATA ARRAYS:
!
!    uprt    perturbation x component of velocity (m/s)
!    vprt    perturbation y component of velocity (m/s)
!    wprt    perturbation z component of velocity (m/s)
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
  INCLUDE 'indtflg.inc'
  
  INTEGER :: nx,ny,nz
  INTEGER :: vnx,vny,vnz
!  PARAMETER (vnx=nx,vny=ny,vnz=nz)
!
!-----------------------------------------------------------------------
!
!  Arrays to be read in:
!
!-----------------------------------------------------------------------
!
  INTEGER :: nstyps              ! Maximum number of soil types in each
                                 ! grid box
!  PARAMETER (nstyps=4)

  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 on the staggered grid.
  REAL, ALLOCATABLE :: zp(:,:,:) ! The physical height coordinate defined at
                                 ! w-point of the staggered grid.

  REAL, ALLOCATABLE :: j1(:,:,:) ! Coordinate transformation Jacobian defined
                                 ! as - d( zp )/d( x )
  REAL, ALLOCATABLE :: j2(:,:,:) ! Coordinate transformation Jacobian defined
                                 ! as - d( zp )/d( y )
  REAL, ALLOCATABLE :: j3(:,:,:) ! Coordinate transformation Jacobian defined
                                 ! as d( zp )/d( z )
  REAL, ALLOCATABLE :: hterain(:,:)      ! Terrain height.

  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
  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

  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 :: 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  (:,:)      ! Cumulus convective rain
  REAL, ALLOCATABLE :: rainc  (:,:)      ! Cumulus convective rain
!
!-----------------------------------------------------------------------
!
!  Verification Arrays
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: vx     (:)         ! The x-coord. of the physical and
                                          ! computational grid. Defined at u-point.
  REAL, ALLOCATABLE :: vy     (:)         ! The y-coord. of the physical and
                                          ! computational grid. Defined at v-point.
  REAL, ALLOCATABLE :: vz     (:)         ! The z-coord. of the computational grid.
                                          ! Defined at w-point on the staggered grid.
  REAL, ALLOCATABLE :: vzp    (:,:,:)     ! The physical height coordinate defined at
                                          ! w-point of the staggered grid.
  REAL, ALLOCATABLE :: vuprt  (:,:,:) ! Perturbation u-velocity (m/s)
  REAL, ALLOCATABLE :: vvprt  (:,:,:) ! Perturbation v-velocity (m/s)
  REAL, ALLOCATABLE :: vwprt  (:,:,:) ! Perturbation w-velocity (m/s)
  REAL, ALLOCATABLE :: vptprt (:,:,:) ! Perturbation potential temperature (K)
  REAL, ALLOCATABLE :: vpprt  (:,:,:) ! Perturbation pressure (Pascal)
  REAL, ALLOCATABLE :: vqvprt (:,:,:) ! Perturbation water vapor specific humidity
  REAL, ALLOCATABLE :: vqc    (:,:,:) ! Cloud water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: vqr    (:,:,:) ! Rain water mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: vqi    (:,:,:) ! Cloud ice mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: vqs    (:,:,:) ! Snow mixing ratio (kg/kg)
  REAL, ALLOCATABLE :: vqh    (:,:,:) ! Hail mixing ratio (kg/kg)

  REAL, ALLOCATABLE :: vtke   (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2)
  REAL, ALLOCATABLE :: vkmh   (:,:,:) ! Horizontal turb. mixing coef. for
                                      ! momentum. ( m**2/s )
  REAL, ALLOCATABLE :: vkmv   (:,:,:) ! Vertical turb. mixing coef. for
                                      ! momentum. ( m**2/s )
  REAL, ALLOCATABLE :: vubar  (:,:,:) ! Base state u-velocity (m/s)
  REAL, ALLOCATABLE :: vvbar  (:,:,:) ! Base state v-velocity (m/s)
  REAL, ALLOCATABLE :: vwbar  (:,:,:) ! Base state w-velocity (m/s)
  REAL, ALLOCATABLE :: vptbar (:,:,:) ! Base state potential temperature (K)
  REAL, ALLOCATABLE :: vpbar  (:,:,:) ! Base state pressure (Pascal)
  REAL, ALLOCATABLE :: vrhobar(:,:,:) ! Base state air density (kg/m**3)
  REAL, ALLOCATABLE :: vqvbar (:,:,:) ! Base state water vapor specific humidity

  INTEGER, ALLOCATABLE :: vsoiltyp (:,:,:)   ! Soil type
  REAL, ALLOCATABLE :: vstypfrct(:,:,:)      ! Soil type
  INTEGER, ALLOCATABLE :: vvegtyp  (:,:)     ! Vegetation type
  REAL, ALLOCATABLE :: vlai     (:,:)        ! Leaf Area Index
  REAL, ALLOCATABLE :: vroufns  (:,:)        ! Surface roughness
  REAL, ALLOCATABLE :: vveg     (:,:)        ! Vegetation fraction

  REAL, ALLOCATABLE :: vtsfc   (:,:,:)       ! Temperature at surface (K)
  REAL, ALLOCATABLE :: vtsoil  (:,:,:)       ! Deep soil temperature (K)
  REAL, ALLOCATABLE :: vwetsfc (:,:,:)       ! Surface soil moisture
  REAL, ALLOCATABLE :: vwetdp  (:,:,:)       ! Deep soil moisture
  REAL, ALLOCATABLE :: vwetcanp(:,:,:)       ! Canopy water amount
  REAL, ALLOCATABLE :: vsnowdpth(:,:)        ! Snow depth (m)

  REAL, ALLOCATABLE :: vraing(:,:)           ! Grid supersaturation rain
  REAL, ALLOCATABLE :: vrainc(:,:)           ! Cumulus convective rain
!
!-----------------------------------------------------------------------
!
!  Work Arrays
!
!-----------------------------------------------------------------------
!
  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(:,:,:)
  REAL, ALLOCATABLE :: tem2(:,:,:)
  REAL, ALLOCATABLE :: tem3(:,:,:)

  REAL, ALLOCATABLE :: vtem1(:,:,:)
  REAL, ALLOCATABLE :: vtem2(:,:,:)
  REAL, ALLOCATABLE :: vtem3(:,:,:)

  REAL, ALLOCATABLE :: xs(:)
  REAL, ALLOCATABLE :: ys(:)
  REAL, ALLOCATABLE :: zps(:,:,:)
  REAL, ALLOCATABLE :: x2d(:,:)
  REAL, ALLOCATABLE :: y2d(:,:)
  REAL, ALLOCATABLE :: lat(:,:),lon(:,:)

  REAL, ALLOCATABLE :: vxs(:)
  REAL, ALLOCATABLE :: vys(:)
  REAL, ALLOCATABLE :: vzps(:,:,:)
  REAL, ALLOCATABLE :: dxfld(:)
  REAL, ALLOCATABLE :: dyfld(:)
  REAL, ALLOCATABLE :: rdxfld(:)
  REAL, ALLOCATABLE :: rdyfld(:)

  INTEGER, ALLOCATABLE :: iloc(:,:),jloc(:,:)
  REAL, ALLOCATABLE :: zpver(:,:,:)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: fcrnam,runnmin
  CHARACTER (LEN=132) :: filename(3),grdbasfn(3)
  INTEGER :: lengbf(3),lenfil(3)
  INTEGER :: ifproj,ivproj
  REAL :: flatnot(2),vlatnot(2)
  REAL :: fscale,ftrulon,fdx,fdy,fx0,fy0
  REAL :: fctrlat,fctrlon
  REAL :: vscale,vtrulon,vdx,vdy,vx0,vy0
  REAL :: vctrlat,vctrlon
  REAL :: time,xctr,yctr
  INTEGER :: i,j,k
  INTEGER :: grdbas
  INTEGER :: iorder,hinfmt,iread,isread,iensopt,inibred,iseed
  INTEGER :: ireturn
  LOGICAL :: comcoord
  INTEGER :: nch
  INTEGER :: tsfcin,tsoilin,wsfcin,wdpin,wcanpin,snowdin
!
  INTEGER :: istatus

  REAL :: utot,utot2,usd,uscl
  REAL :: vtot,vtot2,vsd,vscl
  REAL :: wtot,wtot2,wsd
  REAL :: pttot,pttot2,ptsd,ptscl
  REAL :: qvtot,qvtot2,qvsd,qvscl
  REAL :: ptot,ptot2,psd,pscl
  REAL :: totalpoint

  NAMELIST /prtbpara/ iensopt,inibred,iseed,iorder,isread,soilinfl,iread
  NAMELIST /input_data/ hinfmt,grdbasfn,filename
  NAMELIST /outpt_data/ hdmpfmt,runnmin,grdout,basout,varout,mstout,    &
                     iceout,trbout,sfcout,rainout,snowout,filcmprs
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  WRITE(6,'(/6(/5x,a)/)')                                                 &
      '###############################################################', &
      '#                                                             #', &
      '# Welcome to ARPSENSIC, a program that reads in history files #', &
      '# generated by ARPS and produces perturbation grids variables #', &
      '#                                                             #', &
      '###############################################################'

!------------------------------------------------------------------------
!  Read namelist
!------------------------------------------------------------------------

  READ(5,prtbpara,END=999)
!  PRINT *,iensopt,inibred,iseed,iorder,isread,soilinfl,iread
  write(6,prtbpara)

  READ(5,input_data,END=999)
!  PRINT *,hinfmt, grdbasfn, filename
  write(6,input_data)

  READ(5,outpt_data,END=999)
!  PRINT *,hdmpfmt, runnmin,                                             &
!          grdout,basout,varout,mstout,                                  &
!          iceout,trbout,sfcout,rainout,snowout,filcmprs
  write(6,outpt_data)

  GO TO 1001
  999 PRINT *, 'ERROR in reading from input file,  Program Terminated'
  STOP
  1001 CONTINUE

!  WRITE(6,'(/a,i5/a,i5)')                                               &
!      ' The  scaling factor is :',iorder,                               &
!      ' and the soil field reading flag is :',isread,                   &
!      ' The 3rd set data reading flag is: ',iread
!
!  WRITE(6,'(a,i5)')' The data format flag value is: ',hinfmt

!------------------------------------------------------------------------
!
!  Get the dimensions of the input data files
!
!------------------------------------------------------------------------

  lengbf(1) = len_trim(grdbasfn(1))

  CALL get_dims_from_data(hinfmt,grdbasfn(1)(1:lengbf(1)),                    &
       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

  print*,'nstyps =', nstyps
  
   vnx=nx
   vny=ny
   vnz=nz

  totalpoint = nx*ny*nz

!-----------------------------------------------------------------------
!
!  Allocate the variables and initialize the them to zero
!
!-----------------------------------------------------------------------
  ALLOCATE(  x=0
  ALLOCATE(  y=0
  ALLOCATE(  z=0
  ALLOCATE(zp(nx,ny,nz),STAT=istatus)
  zp=0
  ALLOCATE(j1(nx,ny,nz),STAT=istatus)
  j1=0
  ALLOCATE(j2(nx,ny,nz),STAT=istatus)
  j2=0
  ALLOCATE(j3(nx,ny,nz),STAT=istatus)
  j3=0
  ALLOCATE(hterain(nx,ny),STAT=istatus)
  hterain=0
  ALLOCATE(uprt(nx,ny,nz),STAT=istatus)
  uprt=0
  ALLOCATE(vprt(nx,ny,nz),STAT=istatus)
  vprt=0
  ALLOCATE(wprt(nx,ny,nz),STAT=istatus)
  wprt=0
  ALLOCATE(ptprt(nx,ny,nz),STAT=istatus)
  ptprt=0
  ALLOCATE(pprt(nx,ny,nz),STAT=istatus)
  pprt=0
  ALLOCATE(qvprt(nx,ny,nz),STAT=istatus)
  qvprt=0
  ALLOCATE(qc(nx,ny,nz),STAT=istatus)
  qc=0
  ALLOCATE(qr(nx,ny,nz),STAT=istatus)
  qr=0
  ALLOCATE(qi(nx,ny,nz),STAT=istatus)
  qi=0
  ALLOCATE(qs(nx,ny,nz),STAT=istatus)
  qs=0
  ALLOCATE(qh(nx,ny,nz),STAT=istatus)
  qh=0
  ALLOCATE(tke(nx,ny,nz),STAT=istatus)
  tke=0
  ALLOCATE(kmh(nx,ny,nz),STAT=istatus)
  kmh=0
  ALLOCATE(kmv(nx,ny,nz),STAT=istatus)
  kmv=0
  ALLOCATE(ubar(nx,ny,nz),STAT=istatus)
  ubar=0
  ALLOCATE(vbar(nx,ny,nz),STAT=istatus)
  vbar=0
  ALLOCATE(wbar(nx,ny,nz),STAT=istatus)
  wbar=0
  ALLOCATE(ptbar(nx,ny,nz),STAT=istatus)
  ptbar=0
  ALLOCATE(pbar(nx,ny,nz),STAT=istatus)
  pbar=0
  ALLOCATE(rhobar(nx,ny,nz),STAT=istatus)
  rhobar=0
  ALLOCATE(qvbar(nx,ny,nz),STAT=istatus)
  qvbar=0
  ALLOCATE(soiltyp(nx,ny,nstyps),STAT=istatus)
  soiltyp=0
  ALLOCATE(stypfrct(nx,ny,nstyps),STAT=istatus)
  stypfrct=0
  ALLOCATE(vegtyp(nx,ny),STAT=istatus)
  vegtyp=0
  ALLOCATE(lai(nx,ny),STAT=istatus)
  lai=0
  ALLOCATE(roufns(nx,ny),STAT=istatus)
  roufns=0
  ALLOCATE(veg(nx,ny),STAT=istatus)
  veg=0
  ALLOCATE(tsfc(nx,ny,0:nstyps),STAT=istatus)
  tsfc=0
  ALLOCATE(tsoil(nx,ny,0:nstyps),STAT=istatus)
  tsoil=0
  ALLOCATE(wetsfc(nx,ny,0:nstyps),STAT=istatus)
  wetsfc=0
  ALLOCATE(wetdp(nx,ny,0:nstyps),STAT=istatus)
  wetdp=0
  ALLOCATE(wetcanp(nx,ny,0:nstyps),STAT=istatus)
  wetcanp=0
  ALLOCATE(snowdpth(nx,ny),STAT=istatus)
  snowdpth=0
  ALLOCATE(raing(nx,ny),STAT=istatus)
  raing=0
  ALLOCATE(rainc(nx,ny),STAT=istatus)
  rainc=0
  ALLOCATE(vx(vnx),STAT=istatus)
  vx=0
  ALLOCATE(vy(vny),STAT=istatus)
  vy=0
  ALLOCATE(vz(vnz),STAT=istatus)
  vz=0
  ALLOCATE(vzp(vnx,vny,vnz),STAT=istatus)
  vzp=0
  ALLOCATE(vuprt(vnx,vny,vnz),STAT=istatus)
  vuprt=0
  ALLOCATE(vvprt(vnx,vny,vnz),STAT=istatus)
  vvprt=0
  ALLOCATE(vwprt(vnx,vny,vnz),STAT=istatus)
  vwprt=0
  ALLOCATE(vptprt(vnx,vny,vnz),STAT=istatus)
  vptprt=0
  ALLOCATE(vpprt(vnx,vny,vnz),STAT=istatus)
  vpprt=0
  ALLOCATE(vqvprt(vnx,vny,vnz),STAT=istatus)
  vqvprt=0
  ALLOCATE(vqc(vnx,vny,vnz),STAT=istatus)
  vqc=0
  ALLOCATE(vqr(vnx,vny,vnz),STAT=istatus)
  vqr=0
  ALLOCATE(vqi(vnx,vny,vnz),STAT=istatus)
  vqi=0
  ALLOCATE(vqs(vnx,vny,vnz),STAT=istatus)
  vqs=0
  ALLOCATE(vqh(vnx,vny,vnz),STAT=istatus)
  vqh=0
  ALLOCATE(vtke(nx,ny,nz),STAT=istatus)
  vtke=0
  ALLOCATE(vkmh(nx,ny,nz),STAT=istatus)
  vkmh=0
  ALLOCATE(vkmv(nx,ny,nz),STAT=istatus)
  vkmv=0
  ALLOCATE(vubar(vnx,vny,vnz),STAT=istatus)
  vubar=0
  ALLOCATE(vvbar(vnx,vny,vnz),STAT=istatus)
  vvbar=0
  ALLOCATE(vwbar(vnx,vny,vnz),STAT=istatus)
  vwbar=0
  ALLOCATE(vptbar(vnx,vny,vnz),STAT=istatus)
  vptbar=0
  ALLOCATE(vpbar(vnx,vny,vnz),STAT=istatus)
  vpbar=0
  ALLOCATE(vrhobar(vnx,vny,vnz),STAT=istatus)
  vrhobar=0
  ALLOCATE(vqvbar(vnx,vny,vnz),STAT=istatus)
  vqvbar=0
  ALLOCATE(vsoiltyp(vnx,vny,nstyps),STAT=istatus)
  vsoiltyp=0
  ALLOCATE(vstypfrct(vnx,vny,nstyps),STAT=istatus)
  vstypfrct=0
  ALLOCATE(vvegtyp(vnx,vny),STAT=istatus)
  vvegtyp=0
  ALLOCATE(vlai(vnx,vny),STAT=istatus)
  vlai=0
  ALLOCATE(vroufns(vnx,vny),STAT=istatus)
  vroufns=0
  ALLOCATE(vveg(vnx,vny),STAT=istatus)
  vveg=0
  ALLOCATE(vtsfc(vnx,vny,0:nstyps),STAT=istatus)
  vtsfc=0
  ALLOCATE(vtsoil(vnx,vny,0:nstyps),STAT=istatus)
  vtsoil=0
  ALLOCATE(vwetsfc(vnx,vny,0:nstyps),STAT=istatus)
  vwetsfc=0
  ALLOCATE(vwetdp(vnx,vny,0:nstyps),STAT=istatus)
  vwetdp=0
  ALLOCATE(vwetcanp(vnx,vny,0:nstyps),STAT=istatus)
  vwetcanp=0
  ALLOCATE(vsnowdpth(nx,ny),STAT=istatus)
  vsnowdpth=0
  ALLOCATE(vraing(vnx,vny),STAT=istatus)
  vraing=0
  ALLOCATE(vrainc(vnx,vny),STAT=istatus)
  vrainc=0
  ALLOCATE(prcrate(nx,ny,4),STAT=istatus)
  prcrate=0
  ALLOCATE(radfrc(nx,ny,nz),STAT=istatus)
  radfrc=0
  ALLOCATE(radsw(nx,ny),STAT=istatus)
  radsw=0
  ALLOCATE(rnflx(nx,ny),STAT=istatus)
  rnflx=0
  ALLOCATE(usflx(nx,ny),STAT=istatus)
  usflx=0
  ALLOCATE(vsflx(nx,ny),STAT=istatus)
  vsflx=0
  ALLOCATE(ptsflx(nx,ny),STAT=istatus)
  ptsflx=0
  ALLOCATE(qvsflx(nx,ny),STAT=istatus)
  qvsflx=0
  ALLOCATE(tem1(nx,ny,nz),STAT=istatus)
  tem1=0
  ALLOCATE(tem2(nx,ny,nz),STAT=istatus)
  tem2=0
  ALLOCATE(tem3(nx,ny,nz),STAT=istatus)
  tem3=0
  ALLOCATE(vtem1(vnx,vny,vnz),STAT=istatus)
  vtem1=0
  ALLOCATE(vtem2(vnx,vny,vnz),STAT=istatus)
  vtem2=0
  ALLOCATE(vtem3(vnx,vny,vnz),STAT=istatus)
  vtem3=0
  ALLOCATE(xs(nx),STAT=istatus)
  xs=0
  ALLOCATE(ys(ny),STAT=istatus)
  ys=0
  ALLOCATE(zps(nx,ny,nz),STAT=istatus)
  zps=0
  ALLOCATE(x2d(nx,ny),STAT=istatus)
  x2d=0
  ALLOCATE(y2d(nx,ny),STAT=istatus)
  y2d=0
  ALLOCATE(lat(nx,ny),STAT=istatus)
  ALLOCATE(lon(nx,ny),STAT=istatus)
  lat=0
  lon=0
  ALLOCATE(vxs(nx),STAT=istatus)
  vxs=0
  ALLOCATE(vys(ny),STAT=istatus)
  vys=0
  ALLOCATE(vzps(vnx,vny,vnz),STAT=istatus)
  vzps=0
  ALLOCATE(dxfld(vnx),STAT=istatus)
  dxfld=0
  ALLOCATE(dyfld(vny),STAT=istatus)
  dyfld=0
  ALLOCATE(rdxfld(vnx),STAT=istatus)
  rdxfld=0
  ALLOCATE(rdyfld(vny),STAT=istatus)
  rdyfld=0
  ALLOCATE(iloc(nx,ny),STAT=istatus)
  ALLOCATE(jloc(nx,ny),STAT=istatus)
  iloc=0
  jloc=0
  ALLOCATE(zpver(nx,ny,vnz),STAT=istatus)
  zpver=0

  
  101  CONTINUE
!
!-----------------------------------------------------------------------
!
!  Get the name of the grid/base data set.
!
!-----------------------------------------------------------------------
!

  DO i=1,3

    WRITE(6,'(/a,i5,a)')' For data set ',i,':'
    lengbf(i) = len_trim(grdbasfn(i))
!    lengbf(i)=LEN(grdbasfn(i))
!    CALL strlnth( grdbasfn(i), lengbf(i))
    WRITE(6,'(/a,a)')' The grid/base name is ',                         &
                        grdbasfn(i)(1:lengbf(i))

    lenfil(i) = len_trim(filename(i))
!    lenfil(i) = LEN(filename(i))
!    CALL strlnth( filename(i), lenfil(i))
    WRITE(6,'(/a,/1x,a)')' The data set name is ',                      &
                           filename(i)(1:lenfil(i))
  END DO
!
  WRITE(6,'(/5x,a,a)') 'The output run name is: ', runnmin
  WRITE(6,'(a,i5)')' The  output data format flag  is: ',hdmpfmt
!-----------------------------------------------------------------------
!
!  Read all input data arrays (a - forecast data)
!
!-----------------------------------------------------------------------
!
  CALL dtaread(nx,ny,nz,nstyps,                                         &
               hinfmt,nch,grdbasfn(1)(1:lengbf(1)),lengbf(1),           &
               filename(1)(1:lenfil(1)),lenfil(1),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)

  IF (isread == 1) THEN
    CALL readsoil(nx,ny,nstyps,soilinfl,dx,dy,                          &
               mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,     &
                  tsfcin,tsoilin,wsfcin,wdpin,wcanpin,snowdin,          &
                  tsfc,  tsoil,  wetsfc,wetdp,wetcanp,snowdpth)
  END IF
!
!-----------------------------------------------------------------------
!
!  ireturn = 0 for a successful read
!
!-----------------------------------------------------------------------
!
  IF( ireturn /= 0 ) THEN   ! failed read
    PRINT *,'Problem reading forecast data, forced to STOP'
    STOP
  END IF                      ! successful read

  IF (inibred == 1) THEN

  utot=0.0
  vtot=0.0
  wtot=0.0
  pttot=0.0
  qvtot=0.0
  ptot=0.0
  utot2=0.0
  vtot2=0.0
  wtot2=0.0
  pttot2=0.0
  qvtot2=0.0
  ptot2=0.0
  do k=1,nz
  do i=1,nx
  do j=1,ny
  utot=utot+uprt(i,j,k)
  vtot=vtot+vprt(i,j,k)
  wtot=wtot+wprt(i,j,k)
  pttot=pttot+ptprt(i,j,k)
  qvtot=qvtot+qvprt(i,j,k)
  ptot=ptot+pprt(i,j,k)
  utot2=utot2+uprt(i,j,k)*uprt(i,j,k)
  vtot2=vtot2+vprt(i,j,k)*vprt(i,j,k)
  wtot2=wtot2+wprt(i,j,k)*wprt(i,j,k)
  pttot2=pttot2+ptprt(i,j,k)*ptprt(i,j,k)
  qvtot2=qvtot2+qvprt(i,j,k)*qvprt(i,j,k)
  ptot2=ptot2+pprt(i,j,k)*pprt(i,j,k)
  enddo
  enddo
  enddo
  utot=utot/totalpoint
  vtot=utot/totalpoint
  wtot=wtot/totalpoint
  pttot=pttot/totalpoint
  qvtot=qvtot/totalpoint
  ptot=ptot/totalpoint
  usd=sqrt(utot2/totalpoint-utot*utot)
  vsd=sqrt(vtot2/totalpoint-vtot*vtot)
  wsd=sqrt(wtot2/totalpoint-wtot*wtot)
  ptsd=sqrt(pttot2/totalpoint-pttot*pttot)
  qvsd=sqrt(qvtot2/totalpoint-qvtot*qvtot)
  psd=sqrt(ptot2/totalpoint-ptot*ptot)

  print *,'totalpoint=',totalpoint
  print *,'usd=',usd,' utot=',utot,' utot2=',utot2/totalpoint
  print *,'vsd=',vsd,' vtot=',vtot,' vtot2=',vtot2/totalpoint
  print *,'wsd=',wsd,' wtot=',wtot,' wtot2=',wtot2/totalpoint
  print *,'ptsd=',ptsd,' pttot=',pttot,' pttot2=',pttot2/totalpoint
  print *,'qvsd=',qvsd,' qvtot=',qvtot,' qvtot2=',qvtot2/totalpoint
  print *,'psd=',psd,' ptot=',ptot,' ptot2=',ptot2/totalpoint

  END IF

  curtim=time
  fcrnam=runname
  ifproj=mapproj
  fscale=sclfct
  flatnot(1)=trulat1
  flatnot(2)=trulat2
  ftrulon=trulon
  fdx=x(3)-x(2)
  fdy=y(3)-y(2)
  fctrlat=ctrlat
  fctrlon=ctrlon
  CALL setmapr(ifproj,fscale,flatnot,ftrulon)
  CALL lltoxy(1,1,fctrlat,fctrlon,xctr,yctr)
  fx0=xctr-fdx*((nx-3)/2)
  fy0=yctr-fdy*((ny-3)/2)
  CALL setorig(1,fx0,fy0)
!
!-----------------------------------------------------------------------
!
!  Get verification data (b - analysis (or gridded obs.) data)
!
!-----------------------------------------------------------------------
!
!  Set the gridread parameter to 0 so that the verification
!  grid/base file will be read.
!
!-----------------------------------------------------------------------
!
  CALL setgbrd (0)
!
!-----------------------------------------------------------------------
!
!  Read in the verification data.
!
!-----------------------------------------------------------------------
!
  CALL dtaread(vnx,vny,vnz,nstyps,                                      &
               hinfmt,nch,grdbasfn(2)(1:lengbf(2)),lengbf(2),           &
               filename(2)(1:lenfil(2)),lenfil(2),time,                 &
               vx,vy,vz,vzp, vuprt ,vvprt ,vwprt ,vptprt, vpprt ,       &
               vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv,         &
               vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar,     &
               vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg,            &
               vtsfc,vtsoil,vwetsfc,vwetdp,vwetcanp,vsnowdpth,          &
       raing,rainc,prcrate,                                             &
       radfrc,radsw,rnflx,                                              &
       usflx,vsflx,ptsflx,qvsflx,                                       &
               ireturn, vtem1,vtem2,vtem3)
  IF (isread == 1) THEN
    CALL readsoil(nx,ny,nstyps,soilinfl,dx,dy,                          &
               mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,     &
                  tsfcin,tsoilin,wsfcin,wdpin,wcanpin,snowdin,          &
                  vtsfc, vtsoil, vwetsfc,vwetdp,vwetcanp,vsnowdpth)
  END IF
  IF( ireturn /= 0 ) THEN   ! failed read
    PRINT *,'Problem reading verification data, forced to STOP'
    STOP
  END IF                      ! successful read
!  IF (vnx /= nx.OR.vny /= ny.OR.vnz /= nz) THEN
!    PRINT *,'nx,ny,nz,','=/=','vnx,vny,vnz'
!    PRINT *,nx,ny,nz,'=/=',vnx,vny,vnz
!    PRINT *, 'forced to stop'
!    STOP
!  END IF            ! don't need
  ivproj=mapproj
  vscale=sclfct
  vlatnot(1)=trulat1
  vlatnot(2)=trulat2
  vtrulon=trulon
  vdx=vx(3)-vx(2)
  vdy=vy(3)-vy(2)
  vctrlat=ctrlat
  vctrlon=ctrlon
  CALL setmapr(ivproj,vscale,vlatnot,vtrulon)
  CALL lltoxy(1,1,vctrlat,vctrlon,xctr,yctr)
  vx0=xctr-vdx*((vnx-3)/2)
  vy0=yctr-vdy*((vny-3)/2)
  CALL setorig(1,vx0,vy0)

  IF (fx0 == vx0.AND.fy0 == vy0.AND.                                    &
        flatnot(1) == vlatnot(1).AND.flatnot(2) == vlatnot(2).AND.      &
        ftrulon == vtrulon.AND.ifproj == ivproj .AND.                   &
        fscale == vscale ) THEN
    PRINT *, 'Grids 1 and 2 shares a common coordinate system'
  ELSE
    PRINT *, 'Grids 1/2 are  different, CHECK the PROFRAM or data'
    PRINT *, 'Forced to STOP'
    STOP
  END IF

  IF (inibred == 0) THEN
!
!-----------------------------------------------------------------------
!
!  Find   difference = forecast - verification
!       a=a-b (note the forth parameter is 0)
!  To reduce memory requirements, the difference fields are
!  written to the same arrays as the interpolated fields.
!
!-----------------------------------------------------------------------
!
  CALL prtfield(nx,ny,nz,0,                                             &
                uprt, vprt, wprt, ptprt, pprt,                          &
                qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                 &
                ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,           &
                tsfc,tsoil,wetsfc,wetdp,wetcanp,                        &
                raing,rainc,                                            &
                vuprt, vvprt, vwprt, vptprt, vpprt,                     &
                vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv,        &
                vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar,    &
                vtsfc,vtsoil,vwetsfc,vwetdp,vwetcanp,                   &
                vraing,vrainc,                                          &
                uprt, vprt, wprt, ptprt, pprt,                          &
                qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,tem1,            &
                tsfc,tsoil,wetsfc,wetdp,wetcanp,                        &
                raing,rainc,                                            &
                ireturn )

  IF (iensopt == 1) THEN   ! scaling the perturbation in breeding IC

  utot2=0.0
  vtot2=0.0
  pttot2=0.0
  qvtot2=0.0
  ptot2=0.0
  do k=1,nz
  do i=1,nx
  do j=1,ny
  utot2=utot2+uprt(i,j,k)*uprt(i,j,k)
  vtot2=vtot2+vprt(i,j,k)*vprt(i,j,k)
  pttot2=pttot2+ptprt(i,j,k)*ptprt(i,j,k)
  qvtot2=qvtot2+qvprt(i,j,k)*qvprt(i,j,k)
  ptot2=ptot2+pprt(i,j,k)*pprt(i,j,k)
  enddo
  enddo
  enddo
  usd=sqrt(utot2/totalpoint)
  vsd=sqrt(vtot2/totalpoint)
  ptsd=sqrt(pttot2/totalpoint)
  qvsd=sqrt(qvtot2/totalpoint)
  psd=sqrt(ptot2/totalpoint)

  if(usd > 0.0) uscl = min(0.20*10.0/usd, 1.0)
  vscl = uscl
  if(ptsd > 0.0) ptscl = min(0.20*11.0/ptsd, 1.0)
  if(qvsd > 0.0) qvscl = min(0.20*2e-3/qvsd, 1.0)
  if(psd > 0.0) pscl = min(0.20*1000.0/psd, 1.0)

  print *,'totalpoint=',totalpoint
  print *,'usd=',usd,' uscl=',uscl
  print *,'vsd=',vsd,' vscl=',vscl
  print *,'ptsd=',ptsd,' ptscl=',ptscl
  print *,'qvsd=',qvsd,' qvscl=',qvscl
  print *,'psd=',psd,' pscl=',pscl

  uprt = uscl*uprt
  vprt = vscl*vprt
  wprt = 0.0
  ptprt = ptscl*ptprt
  qvprt = qvscl*qvprt
  pprt = pscl*pprt

  END IF

  ELSE IF (inibred ==1) THEN    ! generate random perturbations

  do k=1,nz-1
  do i=1,nx
  do j=1,ny-1
  iseed=mod(iseed*7141+54773,259200)
  uprt(i,j,k) = 2.*0.20*usd*(iseed/259199.-.5)
  enddo
  enddo
  enddo
  do k=1,nz-1
  do i=1,nx-1
  do j=1,ny
  iseed=mod(iseed*7141+54773,259200)
  vprt(i,j,k) = 2.*0.20*vsd*(iseed/259199.-.5)
  enddo
  enddo
  enddo
  do k=1,nz-1
  do i=1,nx-1
  do j=1,ny-1
  iseed=mod(iseed*7141+54773,259200)
  ptprt(i,j,k) = 2.*0.20*ptsd*(iseed/259199.-.5)
  enddo
  enddo
  enddo
  do k=1,nz-1
  do i=1,nx-1
  do j=1,ny-1
  iseed=mod(iseed*7141+54773,259200)
  qvprt(i,j,k) = 2.*0.20*qvsd*(iseed/259199.-.5)
  enddo
  enddo
  enddo
  do k=1,nz-1
  do i=1,nx-1
  do j=1,ny-1
  iseed=mod(iseed*7141+54773,259200)
  pprt(i,j,k) = 2.*0.20*psd*(iseed/259199.-.5)
  enddo
  enddo
  enddo

  wprt=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
  tsfc=0.0
  tsoil=0.0
  wetsfc=0.0
  wetdp=0.0
  wetcanp=0.0
  raing=0.0
  rainc=0.0

  END IF
!
!-----------------------------------------------------------------------
!
!  Get the name of the CONTROL data set, field c or c=b (if iread=0)
!  (note: It is needed only if control .ne. verification)
!
!-----------------------------------------------------------------------
!
  IF (iread /= 0) THEN
!
!-----------------------------------------------------------------------
!
!  Read all input data arrays
!
!-----------------------------------------------------------------------
!
    CALL setgbrd (0)
    CALL dtaread(vnx,vny,vnz,nstyps,                                    &
                 hinfmt,nch,grdbasfn(3)(1:lengbf(3)),lengbf(3),         &
                 filename(3)(1:lenfil(3)),lenfil(3),time,               &
                 vx,vy,vz,vzp, vuprt ,vvprt ,vwprt ,vptprt, vpprt ,     &
                 vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv,       &
                 vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar,   &
                 vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg,          &
                 vtsfc,vtsoil,vwetsfc,vwetdp,vwetcanp, vsnowdpth,       &
         raing,rainc,prcrate,                                           &
         radfrc,radsw,rnflx,                                            &
         usflx,vsflx,ptsflx,qvsflx,                                     &
                 ireturn, vtem1,vtem2,vtem3)
    IF (isread == 1) THEN
      CALL readsoil(nx,ny,nstyps,soilinfl,dx,dy,                        &
                 mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,   &
                    tsfcin,tsoilin,wsfcin,wdpin,wcanpin,snowdin,        &
                    vtsfc, vtsoil, vwetsfc,vwetdp,vwetcanp,vsnowdpth)
    END IF
  IF( ireturn /= 0 ) THEN   ! failed read
    PRINT *,'Problem reading control data, forced to STOP'
    STOP
  END IF                      ! successful read
!    IF (vnx /= nx.OR.vny /= ny.OR.vnz /= nz) THEN
!      PRINT *,'nx,ny,nz','=/=','vnx,vny,vnz'
!      PRINT *,nx,ny,nz,'=/=',vnx,vny,vnz
!      PRINT *, 'forced to stop'
!      STOP
!    END IF
    ivproj=mapproj
    vscale=sclfct
    vlatnot(1)=trulat1
    vlatnot(2)=trulat2
    vtrulon=trulon
    vdx=vx(3)-vx(2)
    vdy=vy(3)-vy(2)
    vctrlat=ctrlat
    vctrlon=ctrlon
    CALL setmapr(ivproj,vscale,vlatnot,vtrulon)
    CALL lltoxy(1,1,vctrlat,vctrlon,xctr,yctr)
    vx0=xctr-vdx*((vnx-3)/2)
    vy0=yctr-vdy*((vny-3)/2)
    CALL setorig(1,vx0,vy0)
!
    IF (fx0 == vx0.AND.fy0 == vy0.AND.                                  &
          flatnot(1) == vlatnot(1).AND.flatnot(2) == vlatnot(2).AND.    &
          ftrulon == vtrulon.AND.ifproj == ivproj .AND.                 &
          fscale == vscale ) THEN
      PRINT *, 'Grids 2 and 3 shares a common coordinate system'
    ELSE
      PRINT *, 'Grids 2/3 are different, CHECK the PROFRAM or data'
      PRINT *, 'Forced to STOP'
      STOP
    END IF
  END IF
!-----------------------------------------------------------------------
!
!  Set output variables to forecast coordinates
!
!-----------------------------------------------------------------------
!
  curtim=time
  mapproj=ivproj
  sclfct=vscale
  trulat1=vlatnot(1)
  trulat2=vlatnot(2)
  trulon=vtrulon
  ctrlat=vctrlat
  ctrlon=vctrlon
!
!-----------------------------------------------------------------------
!
!  Find   ptb field =  filed  c +  a/n (note iorder=/=0)
!       (b=b+a/n  in the code)
!  To reduce memory requirements, the resulted fields are
!  written to the same arrays as the control fields.
!
!-----------------------------------------------------------------------
!
  IF (iorder /= 0) THEN
    CALL prtfield(nx,ny,nz,iorder,                                      &
                  vuprt, vvprt, vwprt, vptprt, vpprt,                   &
                  vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv,      &
                  vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar,  &
                  vtsfc,vtsoil,vwetsfc,vwetdp,vwetcanp,                 &
                  vraing,vrainc,                                        &
                  uprt, vprt, wprt, ptprt, pprt,                        &
                  qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,               &
                  ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,         &
                  tsfc,tsoil,wetsfc,wetdp,wetcanp,                      &
                  raing,rainc,                                          &
                  vuprt, vvprt, vwprt, vptprt, vpprt,                   &
               vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv,vtem1,   &
                  vtsfc,vtsoil,vwetsfc,vwetdp,vwetcanp,                 &
                  vraing,vrainc,                                        &
                  ireturn )
  END IF
!-----------------------------------------------------------------------
!  Assign the average of soil variables to every soil type
!-----------------------------------------------------------------------
  CALL dhslini(nx,ny,nstyps,                                            &
                tsfc,tsoil,wetsfc,wetdp,wetcanp)
!
!-----------------------------------------------------------------------
!
!  Get output info
!
!-----------------------------------------------------------------------
!
  runname=runnmin
!
!-----------------------------------------------------------------------
!
!  Find out the number of characters to be used to construct file
!  names.
!
!-----------------------------------------------------------------------
!
  CALL gtlfnkey( runname, lfnkey )
!
!-----------------------------------------------------------------------
!
!  Find out the number of characters to be used to construct file
!  names.
!
!-----------------------------------------------------------------------
!
  CALL gtlfnkey( runname, lfnkey )
!
  IF (hdmpfmt == 10.AND.grbpkbit == 0) THEN
    grbpkbit=16
  END IF
!
!-----------------------------------------------------------------------
!
!  Set control parameters for
!  grid, base state, moisture, and ice variable dumping.
!
!-----------------------------------------------------------------------
!
  varout=1
  totout = totin
  grdout = grdin
  basout = basin
  varout = varin
  mstout = mstin
  rainout = rainin
  prcout = prcin
  iceout = icein
  tkeout = tkein
  trbout = trbin
  sfcout = sfcin
  snowout = snowin
  landout = landin
  radout = radin
  flxout = flxin

  CALL gtbasfn(runname(1:lfnkey),'./',2,hdmpfmt,mgrid,nestgrd,          &
               grdbasfn(1), lengbf(1))

  WRITE(6,'(/1x,a,a)')                                                  &
      'Output grid/base state file is ', grdbasfn(1)(1:lengbf(1))
  nchdmp = 80
  grdbas = 1      ! Dump out grd and base state arrays only

  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        vuprt(i,j,k)=vubar(i,j,k)+vuprt(i,j,k)
        vvprt(i,j,k)=vvbar(i,j,k)+vvprt(i,j,k)
        vwprt(i,j,k)=vwbar(i,j,k)+vwprt(i,j,k)
        vqvprt(i,j,k)=vqvbar(i,j,k)+vqvprt(i,j,k)
      END DO
    END DO
  END DO
  IF (iorder /= 0) THEN
    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          IF (k /= 1) THEN
            IF ((vpprt(i,j,k)+vpbar(i,j,k)) >=                          &
                  (vpprt(i,j,k-1)+vpbar(i,j,k-1))  ) THEN
              vpprt(i,j,k)=vpprt(i,j,k-1)+vpbar(i,j,k-1)-vpbar(i,j,k)   &
                  -(vpbar(i,j,k-1)-vpbar(i,j,k))*0.001
            END IF
          END IF
          IF (vqvprt(i,j,k) <= 1.0E-8) THEN
            vqvprt(i,j,k)=1.0E-8
          END IF
        END DO
      END DO
    END DO
  END IF

  CALL dtadump(nx,ny,nz,nstyps,                                         &
       hdmpfmt,nchdmp,grdbasfn (1)(1:lengbf(1)),grdbas,filcmprs,        &
               vuprt,vvprt,vwprt,vptprt,vpprt,                          &
               vqvprt,vqc,vqr,vqi,vqs,vqh,vtke,vkmh,vkmv,               &
               vubar,vvbar,vwbar,vptbar,vpbar,vrhobar,vqvbar,           &
               vx,vy,vz,vzp,hterain,j1,j2,j3,                           &
               vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg,            &
               vtsfc,vtsoil,vwetsfc,vwetdp,vwetcanp,vsnowdpth,          &
             vraing,vrainc,prcrate,                                     &
             radfrc,radsw,rnflx,                                        &
             usflx,vsflx,ptsflx,qvsflx,                                 &
               vtem1,vtem2,vtem3)
!
!-----------------------------------------------------------------------
!
!  Find a unique name hdmpfn(1:ldmpf) for history dump data set
!  at time 'curtim'.
!
!-----------------------------------------------------------------------
!
  CALL gtdmpfn(runname(1:lfnkey),'./',2,                                &
               curtim,hdmpfmt,                                          &
               mgrid,nestgrd, hdmpfn, ldmpf)

  WRITE(6,'(/1x,a,f10.0,a,a)')                                          &
      'Output file at time ',curtim,' (s) is ', hdmpfn(1:ldmpf)

  grdbas = 0      ! Not just dump out grd and base state arrays only

  CALL dtadump(nx,ny,nz,nstyps,                                         &
               hdmpfmt,nchdmp,hdmpfn(1:ldmpf),grdbas,filcmprs,          &
               vuprt,vvprt,vwprt,vptprt,vpprt,                          &
               vqvprt,vqc,vqr,vqi,vqs,vqh,vtke,vkmh,vkmv,               &
               vubar,vvbar,vwbar,vptbar,vpbar,vrhobar,vqvbar,           &
               vx,vy,vz,vzp,hterain, j1,j2,j3,                          &
               vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg,            &
               vtsfc,vtsoil,vwetsfc,vwetdp,vwetcanp,vsnowdpth,          &
             vraing,vrainc,prcrate,                                     &
             radfrc,radsw,rnflx,                                        &
             usflx,vsflx,ptsflx,qvsflx,                                 &
               vtem1,vtem2,vtem3)

  STOP

END PROGRAM arpsensic
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE PRTFIELD                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE prtfield(nx,ny,nz,nscale,                                    & 2,21
           uprt, vprt, wprt, ptprt, pprt,                               &
           qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                      &
           ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,                &
           tsfc,tsoil,wetsfc,wetdp,wetcanp,                             &
           raing,rainc,                                                 &
           auprt, avprt, awprt, aptprt, apprt,                          &
           aqvprt, aqc, aqr, aqi, aqs, aqh, atke,akmh,akmv,             &
           aubar, avbar, awbar, aptbar, apbar, arhobar, aqvbar,         &
           atsfc,atsoil,awetsfc,awetdp,awetcanp,                        &
           araing,arainc,                                               &
           duprt, dvprt, dwprt, dptprt, dpprt,                          &
           dqvprt, dqc, dqr, dqi, dqs, dqh, dtke,dkmh,dkmv,tem1,        &
           dtsfc,dtsoil,dwetsfc,dwetdp,dwetcanp,                        &
           draing,drainc,                                               &
           ireturn)
!
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  (1) nscale=0, Subtract the forecast fields from the verification
!  fields (names beginning with "a") and output to the difference
!  fields (names beginning with "d"). d=( )bar+( )-abar-a
!  (2) nscale=/=0. Generate  perturbation fields d=( )+a/nscale with
!  the bar arraies being neglected.
!     The input (  ( )   ) and output  ( d )
!  fields may share the same storage location.  For this subroutine
!  it is assumed the forecast and corresponding verification
!  data are at the same physical location, however, the physical
!  location may differ between variables. That is uprt and auprt
!  are at the same location, but that may differ from pprt and apprt.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster  Ou School of Meteorology. April 1992
!
!  MODIFICATION HISTORY:
!   14 May 1992  (KB) changed from arps2.5 to arps3.0
!   03 Aug 1992  (KB) updated to account for changes in arps3.0
!
!   09/07/1995  (KB)
!   Added differencing of surface (soil) fields.
!   15/05/1998  (DH)
!   converted from the difference scheme to the current multi-purpose
!   version used in ARPSENS group.
!   15/12/1998  (DH)
!   Added the 3-Dimensionity of surface (soil) fields.
!
!-----------------------------------------------------------------------
!
!  INPUT :
!    nx,ny,nz Array dimensions for forecast field.
!
!    FORECAST FIELDS:
!
!    uprt     perturbation x component of velocity (m/s)
!    vprt     perturbation y component of velocity (m/s)
!    wprt     perturbation vertical component of velocity in Cartesian
!             coordinates (m/s).
!
!    ptprt    perturbation potential temperature (K)
!    pprt     perturbation pressure (Pascal)
!
!    qvprt    perturbation water vapor mixing ratio (kg/kg)
!    qc       Cloud water mixing ratio (kg/kg)
!    qr       Rainwater mixing ratio (kg/kg)
!    qi       Cloud ice mixing ratio (kg/kg)
!    qs       Snow mixing ratio (kg/kg)
!    qh       Hail mixing ratio (kg/kg)
!
!    ubar     Base state x velocity component (m/s)
!    vbar     Base state y velocity component (m/s)
!    wbar     Base state z velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhobar   Base state density (kg/m**3)
!    qvbar    Base state water vapor mixing ratio (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!    kmh      Horizontal turbulent mixing coefficient (m**2/s)
!    kmv      Vertical turbulent mixing coefficient (m**2/s)
!
!    tsfc     Temperature at surface (K)
!    tsoil    Deep soil temperature (K)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!
!    INTERPOLATED VERIFICATION FIELDS:
!
!    auprt     perturbation x component of velocity (m/s)
!    avprt     perturbation y component of velocity (m/s)
!    awprt     perturbation vertical component of velocity in Cartesian
!              coordinates (m/s).
!
!    aptprt    perturbation potential temperature (K)
!    apprt     perturbation pressure (Pascal)
!
!    aqvprt    perturbation water vapor mixing ratio (kg/kg)
!    aqc       Cloud water mixing ratio (kg/kg)
!    aqr       Rainwater mixing ratio (kg/kg)
!    aqi       Cloud ice mixing ratio (kg/kg)
!    aqs       Snow mixing ratio (kg/kg)
!    aqh       Hail mixing ratio (kg/kg)
!
!    aubar     Base state x velocity component (m/s)
!    avbar     Base state y velocity component (m/s)
!    awbar     Base state z velocity component (m/s)
!    aptbar    Base state potential temperature (K)
!    apbar     Base state pressure (Pascal)
!    arhobar   Base state density (kg/m**3)
!    aqvbar    Base state water vapor mixing ratio (kg/kg)
!
!    atsfc     Temperature at surface (K)
!    atsoil    Deep soil temperature (K)
!    awetsfc   Surface soil moisture
!    awetdp    Deep soil moisture
!    awetcanp  Canopy water amount
!
!    araing    Grid supersaturation rain
!    arainc    Cumulus convective rain
!
!  OUTPUT :
!
!    DIFFERENCE FIELDS (may share storage with forecast fields
!                       or interpolated fields in calling program):
!
!    duprt     perturbation x component of velocity (m/s)
!    dvprt     perturbation y component of velocity (m/s)
!    dwprt     perturbation vertical component of velocity in Cartesian
!              coordinates (m/s).
!
!    dptprt    perturbation potential temperature (K)
!    dpprt     perturbation pressure (Pascal)
!
!    dqvprt    perturbation water vapor mixing ratio (kg/kg)
!    dqc       Cloud water mixing ratio (kg/kg)
!    dqr       Rainwater mixing ratio (kg/kg)
!    dqi       Cloud ice mixing ratio (kg/kg)
!    dqs       Snow mixing ratio (kg/kg)
!    dqh       Hail mixing ratio (kg/kg)
!
!    dtsfc     Temperature at surface (K)
!    dtsoil    Deep soil temperature (K)
!    dwetsfc   Surface soil moisture
!    dwetdp    Deep soil moisture
!    dwetcanp  Canopy water amount
!
!    draing    Grid supersaturation rain
!    drainc    Cumulus convective rain
!
!    tem1      Work array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz         ! 3 dimensions of array
  INTEGER :: nscale
!
!-----------------------------------------------------------------------
!
!  Model Arrays
!
!-----------------------------------------------------------------------
!
  REAL :: uprt   (nx,ny,nz)    ! Perturbation u-velocity (m/s)
  REAL :: vprt   (nx,ny,nz)    ! Perturbation v-velocity (m/s)
  REAL :: wprt   (nx,ny,nz)    ! Perturbation w-velocity (m/s)
  REAL :: ptprt  (nx,ny,nz)    ! Perturbation potential temperature (K)
  REAL :: pprt   (nx,ny,nz)    ! Perturbation pressure (Pascal)
  REAL :: qvprt  (nx,ny,nz)    ! Perturbation water vapor specific humidity
  REAL :: qc     (nx,ny,nz)    ! Cloud water mixing ratio (kg/kg)
  REAL :: qr     (nx,ny,nz)    ! Rain water mixing ratio (kg/kg)
  REAL :: qi     (nx,ny,nz)    ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs     (nx,ny,nz)    ! Snow mixing ratio (kg/kg)
  REAL :: qh     (nx,ny,nz)    ! Hail mixing ratio (kg/kg)

  REAL :: tke    (nx,ny,nz)    ! Turbulent Kinetic Energy ((m/s)**2)
  REAL :: kmh    (nx,ny,nz)    ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: kmv    (nx,ny,nz)    ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )

  REAL :: ubar   (nx,ny,nz)    ! Base state u-velocity (m/s)
  REAL :: vbar   (nx,ny,nz)    ! Base state v-velocity (m/s)
  REAL :: wbar   (nx,ny,nz)    ! Base state w-velocity (m/s)
  REAL :: ptbar  (nx,ny,nz)    ! Base state potential temperature (K)
  REAL :: pbar   (nx,ny,nz)    ! Base state pressure (Pascal)
  REAL :: rhobar (nx,ny,nz)    ! Base state air density (kg/m**3)
  REAL :: qvbar  (nx,ny,nz)    ! Base state water vapor specific humidity

  REAL :: tsfc  (nx,ny)       ! Temperature at surface (K)
  REAL :: tsoil (nx,ny)       ! Deep soil temperature (K)
  REAL :: wetsfc(nx,ny)       ! Surface soil moisture
  REAL :: wetdp (nx,ny)       ! Deep soil moisture
  REAL :: wetcanp(nx,ny)      ! Canopy water amount

  REAL :: raing (nx,ny)       ! Grid supersaturation rain
  REAL :: rainc (nx,ny)       ! Cumulus convective rain
!
!-----------------------------------------------------------------------
!
!  Verification data interpolated to model grid
!
!-----------------------------------------------------------------------
!
  REAL :: auprt  (nx,ny,nz)    ! Perturbation u-velocity (m/s)
  REAL :: avprt  (nx,ny,nz)    ! Perturbation v-velocity (m/s)
  REAL :: awprt  (nx,ny,nz)    ! Perturbation w-velocity (m/s)
  REAL :: aptprt (nx,ny,nz)    ! Perturbation potential temperature (K)
  REAL :: apprt  (nx,ny,nz)    ! Perturbation pressure (Pascal)

  REAL :: aqvprt (nx,ny,nz)    ! Perturbation water vapor specific humidity
  REAL :: aqc    (nx,ny,nz)    ! Cloud water mixing ratio (kg/kg)
  REAL :: aqr    (nx,ny,nz)    ! Rain water mixing ratio (kg/kg)
  REAL :: aqi    (nx,ny,nz)    ! Cloud ice mixing ratio (kg/kg)
  REAL :: aqs    (nx,ny,nz)    ! Snow mixing ratio (kg/kg)
  REAL :: aqh    (nx,ny,nz)    ! Hail mixing ratio (kg/kg)

  REAL :: atke   (nx,ny,nz)    ! Turbulent Kinetic Energy ((m/s)**2)
  REAL :: akmh   (nx,ny,nz)    ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: akmv   (nx,ny,nz)    ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )

  REAL :: aubar  (nx,ny,nz)    ! Base state u-velocity (m/s)
  REAL :: avbar  (nx,ny,nz)    ! Base state v-velocity (m/s)
  REAL :: awbar  (nx,ny,nz)    ! Base state w-velocity (m/s)
  REAL :: aptbar (nx,ny,nz)    ! Base state potential temperature (K)
  REAL :: arhobar(nx,ny,nz)    ! Base state density (kg/m**3)
  REAL :: apbar  (nx,ny,nz)    ! Base state pressure (Pascal)
  REAL :: aqvbar (nx,ny,nz)    ! Base state water vapor specific humidity

  REAL :: atsfc  (nx,ny)       ! Temperature at surface (K)
  REAL :: atsoil (nx,ny)       ! Deep soil temperature (K)
  REAL :: awetsfc(nx,ny)       ! Surface soil moisture
  REAL :: awetdp (nx,ny)       ! Deep soil moisture
  REAL :: awetcanp(nx,ny)      ! Canopy water amount

  REAL :: araing (nx,ny)       ! Grid supersaturation rain
  REAL :: arainc (nx,ny)       ! Cumulus convective rain
!
!------------------------------------------------------------------------
!
!  Difference arrays
!
!-----------------------------------------------------------------------
!
  REAL :: duprt  (nx,ny,nz)    ! perturbation x component of velocity (m/s)
  REAL :: dvprt  (nx,ny,nz)    ! perturbation y component of velocity (m/s)
  REAL :: dwprt  (nx,ny,nz)    ! perturbation vertical component of
                               ! velocity in Cartesian coordinates (m/s)
  REAL :: dptprt (nx,ny,nz)    ! perturbation potential temperature (K)
  REAL :: dpprt  (nx,ny,nz)    ! perturbation pressure (Pascal)
  REAL :: dqvprt (nx,ny,nz)    ! perturbation water vapor mixing ratio (kg/kg)
  REAL :: dqc    (nx,ny,nz)    ! Cloud water mixing ratio (kg/kg)
  REAL :: dqr    (nx,ny,nz)    ! Rainwater mixing ratio (kg/kg)
  REAL :: dqi    (nx,ny,nz)    ! Cloud ice mixing ratio (kg/kg)
  REAL :: dqs    (nx,ny,nz)    ! Snow mixing ratio (kg/kg)
  REAL :: dqh    (nx,ny,nz)    ! Hail mixing ratio (kg/kg)
  REAL :: dtke   (nx,ny,nz)    ! Turbulent Kinetic Energy ((m/s)**2)
  REAL :: dkmh   (nx,ny,nz)    ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: dkmv   (nx,ny,nz)    ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )

  REAL :: dtsfc  (nx,ny)       ! Temperature at surface (K)
  REAL :: dtsoil (nx,ny)       ! Deep soil temperature (K)
  REAL :: dwetsfc(nx,ny)       ! Surface soil moisture
  REAL :: dwetdp (nx,ny)       ! Deep soil moisture
  REAL :: dwetcanp(nx,ny)      ! Canopy water amount

  REAL :: draing (nx,ny)       ! Grid supersaturation rain
  REAL :: drainc (nx,ny)       ! Cumulus convective rain

  REAL :: tem1   (nx,ny,nz)    ! A work array
  INTEGER :: ireturn, i,j,k
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: is,js,ks,ie,je,ke
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  is=1
  js=1
  ks=1
  ie=nx-1
  je=ny-1
  ke=nz-1

!-----------------------------------------------------------------------
!
!  Scalars
!
!-----------------------------------------------------------------------

  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        tem1(i,j,k)=0.0
      END DO
    END DO
  END DO

  PRINT *, ' ptprt: '
  CALL subtrprt(nx,ny,nz, ptprt,ptbar,aptprt,aptbar,dptprt,nscale,      &
             is,js,ks,ie,je,ke)
  PRINT *, ' pprt: '
  CALL subtrprt(nx,ny,nz,  pprt, pbar, apprt, apbar, dpprt,nscale,      &
             is,js,ks,ie,je,ke)
  PRINT *, ' qvprt: '
  CALL subtrprt(nx,ny,nz, qvprt,qvbar,aqvprt,aqvbar,dqvprt,nscale,      &
             is,js,ks,ie,je,ke)
  PRINT *, ' qc: '
  CALL subtrprt(nx,ny,nz,    qc, tem1,   aqc,  tem1,   dqc,nscale,      &
             is,js,ks,ie,je,ke)
  PRINT *, ' qr: '
  CALL subtrprt(nx,ny,nz,    qr, tem1,   aqr,  tem1,   dqr,nscale,      &
             is,js,ks,ie,je,ke)
  PRINT *, ' qi: '
  CALL subtrprt(nx,ny,nz,    qi, tem1,   aqi,  tem1,   dqi,nscale,      &
             is,js,ks,ie,je,ke)
  PRINT *, ' qs: '
  CALL subtrprt(nx,ny,nz,    qs, tem1,   aqs,  tem1,   dqs,nscale,      &
             is,js,ks,ie,je,ke)
  PRINT *, ' qh: '
  CALL subtrprt(nx,ny,nz,    qh, tem1,   aqh,  tem1,   dqh,nscale,      &
             is,js,ks,ie,je,ke)
  PRINT *, ' tke: '
  CALL subtrprt(nx,ny,nz,    tke, tem1,   atke,  tem1, dtke,nscale,     &
             is,js,ks,ie,je,ke)
  PRINT *, ' kmh: '
  CALL subtrprt(nx,ny,nz,    kmh, tem1,   akmh,  tem1, dkmh,nscale,     &
             is,js,ks,ie,je,ke)
  PRINT *, ' kmv: '
  CALL subtrprt(nx,ny,nz,    kmv, tem1,   akmv,  tem1, dkmv,nscale,     &
             is,js,ks,ie,je,ke)

!-----------------------------------------------------------------------
!
!  u wind components
!
!-----------------------------------------------------------------------

  ie=nx
  PRINT *, ' uprt: '
  CALL subtrprt(nx,ny,nz,uprt,ubar,auprt,aubar,duprt,nscale,            &
             is,js,ks,ie,je,ke)

!-----------------------------------------------------------------------
!
!  v wind components
!
!-----------------------------------------------------------------------

  ie=nx-1
  je=ny
  PRINT *, ' vprt: '
  CALL subtrprt(nx,ny,nz,vprt,vbar,avprt,avbar,dvprt,nscale,            &
             is,js,ks,ie,je,ke)

!-----------------------------------------------------------------------
!
!  w wind components
!
!-----------------------------------------------------------------------

  je=ny-1
  ke=nz
  PRINT *, ' wprt: '
  CALL subtrprt(nx,ny,nz,wprt,tem1,awprt,tem1,dwprt,nscale,             &
             is,js,ks,ie,je,ke)

!-----------------------------------------------------------------------
!
!  2-d surface (soil) variables
!
!-----------------------------------------------------------------------

  ie=nx-1
  je=ny-1
  ks=1
  ke=1
!
  PRINT *, ' tsfc:'
  CALL subtrprt(nx,ny,1,   tsfc,tem1,   atsfc,tem1,   dtsfc,nscale,     &
             is,js,ks,ie,je,ke)
  PRINT *, ' tsoil:'
  CALL subtrprt(nx,ny,1,  tsoil,tem1,  atsoil,tem1,  dtsoil,nscale,     &
             is,js,ks,ie,je,ke)
  PRINT *, ' wetsfc:'
  CALL subtrprt(nx,ny,1, wetsfc,tem1, awetsfc,tem1, dwetsfc,nscale,     &
             is,js,ks,ie,je,ke)
  PRINT *, ' wetdp:'
  CALL subtrprt(nx,ny,1,  wetdp,tem1,  awetdp,tem1,  dwetdp,nscale,     &
             is,js,ks,ie,je,ke)
  PRINT *, ' wetcanp:'
  CALL subtrprt(nx,ny,1,wetcanp,tem1,awetcanp,tem1,dwetcanp,nscale,     &
             is,js,ks,ie,je,ke)
  PRINT *, ' raing:'
  CALL subtrprt(nx,ny,1,  raing,tem1,  araing,tem1,  draing,nscale,     &
             is,js,ks,ie,je,ke)
  PRINT *, ' rainc:'
  CALL subtrprt(nx,ny,1,  rainc,tem1,  arainc,tem1,  drainc,nscale,     &
             is,js,ks,ie,je,ke)
!
  RETURN
END SUBROUTINE prtfield
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE SUBTRPRT                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE subtrprt(nx,ny,nz, a,abar,b,bbar,c,nscale,                   & 21
           istr,jstr,kstr,iend,jend,kend)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Subtracts 2 three-dimensional arrays, represented by
!  means plus perturbations.
!
!  AUTHOR: Keith Brewster  OU School of Meteorology.  Feb 1992
!
!  MODIFICATION HISTORY:
!   11 Aug 1992  (KB) changed from arps2.5 to arps3.0
!
!
!-----------------------------------------------------------------------
!
!  INPUT:
!    a        perturbation data array
!    abar     mean data array
!    b        perturbation data array to subtract from a
!    bbar     mean data array to subtract from a
!
!  OUTPUT:
!    c        difference array a-b
!             (may share storage in calling program with array a or b)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! 3 dimensions of array
  INTEGER :: nscale

  REAL :: a   (nx,ny,nz)       ! data array
  REAL :: abar(nx,ny,nz)       ! base state of data array a
  REAL :: b   (nx,ny,nz)       ! data array to subtract from a
  REAL :: bbar(nx,ny,nz)       ! base state of data arrya b
  REAL :: c   (nx,ny,nz)       ! difference array a-b
  INTEGER :: istr,jstr,kstr
  INTEGER :: iend,jend,kend
  INTEGER :: i,j,k,imid,jmid,kmid
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  imid=nint(0.5*FLOAT(istr+iend))
  jmid=nint(0.5*FLOAT(jstr+jend))
  kmid=nint(0.5*FLOAT(kstr+kend))
  IF (nz < 10) kmid=kstr
  PRINT *, 'imid,jmid,kmid: ',imid,jmid,kmid
!
!-----------------------------------------------------------------------
!
!  Tell us about a sample input point
!
!-----------------------------------------------------------------------
!
  IF (nscale == 0) THEN
    PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)),   &
                     ' b= ',(b(imid,jmid,kmid)+bbar(imid,jmid,kmid))
  ELSE
    PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)),   &
                     ' b= ',b(imid,jmid,kmid)
  END IF
!
!-----------------------------------------------------------------------
!
!  Subtraction
!
!-----------------------------------------------------------------------
!
  DO k=kstr,kend
    DO j=jstr,jend
      DO i=istr,iend
        IF (nscale == 0) THEN
          c(i,j,k)=a(i,j,k)+abar(i,j,k)-(b(i,j,k)+bbar(i,j,k))
        ELSE
          c(i,j,k)=a(i,j,k)+b(i,j,k)/FLOAT(nscale)
        END IF
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Tell us about a sample output point
!
!-----------------------------------------------------------------------
!
  IF (nscale == 0) THEN
    PRINT *, '         c= ',c(imid,jmid,kmid)
  ELSE
    PRINT *, '         c= ',c(imid,jmid,kmid)+abar(imid,jmid,kmid)
  END IF
  RETURN
END SUBROUTINE subtrprt


SUBROUTINE dhslini(nx,ny,nstyps,                                        & 1
           tsfc,tsoil,wetsfc,wetdp,wetcanp)
  IMPLICIT NONE
  INTEGER :: nx,ny,nstyps
  INTEGER :: i,j,k

  REAL :: tsfc   (nx,ny,0:nstyps)      ! Temperature at surface (K)
  REAL :: tsoil  (nx,ny,0:nstyps)      ! Deep soil temperature (K)
  REAL :: wetsfc (nx,ny,0:nstyps)      ! Surface soil moisture
  REAL :: wetdp  (nx,ny,0:nstyps)      ! Deep soil moisture
  REAL :: wetcanp(nx,ny,0:nstyps)      ! Canopy water amount
  DO k=1,nstyps
    DO j=1,ny
      DO i=1,nx
        tsfc(i,j,k)=tsfc(i,j,0)
        tsoil(i,j,k)=tsoil(i,j,0)
        wetsfc(i,j,k)=wetsfc(i,j,0)
        wetdp(i,j,k)=wetdp(i,j,0)
        wetcanp(i,j,k)=wetcanp(i,j,0)
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE dhslini