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


SUBROUTINE get_ged(dtype,filenm,nx,ny,xdims,ydims,resl,glon,glat,  & 4,5
                   iout,rout)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in the surface data set and determine the region just
!  covering the model domain.
!
!  For soil data, the resolution is   1 x 1   degrees of lat x lon
!  For NDVI data, the resolution is 1/6 x 1/6 degrees of lat x lon
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!
!  1/28/94
!
!  MODIFICATIONS:
!
!  2001/07/26 Gene Bassett
!  Moved calls to getstyp1, getvtyp1, getndvi1 inside this subroutine.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  dtype    Data flag: 1 for soil type data
!                      2 for vegetation type data
!  filenm   File name of usrface data set
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!  xdims    X-dimsion size of the data set
!  ydims    Y-dimsion size of the data set
!
!  glon     Longitude values of model grid points
!  glat     Latitude  values of model grid points
!
!  OUTPUT:
!
!  dtlat    Latitude  values of data grid points
!  dtlon    Longitude values of data grid points
!
!  dtarr    Integer data array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: dtype        ! Data flag: 1 for soil; 2 for ndvi
  CHARACTER (LEN=*  ) :: filenm ! File name of usrface data set
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
  INTEGER :: xdims        ! Columns of/ the data set
  INTEGER :: ydims        ! Rows of the data set

  REAL :: resl         ! Resolution of data
  INTEGER :: dtproj       ! Projection of data

  REAL :: glon(nx,ny)  ! Longitude values of model grid points
  REAL :: glat(nx,ny)  ! Latitude  values of model grid points

  INTEGER :: iout(1)   ! Points to vegtyp with dtype=2
  REAL :: rout(1)     ! Points to soiltyp when dtype=1, ndvi when dtype=3
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  REAL, ALLOCATABLE :: dtlat(:) ! Latitude  values of data grid points
  REAL, ALLOCATABLE :: dtlon(:) ! Longitude values of data grid points
  INTEGER, ALLOCATABLE :: dtarr(:,:)  ! Data array

  INTEGER :: lenfl        ! Length of filenm
  INTEGER :: flunit       ! IO unit of filenm
  CHARACTER (LEN=20) :: title     ! Title of data
  CHARACTER (LEN=20) :: temtit    ! Temporary string of data title
  INTEGER :: colmn        ! Columns of data
  INTEGER :: row          ! Rows of data
  REAL :: resl0        ! Resolution of data
!
  INTEGER :: i,j
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ALLOCATE (dtlat(ydims))
  ALLOCATE (dtlon(xdims))
  ALLOCATE (dtarr(xdims,ydims))

  lenfl = LEN( filenm )
  CALL strlnth( filenm, lenfl )
  CALL getunit( flunit )

  WRITE(6,'(a,a)') '    Opening ',filenm(1:lenfl)
  OPEN (UNIT = flunit, FILE = filenm(1:lenfl), STATUS = 'old')
!    :      form = 'unformatted', access = 'sequential')

  READ (flunit,'(a20,i8,i8,e20.10,i8)')                                 &
       title, colmn, row, resl, dtproj

  IF ( dtype == 1 ) THEN
    temtit = 'W & H-S Soil Type   '
    resl0 = 1.0
  ELSE IF ( dtype == 2 )THEN
    temtit = 'World Ecosystems    '
    resl0 = 1.0/6.0
  ELSE IF ( dtype == 3 )THEN
    temtit = 'NDVI data           '
    resl0 = 1.0/6.0
  ELSE
    WRITE (6, '(a/a)')                                                  &
        'Data type is not correct.',                                    &
        'Program stops here in subroutine GET_GED.'
    STOP
  END IF

  IF ( title /= temtit ) THEN
    WRITE (6, '(a/a/a/a/a)')                                            &
        'Data title is not correct.',                                   &
        'The file opened perhaps was not what you wanted.',             &
        'title required = '//temtit,                                    &
        'title in the file = '//title,                                  &
        'Program stops here in subroutine GET_GED.'
    STOP
  ELSE IF ( colmn /= xdims .OR. row /= ydims ) THEN
    WRITE (6, '(a/a/a,i5,a,i5/a,i5,a,i5/a)')                            &
        'Data dimension size is not correct.',                          &
        'The file opened perhaps was not what you wanted.',             &
        'column required = ',xdims, ', column in the file = ', colmn,   &
        'row required = ',ydims, ', row in the file = ', row,           &
        'Program stops here in subroutine GET_GED.'
    STOP
  ELSE IF ( ABS( ( resl - resl0 ) / resl ) > 0.01 ) THEN
    WRITE (6, '(a/a/a,i5/a,i5/a)')                                      &
        'Data resolution is not correct.',                              &
        'The file opened perhaps was not what you wanted.',             &
        'resolution required = ',resl0,                                 &
        'resolution in the file = ', resl,                              &
        'Program stops here in subroutine GET_GED.'
    STOP
  ELSE IF ( dtproj /= 1 ) THEN
    WRITE (6, '(a/a/a,i5/a,i5/a)')                                      &
        'Data projection is not correct.',                              &
        'The file opened perhaps was not what you wanted.',             &
        'projection required = ', 1,                                    &
        'projection in the file = ', dtproj,                            &
        'Program stops here in subroutine GET_GED.'
    STOP
  END IF

  WRITE(6,'(a)') '    Reading data set for GED '//title
  READ (flunit,'(20i4)') dtarr

  DO i = 1, xdims
    dtlon(i) = resl*(FLOAT(i)-.5)
  END DO

  DO j = 1, ydims
    dtlat(j) = resl*(FLOAT(j)-.5) - 90.0
  END DO

  CLOSE ( flunit )
  CALL retunit ( flunit )

  IF (dtype == 1) THEN

!-----------------------------------------------------------------------
!  Translate the soil classes data (stypdat) into the USDA soil type
!  category and fill them into the model grid array soiltyp (in iout).
!-----------------------------------------------------------------------

    CALL getstyp2(xdims,ydims,nx,ny,dtlon,dtlat,resl,glon,glat,dtarr,iout) 

  ELSE IF (dtype == 2) THEN

!-----------------------------------------------------------------------
!  Transfer the World Ecosystem Classes data into the 12 vegetation
!  type categories and fill into model grid array, vegtyp(nx,ny) (in iout)
!-----------------------------------------------------------------------

    CALL getvtyp2(xdims,ydims,nx,ny,dtlon,dtlat,resl,glon,glat,dtarr,iout)

  ELSE ! dtype = 3

!-----------------------------------------------------------------------
!  Calculate the real-NDVI data from byte-NDVI and fill into model
!  grid array, ndvi(nx,ny) (in rout)
!-----------------------------------------------------------------------

    CALL getndvi2(xdims,ydims,nx,ny,dtlon,dtlat,resl,glon,glat,dtarr,rout)

  ENDIF

!wdt-williams update
  DEALLOCATE (dtlat)
  DEALLOCATE (dtlon)
  DEALLOCATE (dtarr)

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

SUBROUTINE get_statsgo(flag,filenm1,filenm2,nx,ny,nstyps,x1,y1,        & 3,8
           col,row,glat,glon,resl,  &
           iout,rout)                
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in the surface property data set & determine the region just
!  covering the model domain.
!
!  the resolution is   1 * 1 km  of grid
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!
!  Leilei Wang ,Vince Wong
!  3/18/97
!
!  MODIFICATIONS:
!
!  2001/07/26 Gene Bassett
!  Moved calls to getstyp1, getvtyp1, getndvi1 inside this subroutine.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  filenm   File name of usrface data set
!
!  OUTPUT:
!
!  iout     Integer data array (soiltyp when flag=1, vegtyp when flag=2)
!  rout     Real data array (stypfrct when flag=1, ndvi when flag=3)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: flag         ! Flag for what kind of data
  CHARACTER (LEN=* ) :: filenm1 ! File name of usrface data set
  CHARACTER (LEN=* ) :: filenm2 ! File name of GED     data set
  INTEGER :: nx
  INTEGER :: ny
  INTEGER :: nstyps

  INTEGER :: dcol,drow
  INTEGER :: col,row      ! Columns of the data set

  INTEGER :: xdims        ! Actual colmns in data array
  INTEGER :: ydims        ! Actual rows in data array
