!wdt update acct variables & f_walltime to double precision
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE RANARY                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE ranary(nx,ny,nx1,nx2,ny1,ny2,iseed,amplit, rantp) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Generate a 2-D array of machine-independent random numbers
!  between -amplit and +amplit with average value equal to zero.
!  The input parameter, iseed, must be a negative integer number
!  which needs to be defined only once in the calling subroutine
!  program.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: V.Wong and X.Song
!  7/20/1992.
!
!  Modification history:
!  7/25/1992. (MX)
!  Modified to allow assignment on a portion of the array.
!
!  12/11/1992 (MX)
!  Bug fix to the index bounds in loop 10.
!
!  9/1/94 (J. Levit & Y. Lu)
!  Cleaned up documentation
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx         Number of random numbers in the x-direction
!    ny         Number of random numbers in the y-direction
!    nx1,nx2    Begin and end 1st array indicies specifying a
!               subdomain in which the average is set to zero
!    ny1,ny2    Begin and end 2nd array indicies specifying a
!               subdomain in which the average is set to zero
!    iseed      an arbitrary negative integer as a seed for a
!               sequence of random numbers
!    amplit     The generated numbers stay within the range
!               [-amplit,amplit]
!
!  OUTPUT:
!
!    rantp      The array for storing random numbers
!
!  WORK ARRAY:
!
!    rantp      Temporary working array
!
!
!-----------------------------------------------------------------------
!
!
!
!-----------------------------------------------------------------------
!
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE        ! Force explicit declarations

  INTEGER :: i,j
  INTEGER :: nx,ny        ! Number of random numbers in each direction
  INTEGER :: nx1,nx2      ! Begin and end 1st array indecies specifying
                          ! a subdomain in which the average is set to
                          ! zero
  INTEGER :: ny1,ny2      ! Begin and end 2nd array indecies specifying
                          ! a subdomain in which the average is set to
                          ! zero
  INTEGER :: iseed        ! The seed for random number generation
  REAL :: amplit          ! The generated numbers stay within
                          ! [-amplit, amplit]
  REAL :: rantp (nx,ny)   ! Temporary working array for storing the
                          ! numbers
  REAL :: ran3            ! The function to generate random number.
!
!-----------------------------------------------------------------------
!
!  Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!
  REAL :: wv,ranave,rantpmax
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!-----------------------------------------------------------------------
!
!  Generate random numbers between 0 and 1, for a given iseed.
!
!-----------------------------------------------------------------------
!
  ranave=0.0
  DO j=ny1,ny2
    DO i=nx1,nx2
      rantp(i,j)=ran3(iseed)
      ranave=ranave+rantp(i,j)
    END DO
  END DO
  ranave=ranave/(FLOAT( (nx2-nx1+1)*(ny2-ny1+1) ))
!
!-----------------------------------------------------------------------
!
!  Adjust the random numbers so that the new plane average value
!  equals zero and the random numbers stay within the range
!  [-amplit, amplit].
!
!-----------------------------------------------------------------------
!
  rantpmax=0.0
  DO j=ny1,ny2
    DO i=nx1,nx2
      rantpmax=AMAX1(rantpmax,ABS(ranave-rantp(i,j)))
    END DO
  END DO

  wv=amplit/rantpmax

  DO j=ny1,ny2
    DO i=nx1,nx2
      rantp(i,j)=(rantp(i,j)-ranave)*wv
    END DO
  END DO

  RETURN
END SUBROUTINE ranary

!
!##################################################################
!##################################################################
!######                                                      ######
!######                  FUNCTION RAN3                       ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


  FUNCTION ran3(iseed)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Generates a random number between 0 and 1 by feeding
!  a negative integer iseed.
!
!  Reference: "Seminumerical Algorithms" by Donald Knuth
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    iseed      an arbitrary negative integer as a seed for a
!               sequence of random numbers
!
!  OUTPUT:
!
!    ran3       A random number between 0 and 1.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE        ! Force explicit declarations

  INTEGER :: iseed        ! The seed for random number generation
  REAL :: ran3            ! The function to generate random number.