!
  REAL :: resl              ! Resolution of data
  REAL :: resl_ged          ! Resolution of GED soil data

  REAL :: dctrlat           ! LAT. of projection center of data
  REAL :: dctrlon           ! LON. of projection center of data

  REAL :: x1(nx,ny)         ! X_coordinate of arps model grid
  REAL :: y1(nx,ny)         ! Y_coordinate of arps model grid
  REAL :: glat(nx,ny)       ! latitude   of arps model grid
  REAL :: glon(nx,ny)       ! longitude   of arps model grid
  REAL, ALLOCATABLE :: xdat(:)  ! X_coordinate of data grid covering model
                                ! domain
  REAL, ALLOCATABLE :: ydat(:)  ! Y_coordinate of data grid covering model
                                ! domain

  INTEGER, ALLOCATABLE :: idtarr (:,:)       ! Output data array

  INTEGER :: iout(1)   ! Points to soiltyp when flag=1, vegtyp when flag=2
  REAL :: rout(1)      ! Points to stypfrct when flag=1, ndvi when flag=3
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: lenfl
  INTEGER :: flunit
  INTEGER :: xmin,ymin
  INTEGER :: xmax,ymax
  REAL :: dnwx,dnwy,dsex,dsey
  REAL :: xxmin,yymax
  REAL :: xorgmin,xorgmax
  REAL :: yorgmin,yorgmax
  INTEGER :: i,j,m,n
  
  INTEGER, ALLOCATABLE :: item1(:,:) ! add by wyh
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
  INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  
  lenfl = 70
  CALL strlnth( filenm1, lenfl )
  CALL getunit( flunit )

  resl = 1000.
  dctrlon = -100.00

  IF (flag == 1) THEN
    dcol = 4587
    drow = 2889
    dctrlat =    45.0
    dnwx =  - 2050500
    dnwy =    752500
    dsex =    2536500
    dsey =  - 2136500
  ELSE
    dcol = 9223
    drow = 8996
    dctrlat =    50.0
    dnwx =  - 4486550
    dnwy =    4479868
    dsex =    4735450
    dsey =  - 4515132
  END IF

  DO i = 1,nx
    DO j = 1,ny
      CALL lae_lltoxy(1,1,dctrlat,dctrlon,glat(i,j),glon(i,j),          &
                      x1(i,j),y1(i,j))
    END DO
  END DO
!
  xorgmin = x1(1, 1)
  xorgmax = x1(nx,1)
  yorgmin = y1(1, 1)
  yorgmax = y1(1,ny)

  DO i = 1, nx
    DO j = 1, ny
      IF(x1(i,j) < xorgmin) THEN
        xorgmin = x1(i,j)
      END IF
      IF(x1(i,j) > xorgmax) THEN
        xorgmax = x1(i,j)
      END IF
      IF(y1(i,j) < yorgmin) THEN
        yorgmin = y1(i,j)
      END IF
      IF(y1(i,j) > yorgmax) THEN
        yorgmax = y1(i,j)
      END IF
    END DO
  END DO
!
  IF(dx < resl) THEN
    xmin = INT( (xorgmin-dnwx-2*resl)/resl )
    xmax = INT( (xorgmax-dnwx+2*resl)/resl )+1
  ELSE
    xmin = INT( (xorgmin-dnwx-2*dx)/resl )
    xmax = INT( (xorgmax-dnwx+2*dx)/resl )+1
  END IF

  IF(dy < resl) THEN
    ymin = INT( (dnwy-yorgmax-2*resl)/resl )
    ymax = INT( (dnwy-yorgmin+2*resl)/resl )+1
  ELSE
    ymin = INT( (dnwy-yorgmax-2*dy)/resl )
    ymax = INT( (dnwy-yorgmin+2*dy)/resl )+1
  END IF

  CALL get_ged(flag,filenm2,nx,ny,col,row,resl_ged,glon,glat,iout,rout)

  IF ( xmin < 1   ) xmin = 1
  IF ( xmax > dcol) xmax = dcol
  IF ( ymin < 1   ) ymin = 1
  IF ( ymax > drow) ymax = drow

  xdims = xmax - xmin + 1
  ydims = ymax - ymin + 1

  ALLOCATE(xdat(xdims))
  ALLOCATE(ydat(ydims))
  ALLOCATE(idtarr(xdims,ydims))
  ALLOCATE(item1(xdims, ydims))

  xxmin = dnwx + (xmin-1) * resl
  yymax = dnwy - (ymax-1) * resl

  item1 = 0

  CALL readbsq (filenm1, flunit, 1, drow,dcol, 1,                       &
      1, ymin, xmin, 1, ydims , xdims,item1,xdims,ydims)

  DO j=1,ydims
    DO i=1,xdims
      idtarr(i,j) = item1(i,ydims-j+1)
    END DO
  END DO

  DEALLOCATE(item1)

  DO i = 1 ,xdims
    xdat( i ) = xxmin + (i-1)*resl
  END DO

  DO j = 1 ,ydims
    ydat( j ) = yymax + (j-1)*resl
  END DO

  CLOSE ( flunit )
  CALL retunit ( flunit )

  IF (flag == 1) THEN

!-----------------------------------------------------------------------
!  Translate the soil classes data storged in stypdat into the USDA soil 
!  type categories. nstyp soil types with the highest percentages are 
!  then brought to the model grid points, with each type being associated
!  with a value of soil type fraction. (iout is soiltyp, rout is stypfrct)
!-----------------------------------------------------------------------

    CALL getstyp1(xdims,ydims, nx,ny,nstyps,glon,glat,              &
                  x1,y1,idtarr,xdat,ydat,iout,rout)

  ELSE IF (flag == 2) THEN

!-----------------------------------------------------------------------
!  Calculate vegetation type array using statistical method to get
!  vegtyp(nx,ny) (in iout)
!-----------------------------------------------------------------------

    CALL getvtyp1(xdims,ydims, nx,ny,glon,glat,                     &
                  x1,y1,idtarr,xdat,ydat,iout)

  ELSE  ! flag = 3 
  
!-----------------------------------------------------------------------
!  Calculate real NDVI type data using statistical method to get
!  ndvi(nx,ny) (in rout)
!-----------------------------------------------------------------------

    CALL getndvi1(xdims,ydims,nx,ny,                              &
                  x1,y1,idtarr,xdat,ydat,rout)

  ENDIF

!wdt-williams update
  DEALLOCATE(xdat)
  DEALLOCATE(ydat)
  DEALLOCATE(idtarr)

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

SUBROUTINE getstyp1(xdims,ydims,nx,ny,nstyps,glon,glat,                 & 1
           x1,y1,stypdat,xdat,ydat,soiltyp,stypfrct)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Transfer the soil classes data (soiltyp) into the USDA soil type
!  category and fill them into the model grid.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!  Leilei Wang ,Vince Wong
!  3/20/97
!
!-----------------------------------------------------------------------
!  MODIFICATIONS:
!
!  Yunheng Wang
!  8/04/01
!  Rewrite some code in Fortran 90.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!  xdims    X-dimsion size of the data set
!  ydims    Y-dimsion size of the data set
!
!  glon     Longitude values of model grid points
!  glat     Latitude  values of model grid points
!
!  dtlon    Longitude values of data grid points
!  dtlat    Latitude  values of data grid points
!
!  stypdat  Soil data array
!  sfrctdat Soil type fraction data array
!
!  OUTPUT:
!
!  soiltyp  Soil type array
!  stypfrct Soil type fraction array
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
  INTEGER :: xdims        ! X-dimension size of the data set
  INTEGER :: ydims        ! Y-dimension size of the data set
  REAL :: x1  (nx,ny)  ! X_coordinate of model grid
  REAL :: y1  (nx,ny)  ! X_coordinate of model grid

  INTEGER :: nstyps       ! Number of soil types in each grid point
!
  REAL :: xdat(xdims)  ! X_coor. of soil data
  REAL :: ydat(ydims)  ! Y_coor. of soil data
  REAL :: glon(nx,ny)  ! Lon. of  model grid
  REAL :: glat(nx,ny)  ! Lat. of  model grid

  INTEGER :: stypdat (xdims,ydims)        ! Soil type data array
  INTEGER :: soiltyp (nx,ny,nstyps)       ! Soil type array
  REAL :: stypfrct(nx,ny,nstyps)       ! Soil type fraction array
  REAL :: dresl                        ! data resolution
  INTEGER :: item1  (nsoiltyp)            ! Temporary array
  INTEGER :: item2  (nsoiltyp)            ! Temporary array
  REAL :: tem3  (nsoiltyp)             ! Temporary array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,im, iis,is,ii,jj
  INTEGER :: imax
  INTEGER :: count
  INTEGER :: xleft,xright,yleft,yright
  REAL :: tot,temp1
  INTEGER :: itemp2
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
!  Fill the soil types into the model grid using statistics results
!  by chosing soil types of highest weights and their corresponding
!  fractions for each model grid point
!
!-----------------------------------------------------------------------
!
  dresl =1000.

  DO j = 1,ydims
    DO i = 1,xdims

      SELECT CASE (stypdat(i,j))
        CASE (5:12)
          stypdat(i,j) = stypdat(i,j) -1
        CASE (13)
          stypdat(i,j) = 4 
        CASE (14)
          stypdat(i,j) = 13
        CASE (15, 16)
          stypdat(i,j) = 11
      END SELECT

    END DO
  END DO

  DO i = 1,nx
    DO j = 1,ny

      item1(:) = 0
      item2(:) = 0
      tem3(:)  = 0.0

      count = 0
      tot   = 0.0

      IF(dx < dresl) THEN
        xleft = INT( (x1(i,j) - xdat(1)-dresl/2)/dresl )+1
        xright= INT( (x1(i,j) - xdat(1)+dresl/2)/dresl )
      ELSE
        xleft = INT( (x1(i,j) - xdat(1)-dx/2)/dresl )+1
        xright= INT( (x1(i,j) - xdat(1)+dx/2)/dresl )
      END IF

      IF(dy < 1000.0) THEN
        yleft = INT( (y1(i,j) - ydat(1)-dresl/2)/dresl )+1
        yright= INT( (y1(i,j) - ydat(1)+dresl/2)/dresl )
      ELSE
        yleft = INT( (y1(i,j) - ydat(1)-dy/2)/dresl )+1
        yright= INT( (y1(i,j) - ydat(1)+dy/2)/dresl )
      END IF

      DO ii = xleft, xright
        DO jj = yleft, yright
          IF(xleft < 1.OR.yleft < 1) GO TO 800
!wdt update
          IF(xright > xdims.OR.yright > ydims) GO TO 800
          IF(stypdat(ii,jj) /= 0) count=count + 1
          DO is=1,nsoiltyp
            IF(stypdat(ii,jj) == is) THEN
              item1(is) = item1(is) + 1
              item2(is) = is
            END IF
          END DO
800     CONTINUE
        END DO
      END DO

      IF(count == 0) GOTO 100 

      tem3(:) = FLOAT(item1(:))/count

      DO iis=1,nsoiltyp
        imax= iis
        DO im=iis, nsoiltyp
          IF (tem3(im) > tem3(imax)) imax = im
        END DO 
        temp1=tem3(iis)
        tem3(iis)=tem3(imax)
        tem3(imax)=temp1
        itemp2=item2(iis)
        item2(iis)=item2(imax)
        item2(imax)=itemp2
      END DO

      tot = SUM(tem3(1:nstyp))

100   CONTINUE

      DO is=1,nstyp

        IF(tot /= 0.0)THEN
          soiltyp(i,j,is) = item2(is)
          stypfrct(i,j,is)= tem3(is)/tot
        ELSE
          stypfrct(i,j,is) = 0.0
          stypfrct(i,j,1) = 1.0
        END IF

      END DO

    END DO
  END DO

  DO is = 2, nstyp
     WHERE(soiltyp(:,:,is) == 0)
       soiltyp(:,:,is) = soiltyp (:,:,is-1)
       stypfrct(:,:,is) = 0.0
     END WHERE
	  
     WHERE (soiltyp(:,:,is) == soiltyp(:,:,is-1))
       stypfrct(:,:,is) = 0.0
     END WHERE
  END DO

!  open(234, file='soiltyp')
!  write(234,'(10I4)') soiltyp

!
!-----------------------------------------------------------------------
!
!  The following DO loop is just for the Vortex experiment. On the
!  south boundary Vortex domain the vegetation type changes sharply,
!  that will pretty much influence the boundary conditions. Therefore
!  the following DO loop will use the values of inner grids to those
!  boundary grids.
!
!-----------------------------------------------------------------------
!
!  DO 300 j = 1, 7
!  DO 300 i = 1, nx
!    soiltyp(i,j) = soiltyp(i,8)
!300  CONTINUE

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


SUBROUTINE getstyp2(xdims,ydims, nx,ny,                                 & 1
           dtlon,dtlat,resl, glon,glat,                                 &
           stypdat, soiltyp)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Translate the soil classes data (stypdat) into the USDA soil type
!  category and fill them into the model grid array soiltyp.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!
!  1/29/94
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!  xdims    X-dimsion size of the data set
!  ydims    Y-dimsion size of the data set
!
!  glon     Longitude values of model grid points
!  glat     Latitude  values of model grid points
!
!  dtlon    Longitude values of data grid points
!  dtlat    Latitude  values of data grid points
!
!  stypdat  Soil data array
!
!  OUTPUT:
!
!  soiltyp  Soil type array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
  INTEGER :: xdims        ! X-dimsion size of the data set
  INTEGER :: ydims        ! Y-dimsion size of the data set
!
  REAL :: glon(nx,ny)  ! Longitude values of model grid points
  REAL :: glat(nx,ny)  ! Latitude  values of model grid points
!
  REAL :: dtlon(xdims) ! Longitude values of data grid points
  REAL :: dtlat(ydims) ! Latitude  values of data grid points
  REAL :: resl         ! Data resolution
!
  INTEGER :: stypdat(xdims,ydims)  ! Soil data array
!
  INTEGER :: soiltyp(nx,ny)        ! Soil type array

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j, ii, jj, ix, jy
  INTEGER :: ixmin,ixmax,jymin,jymax
  REAL :: r,r1
  
  LOGICAL :: find
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  
  DO j = 1, ydims
    DO i = 1, xdims

      SELECT CASE (stypdat(i,j))

        CASE ( 11, 17, 23 ) 

          stypdat(i,j) = 2

        CASE ( 14, 20, 26, 27 )

          stypdat(i,j) = 3

        CASE ( 12, 18, 24 )

          stypdat(i,j) = 5

        CASE ( 15, 21, 28 )

          stypdat(i,j) = 6

        CASE ( 13 )

          stypdat(i,j) = 8

        CASE ( 19, 25 )

          stypdat(i,j) = 9

        CASE ( 16, 22 )

          stypdat(i,j) = 10

        CASE ( 29, 30, 31 )

          stypdat(i,j) = 11

        CASE ( 34 )

          stypdat(i,j) = 12

        CASE ( 00 )

          stypdat(i,j) = 13

        CASE DEFAULT

          WRITE (6, '(a/a,i4,a,i4,a,i2/a)')                           &
            'The value for soil type is invalid.',                    &
            'stypdat(', i, ',', j, ') = ', stypdat(i,j),              &
            'Program stops in subroutine GETSTYP2.'

          STOP

        END SELECT

    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Fill the soil type into the model grid by chosing the values
!  at the nearest data grid points.
!
!-----------------------------------------------------------------------
!
  ixmin = MIN( xdims-1,                                                 &
               INT( ( glon(1,1)-dtlon(1) ) / resl ) + 1,                &
               INT( ( glon(1,ny)-dtlon(1) ) / resl ) + 1 )

  ixmin = MAX( 1, ixmin - 5 )

  jymin = MIN( ydims-1,                                                 &
               INT( ( glat(1,1)-dtlat(1) ) / resl ) + 1,                &
               INT( ( glat(nx,1)-dtlat(1) ) / resl ) + 1 )

  jymin = MAX( 1, jymin - 5 )

  ixmax = MIN( xdims,                                                   &
               nint( ( glon(nx,1)-dtlon(1) ) / resl ) + 1,              &
               nint( ( glon(nx,ny)-dtlon(1) ) / resl ) + 1 )

  ixmax = MIN( xdims, ixmax + 5 )

  jymax = MIN( ydims,                                                   &
               nint( ( glat(1,ny)-dtlat(1) ) / resl ) + 1,              &
               nint( ( glat(nx,ny)-dtlat(1) ) / resl ) + 1 )

  jymax = MIN( ydims, jymax + 5 )

  DO j = 1, ny
    DO i = 1, nx

      r1 = 1.0E20
      find = .FALSE.
      ii = 0
      jj = 0

      DO jy = jymin, jymax
        DO ix = ixmin, ixmax
          r = ( dtlon(ix) - glon(i,j) ) ** 2                         &
              + ( dtlat(jy) - glat(i,j) ) ** 2
          IF(.NOT. find .AND. r < r1) find = .TRUE.
          IF(.NOT. find) write(6,*) 'Warnning in GETSTYP2 ......'
          IF ( r < r1 ) THEN
            r1 = r
            ii = ix
            jj = jy
          END IF
        END DO
      END DO

      IF(find) THEN
        soiltyp(i,j) = stypdat(ii,jj)
      ELSE
        WRITE(6,*) 'GETSTYP2: Error in getstyp2 for some machine',    &
                   ' Please try again: makearps -opt 1 arpssfc'
        STOP
      END IF

    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  The following DO loop is just for the Vortex experiment. On the
!  south boundary Vortex domain the vegetation type changes sharply,
!  that will pretty much influence the boundary conditions. Therefore
!  the following DO loop will use the values of inner grids to those
!  boundary grids.
!
!-----------------------------------------------------------------------
!
!  DO 300 j = 1, 7
!  DO 300 i = 1, nx
!    soiltyp(i,j) = soiltyp(i,8)
!300  CONTINUE

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


SUBROUTINE getvtyp1(xdims,ydims,nx,ny,glon,glat,                        & 1
           x1,y1,vtypdat,xdat,ydat,vegtyp)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Converse the OWE classes of vegetation type (vegtyp)into ARPS