!
!-----------------------------------------------------------------------
!
!  Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: mbig,mseed,mz,k,ii,inext,inextp,i,iff,mk,mj
  REAL :: fac

  PARAMETER (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
  INTEGER :: ma(55)
  SAVE iff, inext, inextp, ma
  DATA iff /0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!----------------------------------------------------------------------
!
!  Initialize the sequence of random numbers between 0 and 1,
!  using iseed.
!
!----------------------------------------------------------------------
!
  IF (iseed < 0.OR.iff == 0) THEN
    iff=1
    mj=mseed-IABS(iseed)
    mj=MOD(mj,mbig)
    ma(55)=mj
    mk=1
    DO i=1,54
      ii=MOD(21*i,55)
      ma(ii)=mk
      mk=mj-mk
      IF (mk < mz) mk=mk+mbig
      mj=ma(ii)
    END DO

    DO k=1,4
      DO i=1,55
        ma(i)=ma(i)-ma(1+MOD(i+30,55))
        IF (ma(i) < mz) ma(i)=ma(i)+mbig
      END DO
    END DO

    inext=0
    inextp=31
    iseed=1
  END IF
!
!----------------------------------------------------------------------
!
!  Start to generate a random number.
!
!----------------------------------------------------------------------
!
  inext=inext+1
  IF (inext == 56) inext=1
  inextp=inextp+1
  IF (inextp == 56) inextp=1
  mj=ma(inext)-ma(inextp)
  IF (mj < mz) mj=mj+mbig
  ma(inext)=mj
  ran3=mj*fac

  RETURN
  END FUNCTION ran3
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE SWAP4BYTE                  ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE swap4byte(a,n) 2
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Reverse order of bytes in integer*4 words
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: n
  INTEGER*4   a(n)

  INTEGER :: i

  INTEGER*4   itema
  CHARACTER (LEN=1) :: ctema(4)
  CHARACTER (LEN=1) :: ctemb

  EQUIVALENCE ( itema,ctema(1) )
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO i = 1,n
    itema    = a(i)
    ctemb    = ctema(4)
    ctema(4) = ctema(1)
    ctema(1) = ctemb
    ctemb    = ctema(3)
    ctema(3) = ctema(2)
    ctema(2) = ctemb
    a(i)     = itema
  END DO

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


SUBROUTINE findvartype(iendn,itypec,lw) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Find the type of endian, type of character set and length of
!  machine word.  Renamed version of subroutine w3fi04 from nmcdecode.f90.
!
!-----------------------------------------------------------------------
!
!SUBROUTINE w3fi04(iendn,itypec,lw)
!$$$   SUBPROGRAM  DOCUMENTATION  BLOCK
!
! SUBPROGRAM: W3FI04         FIND WORD SIZE, ENDIAN, CHARACTER SET
!   PRGMNR: JONES,R.E.       ORG: W/NMC42       DATE: 94-10-07
!
! ABSTRACT: SUBROUTINE COMPUTES WORD SIZE, THE TYPE OF CHARACTER
!   SET, ASCII OR EBCDIC, AND IF THE COMPUTER IS BIG-ENDIAN, OR
!   LITTLE-ENDIAN.
!
! PROGRAM HISTORY LOG:
!   94-10-07  R.E.JONES
!
! USAGE:  CALL W3FI04 (IENDN, ITYPEC, LW)
!
!   OUTPUT ARGUMENT LIST:
!  IENDN     -  INTEGER FOR BIG-ENDIAN OR LITTLE-ENDIAN
!               = 0   BIG-ENDIAN
!               = 1   LITTLE-ENDIAN
!               = 2   CANNOT COMPUTE
!  ITYPEC    -  INTEGER FOR TYPE OF CHARACTER SET
!               = 0   ASCII  CHARACTER SET
!               = 1   EBCDIC CHARACTER SET
!               = 2   NOT ASCII OR EBCDIC
!  LW        -  INTEGER FOR WORDS SIZE OF COMPUTER IN BYTES
!               = 4   FOR 32 BIT COMPUTERS
!               = 8   FOR 64 BIT COMPUTERS
!
! ATTRIBUTES:
!   LANGUAGE: SiliconGraphics 3.5 FORTRAN 77
!   MACHINE:  SiliconGraphics IRIS-4D/25, 35, INDIGO, INDY
!
!$$$
!
  INTEGER :: itest1
  INTEGER :: itest2
  INTEGER :: itest3
  INTEGER :: iendn
  INTEGER :: itypec
  INTEGER :: lw
!
  CHARACTER (LEN=8) :: ctest1
  CHARACTER (LEN=8) :: ctest2
  CHARACTER (LEN=1) :: ctest3(8)
  CHARACTER (LEN=1) :: BLANK
!
  EQUIVALENCE   (ctest1,itest1),(ctest2,itest2)
!
  EQUIVALENCE   (itest3,ctest3(1))
!
  DATA  ctest1/'12345678'/
  DATA  itest3/z'01020304'/
  DATA  BLANK /' '/
!
  SAVE
!
!  TEST FOR TYPE OF CHARACTER SET
!  BLANK IS 32 (20 HEX) IN ASCII, 64 (40 HEX) IN EBCDEC
!
  IF (ICHAR(BLANK) == 32) THEN
    itypec = 0
  ELSE IF (ICHAR(BLANK) == 64) THEN
!
!  COMPUTER IS PROBABLY AN IBM360, 370, OR 390 WITH
!  A 32 BIT WORD SIZE, AND BIG-ENDIAN.
!
    itypec = 1
  ELSE
    itypec = 2
  END IF
!
!  TEST FOR WORD SIZE, SET LW TO 4 FOR 32 BIT COMPUTER,
!  8 FOR FOR 64 BIT COMPUTERS
!
  itest2 = itest1
  IF (ctest1 == ctest2) THEN
!
!  COMPUTER MAY BE A CRAY, OR COULD BE DEC VAX ALPHA
!  OR SGI WITH R4000, R4400, R8800 AFTER THEY CHANGE
!  FORTRAN COMPILERS FOR 64 BIT INTEGER.
!
    lw = 8
  ELSE
    lw = 4
  END IF
!
!  IF CHARACTER SET IS NOT ASCII OR EBCDIC SET IENDN = 2
!  CAN NOT TEST FOR ENDIAN TYPE
!
  IF (itypec == 2) THEN
    iendn = 2
    RETURN
  END IF
!
!  USING ITEST3 WITH Z'01020304' EQUIVALNCED TO CTEST3
!  ON A 32 BIT BIG-ENDIAN COMPUTER 03 IS IN THE 3RD
!  BYTE OF A 4 BYTE WORD. ON A 32 BIT LITTLE-ENDIAN
!  COMPUTER IT IS IN 2ND BYTE.
!  ON A 64 BIT COMPUTER Z'01020304' IS RIGHT ADJUSTED IN
!  A 64 BIT WORD, 03 IS IN THE 7TH BYTE.  ON A LITTLE-
!  ENDIAN 64 BIT COMPUTER IT IS IN THE 2ND BYTE.
!
  IF (lw == 4) THEN
    IF (ICHAR(ctest3(3)) == 3) THEN
      iendn = 0
    ELSE IF (ICHAR(ctest3(3)) == 2) THEN
      iendn = 1
    ELSE
      iendn = 2
    END IF
  ELSE IF (lw == 8) THEN
    IF (ICHAR(ctest3(7)) == 3) THEN
      iendn = 0
    ELSE IF (ICHAR(ctest3(2)) == 3) THEN
      iendn = 1
    ELSE
      iendn = 2
    END IF
  ELSE
    iendn = 2
  END IF
!
  RETURN
END SUBROUTINE findvartype
!
!##################################################################
!##################################################################
!######                                                      ######
!######               SUBROUTINE FIX_LAKE_ELIV               ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE fix_lake_eliv(h,lat,lon,nx,ny) 2
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Fill in Great Lakes.  The original data set corresponds to the
!  bottom of the lakes.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Richard Carpenter.
!  2000/02/03
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT/OUPUT :
!
!    h        Height of the grid point.
!
!  INPUT :
!
!    lat      Latitude of the grid point.
!    lon      Longitude of the grid point.
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!
!
!  OUTPUT:
!
!  WORK ARRAYS:
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx, ny            ! x- y-Dimensions.

  REAL :: h(nx,ny)             ! Height
  REAL :: lat(nx,ny)           ! Latitude
  REAL :: lon(nx,ny)           ! Longitude

!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
  INTEGER :: nlake, maxlake
  PARAMETER (maxlake=12)
  REAL :: lakelat1(maxlake), lakelat2(maxlake),                         &
      lakelon1(maxlake), lakelon2(maxlake), lakeelev(maxlake)

  REAL :: latpnt,lonpnt
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  PRINT *, 'ARPSTRN: Great Lakes terrain'
!
!-----------------------------------------------------------------------
!
!  Define the lat/lon boxes which contain each lake
!
!-----------------------------------------------------------------------
!
  k = 0

! Lake Superior
  k = k + 1
  lakelat1(k) = 46.4
  lakelat2(k) = 49.0
  lakelon1(k) = -92.5
  lakelon2(k) = -84.1
  lakeelev(k) = 182.9

! Lake Michigan
  k = k + 1
  lakelat1(k) = 41.0
  lakelat2(k) = 46.3
  lakelon1(k) = -88.1
  lakelon2(k) = -84.7
  lakeelev(k) = 176.5

! Lake Huron
  k = k + 1
  lakelat1(k) = 43.0
  lakelat2(k) = 46.4
  lakelon1(k) = -84.7
  lakelon2(k) = -79.5
  lakeelev(k) = 176.5

! Lake St. Clair
  k = k + 1
  lakelat1(k) = 42.2
  lakelat2(k) = 42.8
  lakelon1(k) = -83.0
  lakelon2(k) = -82.2
  lakeelev(k) = 174.6

! Lake Erie
  k = k + 1
  lakelat1(k) = 41.1
  lakelat2(k) = 43.0
  lakelon1(k) = -83.7
  lakelon2(k) = -78.6
  lakeelev(k) = 173.7

! Lake Ontario
  k = k + 1
  lakelat1(k) = 43.0
  lakelat2(k) = 45.6
  lakelon1(k) = -80.0
  lakelon2(k) = -76.0
  lakeelev(k) = 74.1

  nlake = k

  PRINT *, 'ARPSTRN: Number of Great Lakes: ', nlake
!
!-----------------------------------------------------------------------
!
!  Do no allow the elevation in a box to go below the lake's elevation.
!
!-----------------------------------------------------------------------
!
  DO j=1,ny
    DO i=1,nx
      latpnt = lat(i,j)
      lonpnt = lon(i,j)
      IF(lonpnt > 180) lonpnt = lonpnt - 360

      DO k=1,nlake
        IF (latpnt > lakelat1(k) .AND. latpnt < lakelat2(k) .AND.       &
                lonpnt > lakelon1(k) .AND. lonpnt < lakelon2(k) .AND.   &
                h(i,j) < lakeelev(k) ) THEN
!        write (*,'(a,2i4,2f8.1,i4,2f7.1)')
!    :           "FIXELEVL",i,j,lakeelev(k),h(i,j),k,latpnt,lonpnt
          h(i,j) = lakeelev(k)
        END IF
      END DO
!    write (*,'(a,2i4,f8.1,12x,2f7.1)')
!    :           "FIXELEV ",i,j,h(i,j),latpnt,lonpnt

    END DO
  END DO


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


SUBROUTINE extmnsnd(nx,ny,nx_ext,ny_ext,nz_ext,lvlprof,                 & 2,5
           iextmn,iextmx,jextmn,jextmx,kextmn,kextmx,                   &
           xscl,yscl,x_ext,y_ext,                                       &
           hgt_ext,pln_ext,t_ext,                                       &
           rhs_ext,u_ext,v_ext,                                         &
           zsnd,plsnd,psnd,tsnd,ptsnd,rhssnd,qvsnd,                     &
           rhosnd,usnd,vsnd)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Finds a mean sounding of all external data within the ARPS
!  domain.  This is used to define the model base state.
!  The data are interpolated to constant height levels and
!  averaged in the horizontal.   The final sounding is hydrostatic.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster
!  August, 1995  Replaces MNSOUND
!
!  MODIFICATION HISTORY:
!  12/20/95  Keith Brewster Added code to insure at least one
!            point goes into creating the mean sounding.
!
!  2000/04/05 (Gene Bassett)
!  Added qvsnd array.
!
!  04/02/2001 (K. Brewster)
!  Corrected calculations for qvsnd array that were also impacting
!  rhssnd output array.
!
!  Also, added parameter rhmin=0.05 (5 percent) as a minimum relative
!  humidity in place of hardcoded check value of 1.0E-04.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Sizing variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx,ny
  INTEGER :: nx_ext,ny_ext,nz_ext
  INTEGER :: lvlprof
!
!-----------------------------------------------------------------------
!
!  Limits of calculations
!
!-----------------------------------------------------------------------
!
  INTEGER :: iextmn,iextmx
  INTEGER :: jextmn,jextmx
  INTEGER :: kextmn,kextmx
!
!-----------------------------------------------------------------------
!
!  ARPS grid variables
!
!-----------------------------------------------------------------------
!
  REAL :: xscl(nx)
  REAL :: yscl(ny)
!
!-----------------------------------------------------------------------
!
!  External grid variables
!
!-----------------------------------------------------------------------
!
  REAL :: x_ext(nx_ext,ny_ext)
  REAL :: y_ext(nx_ext,ny_ext)
!
  REAL :: hgt_ext(nx_ext,ny_ext,nz_ext)         ! Height MSL
  REAL :: t_ext(nx_ext,ny_ext,nz_ext)           ! Temperature (K)
  REAL :: pln_ext(nx_ext,ny_ext,nz_ext)         ! ln of total pressure
  REAL :: rhs_ext(nx_ext,ny_ext,nz_ext)         ! RHstar=SQRT(1.-RH)
  REAL :: u_ext(nx_ext,ny_ext,nz_ext)           ! u wind component
  REAL :: v_ext(nx_ext,ny_ext,nz_ext)           ! v wind component
!
!-----------------------------------------------------------------------
!
!  Input sounding variables
!
!-----------------------------------------------------------------------
!
  REAL :: zsnd(lvlprof)
!
!-----------------------------------------------------------------------
!
!  Output sounding variables
!
!-----------------------------------------------------------------------
!
  REAL :: plsnd(lvlprof)         ! ln of pressure
  REAL :: psnd(lvlprof)          ! Pressure (Pa)
  REAL :: tsnd(lvlprof)          ! Temperature (K)
  REAL :: ptsnd(lvlprof)         ! Potential temperature (K)
  REAL :: rhssnd(lvlprof)        ! RHstar
  REAL :: qvsnd(lvlprof)         ! Specific humidity (kg/kg)
  REAL :: rhosnd(lvlprof)        ! Density
  REAL :: usnd(lvlprof)          ! u wind component
  REAL :: vsnd(lvlprof)          ! v wind component
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Lapse rate and a constant for hydrostatic integration
!
!-----------------------------------------------------------------------
!
  REAL :: gamma,pconst
  PARAMETER (gamma = 0.0068,      & ! 6.8 degrees per km
         pconst = (g/rd))
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  REAL, PARAMETER :: rhmin=0.05
  REAL, PARAMETER :: rhmax=1.0
  INTEGER :: i,j,k,kext,imid,jmid
  INTEGER :: knt,knttot,ktop,kbot
  REAL :: wlow,whigh,tvbar,xmid,ymid,dist2,distmn,accept
  REAL :: c1,c2,pres,qvsat,qvbot,qvtop,rh,epsilon
!
!-----------------------------------------------------------------------
!
!  Function f_qvsat and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_qvsat

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

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  epsilon=0.01*(zsnd(2)-zsnd(1))
!
!-----------------------------------------------------------------------
!
!  Zero-out all sounding arrays
!  ptsnd array is used temporarily to store a count.
!
!-----------------------------------------------------------------------
!
  DO k=1,lvlprof
    ptsnd(k)=0.
    plsnd(k)=0.
    tsnd(k)=0.
    rhssnd(k)=0.
    qvsnd(k)=0.
    usnd(k)=0.
    vsnd(k)=0.
  END DO
!
!-----------------------------------------------------------------------
!
!  Find the external data grid location that is closest
!  to the middle of the grid.
!  For the case where no external grid point lies within the
!  grid this grid column will be used as the mean sounding.
!
!-----------------------------------------------------------------------
!
  xmid=0.5*(xscl(1)+xscl(nx-1))
  ymid=0.5*(yscl(1)+yscl(ny-1))
  distmn=(x_ext(1,1)-xmid)*(x_ext(1,1)-xmid) +                          &
         (y_ext(1,1)-ymid)*(y_ext(1,1)-ymid)
  imid=1
  jmid=1
  DO j=1,ny_ext
    DO i=1,nx_ext
      dist2=(x_ext(i,j)-xmid)*(x_ext(i,j)-xmid) +                       &
            (y_ext(i,j)-ymid)*(y_ext(i,j)-ymid)
      IF(dist2 < distmn) THEN
        distmn=dist2
        imid=i
        jmid=j
      END IF
    END DO
  END DO
  iextmn=MIN(iextmn,imid)
  iextmx=MAX(iextmx,imid)
  jextmn=MIN(jextmn,jmid)
  jextmx=MAX(jextmx,jmid)
!
!-----------------------------------------------------------------------
!
!  Look at all external points possibly within domain
!
!-----------------------------------------------------------------------
!
  knt=0
  knttot=0
  DO j=jextmn,jextmx
    DO i=iextmn,iextmx
!
!-----------------------------------------------------------------------
!
!  Is this point within the ARPS domain?
!  Since the ARPS grid is Cartesian, need only compare
!  the external grid coordinates to the x and y limits
!  of the ARPS grid.
!
!-----------------------------------------------------------------------
!
      knttot=knttot+1
      IF((x_ext(i,j) >= xscl(1) .AND. x_ext(i,j) <= xscl(nx) .AND.      &
            y_ext(i,j) >= yscl(1) .AND. y_ext(i,j) <= yscl(ny)) .OR.    &
            (i == imid .AND. j == jmid)) THEN
        knt=knt+1
!
!-----------------------------------------------------------------------
!
!  Interpolate external data in vertical onto profile
!  arrays.
!
!-----------------------------------------------------------------------
!
        DO k=1,lvlprof
          IF(zsnd(k) >= hgt_ext(i,j,kextmx)) EXIT
          IF((hgt_ext(i,j,1)-zsnd(k)) < epsilon) THEN
            DO kext=kextmn,kextmx-1
              IF(hgt_ext(i,j,kext) >= zsnd(k)) EXIT
            END DO
            whigh=(zsnd(k)-hgt_ext(i,j,kext-1))/                        &
                  (hgt_ext(i,j,kext)-hgt_ext(i,j,kext-1))
            wlow=1.-whigh
            ptsnd(k)=ptsnd(k)+1.
            plsnd(k)=plsnd(k)+                                          &
                  whigh*pln_ext(i,j,kext)+wlow*pln_ext(i,j,kext-1)
            tsnd(k)=tsnd(k)+                                            &
                  whigh*t_ext(i,j,kext) + wlow*t_ext(i,j,kext-1)
            rhssnd(k)=rhssnd(k)+                                        &
                  whigh*rhs_ext(i,j,kext)+wlow*rhs_ext(i,j,kext-1)
            usnd(k)=usnd(k)+                                            &
                  whigh*u_ext(i,j,kext) + wlow*u_ext(i,j,kext-1)
            vsnd(k)=vsnd(k)+                                            &
                  whigh*v_ext(i,j,kext) + wlow*v_ext(i,j,kext-1)
          END IF
        END DO
      END IF
    END DO
  END DO

  WRITE(6,'(/a,i6,a,/a,i6,a/)')                                         &
       '  extmnsnd found ',knt,' points within ARPS domain',            &
       '  of ',knttot,' checked.'
!
!-----------------------------------------------------------------------
!
!  Find lowest height with data
!
!-----------------------------------------------------------------------
!
  WRITE(6,'(a)') ' Finding range of mean sounding data ...'
  accept=0.3*knt
  DO k=1,lvlprof-1
    WRITE(6,'(a,f10.2,a,f6.0,a,f10.0)') ' z = ',zsnd(k),                &
                  ' knt = ',ptsnd(k),' accept = ',accept
    IF(ptsnd(k) > accept) EXIT
  END DO
  kbot=k
!
!-----------------------------------------------------------------------
!
!  Find highest height with data
!
!-----------------------------------------------------------------------
!
  DO k=lvlprof,2,-1
    WRITE(6,'(a,f10.2,a,f6.0,a,f10.0)') ' z = ',zsnd(k),                &
                  ' knt = ',ptsnd(k),' accept = ',accept
    IF(ptsnd(k) > accept) EXIT
  END DO
  ktop=k
!
  WRITE(6,'(a,f10.2,a,f10.2,a)')                                        &
               ' Height of external data for mean spans from ',         &
                 zsnd(kbot),' to ',zsnd(ktop),' meters.'
!
!-----------------------------------------------------------------------
!
!  Divide through to find average
!
!-----------------------------------------------------------------------
!
  DO k=kbot,ktop
    plsnd(k)=plsnd(k)/ptsnd(k)
    tsnd(k)=tsnd(k)/ptsnd(k)
    rhssnd(k)=rhssnd(k)/ptsnd(k)
    usnd(k)=usnd(k)/ptsnd(k)
    vsnd(k)=vsnd(k)/ptsnd(k)
  END DO
!
!-----------------------------------------------------------------------
!
!  Set variables "below-ground"
!  Use a constant lapse rate, gamma.
!  plsnd is a sort of first-guess log(pressure), needed for qv
!  calculation from rhstar.
!
!-----------------------------------------------------------------------
!
  pres = EXP(plsnd(kbot))
  rh=MAX(rhmin,rhmax-(rhssnd(kbot)*rhssnd(kbot)))
  qvsat = f_qvsat( pres, tsnd(kbot) )
  qvbot=rh*qvsat
  c1=g/(rd*gamma)
  DO k=kbot-1,1,-1
    tsnd(k)=tsnd(kbot)-gamma*(zsnd(k)-zsnd(kbot))
    plsnd(k)=plsnd(kbot)+                                               &
        c1*ALOG((tsnd(kbot)-gamma*(zsnd(k)-zsnd(kbot)))/tsnd(kbot))
    psnd(k) = EXP(plsnd(k))
    qvsat = f_qvsat( psnd(k), tsnd(k) )
    rh=qvbot/qvsat
    rhssnd(k)=SQRT(MAX(0.,(rhmax-rh)))
    usnd(k)=usnd(kbot)
    vsnd(k)=vsnd(kbot)
  END DO
!
!-----------------------------------------------------------------------
!
!  Set variables "above-top"
!  Use a constant temperature.
!  We're assuming stratosphere here.
!
!-----------------------------------------------------------------------
!
  pres = EXP(plsnd(ktop))
  rh=MAX(rhmin,rhmax-(rhssnd(ktop)*rhssnd(ktop)))
  qvsat = f_qvsat( pres, tsnd(ktop) )
  qvtop=rh*qvsat
  c2=g/rd
  DO k=ktop+1,lvlprof
    tsnd(k)=tsnd(ktop)
    plsnd(k)=plsnd(ktop)-c2*(zsnd(k)-zsnd(ktop))/tsnd(ktop)
    psnd(k) = EXP(plsnd(k))
    qvsat = f_qvsat( psnd(k), tsnd(k) )
    rh=qvtop/qvsat
    rhssnd(k)=SQRT(MAX(0.,(rhmax-rh)))
    usnd(k)=usnd(ktop)
    vsnd(k)=vsnd(ktop)
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate qv profile from RH-star.
!  Temporarily use rhosnd to store virtual temperature
!
!-----------------------------------------------------------------------
!
  DO k=1,lvlprof
    pres = EXP(plsnd(k))
    rh=MAX(rhmin,rhmax-(rhssnd(k)*rhssnd(k)))
    qvsat = f_qvsat( pres, tsnd(k) )
    qvsnd(k)=rh*qvsat
    rhosnd(k) = tsnd(k)*(1.0+rvdrd*qvsnd(k))/(1.0+qvsnd(k))
  END DO
!
!-----------------------------------------------------------------------
!
!  Make sure that the sounding is hydrostatic by integrating
!  from kbot. rhosnd is really virtual temperature here.
!
!-----------------------------------------------------------------------
!
  psnd(kbot)=EXP(plsnd(kbot))
  DO k=kbot-1,1,-1
    tvbar=0.5*(rhosnd(k+1)+rhosnd(k))
    psnd(k)=psnd(k+1)*EXP(pconst*(zsnd(k+1)-zsnd(k))/tvbar)
  END DO
!
  DO k=kbot+1,lvlprof
    tvbar=0.5*(rhosnd(k-1)+rhosnd(k))
    psnd(k)=psnd(k-1)*EXP(pconst*(zsnd(k-1)-zsnd(k))/tvbar)
  END DO
!
!-----------------------------------------------------------------------
!
!  Derived variable calculations
!  Compute density from virtual temperature.
!  Compute potential temperature.
!  Compute log of pressure.
!
!-----------------------------------------------------------------------
!
  DO k=1,lvlprof
    rhosnd(k)=psnd(k)/(rd*rhosnd(k))
    ptsnd(k)=tsnd(k)*((p0/psnd(k))**rddcp)
!    pres=exp(plsnd(k))
!    print *,'psnd old, psnd new: ',pres,psnd(k),(psnd(k)-pres)
    plsnd(k)=ALOG(psnd(k))
  END DO
!
  RETURN
END SUBROUTINE extmnsnd

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


SUBROUTINE arpsstop(comment,outerr) 148,1

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Close down any necessary options (e.g. MPI) and stop execution of
!  the program.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/04/21
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    comment    Comment string.
!    outerr     Output error indicator (0-send to standard out,
!                  1-send to standard error)
!
!  OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------

  IMPLICIT NONE        ! Force explicit declarations

  CHARACTER (LEN=*) :: comment     ! comment
  INTEGER :: outerr       ! error code (0-standard out, 1-standard error)

!-----------------------------------------------------------------------
!
!  Miscellaneous local variables:
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'globcst.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  IF (outerr == 0) THEN
    WRITE (6,*) trim(comment)
  ELSE
    WRITE (0,*) trim(comment)
  END IF

  IF (mp_opt > 0) THEN
    CALL mpexit(outerr)
  END IF

  STOP
END SUBROUTINE arpsstop

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


SUBROUTINE wrtcomment(comment,outerr) 3

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write out a message to standard our or standard error.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/04/21
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    comment    Comment string.
!    outerr     Output error indicator (0-send to standard out,
!                  1-send to standard error)
!
!  OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------

  IMPLICIT NONE        ! Force explicit declarations

  CHARACTER (LEN=*) :: comment     ! comment
  INTEGER :: outerr       ! error code (0-standard out, 1-standard error)

!-----------------------------------------------------------------------
!
!  Miscellaneous local variables:
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'globcst.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  IF (outerr == 0) THEN
    IF (myproc == 0) WRITE (6,*) trim(comment)
  ELSE
    WRITE (0,*) trim(comment)
  END IF

  RETURN
END SUBROUTINE wrtcomment

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


SUBROUTINE acct_init 1,3

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Initialize cpu and wall clock timers.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/09/14
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------

  IMPLICIT NONE        ! Force explicit declarations

  REAL :: f_cputime
  DOUBLE PRECISION :: f_walltime

!-----------------------------------------------------------------------
!
!  Miscellaneous local variables:
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'globcst.inc'

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  use_acct = 1

  CALL init_walltime

  wall_init = f_walltime()
  cpu_init = dble(f_cputime())

  acct_cpu_time = 0d0
  acct_wall_time = 0d0

  current_acct = 0
  interrupt_acct = 0

  RETURN
END SUBROUTINE acct_init

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


SUBROUTINE set_acct(account) 75,3

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Start cpu and wall timers for the specified account number after
!  adding the values of the current account to its totals.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/09/14
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    account    Integer pointer to the account to start timing.
!
!  OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------

  IMPLICIT NONE        ! Force explicit declarations

  INTEGER :: account

  REAL :: f_cputime
  DOUBLE PRECISION :: f_walltime

!-----------------------------------------------------------------------
!
!  Miscellaneous local variables:
!
!-----------------------------------------------------------------------

  DOUBLE PRECISION :: cpu0, wall0

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'globcst.inc'

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  IF (use_acct == 0) RETURN

  cpu0 = dble(f_cputime())
  wall0 = f_walltime()

  IF (current_acct /= 0) THEN  ! close off current account

    acct_cpu_time(current_acct) =   &
       acct_cpu_time(current_acct) + cpu0 - cpu_acct
    acct_wall_time(current_acct) =   &
       acct_wall_time(current_acct) + wall0 - wall_acct
    
    IF (interrupt_acct /= 0) THEN  ! reset the interrupt

      CALL acct_interrupt(interrupt_acct)

    ENDIF

  ENDIF

  current_acct = account

  cpu_acct = cpu0
  wall_acct = wall0

  RETURN
END SUBROUTINE set_acct

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


SUBROUTINE acct_interrupt(account) 210,3

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Interrupt accounting for the main timer to start tracting another
!  account.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/09/14
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    account    Integer pointer to the account to start timing.  A value
!               of zero will turn off all but to total cpu/wall timers.
!
!  OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------

  IMPLICIT NONE        ! Force explicit declarations

  INTEGER :: account

  REAL :: f_cputime
  DOUBLE PRECISION :: f_walltime

!-----------------------------------------------------------------------
!
!  Miscellaneous local variables:
!
!-----------------------------------------------------------------------

  DOUBLE PRECISION :: cpu0, wall0

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'globcst.inc'

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  IF (use_acct == 0) RETURN

  cpu0 = dble(f_cputime())
  wall0 = f_walltime()

  IF (interrupt_acct /= 0) THEN  ! stop an existing interrupt
    CALL acct_stop_inter
  ENDIF

  interrupt_acct = account

  cpu_inter = cpu0
  wall_inter = wall0

  RETURN
END SUBROUTINE acct_interrupt

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


SUBROUTINE acct_stop_inter 132,2

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Stop the interrupt timer and adjust the main timer accordingly.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/09/14
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------

  IMPLICIT NONE        ! Force explicit declarations

  REAL :: f_cputime
  DOUBLE PRECISION :: f_walltime

!-----------------------------------------------------------------------
!
!  Miscellaneous local variables:
!
!-----------------------------------------------------------------------

  DOUBLE PRECISION :: cpu0, wall0

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'globcst.inc'

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  IF (use_acct == 0) RETURN

  cpu0 = dble(f_cputime())
  wall0 = f_walltime()

  IF (interrupt_acct /= 0) THEN 

    acct_cpu_time(current_acct) =   &
       acct_cpu_time(current_acct) - (cpu0 - cpu_inter)
    acct_wall_time(current_acct) =   &
       acct_wall_time(current_acct) - (wall0 - wall_inter)

    acct_cpu_time(interrupt_acct) =   &
       acct_cpu_time(interrupt_acct) + cpu0 - cpu_inter
    acct_wall_time(interrupt_acct) =   &
       acct_wall_time(interrupt_acct) + wall0 - wall_inter

  ENDIF

  interrupt_acct = 0

  RETURN
END SUBROUTINE acct_stop_inter

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


SUBROUTINE acct_finish 1,8

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Stop cpu/wall timers and calculate totals.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/09/14
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------

  IMPLICIT NONE        ! Force explicit declarations

  REAL :: f_cputime
  DOUBLE PRECISION :: f_walltime

!-----------------------------------------------------------------------
!
!  Miscellaneous local variables:
!
!-----------------------------------------------------------------------

  DOUBLE PRECISION :: cpu0, wall0
  REAL :: time
  INTEGER :: acct

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'globcst.inc'
  INCLUDE 'mp.inc'

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  IF (use_acct == 0) RETURN

  cpu0 = dble(f_cputime())
  wall0 = f_walltime()

  total_cpu = cpu0 - cpu_init
  total_wall = wall0 - wall_init

  IF (interrupt_acct /= 0) THEN  ! stop an existing interrupt
    CALL acct_stop_inter
  END IF

  IF (current_acct /= 0) THEN  
    CALL set_acct(0)
  END IF

!-----------------------------------------------------------------------
!
!  Sum up the CPU times used by all processors if on a multi-processor
!  system.
!
!-----------------------------------------------------------------------

  IF (mp_opt > 0) THEN
    time = real(total_wall)
    CALL mptotal(time)
    total_wall = time
    time = real(total_cpu)
    CALL mptotal(time)
    total_cpu = time
    DO acct = 1,max_acct
      time = real(acct_cpu_time(acct))
      CALL mptotal(time)
      acct_cpu_time(acct) = dble(time)
      time = real(acct_wall_time(acct))
      CALL mptotal(time)
      acct_wall_time(acct) = dble(time)
    END DO
  END IF

  RETURN
END SUBROUTINE acct_finish

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


SUBROUTINE acct_report_arps 1,4

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Output the timing statistics for the ARPS model.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/09/14
!
!  MODIFICATION HISTORY:
!
!  13 March 2002 (Eric Kemp)
!  Added time accounting for WRF BMJ scheme.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  OUTPUT:
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------

  IMPLICIT NONE        ! Force explicit declarations

!-----------------------------------------------------------------------
!
!  Miscellaneous local variables:
!
!-----------------------------------------------------------------------

  REAL :: tc, tw

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'globcst.inc'
  INCLUDE 'mp.inc'

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
!  Print out CPU usage.
!
!-----------------------------------------------------------------------
!
  CALL flush(6) !wdt
!wdt update
  ! A barrier before and after printing stats avoids spurious output
  ! in the body of the stats output.
  CALL mpbarrier

1000 FORMAT  &
        (a/a,25(/1x,a,2x,1p,e13.6,a,2x,0p,f6.2,a,4x,1p,e13.6,a,2x,0p,f6.2,a))
  IF (myproc == 0) THEN


    IF(total_cpu  == 0d0) total_cpu =1d0
    IF(total_wall == 0d0) total_wall=1d0
    tc = 100.0/total_cpu
    tw = 100.0/total_wall

    WRITE (6,'(/a/)') "ARPS CPU Summary:"
    WRITE  (6,1000) &
    ' Process             CPU time                   WALL CLOCK time',  &
    ' -----------------   ----------------------     ----------------------',&
    'Initialization  :',acct_cpu_time(init_acct),                     's', &
     acct_cpu_time(init_acct)*tc,  '%', acct_wall_time(init_acct),    's', &
     acct_wall_time(init_acct)*tw, '%',                                    &
    'Data output     :',acct_cpu_time(output_acct),                   's', &
     acct_cpu_time(output_acct)*tc,'%', acct_wall_time(output_acct),  's', &
     acct_wall_time(output_acct)*tw,'%',                                   &
    'Wind advection  :',acct_cpu_time(advuvw_acct),                   's', &
     acct_cpu_time(advuvw_acct)*tc,'%', acct_wall_time(advuvw_acct),  's', &
     acct_wall_time(advuvw_acct)*tw,'%',                                   &
    'Scalar advection:',acct_cpu_time(advs_acct),                     's', &
     acct_cpu_time(advs_acct)*tc,  '%', acct_wall_time(advs_acct),    's', &
     acct_wall_time(advs_acct)*tw, '%',                                    &
    'Coriolis force  :',acct_cpu_time(coriol_acct),                   's', &
     acct_cpu_time(coriol_acct)*tc,'%', acct_wall_time(coriol_acct),  's', &
     acct_wall_time(coriol_acct)*tw,'%',                                   &
    'Buoyancy term   :',acct_cpu_time(buoy_acct),                     's', &
     acct_cpu_time(buoy_acct)*tc,  '%', acct_wall_time(buoy_acct),    's', &
     acct_wall_time(buoy_acct)*tw,'%',                                     &
    'Misc Large tstep:',acct_cpu_time(tinteg_acct),                   's', &
     acct_cpu_time(tinteg_acct)*tc,'%', acct_wall_time(tinteg_acct),  's', &
     acct_wall_time(tinteg_acct)*tw,'%',                                   &
    'Small time steps:',acct_cpu_time(smlstp_acct),                   's', &
     acct_cpu_time(smlstp_acct)*tc,'%', acct_wall_time(smlstp_acct),  's', &
     acct_wall_time(smlstp_acct)*tw,'%',                                   &
    'Radiation       :',acct_cpu_time(rad_acct),                      's', &
     acct_cpu_time(rad_acct)*tc,   '%', acct_wall_time(rad_acct),     's', &
     acct_wall_time(rad_acct)*tw, '%',                                     &
    'Soil model      :',acct_cpu_time(soil_acct),                     's', &
     acct_cpu_time(soil_acct)*tc,  '%', acct_wall_time(soil_acct),    's', &
     acct_wall_time(soil_acct)*tw,'%',                                     &
    'Surface physics :',acct_cpu_time(sfcphy_acct),                   's', &
     acct_cpu_time(sfcphy_acct)*tc,'%', acct_wall_time(sfcphy_acct),  's', &
     acct_wall_time(sfcphy_acct)*tw,'%',                                   &
    'Turbulence      :',acct_cpu_time(tmix_acct),                     's', &
     acct_cpu_time(tmix_acct)*tc,  '%', acct_wall_time(tmix_acct),    's', &
     acct_wall_time(tmix_acct)*tw,'%',                                     &
    'Comput. mixing  :',acct_cpu_time(cmix_acct),                     's', &
     acct_cpu_time(cmix_acct)*tc,  '%', acct_wall_time(cmix_acct),    's', &
     acct_wall_time(cmix_acct)*tw,'%',                                     &
    'Rayleigh damping:',acct_cpu_time(raydmp_acct),                   's', &
     acct_cpu_time(raydmp_acct)*tc,'%', acct_wall_time(raydmp_acct),  's', &
     acct_wall_time(raydmp_acct)*tw,'%',                                   &
    'TKE src terms   :',acct_cpu_time(tkesrc_acct),                   's', &
     acct_cpu_time(tkesrc_acct)*tc,'%', acct_wall_time(tkesrc_acct),  's', &
     acct_wall_time(tkesrc_acct)*tw,'%',                                   &
    'Gridscale precp.:',acct_cpu_time(qpfgrd_acct),                   's', &
     acct_cpu_time(qpfgrd_acct)*tc,'%', acct_wall_time(qpfgrd_acct),  's', &
     acct_wall_time(qpfgrd_acct)*tw,'%',                                   &
    'Kuo cumulus     :',acct_cpu_time(kuocum_acct),                   's', &
     acct_cpu_time(kuocum_acct)*tc,'%', acct_wall_time(kuocum_acct),  's', &
     acct_wall_time(kuocum_acct)*tw,'%',                                   &
    'Kain-Fritsch    :',acct_cpu_time(kfcum_acct),                    's', &
     acct_cpu_time(kfcum_acct)*tc, '%', acct_wall_time(kfcum_acct),   's', &
     acct_wall_time(kfcum_acct)*tw, '%',                                   &
    'BMJ cumulus     :',acct_cpu_time(bmjcum_acct),'s',                    &
     acct_cpu_time(bmjcum_acct)*tc, '%', acct_wall_time(bmjcum_acct), 's', &
     acct_wall_time(bmjcum_acct)*tw, '%',                                  &
    'Warmrain microph:',acct_cpu_time(warmrain_acct),                 's', &
     acct_cpu_time(warmrain_acct)*tc,'%',acct_wall_time(warmrain_acct),'s',&
     acct_wall_time(warmrain_acct)*tw,'%',                                 &
    'Lin ice microph :',acct_cpu_time(linice_acct),                   's', &
     acct_cpu_time(linice_acct)*tc,'%', acct_wall_time(linice_acct),  's', &
     acct_wall_time(linice_acct)*tw,'%',                                   &
    'NEM ice microph :',acct_cpu_time(nemice_acct),                   's', &
     acct_cpu_time(nemice_acct)*tc,'%', acct_wall_time(nemice_acct),  's', &
     acct_wall_time(nemice_acct)*tw,'%',                                   &
    'Hydrometero fall:',acct_cpu_time(qfall_acct),                    's', &
     acct_cpu_time(qfall_acct)*tc, '%', acct_wall_time(qfall_acct),   's', &
     acct_wall_time(qfall_acct)*tw, '%',                                   &
    'Bound.conditions:',acct_cpu_time(bc_acct),                       's', &
     acct_cpu_time(bc_acct)*tc,    '%', acct_wall_time(bc_acct),      's', &
     acct_wall_time(bc_acct)*tw,  '%',                                     &
    'Message passing :',acct_cpu_time(mp_acct),                       's', &
     acct_cpu_time(mp_acct)*tc,    '%', acct_wall_time(mp_acct),      's', &
     acct_wall_time(mp_acct)*tw,  '%',                                     &
    'Miscellaneous   :',acct_cpu_time(misc_acct),                     's', &
     acct_cpu_time(misc_acct)*tc,  '%', acct_wall_time(misc_acct),    's', &
     acct_wall_time(misc_acct)*tw,'%'  

    WRITE(6,'(/1x,a,1p,e13.6,a,2x,0p,7x,4x,1p,e13.6,a,2x,0p,6x)')  &
       'Entire model    :',total_cpu,'s',total_wall,'s'

    WRITE(6,'(/1x,a,1p,e13.6,a,2x,0p,7x,4x,1p,e13.6,a,2x,0p,6x)')  &
       'Without Init/IO :',  &
       total_cpu-acct_cpu_time(init_acct)-acct_cpu_time(output_acct),'s',  &
       total_wall-acct_wall_time(init_acct)-acct_wall_time(output_acct),'s'

  END IF
  CALL flush(6) !wdt
!wdt update
  CALL mpbarrier

  RETURN
END SUBROUTINE acct_report_arps

!wdt update

SUBROUTINE init_walltime 1

  INTEGER :: lastval, rollover

  COMMON /walltime/ lastval, rollover

  lastval = 0
  rollover = 0

  RETURN
END SUBROUTINE init_walltime


DOUBLE PRECISION FUNCTION f_walltime() 5
!
!-----------------------------------------------------------------------
!
!  Function to measure wall clock time.
!
!-----------------------------------------------------------------------
!
  INTEGER :: lastval, rollover

  INTEGER :: wall_int
  INTEGER :: ticspersec,max_count

  COMMON /walltime/ lastval, rollover

  CALL system_clock(wall_int,ticspersec,max_count)

  IF (wall_int < lastval .and. lastval > max_count*0.01) THEN
  ! Assume that the clock has rolled over once.  The test for max_count*0.1
  ! is for Absoft which sometimes seems to have wall_int go backwards a
  ! little at times.
    WRITE (6,*) "F_WALLTIME: wall clock assumed to have rolled over once."
    rollover = rollover + 1
  ENDIF
  lastval = wall_int
  f_walltime = dble(wall_int)/dble(ticspersec)  &
           + (dble(max_count)/dble(ticspersec))*dble(rollover)

  RETURN
END FUNCTION f_walltime

!wdt update - subroutine moved from another file to here:
!##################################################################
!##################################################################
!######                                                      ######
!######              SUBROUTINE FIX_SOIL_NSTYP               ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################


SUBROUTINE fix_soil_nstyp(nx,ny,nstypin,nstyp,  & 7
      tsfc,tsoil,wetsfc,wetdp,wetcanp)
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Adjust soil variables if nstyp differs between what is defined in 
!  ARPS and what is in the data file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/08/11
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx
!  ny
!  nstypin      Number of soil types in the data to be imported
!  nstyp        Number of soil types to use in the arps
!
!  INPUT/OUTPUT:
!
!  tsfc         Temperature at surface (K)
!  tsoil        Deep soil temperature (K)
!  wetsfc       Surface soil moisture
!  wetdp        Deep soil moisture
!  wetcanp      Canopy water amount
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER :: nx
  INTEGER :: ny
  INTEGER :: nstypin,nstyp

  REAL :: tsfc   (nx,ny,0:nstyp)      ! Temperature at surface (K)
  REAL :: tsoil  (nx,ny,0:nstyp)      ! Deep soil temperature (K)
  REAL :: wetsfc (nx,ny,0:nstyp)      ! Surface soil moisture
  REAL :: wetdp  (nx,ny,0:nstyp)      ! Deep soil moisture
  REAL :: wetcanp(nx,ny,0:nstyp)      ! Canopy water amount
 
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: is

!-----------------------------------------------------------------------
!
!  Include file:
!
!-----------------------------------------------------------------------

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

!-----------------------------------------------------------------------
!
! Only one soil type, use average values.
!
!-----------------------------------------------------------------------

  IF ( nstypin == 1 ) THEN
    tsfc   (1:nx-1,1:ny-1,1) = tsfc   (1:nx-1,1:ny-1,0)
    tsoil  (1:nx-1,1:ny-1,1) = tsoil  (1:nx-1,1:ny-1,0)
    wetsfc (1:nx-1,1:ny-1,1) = wetsfc (1:nx-1,1:ny-1,0)
    wetdp  (1:nx-1,1:ny-1,1) = wetdp  (1:nx-1,1:ny-1,0)
    wetcanp(1:nx-1,1:ny-1,1) = wetcanp(1:nx-1,1:ny-1,0)
  END IF

!-----------------------------------------------------------------------
!
! Use average values for soil types beyond that provided in the data.
!
!-----------------------------------------------------------------------

  IF ( nstypin < nstyp ) THEN
    DO is=nstypin+1,nstyp
      tsfc   (1:nx-1,1:ny-1,is) = tsfc   (1:nx-1,1:ny-1,0)
      tsoil  (1:nx-1,1:ny-1,is) = tsoil  (1:nx-1,1:ny-1,0)
      wetsfc (1:nx-1,1:ny-1,is) = wetsfc (1:nx-1,1:ny-1,0)
      wetdp  (1:nx-1,1:ny-1,is) = wetdp  (1:nx-1,1:ny-1,0)
      wetcanp(1:nx-1,1:ny-1,is) = wetcanp(1:nx-1,1:ny-1,0)
    ENDDO
  END IF

  RETURN

END SUBROUTINE fix_soil_nstyp

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


SUBROUTINE fix_stypfrct_nstyp(nx,ny,nstypin,nstyp,stypfrct) 6

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Adjust stypfrct if the number of soil types read in differs from
!  the number specified for arps to use, nstyp.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/08/11
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx
!  ny
!  nstypin      Number of soil types in the data file
!  nstyp        Number of soil types to use in the arps
!
!  INPUT/OUTPUT:
!
!  stypfrct     Soil type fraction
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER :: nx
  INTEGER :: ny
  INTEGER :: nstypin,nstyp
  REAL :: stypfrct(nx,ny,nstyp)

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: i,j,is
  REAL :: frctot

!-----------------------------------------------------------------------
!
!  Include file:
!
!-----------------------------------------------------------------------

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

!-----------------------------------------------------------------------
!
! Readjust stypfrct if only part of the soil type were read in.
!
!-----------------------------------------------------------------------

  IF (nstypin > nstyp) THEN
    DO j=1,ny
      DO i=1,nx
        frctot = 0.0
        DO is=1,nstyp
          frctot = frctot + stypfrct(i,j,is)
        END DO
        IF (frctot == 0) THEN
          frctot = 1.
          stypfrct(i,j,1) = 1.
        ENDIF
        DO is=1,nstyp
          stypfrct(i,j,is) = stypfrct(i,j,is) / frctot
        END DO
      END DO
    END DO
  ENDIF

!-----------------------------------------------------------------------
!
! Fill in any unset values.
!
!-----------------------------------------------------------------------

  IF ( nstypin < nstyp ) THEN
    stypfrct(1:nx-1,1:ny-1,nstypin+1:nstyp) = 0.
  END IF

!-----------------------------------------------------------------------
!
! Make sure the the fraction totals up to 1 (or 0)
!
!-----------------------------------------------------------------------

  DO j=1,ny
    DO i=1,nx
      frctot = 0.
      DO is = 1,nstyp
        frctot = frctot + stypfrct(i,j,is)
      END DO
      IF (frctot /= 0) THEN
        DO is = 1,nstyp
          stypfrct(i,j,is) = stypfrct(i,j,is) / frctot
        END DO
      END IF
    END DO
  END DO

  RETURN

END SUBROUTINE fix_stypfrct_nstyp

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!##################################################################
!##################################################################
!######                                                      ######
!######              SUBROUTINE REMAP_SOIL_VARS              ######
!######                                                      ######
!######                    Developed by                      ######
!######            Weather Decision Technologies             ######
!######                                                      ######
!##################################################################
!##################################################################


SUBROUTINE remap_soil_vars(nx,ny,nstypin,nstyp,  & 1
       tsfc_in,tsoil_in,wetsfc_in,wetdp_in,wetcanp_in,soiltyp_in,  &
       tsfcin,tsoilin,wsfcin,wdpin,wcanpin,  &
       tsfc,tsoil,wetsfc,wetdp,wetcanp,soiltyp)
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Rearrange soil variables to be consistent with any difference between
!  the soil types to be used and the ones in the data set.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/10/27
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx
!  ny
!  nstypin      Number of soil types in the data file
!  nstyp        Number of soil types to use in the arps
!  tsfcin
!  tsoilin
!  wsfcin
!  wdpin
!  wcanpin
!  tsfc_in      Temperature at surface in data set (K)
!  tsoil_in     Deep soil temperature in data set (K)
!  wetsfc_in    Surface soil moisture in data set
!  wetdp_in     Deep soil moisture in data set
!  wetcanp_in   Canopy water amount in data set
!  soiltyp_in   Soil type in data set
!
!  OUTPUT:
!
!  tsfc         Temperature at surface (K)
!  tsoil        Deep soil temperature (K)
!  wetsfc       Surface soil moisture
!  wetdp        Deep soil moisture
!  wetcanp      Canopy water amount
!  soiltyp      Soil type
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER :: nx
  INTEGER :: ny
  INTEGER :: nstypin,nstyp
  INTEGER ::  tsfcin,tsoilin,wsfcin,wdpin,wcanpin

  REAL :: tsfc_in   (nx,ny,0:nstypin)   ! Temperature at surface (K)
  REAL :: tsoil_in  (nx,ny,0:nstypin)   ! Deep soil temperature (K)
  REAL :: wetsfc_in (nx,ny,0:nstypin)   ! Surface soil moisture
  REAL :: wetdp_in  (nx,ny,0:nstypin)   ! Deep soil moisture
  REAL :: wetcanp_in(nx,ny,0:nstypin)   ! Canopy water amount
  INTEGER :: soiltyp_in(nx,ny,nstypin)  ! Soil type in model domain

  REAL :: tsfc   (nx,ny,0:nstyp)      ! Temperature at surface (K)
  REAL :: tsoil  (nx,ny,0:nstyp)      ! Deep soil temperature (K)
  REAL :: wetsfc (nx,ny,0:nstyp)      ! Surface soil moisture
  REAL :: wetdp  (nx,ny,0:nstyp)      ! Deep soil moisture
  REAL :: wetcanp(nx,ny,0:nstyp)      ! Canopy water amount
  INTEGER :: soiltyp(nx,ny,nstyp)     ! Soil type in model domain
 
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: i,j,k,is,isin
  INTEGER :: ifound_styp

!-----------------------------------------------------------------------
!
!  Include file:
!
!-----------------------------------------------------------------------

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

!-----------------------------------------------------------------------
!
! Match the soil type in the old data set.  If not found then set to 
! average soil properties (except water).
!
!-----------------------------------------------------------------------

  DO j=1,ny
    DO i=1,nx
      DO is=1,nstyp
        ifound_styp = 0
        DO isin=1,nstypin
          IF (soiltyp(i,j,is) == soiltyp_in(i,j,isin)) THEN
            IF (tsfcin == 1) tsfc(i,j,is) = tsfc_in(i,j,isin)
            IF (tsoilin == 1) tsoil(i,j,is) = tsoil_in(i,j,isin)
            IF (wsfcin == 1) wetsfc(i,j,is) = wetsfc_in(i,j,isin)
            IF (wdpin == 1) wetdp(i,j,is) = wetdp_in(i,j,isin)
            IF (wcanpin == 1) wetcanp(i,j,is) = wetcanp_in(i,j,isin)
            ifound_styp = 1
            EXIT
          END IF
        END DO
        IF (ifound_styp == 0) THEN
          IF (tsfcin == 1) tsfc(i,j,is) = tsfc_in(i,j,0)
          IF (tsoilin == 1) tsoil(i,j,is) = tsoil_in(i,j,0)
          IF (wsfcin == 1) wetsfc(i,j,is) = wetsfc_in(i,j,0)
          IF (wdpin == 1) wetdp(i,j,is) = wetdp_in(i,j,0)
          IF (wcanpin == 1) wetcanp(i,j,is) = wetcanp_in(i,j,0)
          
          IF (tsfcin == 1 .and. soiltyp(i,j,is) == 13) THEN
            tsfc(i,j,is) = tsoil_in(i,j,0)
          END IF
        END IF
      END DO
    END DO
  END DO

END SUBROUTINE remap_soil_vars

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 2000/11/01 GMB soil consistency:
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE INTERP_SOIL               ######
!######                                                      ######
!######                    Developed by                      ######
!######            Weather Decision Technologies             ######
!######                                                      ######
!##################################################################
!##################################################################


SUBROUTINE intrp_soil(nx,ny,nx1,ny1,nstyp,nstyp1,wx,wy,ix,jy,  & 2
                      tsfc,tsoil,wetsfc,wetdp,wetcanp,  &
                      soiltyp,stypfrct,vegtyp,  &
                      tsfc1,tsoil1,wetsfc1,wetdp1,wetcanp1,  &
                      soiltyp1,stypfrct1,vegtyp1)
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Interpolate soil properties onto another grid.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/11/01
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx
!  ny
!  nx1
!  ny1
!  xw           Weight factor in x-direction
!  yw           Weight factor in y-direction
!  ix           Old grid index (lower left) for new grid point
!  jy           Old grid index (lower left) for new grid point
!  nstyp        Number of soil types on original grid
!  nstyp1       Number of soil types to use on new grid
!  tsfc         Temperature at surface in data set (K)
!  tsoil        Deep soil temperature in data set (K)
!  wetsfc       Surface soil moisture in data set
!  wetdp        Deep soil moisture in data set
!  wetcanp      Canopy water amount in data set
!  soiltyp      Soil type in data set
!  stypfrct     Soil type fraction
!  vegtyp       Vegetation type in data set
!
!  OUTPUT:
!
!  tsfc1        Temperature at surface in data set (K)
!  tsoil1       Deep soil temperature in data set (K)
!  wetsfc1      Surface soil moisture in data set
!  wetdp1       Deep soil moisture in data set
!  wetcanp1     Canopy water amount in data set
!  soiltyp1     Soil type in data set
!  stypfrct1    Soil type fraction
!  vegtyp1      Vegetation type in data set
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER :: nx,ny,nx1,ny1
  INTEGER :: nstyp,nstyp1

  REAL :: wx(nx1,ny1)    ! Weight factor in x-direction
  REAL :: wy(nx1,ny1)    ! Weight factor in y-direction
  INTEGER :: ix(nx1,ny1) ! Old grid index (lower left) for new grid point
  INTEGER :: jy(nx1,ny1) ! Old grid index (lower left) for new grid point

  REAL :: tsfc   (nx,ny,0:nstyp)   ! Temperature at surface (K)
  REAL :: tsoil  (nx,ny,0:nstyp)   ! Deep soil temperature (K)
  REAL :: wetsfc (nx,ny,0:nstyp)   ! Surface soil moisture
  REAL :: wetdp  (nx,ny,0:nstyp)   ! Deep soil moisture
  REAL :: wetcanp(nx,ny,0:nstyp)   ! Canopy water amount
  INTEGER :: soiltyp(nx,ny,nstyp)  ! Soil type in model domain
  REAL :: stypfrct(nx,ny,nstyp)
  INTEGER :: vegtyp(nx,ny)

  REAL :: tsfc1   (nx1,ny1,0:nstyp1)      ! Temperature at surface (K)
  REAL :: tsoil1  (nx1,ny1,0:nstyp1)      ! Deep soil temperature (K)
  REAL :: wetsfc1 (nx1,ny1,0:nstyp1)      ! Surface soil moisture
  REAL :: wetdp1  (nx1,ny1,0:nstyp1)      ! Deep soil moisture
  REAL :: wetcanp1(nx1,ny1,0:nstyp1)      ! Canopy water amount
  INTEGER :: soiltyp1(nx1,ny1,nstyp1)     ! Soil type in model domain
  REAL :: stypfrct1(nx1,ny1,nstyp1)
  INTEGER :: vegtyp1(nx1,ny1)
 
!-----------------------------------------------------------------------
!
!  Include file:
!
!-----------------------------------------------------------------------

  INCLUDE 'arpssfc.inc'

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  ! Arrays start at zero in case no soil type defined (i.e. soiltyp=0)
  REAL :: tsfc1sum    (0:nsoiltyp)     ! Ground sfc. temperature (K)
  REAL :: tsoil1sum   (0:nsoiltyp)     ! Deep soil temperature (K)
  REAL :: wetsfc1sum  (0:nsoiltyp)     ! Surface soil moisture
  REAL :: wetdp1sum   (0:nsoiltyp)     ! Deep soil moisture
  REAL :: wetcanp1sum (0:nsoiltyp)     ! Canopy water amount
  REAL :: stypfrct1sum(0:nsoiltyp)     ! Frction of soil type

  INTEGER :: i,j,i1,j1,is,ii
  REAL :: weight,maxweight,frctot
  REAL :: totweight(0:nsoiltyp)
  INTEGER :: maxtype
  INTEGER :: soiltype

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  tsfc1 = 0
  tsoil1 = 0
  wetsfc1 = 0
  wetdp1 = 0
  wetcanp1 = 0
  soiltyp1 = 0
  stypfrct1 = 0
  vegtyp1 = 0

  IF (nstyp == 1) THEN
    stypfrct = 1   ! make sure stypfrct is set for nstyp=1
    ! copy level 0 to level 1 if level 1 is undefined:
    IF (tsfc(1,1,1) <= 0)   tsfc(:,:,1) = tsfc(:,:,0)
    IF (tsoil(1,1,1) <= 0)  tsoil(:,:,1) = tsoil(:,:,0)
    IF (wetsfc(1,1,1) < 0)  wetsfc(:,:,1) = wetsfc(:,:,0)
    IF (wetdp(1,1,1) < 0)   wetdp(:,:,1) = wetdp(:,:,0)
    IF (wetcanp(1,1,1) < 0) wetcanp(:,:,1) = wetcanp(:,:,0)
  ENDIF

  DO j1=1,ny1-1
    DO i1=1,nx1-1
      tsfc1sum = 0.
      tsoil1sum = 0.
      wetsfc1sum = 0.
      wetdp1sum = 0.
      wetcanp1sum = 0.
      stypfrct1sum = 0.
      !vegtyp1sum = 0.
      maxweight = 0.
      totweight = 0.
      vegtyp1(i1,j1) = 0
      DO j=jy(i1,j1),jy(i1,j1)+1
        DO i=ix(i1,j1),ix(i1,j1)+1
          weight = wx(i1,j1)*(1.+ix(i1,j1)-i)*wy(i1,j1)*(1+jy(i1,j1)-j)  &
                 + (1.-wx(i1,j1))*(i-ix(i1,j1))*wy(i1,j1)*(1+jy(i1,j1)-j)  &
                 + wx(i1,j1)*(1.+ix(i1,j1)-i)*(1.-wy(i1,j1))*(j-jy(i1,j1))  &
                 + (1.-wx(i1,j1))*(i-ix(i1,j1))*(1.-wy(i1,j1))*(j-jy(i1,j1))
          DO is=1,nstyp
            IF (stypfrct(i,j,is) > 0) THEN
              soiltype = soiltyp(i,j,is)
              IF (soiltype < 1 .or. soiltype > nsoiltyp) soiltype = 0
              tsfc1sum(soiltyp(i,j,is)) = tsfc1sum(soiltyp(i,j,is))  &
                 + weight*stypfrct(i,j,is)*tsfc(i,j,is)
              tsoil1sum(soiltyp(i,j,is)) = tsoil1sum(soiltyp(i,j,is))  &
                 + weight*stypfrct(i,j,is)*tsoil(i,j,is)
              wetsfc1sum(soiltyp(i,j,is)) = wetsfc1sum(soiltyp(i,j,is))  &
                 + weight*stypfrct(i,j,is)*wetsfc(i,j,is)
              wetdp1sum(soiltyp(i,j,is)) = wetdp1sum(soiltyp(i,j,is))  &
                 + weight*stypfrct(i,j,is)*wetdp(i,j,is)
              wetcanp1sum(soiltyp(i,j,is)) = wetcanp1sum(soiltyp(i,j,is))  &
                 + weight*stypfrct(i,j,is)*wetcanp(i,j,is)
              stypfrct1sum(soiltyp(i,j,is)) = stypfrct1sum(soiltyp(i,j,is))  &
                 + weight*stypfrct(i,j,is)
              totweight(soiltyp(i,j,is)) = totweight(soiltyp(i,j,is)) + weight
            END IF
          END DO
          ! use the vegtyp of the closest point
          IF (weight > maxweight) THEN
            vegtyp1(i1,j1) = vegtyp(i,j)
            maxweight = weight
          END IF
          !!use most popular vegtyp weighted by distance
          !IF (vegtyp(i,j) > 0 .and. vegtyp(i,j) <= nvegtyp) THEN
          !  vegtyp1sum(vegtyp(i,j)) =  vegtyp1sum(vegtyp(i,j)) + weight
          !END IF
        END DO
      END DO

      !!use most popular vegtyp weighted by distance
      !maxweight = 0
      !vegtyp1(i1,j1) = 0
      !DO ii = 1,nvegtyp
      !  IF (vegtyp1sum(ii) > maxweight) THEN
      !    vegtyp1(i1,j1) = ii
      !    maxweight = vegtyp1sum(ii)
      !  END IF
      !END DO

      ! get the nstyp1 largest weighted soil types

      DO is = 1,nstyp1
        maxtype = -1
        maxweight = 0.
        DO ii = 0,nsoiltyp
          IF (stypfrct1sum(ii) > maxweight) THEN
            maxweight = stypfrct1sum(ii)
            maxtype = ii
          END IF
        END DO
        IF (maxtype /= -1) THEN
          soiltyp1(i1,j1,is) = maxtype
          tsfc1(i1,j1,is) = tsfc1sum(maxtype)/stypfrct1sum(maxtype)
          tsoil1(i1,j1,is) = tsoil1sum(maxtype)/stypfrct1sum(maxtype)
          wetsfc1(i1,j1,is) = wetsfc1sum(maxtype)/stypfrct1sum(maxtype)
          wetdp1(i1,j1,is) = wetdp1sum(maxtype)/stypfrct1sum(maxtype)
          wetcanp1(i1,j1,is) = wetcanp1sum(maxtype)/stypfrct1sum(maxtype)
          stypfrct1(i1,j1,is) = stypfrct1sum(maxtype)/totweight(maxtype)
          stypfrct1sum(maxtype) = 0.
        ELSE 
          IF (is == 1) THEN
            WRITE (6,*)  &
             "INTRP_SOIL: WARNING, no soil type found, variables not assigned!"
          ELSE
            stypfrct1(i1,j1,is) = 0.
          ENDIF
        ENDIF
      END DO

      ! renormalize 

      frctot = 0.
      DO is = 1,nstyp1
        frctot = frctot + stypfrct1(i1,j1,is)
      END DO
      IF (frctot /= 0) THEN
        DO is = 1,nstyp1
          stypfrct1(i1,j1,is) = stypfrct1(i1,j1,is) / frctot
        END DO
      ELSE
        stypfrct1(i1,j1,1) = 1.
        IF (nstyp1 .gt. 1) stypfrct1(i1,j1,2:nstyp1) = 0.
      END IF
      tsfc1(i1,j1,0) = 0.
      tsoil1(i1,j1,0) = 0.
      wetsfc1(i1,j1,0) = 0.
      wetdp1(i1,j1,0) = 0.
      wetcanp1(i1,j1,0) = 0.
      DO is = 1,nstyp1
        tsfc1(i1,j1,0) = tsfc1(i1,j1,0)  &
                         + stypfrct1(i1,j1,is)*tsfc1(i1,j1,is)
        tsoil1(i1,j1,0) = tsoil1(i1,j1,0)  &
                         + stypfrct1(i1,j1,is)*tsoil1(i1,j1,is)
        wetsfc1(i1,j1,0) = wetsfc1(i1,j1,0)  &
                         + stypfrct1(i1,j1,is)*wetsfc1(i1,j1,is)
        wetdp1(i1,j1,0) = wetdp1(i1,j1,0)  &
                         + stypfrct1(i1,j1,is)*wetdp1(i1,j1,is)
        wetcanp1(i1,j1,0) = wetcanp1(i1,j1,0)  &
                         + stypfrct1(i1,j1,is)*wetcanp1(i1,j1,is)
      END DO
!FIXME: remove? GMB
!     tsfc1(i1,j1,0) = tsfc1(i1,j1,0)
!     tsoil1(i1,j1,0) = tsoil1(i1,j1,0)
!     wetsfc1(i1,j1,0) = wetsfc1(i1,j1,0)
!     wetdp1(i1,j1,0) = wetdp1(i1,j1,0)
!     wetcanp1(i1,j1,0) = wetcanp1(i1,j1,0)
    END DO
  END DO

END SUBROUTINE intrp_soil

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc. 
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE LD2DGRID                   ######
!######                                                      ######
!######                    Developed by                      ######
!######            Weather Decision Technologies             ######
!######                                                      ######
!##################################################################
!##################################################################


SUBROUTINE ld2dgrid(dxin,dyin,ctrlatin,ctrlonin,  &
           mprojin,trlat1in,trlat2in,trlonin,sclfctin)

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Load 2-D grid parameters into the ARPS common blocks (useful for
!  C programs using ARPS subroutines).
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2001/08/03
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  OUTPUT:
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE

  REAL :: dxin,dyin,ctrlatin,ctrlonin,trlat1in,trlat2in,trlonin,sclfctin
  INTEGER :: mprojin

!-----------------------------------------------------------------------
!
!  Misc local variables
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  dx = dxin
  dy = dyin
  ctrlat = ctrlatin
  ctrlon = ctrlonin
  mapproj = mprojin
  trulat1 = trlat1in
!wdt update
  trulat2 = trlat2in
  trulon = trlonin
  sclfct = sclfctin

  RETURN
END SUBROUTINE ld2dgrid

!wdt testing:     

      SUBROUTINE test_dump(data,filename,nx,ny,nz,stagdim,incfake),1

      ! write out a unformatted dump to debug MP version

      implicit none

      include 'mp.inc'

      integer nx,ny,nz
      integer stagdim,incfake
      real data(nx,ny,nz)
      character*(*) filename
      character*80 tempname
      integer lenfl

      integer count
      common /save_count/ count

      !return
      count = count + 1

      write(tempname, '(a,a,i4.4)') trim(filename),'_c',count

      if (mp_opt > 0) then
        write(tempname, '(a,a,2i2.2)')  &
           trim(adjustl(tempname)),'_',loc_x,loc_y
      endif

      open (11,file=tempname,form='unformatted')
      write (11) nx,ny,nz,stagdim,incfake
      write (11) data
      close (11)

      CALL test_check_in(tempname)

      RETURN
      END


      SUBROUTINE test_dump2(data,filename,nx,ny,nz,ib,ie,jb,je,kb,ke)

      ! write out a unformatted dump to debug MP version

      implicit none

      include 'mp.inc'

      integer nx,ny,nz
      integer ib,ie,jb,je,kb,ke   ! compare for i=ib,ie ...
      real data(nx,ny,nz)
      character*(*) filename
      character*80 tempname
      integer lenfl

      integer count
      common /save_count/ count

      !return
      count = count + 1

      write(tempname, '(a,a,i4.4)') trim(filename),'_c',count

      if (mp_opt > 0) then
        write(tempname, '(a,a,2i2.2)')  &
           trim(adjustl(tempname)),'_',loc_x,loc_y
      endif

      open (11,file=tempname,form='unformatted')
      write (11) nx,ny,nz,-1,-1
      write (11) ib,ie,jb,je,kb,ke
      write (11) data
      close (11)

      RETURN
      END


      subroutine test_check_in(comment) 1,4
      implicit none
      include 'mp.inc'
      character*(*) comment

      !return
      CALL flush(6)
      CALL flush(0)
      write (*,*) "XXX checking in:",myproc," ",trim(comment)
      CALL flush(6)
      CALL mpbarrier

      return
      end