!  type category and fill them into the model grid.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!  Leilei Wang ,Vince Wong
!  8/11/97
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!  xdims    X-dimsion size of the data set
!  ydims    Y-dimsion size of the data set
!
!  glon     Longitude values of model grid points
!  glat     Latitude  values of model grid points
!
!  dtlon    Longitude values of data grid points
!  dtlat    Latitude  values of data grid points
!
!  nvegtyp  Number of soil types in STATESGO class
!  vtypdat  Vegetation data array
!
!  OUTPUT:
!
!  vegtyp   Vegetation type array
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
  INTEGER :: xdims        ! X-dimension size of the data set
  INTEGER :: ydims        ! Y-dimension size of the data set
  REAL :: x1  (nx,ny)  ! X_coordinate of model grid
  REAL :: y1  (nx,ny)  ! X_coordinate of model grid

  REAL :: xdat(xdims)  ! X_coor. of soil data
  REAL :: ydat(ydims)  ! Y_coor. of soil data
  REAL :: glon(nx,ny)  ! Lon. of  model grid
  REAL :: glat(nx,ny)  ! Lat. of  model grid

  REAL :: dresl                        ! data resolution
  INTEGER :: vtypdat (xdims,ydims)        ! Veg. type data array
  INTEGER :: vegtyp  (nx,ny)              ! Veg. type array

  INTEGER :: item2(nvegtyp)               ! Temporary array
  INTEGER :: item3(nvegtyp)               ! Temporary array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j, iis,is,ii,jj
  INTEGER :: imax,count
  INTEGER :: xleft,xright,yleft,yright
  INTEGER, ALLOCATABLE :: item1(:,:)     ! add by wyh
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-----------------------------------------------------------------------
!
!  Fill the veg. types into the model grid using statistics results
!  by chosing veg. types of highest weights and their corresponding
!  fractions for each model grid point
!
!-----------------------------------------------------------------------
!
  ALLOCATE(item1(xdims, ydims))

  dresl =1000.

  DO j = 1,ydims
    DO i = 1,xdims

      SELECT CASE (vtypdat(i,j))

        CASE ( 01, 08, 50, 69, 71 )

          item1(i,j) = 1

        CASE ( 42, 09, 53 )

          item1(i,j) = 2

        CASE ( 16, 30, 37, 40, 10, 76, 93, 52 )

          item1(i,j) = 3

        CASE ( 02, 41, 43, 49, 07, 82, 83, 85, 94, 87 )

          item1(i,j) = 4

        CASE ( 18, 91, 58 )

          item1(i,j) = 5

        CASE ( 24:27, 29, 56, 61, 03:05, 19, 78, 88, 89 )

          item1(i,j) = 6

        CASE ( 06, 20:23, 46:48, 57, 60, 62, 34, 77 )

          item1(i,j) = 7

        CASE ( 32:33, 54, 79, 86, 90 )

          item1(i,j) = 8

        CASE ( 17, 70 )

          item1(i,j) = 9

        CASE ( 28, 31, 36, 38, 39, 55, 35, 92 )

          item1(i,j) = 10

        CASE ( 44, 45, 13, 74, 75, 80 )

          item1(i,j) = 11

        CASE ( 59, 63, 12, 64 )

          item1(i,j) = 12

        CASE ( 51, 81, 84, 11 )

          item1(i,j) = 13

        CASE ( 00, 65:68, 72, 73, 14, 15 )

          item1(i,j) = 14

        CASE DEFAULT

          WRITE (6, '(a/a,i4,a,i4,a,i2/a)')                               &
            'The value for vegetation type is invalid.',                &
            'vtypdat(', i, ',', j, ') = ', vtypdat(i,j),                &
            'Program stops in subroutine GETVTYP1.'

          STOP

      END SELECT

    END DO
  END DO

  DO j = 1,ny
    DO i = 1,nx

      item2=0
      item3=0

      count = 0
      imax = 0

      IF(dx < dresl) THEN
        xleft = INT( (x1(i,j) - xdat(1)-dresl/2)/dresl )+1
        xright= INT( (x1(i,j) - xdat(1)+dresl/2)/dresl )
      ELSE
        xleft = INT( (x1(i,j) - xdat(1)-dx/2)/dresl )+1
        xright= INT( (x1(i,j) - xdat(1)+dx/2)/dresl )
      END IF

      IF(dy < 1000.0) THEN
        yleft = INT( (y1(i,j) - ydat(1)-dresl/2)/dresl )+1
        yright= INT( (y1(i,j) - ydat(1)+dresl/2)/dresl )
      ELSE
        yleft = INT( (y1(i,j) - ydat(1)-dy/2)/dresl )+1
        yright= INT( (y1(i,j) - ydat(1)+dy/2)/dresl )
      END IF

      DO ii = xleft, xright
        DO jj = yleft, yright

          IF(xleft < 1.OR.yleft < 1) GO TO 800
!wdt update
          IF(xright > xdims.OR.yright > ydims) GO TO 800

          IF(item1(ii,jj) /= 0) count=count + 1
          DO is=1,nvegtyp
            IF(item1(ii,jj) == is) THEN
              item2  (is) = item2(is) + 1
              item3  (is) = is
            END IF
          END DO
800     CONTINUE
        END DO
      END DO

      IF(count /= 0)  THEN
        imax = 1
        DO iis = 1, nvegtyp
          IF(item2(iis) >= item2(imax))  imax=iis
        END DO
      END IF

      IF (imax /= 0) vegtyp(i,j) = item3(imax)

    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  The following DO loop is just for the Vortex experiment. On the
!  south boundary Vortex domain the vegetation type changes sharply,
!  that will pretty much influence the boundary conditions. Therefore
!  the following DO loop will use the values of inner grids to those
!  boundary grids.
!
!-----------------------------------------------------------------------
!
!  DO 300 j = 1, 7
!  DO 300 i = 1, nx
!    vegtyp(i,j) = vegtyp(i,8)
!300  CONTINUE

  DEALLOCATE(item1)

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


SUBROUTINE getvtyp2(xdims,ydims, nx,ny,                                 & 1
           dtlon,dtlat,resl, glon,glat,                                 &
           vtypdat, vegtyp)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Transit the Olson World Ecosystem Classes of vegetation data into
!  the simple vegetation type categories and fill them into the model
!  domain by chosing the values at the nearest data grid points.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!
!  2/20/94
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!  xdims    X-dimsion size of the data set
!  ydims    Y-dimsion size of the data set
!
!  glon     Longitude values of model grid points
!  glat     Latitude  values of model grid points
!
!  dtlon    Longitude values of data grid points
!  dtlat    Latitude  values of data grid points
!
!  vtyp     Vegetation data array
!
!  OUTPUT:
!
!  vegtyp   Vegetation type array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
  INTEGER :: xdims        ! X-dimsion size of the data set
  INTEGER :: ydims        ! Y-dimsion size of the data set
!
  REAL :: glon(nx,ny)  ! Longitude values of model grid points
  REAL :: glat(nx,ny)  ! Latitude  values of model grid points
!
  REAL :: dtlon(xdims) ! Longitude values of data grid points
  REAL :: dtlat(ydims) ! Latitude  values of data grid points
  REAL :: resl         ! Data resolution
!
  INTEGER :: vtypdat(xdims,ydims)  ! Soil data array
!
  INTEGER :: vegtyp(nx,ny)      ! Soil type array

!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j, ii, jj, ix, jy
  INTEGER :: ixmin,ixmax,jymin,jymax
  REAL :: r,r1

!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  DO j = 1, ydims
    DO i = 1, xdims

      IF ( vtypdat(i,j) == 01 .OR.  vtypdat(i,j) == 08 .OR.             &
             vtypdat(i,j) == 50 .OR.                                    &
             vtypdat(i,j) == 69 .OR.                                    &
             vtypdat(i,j) == 71 ) THEN

        vtypdat(i,j) = 1

      ELSE IF ( vtypdat(i,j) == 42 .OR.                                 &
                 vtypdat(i,j) == 53 ) THEN

        vtypdat(i,j) = 2

      ELSE IF ( vtypdat(i,j) == 16 .OR.                                 &
                 vtypdat(i,j) == 30 .OR.                                &
                 vtypdat(i,j) == 37 .OR.                                &
                 vtypdat(i,j) == 40 .OR.                                &
                 vtypdat(i,j) == 52 ) THEN

        vtypdat(i,j) = 3

      ELSE IF ( vtypdat(i,j) == 02 .OR.                                 &
                 vtypdat(i,j) == 41 .OR.                                &
                 vtypdat(i,j) == 43 .OR.                                &
                 vtypdat(i,j) == 49 ) THEN

        vtypdat(i,j) = 4

      ELSE IF ( vtypdat(i,j) == 58 .OR.                                 &
                 vtypdat(i,j) == 47 ) THEN

        vtypdat(i,j) = 5

      ELSE IF ( vtypdat(i,j) == 24 .OR.                                 &
                 vtypdat(i,j) == 25 .OR.                                &
                 vtypdat(i,j) == 26 .OR.                                &
                 vtypdat(i,j) == 27 .OR.                                &
                 vtypdat(i,j) == 29 .OR.                                &
                 vtypdat(i,j) == 56 .OR.                                &
                 vtypdat(i,j) == 61 ) THEN

        vtypdat(i,j) = 6

      ELSE IF ( vtypdat(i,j) == 06 .OR.                                 &
                 vtypdat(i,j) == 20 .OR.                                &
                 vtypdat(i,j) == 21 .OR.                                &
                 vtypdat(i,j) == 22 .OR.                                &
                 vtypdat(i,j) == 23 .OR.                                &
                 vtypdat(i,j) == 46 .OR.                                &
                 vtypdat(i,j) == 48 .OR.                                &
                 vtypdat(i,j) == 57 .OR.                                &
                 vtypdat(i,j) == 60 .OR.                                &
                 vtypdat(i,j) == 62 ) THEN

        vtypdat(i,j) = 7

      ELSE IF ( vtypdat(i,j) == 32 .OR.                                 &
                 vtypdat(i,j) == 33 .OR.                                &
                 vtypdat(i,j) == 54 ) THEN

        vtypdat(i,j) = 8

      ELSE IF ( vtypdat(i,j) == 17 .OR.                                 &
                 vtypdat(i,j) == 70 ) THEN

        vtypdat(i,j) = 9

      ELSE IF ( vtypdat(i,j) == 28 .OR.                                 &
                 vtypdat(i,j) == 31 .OR.                                &
                 vtypdat(i,j) == 36 .OR.                                &
                 vtypdat(i,j) == 38 .OR.                                &
                 vtypdat(i,j) == 39 .OR.                                &
                 vtypdat(i,j) == 55 ) THEN

        vtypdat(i,j) = 10

      ELSE IF ( vtypdat(i,j) == 44 .OR.                                 &
                 vtypdat(i,j) == 45 ) THEN

        vtypdat(i,j) = 11

      ELSE IF ( vtypdat(i,j) == 59 .OR.                                 &
                 vtypdat(i,j) == 63 .OR.                                &
                 vtypdat(i,j) == 64 ) THEN

        vtypdat(i,j) = 12

      ELSE IF ( vtypdat(i,j) == 51 ) THEN

        vtypdat(i,j) = 13

      ELSE IF ( vtypdat(i,j) == 00 .OR.                                 &
                 vtypdat(i,j) == 65 .OR.                                &
                 vtypdat(i,j) == 66 .OR.                                &
                 vtypdat(i,j) == 67 .OR.                                &
                 vtypdat(i,j) == 68 .OR.                                &
                 vtypdat(i,j) == 72 .OR.                                &
                 vtypdat(i,j) == 73 ) THEN

        vtypdat(i,j) = 14
      ELSE

        WRITE (6, '(a/a,i4,a,i4,a,i2/a)')                               &
            'The value for vegetation type is invalid.',                &
            'vtypdat(', i, ',', j, ') = ', vtypdat(i,j),                &
            'Program stops in subroutine GETVTYP.'

        STOP

      END IF

    END DO
  END DO


!
!-----------------------------------------------------------------------
!
!  Fill the vegetation type into the model grid by chosing the values
!  at the nearest data grid points.
!
!-----------------------------------------------------------------------
!
  ixmin = MIN( xdims-1,                                                 &
               INT( ( glon(1,1)-dtlon(1) ) / resl ) + 1,                &
               INT( ( glon(1,ny)-dtlon(1) ) / resl ) + 1 )

  ixmin = MAX( 1, ixmin - 5 )

  jymin = MIN( ydims-1,                                                 &
               INT( ( glat(1,1)-dtlat(1) ) / resl ) + 1,                &
               INT( ( glat(nx,1)-dtlat(1) ) / resl ) + 1 )

  jymin = MAX( 1, jymin - 5 )

  ixmax = MIN( xdims,                                                   &
               nint( ( glon(nx,1)-dtlon(1) ) / resl ) + 1,              &
               nint( ( glon(nx,ny)-dtlon(1) ) / resl ) + 1 )

  ixmax = MIN( xdims, ixmax + 5 )

  jymax = MIN( ydims,                                                   &
               nint( ( glat(1,ny)-dtlat(1) ) / resl ) + 1,              &
               nint( ( glat(nx,ny)-dtlat(1) ) / resl ) + 1 )

  jymax = MIN( ydims, jymax + 5 )

  DO j = 1, ny
    DO i = 1, nx

      r1 = 1.0E20

      DO jy = jymin, jymax
        DO ix = ixmin, ixmax
          r = ( dtlon(ix) - glon(i,j) ) ** 2                            &
              + ( dtlat(jy) - glat(i,j) ) ** 2
          IF ( r < r1 ) THEN
            r1 = r
            ii = ix
            jj = jy
          END IF
        END DO
      END DO

      vegtyp(i,j) = vtypdat(ii,jj)

    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  The following DO loop is just for the Vortex experiment. On the
!  south boundary Vortex domain the vegetation type changes sharply,
!  that will pretty much influence the boundary conditions. Therefore
!  the following DO loop will use the values of inner grids to those
!  boundary grids.
!
!-----------------------------------------------------------------------
!
!  DO 300 j = 1, 7
!  DO 300 i = 1, nx
!    vegtyp(i,j) = vegtyp(i,8)
!300  CONTINUE

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


SUBROUTINE getlai( nx,ny, vtypdat, ndvi, lai ) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the Leaf Area Index (lai) from NDVI data and interpolate
!  them into the model domain.
!
!  The calculation of LAI depends on vegetation type: grass or tree.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!
!  3/22/94
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!
!  vtyp     Vegetation type;
!  ndvi     NDVI data array
!
!  OUTPUT:
!
!  lai      LAI array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
!
  INTEGER :: vtypdat(nx,ny)  ! Vegetation type: 1 -- grass; 2 -- tree
  REAL :: ndvi(nx,ny)  ! Soil data array
!
  REAL :: lai (nx,ny)  ! Soil type array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j = 1, ny
    DO i = 1, nx
      IF ( vtypdat(i,j) == 09 .OR.             & ! Water & ice
       vtypdat(i,j) == 14  ) THEN

      lai(i,j) = 0.

    ELSE IF ( vtypdat(i,j) == 01 .OR.                                   &
               vtypdat(i,j) == 02 .OR.                                  &
               vtypdat(i,j) == 03 .OR.                                  &
               vtypdat(i,j) == 04 .OR.                                  &
               vtypdat(i,j) == 05 .OR.                                  &
               vtypdat(i,j) == 09 .OR.                                  &
               vtypdat(i,j) == 10 .OR.                                  &
               vtypdat(i,j) == 11 .OR.                                  &
               vtypdat(i,j) == 12 .OR.                                  &
               vtypdat(i,j) == 13 .OR.                                  &
               vtypdat(i,j) == 14 ) THEN

      lai(i,j) = - LOG( ( 1. - ndvi(i,j)/.915 ) / .83 ) / 0.96

!    ELSEIF ( vtypdat(i,j) .eq. 04 .or.  ! uncomment this block for 4 & 5
!    :           vtypdat(i,j) .eq. 05 ) THEN
!
!      lai(i,j) = 0.5
!    :             * ( - log( ( 1. - ndvi(i,j)/.915 ) / .83 ) / 0.96
!    :               + 1.625 * exp( ndvi(i,j) / 0.34 ) )

    ELSE

      lai(i,j) = 1.625 * EXP( ndvi(i,j) / 0.34 )

    END IF

    lai(i,j) = MAX( lai(i,j), 0. )

  END DO
  END DO

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

SUBROUTINE getndvi1(xdims,ydims,nx,ny,                                  & 1
           x1,y1,ndvibyt,xdat,ydat,ndvi)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the NDVI data (ndvi) from byte-NDVI data and
!  fill them into the model domain by chosing the average
!  value of most frequently appeared in a grid.
!
!-----------------------------------------------------------------------
!
!  AUTHOR:
!  Leilei Wang ,Vince Wong
!  8/11/97
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!  xdims    X-dimsion size of the data set
!  ydims    Y-dimsion size of the data set
!
!  glon     Longitude values of model grid points
!  glat     Latitude  values of model grid points
!
!  dtlon    Longitude values of data grid points
!  dtlat    Latitude  values of data grid points
!
!  ndvibyt  Ndvi data array
!
!  OUTPUT:
!
!  ndvi     Ndvi type array
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
  INTEGER :: xdims        ! X-dimension size of the data set
  INTEGER :: ydims        ! Y-dimension size of the data set
  REAL :: x1  (nx,ny)  ! X_coordinate of model grid
  REAL :: y1  (nx,ny)  ! X_coordinate of model grid

  REAL :: xdat(xdims)  ! X_coor. of soil data
  REAL :: ydat(ydims)  ! Y_coor. of soil data

  INTEGER :: ndvibyt (xdims,ydims)        ! byte_NDVI data array
  REAL :: ndvi    (nx,ny)              ! NDVI  array
  REAL :: dresl                        ! data resolution

  REAL :: tem2    (20)                 ! Temporary array
  REAL :: tem3    (20)                 ! Temporary array
  INTEGER :: item4(20)
  REAL, ALLOCATABLE :: tem1(:,:)                         ! add by wyh
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j, iis,is,ii,jj
  INTEGER :: imax,count
  INTEGER :: xleft,xright,yleft,yright
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  ALLOCATE(tem1(xdims, ydims))                        ! add by wyh

  DO j = 1, ydims
    DO i = 1, xdims

      tem1(i,j) = MAX( 0.0, ( ndvibyt(i,j) - 100. ) / 100. )

    END DO
  END DO

!-----------------------------------------------------------------------
!
!  Fill the NDVI into the model grid using statistics results
!  by chosing NDVI  of highest weights  for each model grid point
!
!-----------------------------------------------------------------------
!
  dresl =1000.

  DO j = 1,ny
    DO i = 1,nx

      DO is = 1,20
        tem2(is) = -1.0+(is-1)*0.1
        tem3(is) = 0.0
        item4(is) = 0
      END DO

      count = 0
      imax = 0

      IF(dx < dresl) THEN
        xleft = INT( (x1(i,j) - xdat(1)-dresl/2)/dresl )+1
        xright= INT( (x1(i,j) - xdat(1)+dresl/2)/dresl )
      ELSE
        xleft = INT( (x1(i,j) - xdat(1)-dx/2)/dresl )+1
        xright= INT( (x1(i,j) - xdat(1)+dx/2)/dresl )
      END IF

      IF(dy < 1000.0) THEN
        yleft = INT( (y1(i,j) - ydat(1)-dresl/2)/dresl )+1
        yright= INT( (y1(i,j) - ydat(1)+dresl/2)/dresl )
      ELSE
        yleft = INT( (y1(i,j) - ydat(1)-dy/2)/dresl )+1
        yright= INT( (y1(i,j) - ydat(1)+dy/2)/dresl )
      END IF

      DO ii = xleft, xright
        DO jj = yleft, yright

          IF(xleft < 1.OR.yleft < 1) GO TO 800
          IF(xright > xdims.OR.yright > ydims) GO TO 800

          IF(tem1(ii,jj) >= -1.0.AND.tem1(ii,jj) <= 1.0) count=count + 1
          DO is=1,20
            IF(tem1(ii,jj) >= tem2(is).AND. tem1(ii,jj) < tem2(is+1)) THEN
              tem3(is) = tem3(is) + tem1(ii,jj)
              item4(is) = item4(is) + 1
            END IF
          END DO
          800     CONTINUE
        END DO
      END DO

      IF(count == 0) GO TO  100

      imax = 1
      DO iis=1,20
        IF(item4(iis) >= item4(imax)) THEN
          imax=iis
        END IF
      END DO

      100     IF(imax /= 0)THEN
        ndvi(i,j) = tem3(imax)/item4(imax)
      END IF

    END DO
  END DO

  DEALLOCATE(tem1)                    ! add by wyh

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


SUBROUTINE getndvi2(xdims,ydims, nx,ny,                                 & 1
           dtlon,dtlat,resl, glon,glat,                                 &
           ndvibyt, ndvi)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the NDVI data (ndvi) from byte-NDVI data and
!  fill them into the model domain by chosing the values at the
!  nearest data grid points.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!
!  2/20/94
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!  xdims    X-dimsion size of the data set
!  ydims    Y-dimsion size of the data set
!
!  glon     Longitude values of model grid points
!  glat     Latitude  values of model grid points
!
!  dtlon    Longitude values of data grid points
!  dtlat    Latitude  values of data grid points
!
!  ndvibyt  Byte-NDVI data
!
!  OUTPUT:
!
!  ndvi     Vegetation type array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
  INTEGER :: xdims        ! X-dimsion size of the data set
  INTEGER :: ydims        ! Y-dimsion size of the data set
!
  REAL :: glon(nx,ny)  ! Longitude values of model grid points
  REAL :: glat(nx,ny)  ! Latitude  values of model grid points
!
  REAL :: dtlon(xdims) ! Longitude values of data grid points
  REAL :: dtlat(ydims) ! Latitude  values of data grid points
  REAL :: resl         ! Data resolution
!
  INTEGER :: ndvibyt(xdims,ydims)  ! Byte-NDVI data array
!
  REAL :: ndvi(nx,ny)  ! NDVI data array

!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j, ii, jj, ix,jy
  INTEGER :: ixmin,ixmax,jymin,jymax
  REAL :: r,r1

  REAL, ALLOCATABLE :: tem1(:,:)                ! add by wyh
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  ALLOCATE(tem1(xdims, ydims))                     ! add by wyh

  WHERE (ndvibyt > 100)
    tem1 = (ndvibyt - 100. ) / 100. 
  ELSEWHERE
    tem1 = 0.0
  END WHERE

!
!-----------------------------------------------------------------------
!
!  Fill the vegetation type into the model grid by chosing the values
!  at the nearest data grid points.
!
!-----------------------------------------------------------------------
!
  ixmin = MIN( xdims-1,                                                 &
               INT( ( glon(1,1)-dtlon(1) ) / resl ) + 1,                &
               INT( ( glon(1,ny)-dtlon(1) ) / resl ) + 1 )

  ixmin = MAX( 1, ixmin - 5 )

  jymin = MIN( ydims-1,                                                 &
               INT( ( glat(1,1)-dtlat(1) ) / resl ) + 1,                &
               INT( ( glat(nx,1)-dtlat(1) ) / resl ) + 1 )

  jymin = MAX( 1, jymin - 5 )

  ixmax = MIN( xdims,                                                   &
               nint( ( glon(nx,1)-dtlon(1) ) / resl ) + 1,              &
               nint( ( glon(nx,ny)-dtlon(1) ) / resl ) + 1 )

  ixmax = MIN( xdims, ixmax + 5 )

  jymax = MIN( ydims,                                                   &
               nint( ( glat(1,ny)-dtlat(1) ) / resl ) + 1,              &
               nint( ( glat(nx,ny)-dtlat(1) ) / resl ) + 1 )

  jymax = MIN( ydims, jymax + 5 )

  DO j = 1, ny
    DO i = 1, nx

      r1 = 1.0E20

      DO jy = jymin, jymax
        DO ix = ixmin, ixmax
          r = ( dtlon(ix) - glon(i,j) ) ** 2                            &
              + ( dtlat(jy) - glat(i,j) ) ** 2
          IF ( r < r1 ) THEN
            r1 = r
            ii = ix
            jj = jy
          END IF
        END DO
      END DO

      ndvi(i,j) = tem1(ii,jj)

    END DO
  END DO

  DEALLOCATE(tem1)
  RETURN
END SUBROUTINE getndvi2
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GETRFNS                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE getrfns( nx,ny, vtypdat, roufns ) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Generate the surface roughness (roufns) for ARPS model from the
!  table vs vegetation type.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!
!  4/19/94
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!
!  vtyp     Vegetation type;
!
!  OUTPUT:
!
!  roufns   LAI array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
!
  INTEGER :: vtypdat(nx,ny)  ! Vegetation type: 1 -- grass; 2 -- tree
!
  REAL :: roufns(nx,ny)  ! Surface roughness

  REAL :: rfns(14)       ! Table data for surface roughness vs vtyp
  DATA rfns/ 0.002, 0.030, 0.030, 0.100, 0.200, 0.500, 0.750,           &
             1.000, 0.005, 0.150, 0.100, 0.100, 0.050, 0.002/
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j = 1, ny
    DO i = 1, nx

      roufns(i,j) = rfns(vtypdat(i,j))

    END DO
  END DO

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


SUBROUTINE getveg( nx,ny, vtypdat, veg ) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Generate the vegetation fraction (roufns) for ARPS model from the
!  table vs vegetation type.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Yuhe Liu
!
!  02/07/95
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!
!  vtypdat  Vegetation type;
!
!  OUTPUT:
!
!  veg      veg array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
!
  INTEGER :: vtypdat(nx,ny)  ! Vegetation type: 1 -- grass; 2 -- tree
!
  REAL :: veg(nx,ny)      ! Vegetation fraction

  REAL :: vegtab(14)      ! Table data for surface roughness vs vtyp
  DATA vegtab/ 0.10, 0.10, 0.60, 0.40, 0.40, 0.90, 0.99,                &
               0.99, 0.01, 0.30, 0.99, 0.40, 0.20, 0.00/
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO j = 1, ny
    DO i = 1, nx

      veg(i,j) = vegtab(vtypdat(i,j))

    END DO
  END DO

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


SUBROUTINE gtlkuptbl( lkupfl, vegtbl, laitbl, z0tbl ) 1,2
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in vegtbl, a table of veg, roufns, & lai as a funtion of
!  vegitation type.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!
!  02/20/97
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  lkupfl   File name of the vegitation table.
!
!  OUTPUT:
!
!  vegtbl   Mapping of vegitation fraction as a fuction of
!           vegitation type.
!
!  laitbl   Mapping of leaf area index as a fuction of
!           vegitation type.
!
!  vegtbl   Mapping of vegetation fraction as a fuction of
!           vegetation type.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
!
  REAL :: vegtbl(nvegtyp) ! veg as a fuction of vtype.
  REAL :: laitbl(nvegtyp) ! lai as a fuction of vtype.
  REAL :: z0tbl(nvegtyp)  ! z0 (roufns) as a fuction of vtype.

  CHARACTER (LEN=80) :: lkupfl ! File name of vegetation table.
!
  INTEGER :: iv, vin
  REAL :: vegin, laiin, z0in
  CHARACTER (LEN=80) :: instr
  INTEGER :: n, nskip
  INTEGER :: lenfl, flunit
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!-----------------------------------------------------------------------
!
!  Initialize variable.
!
!-----------------------------------------------------------------------
!
  DO iv = 1,nvegtyp

    vegtbl(iv) = 0.0
    laitbl(iv) = 0.0
    z0tbl(iv) = 0.000001

  END DO

  nskip = 0
!
!-----------------------------------------------------------------------
!
!  Read in the vegetation table.
!
!-----------------------------------------------------------------------
!
  lenfl = 70
  CALL strlnth( lkupfl, lenfl )
  CALL getunit( flunit )
  WRITE(6,'(a)') '    Opening vegetation lookup table file '            &
                 //lkupfl(1:lenfl)
  OPEN (UNIT = flunit, FILE = lkupfl(1:lenfl), STATUS = 'old')

  100   CONTINUE

  READ (flunit,'(a80)') instr

  IF ((instr(1:1) /= "c") .AND. (instr(1:1) /= "c") .AND.               &
      (instr(1:1) /= "&") .AND. (instr(1:1) /= "!") .AND.               &
      (instr(1:1) /= "#")) GO TO 700

  nskip = nskip + 1
  GO TO 100

  700   CONTINUE

  CLOSE ( flunit )
  OPEN (UNIT = flunit, FILE = lkupfl(1:lenfl), STATUS = 'old')

  DO n = 1,nskip

    READ (flunit,'(a80)') instr

  END DO

  500   READ (flunit,*,END=200) vin, vegin, laiin, z0in

  IF ( vin > nvegtyp .OR. vin < 1 ) THEN
    WRITE (*,*) "ERROR, entry skipped.  Veg type out of bounds:",vin
  END IF

  vegtbl(vin) = vegin
  laitbl(vin) = laiin
  z0tbl(vin) = z0in

  GO TO 500

  200   CONTINUE

  CLOSE ( flunit )
  CALL retunit ( flunit )

  WRITE (*,*) "Vegataion table used:"
  WRITE (*,*) "vtype   veg     lai     z0"

  DO n = 1,nvegtyp

    WRITE (*,1000) n,vegtbl(n),laitbl(n),z0tbl(n)

  END DO

  1000  FORMAT(i2,7X,3F8.5)

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


SUBROUTINE gtlaitbl( nx,ny, vtyp, laitbl, lai ) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Assign the Leaf Area Index (lai) as a function of vegetation type.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett & Richard Carpenter
!
!  02/20/97
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!
!  vtyp     Vegetation type;
!  laitbl
!
!  OUTPUT:
!
!  lai      LAI array
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
!
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
!
  INTEGER :: vtyp(nx,ny)  ! Vegetation type
  REAL :: laitbl(nvegtyp)
!
  REAL :: lai (nx,ny)  ! Soil type array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  DO j = 1, ny
    DO i = 1, nx

      lai(i,j) = laitbl(vtyp(i,j))

    END DO
  END DO

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


SUBROUTINE gtvegtbl( nx,ny, vtyp, vegtbl, veg ) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Assign the vegetation fraction (veg) as a function of vegetation type.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett & Richard Carpenter
!
!  02/20/97
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!
!  vtyp     Vegetation type;
!  vegtbl
!
!  OUTPUT:
!
!  veg      vegetation fraction
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
!
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
!
  INTEGER :: vtyp(nx,ny)  ! Vegetation type
  REAL :: vegtbl(nvegtyp)
!
  REAL :: veg (nx,ny)  ! Soil type array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  DO j = 1, ny
    DO i = 1, nx

      veg(i,j) = vegtbl(vtyp(i,j))

    END DO
  END DO

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


SUBROUTINE gtrfnstbl( nx,ny, vtyp, z0tbl, roufns ) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Assign the roughness length (roufns) as a function of vegetation type.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett & Richard Carpenter
!
!  02/20/97
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  nx       Number of grid points in the x-direction
!  ny       Number of grid points in the y-direction
!
!  vtyp     Vegetation type;
!  z0tbl
!
!  OUTPUT:
!
!  roufns      roughness length
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'arpssfc.inc'
!
  INTEGER :: nx           ! Number of grid points in the x-direction
  INTEGER :: ny           ! Number of grid points in the y-direction
!
  INTEGER :: vtyp(nx,ny)  ! Vegetation type
  REAL :: z0tbl(nvegtyp)
!
  REAL :: roufns (nx,ny)  ! Soil type array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  DO j = 1, ny
    DO i = 1, nx

      roufns(i,j) = z0tbl(vtyp(i,j))

    END DO
  END DO

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


SUBROUTINE readbsq (filename, lun, inband, inrow, incol, inbyte,        & 1,1
           iband, irow, icol, nband, nrow, ncol, outarray,colomn,roww)
!
!-----------------------------------------------------------------------
!
!  Purpose:
!  This subroutine extracts a specified window (rectangular subset)
!  from a band-sequential (BSQ) binary data array and places it into
!  an INTEGER array specified by the calling program. 
!
!-----------------------------------------------------------------------
!
!  01/96   Initial version, R. A. White, PSU/ESSC,
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  INTEGER     inrow     Number of rows in input data array.
!  INTEGER     incol     Number of columns in input data array.
!  INTEGER     inbyte    Number of bytes per input data value.
!  INTEGER     iband     Initial band for which data are desired.
!  INTEGER     irow      Initial row for which data are desired:
!  INTEGER     icol      Initial column for which data are
!  INTEGER     nband     Number of consecutive bands to be
!                          extracted from input data array
!  INTEGER     nrow      Number of consecutive rows to be
!                          extracted from each band of data array.
!  INTEGER     ncol      Number of consecutive columns to be
!                          extracted from each row of data array.
!  OUTPUT:
!
!  INTEGER     outarray(ncol,nrow,nband)
!                            Array to receive the data from input
!                            file.  If nband = 1, the array may be
!                            declared as outarray(ncol,nrow).
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
!  PARAMETER (maxinbyte =16384)              !  comment out by wyh

  CHARACTER (LEN=*) :: filename
  INTEGER :: lun, inband, inrow, incol, inbyte, iband, irow, icol
  INTEGER :: nband, nrow, ncol
  INTEGER :: colomn,roww
  INTEGER :: outarray(colomn,roww)

  INTEGER :: band, row, col
  CHARACTER (LEN=1) :: inbuf(incol*inbyte)    !  add by wyh

  INTEGER :: lenfl
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  lenfl = LEN( filename )
  CALL strlnth( filename,lenfl )

  WRITE(6,'(a,a)') '    Opening ',filename(1:lenfl)
  OPEN (UNIT=lun, FILE=filename(1:lenfl), ACCESS='DIRECT',              &
                  FORM='UNFORMATTED',     STATUS='OLD',                 &
                  RECL=incol*inbyte)

  DO band=iband,iband+nband-1
    DO row=irow,irow+nrow-1
      READ (lun, REC=row+(band-1)*inrow)                                &
          (inbuf(i), i=1,incol*inbyte)
      DO col=icol,icol+ncol-1
        outarray(col-icol+1, row-irow+1) = ICHAR(inbuf(col))
      END DO
    END DO
  END DO
  CLOSE (lun)

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


SUBROUTINE lae_lltoxy(m,n,ctrlat,ctrlon,lat,lon,x,y) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Determine x,y from latitude, longitude coordinates in
!  Lambert Azimuthal equal area map projection
!
!-----------------------------------------------------------------------
!  AUTHOR: Leilei Wang , Vince Wong and Ming Xue
!
!   3/25/97.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  REAL :: lat,lon
  REAL :: x,y
  INTEGER :: m,n
  REAL :: d2rad,eradius
  PARAMETER (d2rad=3.141592654/180.,                                    &
             eradius = 6371000. )      ! mean earth radius in m

  REAL :: ctrlat,ctrlon                   ! Lat. & Lon. of projection
                                          ! center
  REAL :: const1,k

  const1=SIN(ctrlat*d2rad)*SIN(lat*d2rad)+COS(ctrlat*d2rad)*            &
         COS(lat*d2rad)*COS((lon-ctrlon)*d2rad)
  k = SQRT(2.0/(1+const1))
  x = eradius*k*COS(lat*d2rad)*SIN((lon-ctrlon)*d2rad)
  y = eradius*k*( COS(ctrlat*d2rad)*SIN(lat*d2rad)-                     &
      SIN(ctrlat*d2rad)*COS(lat*d2rad)*COS((lon-ctrlon)*d2rad) )

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


SUBROUTINE lae_xytoll(m,n,x,y,lat,lon)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Determine latitude, longitude coordinates from  the given x,y in
!  Lambert Azimuthal equal area map projection
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Leilei Wang , Vince Wong and Ming Xue
!
!   3/25/97.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  REAL :: lat,lon
  REAL :: x,y
  INTEGER :: m,n
  REAL :: d2rad,eradius
  REAL :: r, c
  PARAMETER (d2rad=3.141592654/180.,                                    &
             eradius = 6371000. )        ! mean earth radius in m
  REAL :: ctrlat,ctrlon                     ! Lat. & Lon. of projection
  PARAMETER (ctrlat=45.0,ctrlon=-100.0)  ! center
!
  r = SQRT(x**2+y**2)
  c = 2* ASIN( r/(2.0*eradius) )

  lat=ASIN( COS(c)*SIN(ctrlat*d2rad)+y*SIN(c)*COS(ctrlat*d2rad)/r)
  lon=ctrlon*d2rad + ATAN (x*SIN(c)/(r*COS(ctrlat*d2rad)*COS(c)         &
       -y*SIN(ctrlat*d2rad)*SIN(c)) )
  lat = lat/d2rad
  lon = lon/d2rad
!
  RETURN
END SUBROUTINE lae_xytoll
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 PROGRAM GVEGFRAC                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE gvegfrac(nx,ny,vfrcdr,glat,glon,vegin1,vegin2,veg) 1,1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Prepare the green vegetation fraction array veg for use in the
!  ARPS soil model.  Calling this subroutine will overwrite and
!  previously defined values of veg that are within the data rich
!  regions.  The pre-existing data is generated using the veg
!  type/table conversion.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Dan Weber
!  04/15/98
!
!  MODIFICATIONS:
!
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  month    month of the forecast initialization time.
!
!  day      day of the initialization time.
!
!  OUTPUT:
!
!  veg      vegetation fraction.
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT   NONE

  INCLUDE   'arpssfc.inc'      ! Includes dimension size for data

  INTEGER :: nx                   ! Number of model grid points in the
                                  ! x direction.
  INTEGER :: ny                   ! Number of model grid points in the
                                  ! y direction.

  REAL :: glat (nx,ny)       ! Lat. of model grid
  REAL :: glon (nx,ny)       ! Lon. of model grid

  REAL :: veg    (nx,ny)     ! Vegetation fraction in model domain
  REAL :: vegin1 (gvnx,gvny) ! Green Vegetation fraction data.
                             ! for use in reading in the NESDIS
                             ! green vegetation data set. month#1
  REAL :: vegin2 (gvnx,gvny) ! Green Vegetation fraction data.
                             ! for use in reading in the NESDIS
                             ! green vegetation data set. month#2
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j
  REAL :: tema, temb, temc, temd, t1, t2, temx, temy
  INTEGER :: itema,itemb,m1,m2,ltdir
  CHARACTER (LEN=80) :: nfile1,nfile2 ! for reading in the green vegetation
                                      ! fraction data set.
  CHARACTER (LEN=70) :: vfrcdr        ! vegetation fraction directory name.
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'      ! Includes month and day information
!
!-----------------------------------------------------------------------
!
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!  determine the monthly data sets to use....
!

  ltdir=LEN(vfrcdr)
  CALL strlnth( vfrcdr , ltdir)
  IF( ltdir  == 0 ) THEN
    vfrcdr = '.'
    ltdir  =1
  END IF

  IF( vfrcdr(ltdir:ltdir) /= '/') THEN
    ltdir=ltdir+1
    vfrcdr(ltdir:ltdir)='/'
  END IF


  IF(month == 2)THEN  ! february...

    IF(day <= 14)THEN
      m1 = month-1
      m2 = month
    ELSE
      m1 = month
      m2 = month+1
    END IF

  ELSE                ! all other months...

    IF(day < 15)THEN
      m1 = month-1
      m2 = month
    ELSE
      m1 = month
      m2 = month+1
    END IF

  END IF               ! end of month definition if block

  IF (m1 == 0) m1 = 12
  IF (m2 == 13) m2 = 1

!
!  approximate the green vegetation data set for a particular day
!  from monthly averaged data sets.  (not all possible combinations
!  are included in the logic, just the majority)....
!

  IF(month == 2)THEN   !  28 day months...
    IF(day == 29)THEN  !  leap year...
      t1 = 0.5
      t2 = 1.0-t1
    ELSE IF(day > 14)THEN
      t1 = 1.0-FLOAT(day-14)/29.
      t2 = 1.0-t1
    ELSE IF(day <= 14)THEN
      t1 = FLOAT(day+16)/30.
      t2 = 1.0-t1
    END IF
  ELSE IF(month == 4.OR.month == 6.OR.month == 9.OR.                    &
           month == 11)THEN       ! 30 day months...
    IF(day > 15)THEN
      t1 = 1.0-FLOAT(day-15)/30.
      t2 = 1.0-t1
    ELSE IF(day <= 15)THEN
      t2 = FLOAT(day+16)/31.
      t1 = 1.0-t2
    END IF

  ELSE                 !  31 day months....
    IF(day > 15)THEN
      t1 = 1.0-FLOAT(day-15)/31.
      t2 = 1.0-t1
    ELSE IF(day <= 15)THEN
      t2 = FLOAT(day+15)/30.
      t1 = 1.0-t2
    END IF
  END IF    ! end of month if block.....
  PRINT *,' t1 and t1 in gvegfract are ',t1,t2
!
!  read in the green vegetation fraction (range 0.0 to 1.0)
!

!  finish writing the data set filename...

  nfile1(1:ltdir) = vfrcdr(1:ltdir)
  WRITE(nfile1(ltdir+1:ltdir+2),'(i2.2)') m1
  WRITE(nfile1(ltdir+3:ltdir+12),'(a10)') '.gvegfract'

  nfile2(1:ltdir) = vfrcdr(1:ltdir)
  WRITE(nfile2(ltdir+1:ltdir+2),'(i2.2)') m2
  WRITE(nfile2(ltdir+3:ltdir+12),'(a10)') '.gvegfract'

  WRITE (*,*) "Opening green veg file 1: ",nfile1(1:ltdir+12)," ..."
  OPEN (12,FILE=nfile1(1:ltdir+12),STATUS='old',FORM='formatted')
  DO j=1,gvny
    DO i=1,gvnx
      READ (12,'(i3)') itema
      vegin1(i,j) = FLOAT(itema)/100.0
    END DO
  END DO
  CLOSE(12)

  WRITE (*,*) "Opening green veg file 2: ",nfile2(1:ltdir+12)," ..."
  OPEN (13,FILE=nfile2(1:ltdir+12),STATUS='old',FORM='formatted')
  DO j=1,gvny
    DO i=1,gvnx
      READ (13,'(i3)') itema
      vegin2(i,j) = FLOAT(itema)/100.0
    END DO
  END DO
  CLOSE(13)

!
!  Compute the veg field by using a linear interpolation (x,y and t).
!  Note the lower left data point is -55.0N and 0W.
!  The j-0.5 is correct (j is an integer).
!  The first box is at lat=75.016, which is considered the center.
!

  temx = 75.016-0.5*0.144
  temy = -55.00
  DO j=1,ny
    DO i=1,nx
      IF (glat(i,j) <= temx.AND.glat(i,j) >= temy)THEN
                                            ! use green veg
! data base
        itema  = glon(i,j)/0.144
        temb   = glon(i,j)-itema*0.144
        tema   = 1.0-temb

        itemb  =(glat(i,j)+55.00)/0.144
        temd   = glat(i,j)-(-55.00+itemb*0.144)
        temc   = 1.0-temd

        veg(i,j) = t1*(temc*(tema*vegin1(itema,itemb)+                  &
                             temb*vegin1(itema+1,itemb))+               &
                       temd*(tema*vegin1(itema,itemb+1)+                &
                             temb*vegin1(itema+1,itemb+1)))             &
                 + t2*(temc*(tema*vegin2(itema,itemb)+                  &
                             temb*vegin2(itema+1,itemb))+               &
                       temd*(tema*vegin2(itema,itemb+1)+                &
                             temb*vegin2(itema+1,itemb+1)))

      END IF
    END DO
  END DO
  PRINT *,'interpolated veg fraction at ',2,2,' is ',veg(2,2)

  RETURN
END SUBROUTINE gvegfrac