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


SUBROUTINE mksoilvar(nx,ny,nz,nzsoil,nstyps,                            & 1,41
           mxdays,mxstns,                                               &
           xs,ys,                                                       &
           hterain,zsc,                                                 &
           ptbar,ptprt,                                                 &
           pbar,pprt,                                                   &
           zpsoil,tsoil,qsoil,wetcanp,snowdpth,                         &
           soiltyp,stypfrct,                                            &
           precfile,staid,                                              &
           iptstn,jptstn,iptapi0,jptapi0,                               &
           api0,obsprec,difprec,totpreca,prec,                          &
           itema,initapi,api1,api2,kk2dep,                              &
           totwt,dprec,tem1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Set the surface variables for ARPS model.
!  Includes calculation of API from precip data to initialize
!  soil moisture.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: John Mewes (precip and API) and Keith Brewster (temps)
!  March, 1997
!
!  MODIFICATION HISTORY:
!
!  04/11/1997 (Keith Brewster)
!  Added processing of NCEP real-time precip data to API calculation.
!  at the same time delete the table data array veg(14).
!
!  04/15/1997 (Keith Brewster)
!  Added option of reading in pre-calculated soil moisture data,
!  since the real-time precip data are typically only once per day.
!  Minor reorganization of parameters. Bells and whistles added.
!  Renamed to version 2.0
!
!  01/05/1998 (Donghai Wang)
!  Fixed a problem in the case of k1< 0.0.
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  2000/01/03 (Gene Bassett)
!  Renamed mstinit to soilinit2 and added capability to update soil
!  data over water and snow cover with data in a soil data file.
!
!  05/25/2002  (J. Brotzge)
!  Added additional soil variables for new soil scheme.
!
!  1 June 2002 Eric Kemp
!  Soil variable updates.
!
!  09/2004 (J.Case, ENSCO Inc.)
!  Included provisions for interpolating RUC SSTs to the ARPS
!  grid at water grid points.
!
!-----------------------------------------------------------------------
!
!  Set various parameters used in API calculation
!
!-----------------------------------------------------------------------
!

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

  INTEGER :: nx,ny,nz,nstyps
  INTEGER :: nzsoil 
  INTEGER :: mxdays,mxstns
!
  REAL :: xs    (nx)        ! The x-coord. of the physical and
                            ! computational grid. Defined at u-point.
  REAL :: ys    (ny)        ! The y-coord. of the physical and
                            ! computational grid. Defined at v-point.

  REAL :: ptprt  (nx,ny,nz)    ! Perturbation potential temperature (K)
  REAL :: pprt   (nx,ny,nz)    ! Perturbation pressure (Pascal)
                               ! humidity (kg/kg)
  REAL :: ptbar  (nx,ny,nz)    ! Base state potential temperature (K)
  REAL :: pbar   (nx,ny,nz)    ! Base state pressure (Pascal)

  INTEGER :: soiltyp(nx,ny,nstyps)    ! Soil type
  REAL :: stypfrct  (nx,ny,nstyps)

  INTEGER :: soilinit2         ! Moisture initialization option
  INTEGER :: prdata            ! Precip data option
  CHARACTER (LEN=80) :: apifile      ! Name of file containing names of the
!                              precip data files.
  CHARACTER (LEN=80) :: apiinit      ! File containing initial API's and
!                              their corresponding locations
  CHARACTER (LEN=80) :: prcpdir      ! Disk directory containing NCEP precip files
  CHARACTER (LEN=80) :: prcplst      ! File containing list of NCEP precip stns

! J.Case, ENSCO Inc, 9/16/2004 -- Changed to include hour in API dates.
  CHARACTER (LEN=13) :: inidate      ! Day for start of precip files
  CHARACTER (LEN=13) :: fnldate      ! Day for stop of precip files

  REAL :: respapi              ! Desired response of the initial API
!                              analysis to wavelengths equal to
!                              2 x (mean initAPI stn sep).
  INTEGER :: ndays             ! Number of days between initial date
!                              and history file date
  REAL :: k1,k2                ! API depletion minimum coefficients for
!                              the ground surface layer and the root
!                              zone layer
  REAL :: range                ! Range parameter in Barnes' wt. fcn.
!                              denominator
  REAL :: respprec             ! Desired response of the final daily
!                              precip analyses at wavelengths equal
!                              to 2 x (mean stn sep).
  REAL :: gamma(2)             ! Range multipliers for 1st and 2nd
!                              passes
!
!-----------------------------------------------------------------------
!
!  API generation arrays:
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=80) :: precfile(mxstns)
  CHARACTER (LEN=7) :: staid(mxstns)
  INTEGER :: iptstn(mxstns)
  INTEGER :: jptstn(mxstns)
  INTEGER :: iptapi0(mxstns)
  INTEGER :: jptapi0(mxstns)
  REAL :: api0(mxstns)
  REAL :: obsprec(mxdays,mxstns)
  REAL :: difprec(mxdays,mxstns)
  REAL :: totpreca(mxstns)
  REAL :: prec(nx,ny)
!
!-----------------------------------------------------------------------
!
!  Computed variables
!
!-----------------------------------------------------------------------
!
  REAL :: zpsoil(nx,ny,nzsoil)          ! Soil layer height (m)
  REAL :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K)
  REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3)
  REAL :: wetcanp(nx,ny,0:nstyps)       ! Canopy water amount
  REAL :: snowdpth(nx,ny)               ! Snow depth (m)

  REAL :: hterain(nx,ny)
  REAL :: zsc(nx,ny,nz)
!
!-----------------------------------------------------------------------
!
!  Work Arrays
!
!-----------------------------------------------------------------------
!
  INTEGER :: itema(nx,ny)
  REAL :: initapi(nx,ny,nstyps)
  REAL :: api1(nx,ny,nstyps)
  REAL :: api2(nx,ny,nstyps)
  REAL :: kk2dep(nx,ny)
  REAL :: totwt(nx,ny)
  REAL :: dprec(nx,ny)
  REAL :: tem1(nx,ny,nz)

  ! temparary arrays for soilinit2 option 5:
  REAL, ALLOCATABLE :: tsoil2 (:,:,:,:)
  REAL, ALLOCATABLE :: qsoil2 (:,:,:,:)
  REAL, ALLOCATABLE :: wetcanp2(:,:,:)
  INTEGER, ALLOCATABLE :: soiltyp2(:,:,:)
!
!-----------------------------------------------------------------------
!
!  API parameters
!
!-----------------------------------------------------------------------
!
  COMMON /apicom1/ apiinit,apifile,prcpdir,prcplst,inidate,fnldate
  COMMON /apicom2/ soilinit2,prdata
  COMMON /apicom3/ respapi,k1,k2,range,respprec,gamma
!
!-----------------------------------------------------------------------
!
!  Map projection variables
!
!-----------------------------------------------------------------------
!
  REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'
  INCLUDE 'phycst.inc'
  INCLUDE 'soilcst.inc'
  INCLUDE 'indtflg.inc'
!
!-----------------------------------------------------------------------
!
!  Miscellaneous local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: sec24hr
  PARAMETER( sec24hr = 86400 )
!
  INTEGER :: i,j,ii,ij,is,is2,istn,jstn,iday,ix,iy,ireturn
  INTEGER :: ipt,jpt,out
  INTEGER :: tstns,totapi0,minday,maxday,length

! J.Case, ENSCO Inc., 9/16/2004 -- Changed variable name for times extracted from fnldate.
  INTEGER :: iniyr,inimo,inidy,inihr,initi,nowti,jlday
  INTEGER :: fnlyr,fnlmo,fnldy,fnlhr,fnlmin,fnlsec
  INTEGER :: itime,iyr,imon,idy,ihr,imin,isec
  INTEGER :: idel,jdel,istr,iend,jstr,jend
  INTEGER :: i1,i2,i3,i4,i5,j1,j2,j3,j4,j5
!
  REAL :: tmp2,tmp3,prs2,prs3,w2,w3,ptk
  REAL :: p0inv,tmpk,pres
  REAL :: xpt,ypt,xctr,yctr,xll,yll
  REAL :: lat,lon,ddx,ddy
  REAL :: wmin,wmax,kk1,kk2,adj1,adj2
  REAL :: apirng,apir2,r2,dist,const,thrdst
  REAL :: sumwt,sumapi,wgt,precip,maxprec
  REAL :: bias,rms,wsfcscl,wdpscl
  REAL :: mindist,totdist,total
  REAL :: d_0,d_1,d_1star,depth1
!
!  INTEGER :: tsfcout,tsoilout,wsfcout,wdpout,wcanpout,snowdout
!  INTEGER :: tsfcin,tsoilin,wsfcin,wdpin,wcanpin,snowdin
!  COMMON /intgcom/tsfcout,tsoilout,wsfcout,wdpout,wcanpout,snowdout

  INTEGER :: zpsoilout,tsoilout,qsoilout,wcanpout,snowdout
  INTEGER :: zpsoilin,tsoilin,qsoilin,wcanpin,snowdin
  COMMON /intgcom/ zpsoilin,tsoilout,qsoilout,wcanpout,snowdout

  REAL :: twater
  COMMON /realcom/twater

  REAL :: amax, amin

  INTEGER :: nbwater,nbw_max
  PARAMETER (nbw_max=128)
  REAL :: tbwater(nbw_max),blat1(nbw_max),blat2(nbw_max)
  REAL :: blon1(nbw_max),blon2(nbw_max)
  COMMON /bwcomi/ nbwater
  COMMON /bwcomr/ tbwater,blat1,blat2,blon1,blon2

  INTEGER :: ib
  INTEGER :: k

!
! J.Manobianco (3/15/2002) Add capability to read/set SST
! J.Manobianco/J.Case (2/6/2003)  Add capability to read RUC SST, soil moisture,
!               and soil temperature on LCC RUCH grid
!
! J.Case
! Other variables
!
  CHARACTER (LEN=10) :: infilename
  INTEGER :: nx_ext,ny_ext,ios
  INTEGER :: iscl(nx,ny),jscl(nx,ny)
  INTEGER :: iproj_ext
  INTEGER :: nisst, njsst, isst, jsst

!  40-km NAM model

  PARAMETER (nx_ext=185,ny_ext=129)

  REAL :: lat_arps(nx,ny),lon_arps(nx,ny)
  REAL :: lat_ruc(nx_ext,ny_ext),lon_ruc(nx_ext,ny_ext)
  REAL :: xs2d(nx,ny),ys2d(nx,ny)
  REAL :: dxfld(nx_ext),dyfld(ny_ext)
  REAL :: rdxfld(nx_ext),rdyfld(ny_ext)
  REAL :: slopey(nx_ext,ny_ext)
  REAL :: alphay(nx_ext,ny_ext)
  REAL :: betay(nx_ext,ny_ext)
  REAL :: blats, blons, res, sstlat, sstlon
  REAL :: sltk_sfc(nx_ext,ny_ext)
  REAL :: sltk_deep(nx_ext,ny_ext)
  REAL :: soim_sfc(nx_ext,ny_ext)
  REAL :: soim_deep(nx_ext,ny_ext)

  REAL :: lat_ext(nx_ext,ny_ext),lon_ext(nx_ext,ny_ext)
  REAL :: x0,y0,x0_ext,y0_ext
  REAL :: x_ext(nx_ext),y_ext(ny_ext)
  REAL :: swlat_ext, swlon_ext
  REAL :: dx_ext, dy_ext

  REAL :: sltk_sfc_interp(nx,ny)
  REAL :: sltk_deep_interp(nx,ny)
  REAL :: soim_sfc_interp(nx,ny)
  REAL :: soim_deep_interp(nx,ny)
  REAL :: goes_sst(nx,ny)
  REAL :: scale_ext, trlon, latnot_ext(2), trlon_ext

  LOGICAL :: sst_flag

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF (soilinit2 == 5) THEN
    ALLOCATE(tsoil2  (nx,ny,nzsoil,0:nstyps))
    ALLOCATE(qsoil2  (nx,ny,nzsoil,0:nstyps))
    ALLOCATE(wetcanp2(nx,ny,0:nstyps))
    ALLOCATE(soiltyp2(nx,ny,nstyps))
  END IF

! J.Case, ENSCO Inc., 9/16/2004 -- corrected this bug in assigning API date to soil output.
!  as a temporary fix, read year, month, day, AND hour from the arpssoil.input file.
!
  READ(inidate,'(i4,1x,i2,1x,i2,1x,i2)') iniyr,inimo,inidy,inihr
  CALL ctim2abss(iniyr,inimo,inidy,inihr,0,0, initi)
  READ(fnldate,'(i4,1x,i2,1x,i2,1x,i2)') fnlyr,fnlmo,fnldy,fnlhr
  fnlmin = 0
  fnlsec = 0
  curtim = 0

  CALL ctim2abss(fnlyr, fnlmo, fnldy, fnlhr, fnlmin, fnlsec, nowti)
  nowti=nowti+curtim

  ndays=((nowti-initi)/sec24hr)+1

  WRITE(6,'(a,i6)') ' Number of days of data for API : ',ndays

  IF(ndays > mxdays.AND.(soilinit2 == 2 .OR. soilinit2 == 3)) THEN
    WRITE(6,'(a,i6,/,a)')                                               &
        ' Number of days of data for API exceeds mxdays dimension: ',   &
        mxdays,'   Increase mxdays or change inidate.'
    STOP
  END IF

!
!-----------------------------------------------------------------------
!
!  Find the center lat,lon and the location (in map coord.) of the
!  lower left corner (xll,yll).
!
!-----------------------------------------------------------------------
!
  latnot(1)=trulat1
  latnot(2)=trulat2
  CALL setmapr(mapproj,sclfct,latnot,trulon)
  CALL lltoxy(1,1,ctrlat,ctrlon,xctr,yctr)
  WRITE(6,*) ' dx= ',dx,' dy= ',dy
  xll=xctr-(0.5*(nx-3)*dx)
  yll=yctr-(0.5*(ny-3)*dy)
!
  WRITE(6,'(2x,a,f10.2,a,f10.2//)')                                     &
      ' Grid center:  lat= ',ctrlat,' lon= ',ctrlon

!
  IF( sfcdat == 1 ) THEN

    DO j=1,ny
      DO i=1,nx
        soiltyp(i,j,1) = styp
        stypfrct(i,j,1) = 1.0
      END DO
    END DO

    DO is=2,nstyps
      DO j=1,ny
        DO i=1,nx
          soiltyp (i,j,is) = 0
          stypfrct(i,j,is) = 0.0
        END DO
      END DO
    END DO

  ELSE IF (sfcdat == 2 .OR. (sfcdat == 3.AND.landin /= 1) ) THEN

    DO is=1,nstyps
      DO j=1,ny
        DO i=1,nx
          soiltyp (i,j,is) = 0
          stypfrct(i,j,is) = 0.0
        END DO
      END DO
    END DO

    CALL readsfcdt(nx,ny,nstyps,sfcdtfl,dx,dy,                          &
             mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,       &
                   soiltyp,stypfrct,itema,tem1,tem1,tem1,tem1)
    WRITE(6,*) ' nstyp = ',nstyp

    IF( nstyp == 1 ) THEN

      DO j=1,ny-1
        DO i=1,nx-1
          stypfrct(i,j,1) = 1.0
        END DO
      END DO

    END IF

  ELSE IF (sfcdat == 3 .AND. landin == 1) THEN

    WRITE(6,'(1x,a/)')                                                  &
        'Surface property data in the history file was used.'

  ELSE

    WRITE(6,'(1x,a,i3,a/)')                                             &
        'Invalid surface data input option. sfcdat =',sfcdat,           &
        '. Program stopped in MKSOILVAR.'
    STOP

  END IF

  DO is=0,nstyp

    PRINT*,'In MKSOILVAR for is =', is

    CALL a3dmax0(tsoil(1,1,1,is),                                       &
         1,nx,1,nx-1,1,ny,1,ny-1,1,nzsoil,1,nzsoil, amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'tsoilmin= ', amin,', tsoilmax =',amax

    CALL a3dmax0(qsoil(1,1,1,is),                                       &
         1,nx,1,nx-1,1,ny,1,ny-1,1,nzsoil,1,nzsoil, amax,amin)
    WRITE(6,'(1x,2(a,e13.6))') 'qsoilmin= ', amin,', qsoilmax =',amax


  END DO


  !
  ! ENSCO inc. J. Case begin
  !
  sst_flag = .true.

  IF(sst_flag) THEN

!
! J.Manobianco (3/15/02) Add setorig call to use xytoll for SST interpolation

    CALL setorig (1,xll,yll)

!   Read in external data

    infilename='sltk.5cm'
    CALL RDEXTSST (nx_ext,ny_ext,infilename,                            &
                   sltk_sfc,blats,blons,res,nisst,njsst,sst_flag)
    print *,'SLTK5 VARS after RDEXT',nisst,njsst,sst_flag
    print *,sltk_sfc(1,1),sltk_sfc(nx_ext,ny_ext),                      &
            sltk_sfc(1,ny_ext),sltk_sfc(nx_ext,1)

    infilename='sltk.40cm'
    CALL RDEXTSST (nx_ext,ny_ext,infilename,                            &
                   sltk_deep,blats,blons,res,nisst,njsst,sst_flag)
    print *,'SLTK40 VARS after RDEXT',nisst,njsst,sst_flag
    print *,sltk_deep(1,1),sltk_deep(nx_ext,ny_ext),                    &
            sltk_deep(1,ny_ext),sltk_deep(nx_ext,1)

!    infilename='soim.1cm'
!    CALL RDEXTSST (nx_ext,ny_ext,infilename,                            &
!                  soim_sfc,blats,blons,res,nisst,njsst,sst_flag)
!    print *,'SOIM1 VARS after RDEXT',nisst,njsst,sst_flag
!    print *,soim_sfc(1,1),soim_sfc(nx_ext,ny_ext),                      &
!           soim_sfc(1,ny_ext),soim_sfc(nx_ext,1)
!
!    infilename='soim.100cm'
!    CALL RDEXTSST (nx_ext,ny_ext,infilename,                            &
!                 soim_deep,blats,blons,res,nisst,njsst,sst_flag)
!    print *,'SOIM100 VARS after RDEXT',nisst,njsst,sst_flag
!    print *,soim_deep(1,1),soim_deep(nx_ext,ny_ext),                    &
!           soim_deep(1,ny_ext),soim_deep(nx_ext,1)

! J.Case, ENSCO Inc. (2/11/2005) -- Reading in GOES SST data from FIT.

     goes_sst=-9999.0
     print *,'Reading GOES SST Product.'
     infilename='sst.data'
     OPEN (unit=45,file=infilename,form='FORMATTED',                     &
           status='OLD',iostat=ios,err=499)
     DO i=1,nx
     DO j=1,ny
       read (45,'(29x,f9.2)') goes_sst(i,j)
       if (goes_sst(i,j) > 0) then
         goes_sst(i,j)=goes_sst(i,j)+273.15
       endif
     ENDDO
     ENDDO
  
     print *,'Normal read of GOES SST data.'
     GOTO 498 

499  print *,'Error reading GOES SST file. IOSTAT = ',ios
498  CONTINUE
                                                                                                                                   
!
! =================================================================
!
!       Compute latitude and longitude of ARPS grid so that x and y
!       locations can be computed relative to the external grid.
!
                                                                                                                                   
    CALL xytoll (nx,ny,xs,ys,lat_arps,lon_arps)
    print *,lat_arps(1,1),lon_arps(1,1),                                 &
           lat_arps(nx-1,ny-1),lon_arps(nx-1,ny-1)
                                                                                                                                   
! J.Case, ENSCO Inc. (9/22/2004)
! Set up external grid information
!
                                                                                                                                   
    iproj_ext = 2        ! LCC
    scale_ext = 1.
    trlon_ext = 265.     ! NAM-40 grid projection
    latnot_ext(1) = 25.  ! NAM-40 grid projection
    latnot_ext(2) = 25.  ! NAM-40 grid projection

!  -- NAM 40 km model
    swlat_ext = 12.190   ! NAM-40 grid projection
    swlon_ext = -133.459 ! NAM-40 grid projection
    dx_ext = 40.63525E+03    ! NAM-40 grid projection
    dy_ext = 40.63525E+03    ! NAM-40 grid projection
                                                                                                                                   
! J.Case, ENSCO Inc. (9/22/2004)
! Set map projection to RUC.
!
                                                                                                                                   
    call setmapr (iproj_ext,scale_ext,latnot_ext,trlon_ext)
    call lltoxy (1,1,swlat_ext,swlon_ext,x0_ext,y0_ext)
                                                                                                                                   
! J.Case, ENSCO Inc.
! x,y locations of NAM:
!
                                                                                                                                   
    x_ext(1)=x0_ext
    y_ext(1)=y0_ext
                                                                                                                                   
    do i=2,nx_ext
      x_ext(i)=x_ext(i-1) + dx_ext
    enddo
    do j=2,ny_ext
      y_ext(j)=y_ext(j-1) + dy_ext
    enddo
                                                                                                                                   
! J.Case, ENSCO Inc.
! lat,lon of external grid:
!
                                                                                                                                   
    call xytoll(nx_ext,ny_ext,x_ext,y_ext,lat_ruc,lon_ruc)
                                                                                                                                   
! J.Case, cont.
! Find x,y locations of external grid.
!
                                                                                                                                   
    call setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
    call setorig(1,x0_ext,y0_ext)
    DO j=1,ny_ext
      call lltoxy(1,1,lat_ruc(1,j),lon_ruc(1,j),                         &
                        x_ext(1),y_ext(j))
    ENDDO
    DO i=1,nx_ext
      call lltoxy(1,1,lat_ruc(i,1),lon_ruc(i,1),                         &
                        x_ext(i),y_ext(1))
    ENDDO
                                                                                                                                   
! J.Case, cont.
!  Find x,y locations of ARPS grid points in terms of external grid.
!
                                                                                                                                   
    CALL lltoxy   (nx,ny,lat_arps,lon_arps,xs2d,ys2d)
    CALL setijloc (nx,ny,nx_ext,ny_ext,xs2d,ys2d,                        &
                     x_ext,y_ext,iscl,jscl)
                                                                                                                                   
!    write (22,*) xs2d
!    write (23,*) ys2d
!    write (24,*) iscl
!    write (25,*) jscl
!    write (26,*) lat_arps
!    write (27,*) lon_arps
                                                                                                                                   
! J.Case, cont.
! Restore map projection to ARPS.
                                                                                                                                   
    CALL setmapr(mapproj,sclfct,latnot,trulon)
    CALL setorig (1,xll,yll)
                                                                                                                                   
! J.Case, cont.
! Interpolate external grid data to ARPS grid points.
                                                                                                                                   
    IF (sst_flag) THEN
      CALL setdxdy(nx_ext,ny_ext,1,nx_ext,1,ny_ext,                      &
                   x_ext,y_ext,dxfld,dyfld,rdxfld,rdyfld)
! J.Case, cont.
! Interpolate shallow soil temperature.
!
      CALL fldint2d(nx,ny,nx_ext,ny_ext,                                 &
                    1,nx,1,ny,1,nx_ext,1,ny_ext,                         &
                    2,xs2d,ys2d,sltk_sfc,                                &
                    x_ext,y_ext,iscl,jscl,                               &
                    dxfld,dyfld,rdxfld,rdyfld,                           &
                    slopey,alphay,betay,sltk_sfc_interp)
! J.Case, cont.
! Interpolate deep soil temperature.
!
      CALL fldint2d(nx,ny,nx_ext,ny_ext,                                 &
                    1,nx,1,ny,1,nx_ext,1,ny_ext,                         &
                    2,xs2d,ys2d,sltk_deep,                               &
                    x_ext,y_ext,iscl,jscl,                               &
                    dxfld,dyfld,rdxfld,rdyfld,                           &
                    slopey,alphay,betay,sltk_deep_interp)
! J.Case, cont.
! Interpolate 1-cm soil moisture.
!
!      CALL fldint2d(nx,ny,nx_ext,ny_ext,                                 &
!                    1,nx,1,ny,1,nx_ext,1,ny_ext,                         &
!                    2,xs2d,ys2d,soim_sfc,                                &
!                    x_ext,y_ext,iscl,jscl,                               &
!                    dxfld,dyfld,rdxfld,rdyfld,                           &
!                    slopey,alphay,betay,soim_sfc_interp)
! J.Case, cont.
! Interpolate 100-cm soil moisture.
!
!      CALL fldint2d(nx,ny,nx_ext,ny_ext,                                 &
!                    1,nx,1,ny,1,nx_ext,1,ny_ext,                         &
!                    2,xs2d,ys2d,soim_deep,                               &
!                    x_ext,y_ext,iscl,jscl,                               &
!                    dxfld,dyfld,rdxfld,rdyfld,                           &
!                    slopey,alphay,betay,soim_deep_interp)                  
        ENDIF   ! sst_flag
                                                                                                                                   
!    write (20,*) sltk_sfc
!    write (21,*) sltk_sfc_interp
                                                                                                                                   
! J.Case, cont.
!  Restore origin back to North Pole (0,0).
!
                                                                                                                                   
    CALL setorig (1,0.,0.)
                                                                                                                               
  ENDIF  ! sst_flag .eq. .true.
                                                                                                                                   
! =================================================================
!

  IF(soilinit2 /= 5) THEN
    p0inv=1./p0
    DO j=1,ny-1
      DO i=1,nx-1
        tmp3 = ptbar(i,j,3) + ptprt(i,j,3)
        prs3 = ALOG(pbar (i,j,3) + pprt (i,j,3))
        tmp2 = ptbar(i,j,2) + ptprt(i,j,2)
        prs2 = ALOG(pbar (i,j,2) + pprt (i,j,2))

        w2=1.+(zsc(i,j,2)-hterain(i,j))/(zsc(i,j,3)-zsc(i,j,2))
        w3=1.-w2

        ptk=w2*tmp2 + w3*tmp3
        pres=EXP((w2*prs2 + w3*prs3))

        tmpk = ptk * (p0inv*pres)**rddcp

        IF (nbwater > 0) CALL xytoll(1,1,xll+dx*(i+0.5),yll+dy*(j+0.5),lat,lon)
!
        DO is=1,nstyp
          IF(soiltyp(i,j,is) /= 0) THEN
            IF(  soiltyp(i,j,is) == 13 ) THEN
              DO k=1,nzsoil 
! J.Case, ENSCO (9/22/2004)
                IF (sst_flag) THEN                
! J.Case, ENSCO (2/11/2005) -- FIT SST product
                  IF (goes_sst(i,j) > 0.0) THEN
                    print *,'Using GOES SST at grid point I/J = ',i,j
                    print *,'Original = ',tsoil(i,j,k,is)
                    print *,'GOES SST = ',goes_sst(i,j)
                    print *,'External = ',sltk_deep_interp(i,j)
                    tsoil (i,j,k,is) = goes_sst(i,j)
                  ELSE
                    print *,'Using external model SST at grid point I/J = ',i,j
                    print *,'Original = ',tsoil(i,j,k,is)
                    print *,'External = ',sltk_deep_interp(i,j)
                    tsoil (i,j,k,is) = sltk_deep_interp(i,j)
                  END IF
                ELSE
                  tsoil (i,j,k,is) = twater + 273.15
                END IF
              END DO 
              DO ib = 1,nbwater
                IF ((lat >= blat1(ib)) .AND. (lat <= blat2(ib))         &
                      .AND. (lon >= blon1(ib)) .AND. (lon <= blon2(ib))) THEN
                  DO k=1,nzsoil 
                    tsoil  (i,j,k,is) = tbwater(ib) + 273.15
                  END DO 
                END IF
              END DO
            ELSE
!TODO:  Replace tsrat and t2rat with array of values for each soil level.
              tsoil (i,j,1,is) = tmpk + ttprt
              DO k=2,nzsoil 
                tsoil  (i,j,k,is) = tmpk + tbprt
              END DO 
            END IF

          END IF
        END DO
      END DO
    END DO

  END IF

  IF ( snowin /= 1 ) THEN
    DO j=1,ny
      DO i=1,nx
        snowdpth(i,j)=snowdpth0
      END DO
    END DO
  END IF

!
!-----------------------------------------------------------------------
!
!  Modify or initialize soil properties using soilinit2 option
!
!-----------------------------------------------------------------------
!
  IF( soilinit2 == 1 ) THEN
!
    DO is=1,nstyp
      DO j=1,ny-1
        DO i=1,nx-1
!
!TODO:  Replace wgrat and w2rat with array of values for each soil level.
          IF( soiltyp(i,j,is) /= 0 ) THEN
            qsoil (i,j,1,is) = wwlt(soiltyp(i,j,is)) + wgrat             &
                 *( wsat(soiltyp(i,j,is)) - wwlt(soiltyp(i,j,is)) )
            DO k=1,nzsoil 
              qsoil (i,j,k,is) = wwlt(soiltyp(i,j,is)) + w2rat           &
                 *( wsat(soiltyp(i,j,is)) - wwlt(soiltyp(i,j,is)) )
            END DO 
            wetcanp(i,j,is) = wetcanp0
          END IF

        END DO
      END DO
    END DO

!
  ELSE IF( soilinit2 == 2 .OR. soilinit2 == 3 ) THEN
!
!-----------------------------------------------------------------------
!
!  Read the initial API data and find its location in the grid
!
!-----------------------------------------------------------------------
!
    totapi0=0
    OPEN(UNIT=3,ERR=205,FILE=apiinit,STATUS='unknown')
!
    out=0
    WRITE(6,*)
    DO istn=1,mxstns
      jstn=istn-out
      READ(3,*,END=201) lat,lon,api0(jstn)
      CALL lltoxy(1,1,lat,lon,xpt,ypt)
      xpt=xpt-xll
      ypt=ypt-yll
      CALL findlc(nx,ny,xs,ys,xpt,ypt,ipt,jpt,ireturn)
      IF( ireturn == 0 ) THEN
        iptapi0(jstn)=ipt
        jptapi0(jstn)=jpt
        WRITE(6,*) 'API obs:',jstn,' found near ',ipt,jpt,              &
            ' with initial API of ',api0(jstn)
        WRITE(6,*) ' --> main soil type:',soiltyp(ipt,jpt,1),           &
                       ' wilting:',wwlt(soiltyp(ipt,jpt,1)),            &
                       ' saturation:',wsat(soiltyp(ipt,jpt,1))
        WRITE(6,*)
      ELSE
        WRITE(6,*)'Initial API station #',jstn,                         &
                 ' is outside the grid..'
        WRITE(6,*)
        out=out+1
      END IF
    END DO
    201   CONTINUE
    totapi0=istn-1-out
    CLOSE(3)
    205   CONTINUE
!
!-----------------------------------------------------------------------
!
!  Read the observed precip data and find its location in the grid
!
!-----------------------------------------------------------------------
!
    IF( prdata == 1 ) THEN

      out=0
      WRITE(6,*)
      OPEN(UNIT=1,ERR=301,FILE=apifile,STATUS='old')
      DO istn=1,mxstns
        jstn=istn-out
        READ(1,1010,END=301) precfile(jstn)
        1010       FORMAT(a80)
        length=LEN(precfile(jstn))
        CALL strlnth(precfile(jstn),length)
        WRITE(6,*) 'File:',precfile(jstn)(1:length)
        OPEN(UNIT=2,FILE=precfile(jstn),STATUS='old')
        READ(2,*) lat,lon
        CALL lltoxy(1,1,lat,lon,xpt,ypt)
        xpt=xpt-xll
        ypt=ypt-yll
        CALL findlc(nx,ny,xs,ys,xpt,ypt,ipt,jpt,ireturn)
        IF( ireturn == 0 ) THEN
!
          iptstn(jstn)=ipt
          jptstn(jstn)=jpt
          totpreca(jstn)=0.0
!
          WRITE(6,*)
          WRITE(6,*) 'Precip stn # ',jstn,' was found near i,j:',       &
              ipt,jpt
          IF( jpt >= nint(ny*0.5) ) THEN
            WRITE(6,*) 'North of Center'
          ELSE
            WRITE(6,*) 'South of Center'
          END IF
          IF( ipt >= nint(nx*0.5)) THEN
            WRITE(6,*) 'East of Center'
          ELSE
            WRITE(6,*) 'West of Center'
          END IF
!
          WRITE(6,*)
          WRITE(6,*)
          DO iday=1,ndays
            READ(2,*) obsprec(iday,jstn) ! in inches...
            obsprec(iday,jstn)=obsprec(iday,jstn)*2.54
!
!  Assume (arbitrarily) that the ground can only absorb 4.0 cm of
!  rainfall per day, with the rest being runoff.
!
            obsprec(iday,jstn)=AMIN1(obsprec(iday,jstn),4.0 )
          END DO
          WRITE(6,*)
        ELSE
          PRINT *,                                                      &
              'Precip station #',jstn,' is outside of the grid..'
          PRINT *
          out=out+1
        END IF
        CLOSE(2)
      END DO
      301     CONTINUE
      CLOSE(1)
      CLOSE(2)
      tstns=istn-1-out

    ELSE IF (prdata == 2) THEN

! J.Case, ENSCO Inc. (9/22/2004)

      CALL rdnceppr(nx,ny,mxstns,mxdays,xs,ys,xll,yll,                  &
                      fnlyr,fnlmo,fnldy,ndays,                          &
                      prcpdir,prcplst,                                  &
                      staid,iptstn,jptstn,obsprec,totpreca,             &
                      tstns)

    ELSE

      WRITE(6,'(i6,a)') prdata,' is not a valid prdata option!'
      STOP


    END IF
!
!---------------------------------------------------------------------
!
!  Analyze the initial root layer API to the grid.  Determine the
!  range parameter using respapi.
!
!---------------------------------------------------------------------
!
    i1=1
    i2=nint(nx*.25)
    i3=nint(nx*.50)
    i4=nint(nx*.75)
    i5=nx-1
!
    j1=1
    j2=nint(ny*.25)
    j3=nint(ny*.50)
    j4=nint(ny*.75)
    j5=ny-1
!
    WRITE(6,*)
    WRITE(6,*) 'Analyzing initial API to grid...'
    WRITE(6,*)
!
    IF( totapi0 > 1 ) THEN
!
      totdist=SQRT((nx*dx*0.001)*(ny*dy*0.001))*                        &
          (1.+SQRT(FLOAT(totapi0)))/FLOAT(totapi0-1) ! using separation that
!        would result if obs were randomly scattered
      apirng=                                                           &
          SQRT(-((2.*totdist/3.14159)**2)*LOG(respapi))
      apir2=apirng*apirng
      const=1.0E-06/apir2
!
      WRITE(6,*)
      WRITE(6,*) 'Using range:',apirng,' km, for initializing API'
      WRITE(6,*)

      DO iy=1,ny-1
        DO ix=1,nx-1
          sumwt=0.
          sumapi=0.
          DO istn=1,totapi0
            ddx=xs(iptapi0(istn))-xs(ix)
            ddy=ys(jptapi0(istn))-ys(iy)
            dist=(ddx*ddx+ddy*ddy)
            wgt=EXP(-dist*const)
            sumwt=sumwt+wgt
            sumapi=sumapi+wgt*api0(istn)
          END DO
          IF( sumwt > 0. ) sumapi=sumapi/sumwt

          DO is=1,nstyp
            IF( soiltyp(ix,iy,is) /= 0 ) initapi(ix,iy,is)=             &
                AMIN1(sumapi,(wsat(soiltyp(ix,iy,is))*d2*100.))
          END DO

        END DO
      END DO

    ELSE ! if no initial values available, set to 0.5*wsat

      DO is=1,nstyp
        DO iy=1,ny-1
          DO ix=1,nx-1
            IF( soiltyp(ix,iy,is) /= 0)                                 &
                initapi(ix,iy,is)=0.5*wsat(soiltyp(ix,iy,is))*d2*100.0
          END DO
        END DO
      END DO

    END IF
!
    WRITE(6,*)
    WRITE(6,*) 'Initial API analysis:'
    WRITE(6,*) '---------------------'
    WRITE(6,*)
!
    WRITE(6,1040) j5,initapi(i1,j5,1),initapi(i2,j5,1),                 &
                     initapi(i3,j5,1),initapi(i4,j5,1),                 &
                     initapi(i5,j5,1)

    WRITE(6,1040) j4,initapi(i1,j4,1),initapi(i2,j4,1),                 &
                     initapi(i3,j4,1),initapi(i4,j4,1),                 &
                     initapi(i5,j4,1)

    WRITE(6,1040) j3,initapi(i1,j3,1),initapi(i2,j3,1),                 &
                     initapi(i3,j3,1),initapi(i4,j3,1),                 &
                     initapi(i5,j3,1)

    WRITE(6,1040) j2,initapi(i1,j2,1),initapi(i2,j2,1),                 &
                     initapi(i3,j2,1),initapi(i4,j2,1),                 &
                     initapi(i5,j2,1)

    WRITE(6,1040) j1,initapi(i1,j1,1),initapi(i2,j1,1),                 &
                     initapi(i3,j1,1),initapi(i4,j1,1),                 &
                     initapi(i5,j1,1)

    WRITE(6,*)
    WRITE(6,1035) i1,i2,i3,i4,i5
    WRITE(6,*)
!
    minday=15
    depth1=0.01
!
    DO iy=1,ny-1
      DO ix=1,nx-1
        kk2dep(ix,iy)=1.0
      END DO
    END DO

    DO is=1,nstyp
      DO iy=1,ny-1
        DO ix=1,nx-1
          api1(ix,iy,is)=initapi(ix,iy,is)*depth1/d2
          api2(ix,iy,is)=initapi(ix,iy,is)
        END DO
      END DO
    END DO
!
    maxprec=0.0
    DO iday=1,ndays
      precip=0.0
      DO istn=1,tstns
        IF( obsprec(iday,istn) >= 0. ) precip=precip+obsprec(iday,istn)
      END DO
      precip=precip/tstns
      IF( precip > maxprec ) THEN
        maxprec=precip
        maxday=iday
      END IF
    END DO
!
!----------------------------------------------------------------------
!
!  If the range is not specified in the input file (i.e. = 0.0)
!  determine the own range based on the input desired response.
!
!----------------------------------------------------------------------
!
    totdist=0.0
    DO ii=1,tstns
      mindist=1.0E30
      DO ij=1,tstns
        IF( ii /= ij ) THEN
          ddx=(xs(iptstn(ii))-xs(iptstn(ij)))
          ddy=(ys(jptstn(ii))-ys(jptstn(ij)))
          mindist=AMIN1(mindist,(ddx*ddx+ddy*ddy))
        END IF
      END DO
      totdist=totdist+0.001*SQRT(mindist)
    END DO
!
    WRITE(6,*) 'Average station separation (km):',totdist/tstns
    WRITE(6,*)
!
    IF( range == 0. ) THEN
      DO ii=1,1000000
        d_0=.000001*ii
        d_1=d_0**(gamma(2)**1.0)
        d_1star=d_0+(1-d_0)*d_1
        IF( d_1star >= respprec ) EXIT
      END DO
      396     CONTINUE
!
      range=SQRT(-((2*(totdist/(tstns*1.0))/3.14159)**2)*LOG(d_0))
      WRITE(6,*)
      WRITE(6,*) 'Range not specified, using range:',range
      WRITE(6,1026) ' -- 1st pass response at 2 x (stn sep):',d_0
      WRITE(6,1026) ' -- 2nd pass response at 2 x (stn sep):',          &
           d_1star
      1026     FORMAT(a40,f10.8)
      WRITE(6,*)
    END IF
!
    r2=range*range
!
!----------------------------------------------------------------------
!
!  Begin main day loop
!  Objectively analyze the precipitation data on 'ndays' different
!  days using a 2-pass analysis.
!
!----------------------------------------------------------------------
!
    DO iday=1,ndays
!
      DO iy=1,ny-1
        DO ix=1,nx-1
          prec(ix,iy)=0.
          totwt(ix,iy)=0.
        END DO
      END DO
      itime=initi+(iday-1)*sec24hr
      CALL abss2ctim( itime,iyr,imon,idy,ihr,imin,isec )
!
!    1st pass....
!
      const=1.0E-06/(r2*gamma(1))
      thrdst=1.0E06*(6.*totdist/tstns)*(6.*totdist/tstns)
      idel=INT(SQRT(thrdst)/dx)+1
      jdel=INT(SQRT(thrdst)/dy)+1
!
      DO istn=1,tstns
        IF( obsprec(iday,istn) >= 0. ) THEN
          istr =MAX(1,(iptstn(istn)-idel))
          iend =MIN((nx-1),(iptstn(istn)+idel))
          jstr =MAX(1,(jptstn(istn)-jdel))
          jend =MIN((ny-1),(jptstn(istn)+jdel))
          DO iy=jstr,jend
            DO ix=istr,iend
              ddx=(xs(iptstn(istn))-xs(ix))
              ddy=(ys(jptstn(istn))-ys(iy))
              dist=(ddx*ddx+ddy*ddy)
              IF( dist < thrdst ) THEN
                wgt=EXP(-dist*const)
                prec(ix,iy)=prec(ix,iy)+wgt*obsprec(iday,istn)
                totwt(ix,iy)=totwt(ix,iy)+wgt
              END IF
            END DO
          END DO
        END IF
      END DO

      DO iy=1,ny-1
        DO ix=1,nx-1
          IF( totwt(ix,iy) > 0. ) prec(ix,iy)=prec(ix,iy)/totwt(ix,iy)
        END DO
      END DO
!
      IF( iday == maxday ) THEN
!
        WRITE(6,'(a,i4,a,i2.2,a,i2.2)')                                 &
                'Example precip analysis on ',                          &
                iyr,'-',imon,'-',idy
        WRITE(6,'(a)') '-------------------------------------'
        WRITE(6,*)
!
        1035       FORMAT(7X,i3,5X,i3,5X,i3,5X,i3,5X,i3)
        1040       FORMAT(i3,3X,f5.2,3X,f5.2,3X,f5.2,3X,f5.2,3X,f5.2)
        WRITE(6,1040) j5,prec(i1,j5),prec(i2,j5),                       &
                         prec(i3,j5),prec(i4,j5),                       &
                         prec(i5,j5)

        WRITE(6,1040) j4,prec(i1,j4),prec(i2,j4),                       &
                         prec(i3,j4),prec(i4,j4),                       &
                         prec(i5,j4)

        WRITE(6,1040) j3,prec(i1,j3),prec(i2,j3),                       &
                         prec(i3,j3),prec(i4,j3),                       &
                         prec(i5,j3)

        WRITE(6,1040) j2,prec(i1,j2),prec(i2,j2),                       &
                         prec(i3,j2),prec(i4,j2),                       &
                         prec(i5,j2)

        WRITE(6,1040) j1,prec(i1,j1),prec(i2,j1),                       &
                         prec(i3,j1),prec(i4,j1),                       &
                         prec(i5,j1)
        WRITE(6,*)
        WRITE(6,1035) i1,i2,i3,i4,i5
        WRITE(6,*)
!
        rms=0.0
        bias=0.0
        total=0.0
        DO istn=1,tstns
          IF( obsprec(iday,istn) >= 0. ) THEN
            1045           FORMAT(a4,1X,f6.3,1X,a6,1X,f6.3)
            WRITE(6,1045) 'obs:',obsprec(iday,istn),' anal:',           &
                    prec(iptstn(istn),jptstn(istn))
            rms=rms+(obsprec(iday,istn)-                                &
                    prec(iptstn(istn),jptstn(istn)))**2
            bias=bias+obsprec(iday,istn)-                               &
                prec(iptstn(istn),jptstn(istn))
            total=total+1.
          END IF
        END DO
        IF( total > 0. ) THEN
          rms=SQRT(rms/total)
          bias=-bias/total
          WRITE(6,*)
          WRITE(6,*) 'RMS Error:',rms,' Bias:',bias
          WRITE(6,*)
        END IF

      END IF
!
!  2nd pass....
!
      DO istn=1,tstns
        IF( obsprec(iday,istn) >= 0. ) THEN
          difprec(iday,istn)=obsprec(iday,istn)-                        &
                      prec(iptstn(istn),jptstn(istn))
        ELSE
          difprec(iday,istn)=-900.
        END IF
      END DO

      DO iy=1,ny-1
        DO ix=1,nx-1
          dprec(ix,iy)=0.
          totwt(ix,iy)=0.
        END DO
      END DO
!
      const=1.0E-06/(r2*gamma(2))
      thrdst=1.0E06*(2.5*totdist/tstns)*(2.5*totdist/tstns)
      idel=INT(SQRT(thrdst)/dx)+1
      jdel=INT(SQRT(thrdst)/dy)+1
      DO istn=1,tstns
        IF( difprec(iday,istn) > -100. ) THEN
          istr =MAX(1,(iptstn(istn)-idel))
          iend =MIN((nx-1),(iptstn(istn)+idel))
          jstr =MAX(1,(jptstn(istn)-jdel))
          jend =MIN((ny-1),(jptstn(istn)+jdel))
          DO iy=jstr,jend
            DO ix=istr,iend
              ddx=xs(iptstn(istn))-xs(ix)
              ddy=ys(jptstn(istn))-ys(iy)
              dist=(ddx*ddx+ddy*ddy)
              IF( dist < thrdst ) THEN
                wgt=EXP(-dist*const)
                dprec(ix,iy)=dprec(ix,iy)+wgt*difprec(iday,istn)
                totwt(ix,iy)=totwt(ix,iy)+wgt
              END IF
            END DO
          END DO
        END IF
      END DO
!
!
!  Sum difference and impose limits on daily precip values...
!
      DO iy=1,ny-1
        DO ix=1,nx-1
          IF( totwt(ix,iy) > 0. )                                       &
              prec(ix,iy)=prec(ix,iy)+(dprec(ix,iy)/totwt(ix,iy))
          prec(ix,iy)=AMAX1(prec(ix,iy),0.0)
          prec(ix,iy)=AMIN1(prec(ix,iy),4.0)
        END DO
      END DO
!
      IF( iday == maxday) THEN
!
        itime=initi+(iday-1)*sec24hr
        CALL abss2ctim( itime,iyr,imon,idy,ihr,imin,isec )
        WRITE(6,'(a,i4,a,i2.2,a,i2.2)')                                 &
                'Example precip analysis on ',                          &
                iyr,'-',imon,'-',idy
        WRITE(6,*) '-------------------------------------'
        WRITE(6,*)
!
        WRITE(6,1040) j5,prec(i1,j5),prec(i2,j5),                       &
                         prec(i3,j5),prec(i4,j5),                       &
                         prec(i5,j5)

        WRITE(6,1040) j4,prec(i1,j4),prec(i2,j4),                       &
                         prec(i3,j4),prec(i4,j4),                       &
                         prec(i5,j4)

        WRITE(6,1040) j3,prec(i1,j3),prec(i2,j3),                       &
                         prec(i3,j3),prec(i4,j3),                       &
                         prec(i5,j3)

        WRITE(6,1040) j2,prec(i1,j2),prec(i2,j2),                       &
                         prec(i3,j2),prec(i4,j2),                       &
                         prec(i5,j2)

        WRITE(6,1040) j1,prec(i1,j1),prec(i2,j1),                       &
                         prec(i3,j1),prec(i4,j1),                       &
                         prec(i5,j1)
        WRITE(6,*)
        WRITE(6,1035) i1,i2,i3,i4,i5
        WRITE(6,*)
        rms=0.0
        bias=0.0
        total=0.0
        DO istn=1,tstns
          IF(obsprec(iday,istn) >= 0.) THEN
            WRITE(6,1045) 'obs:',obsprec(iday,istn),' anal:',           &
                prec(iptstn(istn),jptstn(istn))
            rms=rms+(obsprec(iday,istn)-                                &
                prec(iptstn(istn),jptstn(istn)))**2
            bias=bias+obsprec(iday,istn)-                               &
                prec(iptstn(istn),jptstn(istn))
            total=total+1.
          END IF
        END DO
        IF( total > 0. ) THEN
          rms=SQRT(rms/total)
          bias=-bias/total
          WRITE(6,*)
          WRITE(6,*) 'RMS Error:',rms,' Bias:',bias
          WRITE(6,*)
        END IF
      END IF
!
!---------------------------------------------------------------------
!
!  Calculate the API - assume that API will be relatively
!  insensitive to the initial API after about 90 days.
!  Recall d2=100.0cm and depth1 (m) is set below.  Note that
!  depth1 is not necessarilly equal to d1 in soilcst.inc, which
!  is only a normalizing depth.  In Noilhan-Planton (1989) they
!  suggest that depth1 is only a few millimeters.
!
!  k1 and k2 are minimum depletion coefficients for the topsoil
!  and root layer.  kk1 and kk2 are the time-dependent depletion
!  coefficients where kk? is at a maximum on minday (for minimum
!  moisture depletion - set to January 15, or Julian day 15) and
!  kk? is at a minimum 6 months later (July 15)
!
!  The time-dependent depletion coefficients can THEN be modified
!  up or down according to the apparent volumetric soil water
!  content (i.e. the API) relative to wilting and saturation using
!  the adj1 and adj2 terms, the values of which can be derived
!  empirically using ARM data.  With the limited data available
!  at this time, the adj1 and adj2 terms were simply set to 1.0
!  due to inability to show value in adding this term.  However,
!  intuitively there should be some relationship.
!
!  The API depletion coefficient multiplies the difference
!  between the calculated API and the wilting point (or the
!  initial API, whichever is lower) instead of the multiplying
!  the prior API as in the common formulation.  The new method
!  was derived empirically and appeared to give better correlation
!  with observed data as well as a lower dependence on the initial
!  API using the depletion coefficient formula below.
!
!---------------------------------------------------------------------
!
!
      adj1=1.0
      adj2=1.0
      CALL julday( iyr, imon, idy, jlday )

      kk1=1.0-adj1*                                                     &
             (1.0-k1)*SIN((3.14159*(jlday-minday))                      &
             /365.0)
      kk2=1.0-adj2*                                                     &
             (1.0-k2)*SIN((3.14159*(jlday-minday))                      &
             /365.0)

      DO iy=1,ny-1
        DO ix=1,nx-1
          kk2dep(ix,iy)=kk2dep(ix,iy)*kk2
        END DO
      END DO

      DO is=1,nstyp
        DO iy=1,ny-1
          DO ix=1,nx-1
!
            IF(soiltyp(ix,iy,is) > 0) THEN

              IF( k1 > 0. ) THEN
                api1(ix,iy,is)=prec(ix,iy)+kk1*api1(ix,iy,is)
              ELSE
                api1(ix,iy,is)=0.0
              END IF
              api1(ix,iy,is)=                                           &
                  AMIN1(api1(ix,iy,is),                                 &
                      (wsat(soiltyp(ix,iy,is))*depth1*100.))
!
              api2(ix,iy,is)=prec(ix,iy)+                               &
                  kk2*(api2(ix,iy,is)-AMIN1(initapi(ix,iy,is),          &
                    wwlt(soiltyp(ix,iy,is))*(d2*100.)))+                &
                  AMIN1(initapi(ix,iy,is),                              &
                   wwlt(soiltyp(ix,iy,is))*(d2*100.))
              api2(ix,iy,is)=                                           &
                  AMIN1(api2(ix,iy,is),(wsat(soiltyp(ix,iy,is))*d2*100.))
!
!  Carry out 1/2 of the daily depletion on the final day (the
!  day of the run) without adding in a precip term.
!
              IF( iday == ndays) THEN
                IF( k1 > 0. )                                           &
                    api1(ix,iy,is)=(api1(ix,iy,is)+api1(ix,iy,is)*kk1)/2.0
                api2(ix,iy,is)=(api2(ix,iy,is)+api2(ix,iy,is)*kk2)/2.0
              END IF

            END IF
!
          END DO
        END DO
      END DO

      DO istn=1,tstns
!      write(6,1020) 'Stn #:',istn,'Day:',iyr,'-',imon,'-',idy,
!    :       'Obs. precip (cm):',obsprec(iday,istn),
!    :       'Anal. prec (cm):',prec(iptstn(istn),jptstn(istn))
        totpreca(istn)=totpreca(istn)+                                  &
                    prec(iptstn(istn),jptstn(istn))
      END DO
      1020   FORMAT(a6,i4,3X,a5,i4,a,i2.2,a,i2.2,3X,a17,f5.2,3X,a16,f5.2)

    END DO
!
    WRITE(6,*)
    WRITE(6,*) 'Top layer API analysis:'
    WRITE(6,*) '-----------------------'
    WRITE(6,*)
!
    WRITE(6,1040) j5,api1(i1,j5,1),api1(i2,j5,1),                       &
                     api1(i3,j5,1),api1(i4,j5,1),                       &
                     api1(i5,j5,1)

    WRITE(6,1040) j4,api1(i1,j4,1),api1(i2,j4,1),                       &
                     api1(i3,j4,1),api1(i4,j4,1),                       &
                     api1(i5,j4,1)

    WRITE(6,1040) j3,api1(i1,j3,1),api1(i2,j3,1),                       &
                     api1(i3,j3,1),api1(i4,j3,1),                       &
                     api1(i5,j3,1)

    WRITE(6,1040) j2,api1(i1,j2,1),api1(i2,j2,1),                       &
                     api1(i3,j2,1),api1(i4,j2,1),                       &
                     api1(i5,j2,1)

    WRITE(6,1040) j1,api1(i1,j1,1),api1(i2,j1,1),                       &
                     api1(i3,j1,1),api1(i4,j1,1),                       &
                     api1(i5,j1,1)

    WRITE(6,*)
    WRITE(6,1035) i1,i2,i3,i4,i5
    WRITE(6,*)
!
    WRITE(6,*)
    WRITE(6,*) 'Root zone API analysis:'
    WRITE(6,*) '-----------------------'
    WRITE(6,*)
!
    WRITE(6,1040) j5,api2(i1,j5,1),api2(i2,j5,1),                       &
                     api2(i3,j5,1),api2(i4,j5,1),                       &
                     api2(i5,j5,1)

    WRITE(6,1040) j4,api2(i1,j4,1),api2(i2,j4,1),                       &
                     api2(i3,j4,1),api2(i4,j4,1),                       &
                     api2(i5,j4,1)

    WRITE(6,1040) j3,api2(i1,j3,1),api2(i2,j3,1),                       &
                     api2(i3,j3,1),api2(i4,j3,1),                       &
                     api2(i5,j3,1)

    WRITE(6,1040) j2,api2(i1,j2,1),api2(i2,j2,1),                       &
                     api2(i3,j2,1),api2(i4,j2,1),                       &
                     api2(i5,j2,1)

    WRITE(6,1040) j1,api2(i1,j1,1),api2(i2,j1,1),                       &
                     api2(i3,j1,1),api2(i4,j1,1),                       &
                     api2(i5,j1,1)

    WRITE(6,*)
    WRITE(6,1035) i1,i2,i3,i4,i5
    WRITE(6,*)
!
    WRITE(6,*)
    WRITE(6,*) 'Dependence of root zone API on initial API: '
    WRITE(6,*) '(most useful when adj2 .ne. 1.0 -- see code)'
    WRITE(6,*) '--------------------------------------------'
    WRITE(6,*)
!
    WRITE(6,1040) j5,kk2dep(i1,j5),kk2dep(i2,j5),                       &
                     kk2dep(i3,j5),kk2dep(i4,j5),                       &
                     kk2dep(i5,j5)

    WRITE(6,1040) j4,kk2dep(i1,j4),kk2dep(i2,j4),                       &
                     kk2dep(i3,j4),kk2dep(i4,j4),                       &
                     kk2dep(i5,j4)

    WRITE(6,1040) j3,kk2dep(i1,j3),kk2dep(i2,j3),                       &
                     kk2dep(i3,j3),kk2dep(i4,j3),                       &
                     kk2dep(i5,j3)

    WRITE(6,1040) j2,kk2dep(i1,j2),kk2dep(i2,j2),                       &
                     kk2dep(i3,j2),kk2dep(i4,j2),                       &
                     kk2dep(i5,j2)

    WRITE(6,1040) j1,kk2dep(i1,j1),kk2dep(i2,j1),                       &
                     kk2dep(i3,j1),kk2dep(i4,j1),                       &
                     kk2dep(i5,j1)

    WRITE(6,*)
    WRITE(6,1035) i1,i2,i3,i4,i5
    WRITE(6,*)
!
!--------------------------------------------------------------------
!
!  Convert the API's to wetsfc and wetdp.  Limit both 'wetsfc' and
!  'wetdp' to ~'wsat' as a maximum.  Limit 'wetsfc' to 'wgeq' as
!  a minimum and limit 'wetdp' to ~'wwlt' as a minimum.  If k1=0.0
!  set 'wetsfc' equal to 'wgeq' and if k1<0.0 set 'wetsfc' ~ equal
!  to 0.8*'wgeq'.  See input file for more details.
!
!  Also, convert API from scalar points back to x,y grid.
!
!--------------------------------------------------------------------
!
    IF(soilinit2 == 3) THEN
      wsfcscl=wgrat
      wdpscl=w2rat
    ELSE
      wsfcscl=1.0
      wdpscl=1.0
    END IF
!
    DO is=1,nstyp
      DO iy=1,ny-1
        DO ix=1,nx-1
          IF( soiltyp(ix,iy,is) /= 0 ) THEN

            qsoil(ix,iy,1,is) =wsfcscl*api1(ix,iy,is)/(depth1*100.0)
            DO k = 2,nzsoil
              qsoil(ix,iy,k,is) =wdpscl*api2(ix,iy,is)/(d2*100.0)
            END DO

            wetcanp(ix,iy,is)=wetcanp0
!
            IF( soiltyp(ix,iy,is) /= 13 ) THEN
              wmin=0.95*wwlt(soiltyp(ix,iy,is)) ! arbitrary min,max values
              wmax=0.95*wsat(soiltyp(ix,iy,is)) ! for qsoil(nzsoil)...
            ELSE  ! if soiltyp is water
              wmin=wwlt(soiltyp(ix,iy,is))
              wmax=wsat(soiltyp(ix,iy,is))
            END IF
            DO k = 2,nzsoil
              IF( qsoil(ix,iy,k,is) < wmin ) THEN
                qsoil(ix,iy,k,is)=wmin ! keep qsoil(k) close to wwlt
              ELSE IF( qsoil(ix,iy,k,is) > wmax ) THEN
                qsoil(ix,iy,k,is)=wmax ! keep qsoil(k) slightly less than saturated
              END IF
            END DO
!            wmin=qsoil(ix,iy,nzsoil,is)/wsat(soiltyp(ix,iy,is))
            wmin=qsoil(ix,iy,2,is)/wsat(soiltyp(ix,iy,is))
            wmin=wmin - ( awgeq(soiltyp(ix,iy,is)) *                  &
              wmin**pwgeq(soiltyp(ix,iy,is)) )  *                     &
              ( 1.0 - wmin**(8.0*pwgeq(soiltyp(ix,iy,is))) )
            wmin=wmin*wsat(soiltyp(ix,iy,is))
            IF( qsoil(ix,iy,1,is) > wmax ) THEN
              qsoil(ix,iy,1,is)=wsat(soiltyp(ix,iy,is))
            ELSE IF( qsoil(ix,iy,1,is) < wmin ) THEN
              qsoil(ix,iy,1,is)=wmin
            END IF
!
            IF( k1 == 0. ) THEN
              qsoil(ix,iy,1,is)=wmin
            ELSE IF(  k1 < 0. .AND. soiltyp(ix,iy,is) /= 13) THEN
              qsoil(ix,iy,1,is)=wmin-AMIN1(0.2*wmin,0.05)
            END IF

          END IF

        END DO
      END DO
    END DO

    DO j=1,ny-1
      DO i=1,nx-1
        DO k=1,nzsoil 
          tsoil (i,j,k,0) =  0.0
          qsoil (i,j,k,0) =  0.0
        END DO 
        wetcanp(i,j,0) =  0.0
      END DO
    END DO

    DO is = 1,nstyp
      DO j=1,ny-1
        DO i=1,nx-1
          DO k=1,nzsoil 
            tsoil (i,j,k,0) = tsoil  (i,j,k,0)                           &
                         + tsoil (i,j,k,is) * stypfrct(i,j,is)
            qsoil (i,j,k,0) = qsoil (i,j,k,0)                            &
                         + qsoil (i,j,k,is) * stypfrct(i,j,is)

          END DO
 
          wetcanp(i,j,0) = wetcanp(i,j,0)                               &
                         + wetcanp(i,j,is) * stypfrct(i,j,is)
        END DO
      END DO
    END DO

    DO is=0,nstyp
      CALL edgfill(tsoil (1,1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,nzsoil, &
               1,nzsoil)
      CALL edgfill(qsoil (1,1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,nzsoil,1, &
               1,nzsoil)
      CALL edgfill(wetcanp(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
    END DO
!
!-------------------------------------------------------------------
!
!    Print out the calculated values of wetsfc, wetdp, and wetcanp.
!
!-------------------------------------------------------------------
!
    i5=nx
    j5=ny
!
    WRITE(6,*)
    WRITE(6,*) 'wetsfc analysis:'
    WRITE(6,*) '-----------------------'
    WRITE(6,*)
!
    WRITE(6,1040) j5,qsoil(i1,j5,1,0),qsoil(i2,j5,1,0),                 &
                 qsoil(i3,j5,1,0),qsoil(i4,j5,1,0),                     &
                 qsoil(i5,j5,1,0)

    WRITE(6,1040) j4,qsoil(i1,j4,1,0),qsoil(i2,j4,1,0),                 &
                 qsoil(i3,j4,1,0),qsoil(i4,j4,1,0),                     &
                 qsoil(i5,j4,1,0)

    WRITE(6,1040) j3,qsoil(i1,j3,1,0),qsoil(i2,j3,1,0),                 &
                 qsoil(i3,j3,1,0),qsoil(i4,j3,1,0),                     &
                 qsoil(i5,j3,1,0)

    WRITE(6,1040) j2,qsoil(i1,j2,1,0),qsoil(i2,j2,1,0),                 &
                 qsoil(i3,j2,1,0),qsoil(i4,j2,1,0),                     &
                 qsoil(i5,j2,1,0)

    WRITE(6,1040) j1,qsoil(i1,j1,1,0),qsoil(i2,j1,1,0),                 &
                 qsoil(i3,j1,1,0),qsoil(i4,j1,1,0),                     &
                 qsoil(i5,j1,1,0)

    WRITE(6,*)
    WRITE(6,1035) i1,i2,i3,i4,i5
    WRITE(6,*)
!
    WRITE(6,*)
    WRITE(6,*) 'wetdp analysis:'
    WRITE(6,*) '-----------------------'
    WRITE(6,*)
!
    WRITE(6,1040) j5,qsoil(i1,j5,nzsoil,0),qsoil(i2,j5,nzsoil,0),       &
                 qsoil(i3,j5,nzsoil,0),qsoil(i4,j5,nzsoil,0),           &
                 qsoil(i5,j5,nzsoil,0)

    WRITE(6,1040) j4,qsoil(i1,j4,nzsoil,0),qsoil(i2,j4,nzsoil,0),       &
                 qsoil(i3,j4,nzsoil,0),qsoil(i4,j4,nzsoil,0),           &
                 qsoil(i5,j4,nzsoil,0)

    WRITE(6,1040) j3,qsoil(i1,j3,nzsoil,0),qsoil(i2,j3,nzsoil,0),       &
                 qsoil(i3,j3,nzsoil,0),qsoil(i4,j3,nzsoil,0),           &
                 qsoil(i5,j3,nzsoil,0)

    WRITE(6,1040) j2,qsoil(i1,j2,nzsoil,0),qsoil(i2,j2,nzsoil,0),       &
                 qsoil(i3,j2,nzsoil,0),qsoil(i4,j2,nzsoil,0),           &
                 qsoil(i5,j2,nzsoil,0)

    WRITE(6,1040) j1,qsoil(i1,j1,nzsoil,0),qsoil(i2,j1,nzsoil,0),       &
                 qsoil(i3,j1,nzsoil,0),qsoil(i4,j1,nzsoil,0),           &
                 qsoil(i5,j1,nzsoil,0)

    WRITE(6,*)
    WRITE(6,1035) i1,i2,i3,i4,i5
    WRITE(6,*)
!
!  DO 720 istn=1,tstns
!    totpreco=0.0
!    DO 710 iday=1,ndays
!      IF(obsprec(iday,istn).gt.0.)
!    :    totpreco=totpreco+obsprec(iday,istn)
! 710   CONTINUE
!    write(6,*)
!    write(6,1021) 'Total recorded precip:',totpreco
!    write(6,1021) 'Total analyzed precip:',totpreca(istn)
!1021   format(a22,f6.2)
!    write(6,1022) 'API:',api2(iptstn(istn),jptstn(istn),1)
!1022   format(a4,18x,f6.2)
!    write(6,1023) 'Qsoilsfc:',qsoil(iptstn(istn),jptstn(istn),1,0)
!1023   format(a9,15x,f6.4)
!    write(6,1024) 'Qsoildp:',qsoil(iptstn(istn),jptstn(istn),nzsoil,0)
!1024   format(a8,16x,f6.4)
!    write(6,*)
!    write(6,*)
! 720 CONTINUE
!
    IF( totapi0 > 0 ) THEN
      WRITE(6,*)
      WRITE(6,*) 'Root zone volumetric water contents'                  &
          //' at initial API stations:'
      WRITE(6,*)
      DO i=1,totapi0
        WRITE(6,1042) 'Stn #:',i,' has wetdp of:',                      &
            qsoil(iptapi0(i),jptapi0(i),nzsoil,0)
        1042     FORMAT(a6,1X,i2,a15,1X,f5.3)
      END DO
      WRITE(6,*)
    END IF
!
    WRITE(6,*)
    WRITE(6,*) 'wetcanp analysis:'
    WRITE(6,*) '-----------------------'
    WRITE(6,*)
!
    WRITE(6,1040) j5,wetcanp(i1,j5,0),wetcanp(i2,j5,0),                 &
                 wetcanp(i3,j5,0),wetcanp(i4,j5,0),                     &
                 wetcanp(i5,j5,0)

    WRITE(6,1040) j4,wetcanp(i1,j4,0),wetcanp(i2,j4,0),                 &
                 wetcanp(i3,j4,0),wetcanp(i4,j4,0),                     &
                 wetcanp(i5,j4,0)

    WRITE(6,1040) j3,wetcanp(i1,j3,0),wetcanp(i2,j3,0),                 &
                 wetcanp(i3,j3,0),wetcanp(i4,j3,0),                     &
                 wetcanp(i5,j3,0)

    WRITE(6,1040) j2,wetcanp(i1,j2,0),wetcanp(i2,j2,0),                 &
                 wetcanp(i3,j2,0),wetcanp(i4,j2,0),                     &
                 wetcanp(i5,j2,0)

    WRITE(6,1040) j1,wetcanp(i1,j1,0),wetcanp(i2,j1,0),                 &
                 wetcanp(i3,j1,0),wetcanp(i4,j1,0),                     &
                 wetcanp(i5,j1,0)

    WRITE(6,*)
    WRITE(6,1035) i1,i2,i3,i4,i5
    WRITE(6,*)
!
    1041 FORMAT(i3,3X,f5.1,3X,f5.1,3X,f5.1,3X,f5.1,3X,f5.1)
!
    WRITE(6,*)
    WRITE(6,*) 'tsoil analysis:'
    WRITE(6,*) '-----------------------'
    WRITE(6,*)
!
    DO k=1,nzsoil

    WRITE(6,1041) j5,tsoil(i1,j5,k,0),tsoil(i2,j5,k,0),                   &
                 tsoil(i3,j5,k,0),tsoil(i4,j5,k,0),                       &
                 tsoil(i5,j5,k,0)

    WRITE(6,1041) j4,tsoil(i1,j4,k,0),tsoil(i2,j4,k,0),                   &
                 tsoil(i3,j4,k,0),tsoil(i4,j4,k,0),                       &
                 tsoil(i5,j4,k,0)

    WRITE(6,1041) j3,tsoil(i1,j3,k,0),tsoil(i2,j3,k,0),                   &
                 tsoil(i3,j3,k,0),tsoil(i4,j3,k,0),                       &
                 tsoil(i5,j3,k,0)

    WRITE(6,1041) j2,tsoil(i1,j2,k,0),tsoil(i2,j2,k,0),                   &
                 tsoil(i3,j2,k,0),tsoil(i4,j2,k,0),                       &
                 tsoil(i5,j2,k,0)

    WRITE(6,1041) j1,tsoil(i1,j1,k,0),tsoil(i2,j1,k,0),                   &
                 tsoil(i3,j1,k,0),tsoil(i4,j1,k,0),                       &
                 tsoil(i5,j1,k,0)
    END DO 

    WRITE(6,*)
    WRITE(6,1035) i1,i2,i3,i4,i5
    WRITE(6,*)
!
!--------------------------------------------------------------------
!
!    Read in soil moisture data from a file.
!
!--------------------------------------------------------------------
!
  ELSE IF( soilinit2 == 4 ) THEN

    WRITE(6,'(a,a)') 'Reading soil moisture from ',soilinfl
    CALL readsoil(nx,ny,nzsoil,nstyps,soilinfl,dx,dy,zpsoil,            &
             mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,       &
             zpsoilin,tsoilin,qsoilin,wcanpin,snowdin,          &
                  tem1,qsoil,wetcanp,snowdpth,soiltyp)

    DO is = 1,nstyp
      DO j=1,ny-1
        DO i=1,nx-1
          wetcanp(i,j,is) = wetcanp(i,j,0)

          DO k=1,nzsoil

          qsoil (i,j,k,is) = qsoil (i,j,k,0)
          tsoil (i,j,k,0) = tsoil  (i,j,k,0)                          &
                       + tsoil  (i,j,k,is) * stypfrct(i,j,is)
          END DO 

        END DO
      END DO
    END DO

!
!--------------------------------------------------------------------
!
!    Merge in soil data from soilinfl into soil data provided by history
!    dump where ever there is water, also use snowdpth in soilinfl
!    if present.
!
!--------------------------------------------------------------------
!
  ELSE IF( soilinit2 == 5 ) THEN

    WRITE(6,'(a,a)') 'Merging water soil data from ',soilinfl

    ! Note that nstyp=1 is assumed (or at least only the composite
    ! type is used) for data in soilinfl
    soiltyp2 = soiltyp
    CALL readsoil(nx,ny,nstyp,soilinfl,dx,dy,                          &
             mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon,       &
                  zpsoilin,tsoilin,qsoilin,wcanpin,snowdin,          &
                  tsoil2,qsoil2,wetcanp2,snowdpth,soiltyp2)

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
!   DO is=1,nstyp
!    DO k=1,nzsoil 
!     DO j=1,ny-1
!       DO i=1,nx-1
!         IF (soiltyp(i,j,is) == 13) THEN
!           IF (tsoilin /= 0)  tsoil(i,j,k,is)  = tsoil2(i,j,k,0)
!           IF (qsoilin /= 0)  qsoil(i,j,k,is)  = qsoil2(i,j,k,0)
!           IF (wcanpin /= 0)  wetcanp(i,j,is) = wetcanp2(i,j,0)
!         END IF
!       END DO
!      END DO 
!     END DO
!   END DO

    DO is=1,nstyp
      DO j=1,ny-1
        DO i=1,nx-1
          IF (soiltyp(i,j,is) == 13) THEN
            DO is2=1,nstyp
              IF (soiltyp2(i,j,is2) == 13) THEN
                DO k=1,nzsoil
                  IF (tsoilin /= 0) tsoil(i,j,k,is) = tsoil2(i,j,k,is2)
                  IF (qsoilin /= 0) qsoil(i,j,k,is) = qsoil2(i,j,k,is2)
                END DO 
                IF (wcanpin /= 0)  wetcanp(i,j,is) = wetcanp2(i,j,is2)
                EXIT
              END IF
            END DO
          END IF
        END DO
      END DO
    END DO

        ! is=0 filled in when program returns to main program

!  ELSE IF( soilinit2 == 6 ) THEN
!
!    DO is=1,nstyps
!    DO j=1,ny-1
!    DO i=1,nx-1
!
!      IF( soiltyp(i,j,is) /= 0 ) THEN
!        IF ( soiltyp(i,j,is).eq.13 ) THEN
!          wetsfc (i,j,is) = 1.0
!          wetdp  (i,j,is) = 1.0
!          wetcanp(i,j,is) = 0.0
!        ELSE
!         wetsfc (i,j,is) = soim1_interp(i,j) + 0.05  ! RUC offset
!!          wetsfc (i,j,is) = soim1_interp(i,j)
!          wmin=0.95*(wwlt(soiltyp(i,j,is)))
!          wmax=0.95*(wsat(soiltyp(i,j,is)))
!                                                                                                                                   
!          IF (wetsfc(i,j,is).gt.wmax) THEN
!      print *,'Resetting wetsfc to 0.95*wsat at I/J/IS= ',i,j,is
!      print *,'WETSFC/0.95*WSAT = ',wetsfc(i,j,is),wmax
!            wetsfc(i,j,is)= wmax
!          ENDIF
!          IF (wetsfc(i,j,is).lt.(wmin-amin1(0.2*wmin,0.05))) THEN
!      print *,'Resetting wetsfc to wmin at I/J/IS= ',i,j,is
!      print *,'WETSFC/WMIN = ',wetsfc(i,j,is),(wmin-amin1(0.2*wmin,0.05))
!            wetsfc(i,j,is)= wmin-amin1(0.2*wmin,0.05)
!          ENDIF
!                                                                                                                                   
!          wetdp  (i,j,is) = soim100_interp(i,j) + 0.05  ! RUC offset
!!          wetdp  (i,j,is) = soim100_interp(i,j)
!                                                                                                                                   
!          IF (wetdp(i,j,is).gt.wmax) THEN
!      print *,'Resetting wetdp to 0.95*wsat at I/J/IS= ',i,j,is
!      print *,'WETDP/0.95*WSAT = ',wetdp(i,j,is),wmax
!            wetdp(i,j,is)= wmax
!          ENDIF
!          IF (wetdp(i,j,is).lt.wmin) THEN
!      print *,'Resetting wetdp to wmin at I/J/IS= ',i,j,is
!      print *,'WETDP/WMIN = ',wetdp(i,j,is),wmin
!            wetdp(i,j,is)= wmin
!          ENDIF
!          wetcanp(i,j,is) = wetcanp0
!        ENDIF
!      ENDIF
!
!    ENDDO
!    ENDDO
!    ENDDO
  
  ELSE

    WRITE(6,'(i6,a)') soilinit2,' is not a valid soilinit2 option!'
    STOP

  END IF

  !wdt update
  IF (soilinit2 == 5) THEN
    DEALLOCATE(tsoil2)
    DEALLOCATE(qsoil2)
    DEALLOCATE(wetcanp2)
    DEALLOCATE(soiltyp2)
  END IF
!
!-------------------------------------------------------------------
!
!    End of subroutine
!
!-------------------------------------------------------------------
!
  PRINT *,'Finishes mksoilvar OK!!'
!
  RETURN
!
END SUBROUTINE mksoilvar
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE FINDLC                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE findlc(nx,ny,xs,ys,xpt,ypt,ipt,jpt,ireturn) 10
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Searches in x and y to find i,j location of xpt, ypt.
!
!  X and Y do not have to be on a regular grid, however it is
!  assumed that x and y are monotonically increasing as i and j
!  indices increase.
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster
!  April 1992.
!
!  MODIFICATION HISTORY:
!
!  February, 1993 (K. Brewster)
!  Additional documentation for ARPS 3.1 release
!
!  October, 1994 (K. Brewster)
!  Changed to reference scalar points.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    xs       x coordinate of scalar points in physical/comp. space (m)
!    ys       y coordinate of scalar points in physical/comp. space (m)
!
!    xpt      location to find in x coordinate (m)
!    ypt      location to find in y coordinate (m)
!
!  OUTPUT:
!
!    ipt      i index to the west of desired location
!    jpt      j index to the south of desired location
!
!-----------------------------------------------------------------------
!
!  Arguments
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny          ! Dimensions of ARPS grids
  REAL :: xs(nx)            ! x coordinate of scalar grid points in
                            ! physical/comp. space (m)
  REAL :: ys(ny)            ! y coordinate of grid points in
                            ! physical/comp. space (m)

  REAL :: xpt               ! location to find in x coordinate
  REAL :: ypt               ! location to find in y coordinate
  INTEGER :: ipt            ! i index to the west of desired
                            ! location
  INTEGER :: jpt            ! j index to the south of desired
                            ! location
  INTEGER :: ireturn
!
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ireturn=0
!
  DO i=2,nx
    IF(xpt < xs(i)) EXIT
  END DO
  101 CONTINUE
  ipt=i-1
  IF( xpt > xs(nx-1) .OR. xpt < xs(1) ) ireturn=-1
  DO j=2,ny
    IF( ypt < ys(j)) EXIT
  END DO
  201 CONTINUE
  jpt=j-1
  IF( ypt > ys(ny-1) .OR. ypt < ys(1) ) ireturn=-2

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


SUBROUTINE rdnceppr(nx,ny,mxstns,mxdays,xs,ys,xll,yll,                  & 1,5
           iendyr,iendmo,ienddy,ndays,                                  &
           prcpdir,prcplst,                                             &
           staid,iptstn,jptstn,obsprec,totpreca,                        &
           tstns)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Program to use history data dump and surface property data to
!  create the soil variables
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster
!    04/11/1995
!
!  MODIFICATION HISTORY:
!
!  05/07/1998 (Yuhe Liu)
!  Hub-CAPS's modification to the precipitation format change in
!  extracting the information from each line.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Arguments
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx,ny,mxstns,mxdays
!
  REAL :: xs(nx)
  REAL :: ys(ny)
  REAL :: xll
  REAL :: yll
  INTEGER :: iendyr,iendmo,ienddy,ndays
  CHARACTER (LEN=80) :: prcpdir
  CHARACTER (LEN=80) :: prcplst
  CHARACTER (LEN=7) :: staid(mxstns)
  INTEGER :: iptstn(mxstns)
  INTEGER :: jptstn(mxstns)
  REAL :: obsprec(mxdays,mxstns)
  REAL :: totpreca(mxstns)
  INTEGER :: tstns
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: mxstafile
  PARAMETER( mxstafile = 99000)
  INTEGER :: sec24hr
  PARAMETER( sec24hr = 86400 )

  REAL :: in2cm,maxprcm
! J.Case, ENSCO Inc. (10/26/2004)
!  PARAMETER (in2cm=2.54,maxprcm=4.00)
  PARAMETER (in2cm=2.54,maxprcm=8.00)
!
  CHARACTER (LEN=7) :: rdstaid
  CHARACTER (LEN=7) :: rdst(3)
  CHARACTER (LEN=7) :: dummy
  CHARACTER (LEN=80) :: rfcname
  CHARACTER (LEN=80) :: line
  INTEGER :: iunit,istn,jstn,kstn,iday,out
  INTEGER :: ilat(3),ilon(3)
  INTEGER :: ipt,jpt,ios,ireturn,lendir
  INTEGER :: istart,itime,iyr,imon,idy,ihr,imin,isec
  INTEGER :: nread,nmatch
  REAL :: lat,lon,xpt,ypt
  REAL :: rdlat,rdlon,rdobpr
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!-------------------------------------------------------------------
!
!  Initialize precipitation observations to missing
!
!-------------------------------------------------------------------
!
  DO istn=1,mxstns
    staid(istn)='     '
    DO iday=1,mxdays
      obsprec(iday,istn)=-9.0
    END DO
  END DO
!
!-------------------------------------------------------------------
!
!  Get list of precip stations
!  and keep those that are within the domain.
!
!-------------------------------------------------------------------
!
  iunit=38
  OPEN(iunit,ERR=910,FILE=prcplst,IOSTAT=ios,                           &
             STATUS='old',FORM='formatted')
! J.Case, ENSCO Inc. (10/26/2004)
!  READ(iunit,'(a5)') dummy
  out=0
  nread=0

! J.Case -- commented out read statement -- rewrote.
! J.Case, ENSCO Inc. (10/26/2004)

!  DO istn=1,mxstns,3
!    READ(iunit,'(a7,i5,i6,8x,a7,i5,i6,8x,a7,i5,i6)',                    &
!              ERR=101,END=101)                                          &
!              rdst(1),ilat(1),ilon(1),                                  &
!              rdst(2),ilat(2),ilon(2),                                  &
!              rdst(3),ilat(3),ilon(3)
!    nread=nread+3
!    DO kstn=1,3
!      jstn=istn-out+(kstn-1)
!      lat = 0.01*FLOAT(ilat(kstn))
!      lon =-0.01*FLOAT(ilon(kstn))
!

  DO istn=1,mxstns
    READ(iunit,'(1x,f5.2,3x,f5.2,8x,a7)',ERR=101,END=101)               &
              lat,lon,rdst(1)                                           
    nread=nread+1
    DO kstn=1,1
      jstn=istn-out
      lon = -lon
      CALL lltoxy(1,1,lat,lon,xpt,ypt)
      xpt=xpt-xll
      ypt=ypt-yll
      CALL findlc(nx,ny,xs,ys,xpt,ypt,ipt,jpt,ireturn)
      IF( ireturn == 0 ) THEN

        staid(jstn)=rdst(kstn)
        iptstn(jstn)=ipt
        jptstn(jstn)=jpt
        totpreca(jstn)=0.0
        WRITE(6,*) ' Precip stn ',rdst(kstn),' was found near i,j:',    &
                   ipt,jpt

      ELSE

        write(6,*) ' Precip stn ',rdst(kstn),' is outside of the grid.'
        out=out+1

      END IF
    END DO
  END DO

  101 CONTINUE
  tstns=istn-out-1
  CLOSE(iunit)
  WRITE(6,'(a,i6,a)') ' Read ',nread,' precip station locations.'
  WRITE(6,'(i12,a)') tstns,' inside grid.'
!
!-------------------------------------------------------------------
!
!  Begin main day loop, reading precip file for each day.
!
!-------------------------------------------------------------------
!
  lendir = LEN( prcpdir )
  CALL strlnth( prcpdir, lendir )

  CALL ctim2abss(iendyr,iendmo,ienddy,12,0,0, istart)
  istart=istart-((ndays-1)*sec24hr)
!
  DO iday=1,ndays
    itime = istart + ((iday-1)*sec24hr)
    CALL abss2ctim( itime, iyr, imon, idy, ihr, imin, isec )
    WRITE(rfcname,'(a,a,i4.4,i2.2,i2.2)')                               &
        prcpdir(1:lendir),'/usa-dlyprcp-',iyr,imon,idy
    WRITE(6,'(/a,a)') ' Opening: ',rfcname
! J.Case, ENSCO Inc. (10/26/2004) -- Fixed bug that exited out of loop entirely if 
! one day of data is missing.
    OPEN(iunit,ERR=351,FILE=rfcname,                                    &
         STATUS='old',FORM='formatted')
    nread=0
    nmatch=0
    READ(iunit,'(a5)',ERR=351,END=351) dummy
    DO kstn=1,mxstafile
      READ(iunit,'(a80)',ERR=350,END=351) line
      READ(line,*,ERR=350,END=350) rdlat,rdlon,rdobpr
      IF (line(23:23) == ' ') THEN
        rdstaid = line(24:31)
      ELSE
        rdstaid = line(23:30)
      END IF

      nread=nread+1
      DO istn=1,tstns
        IF( rdstaid == staid(istn) ) THEN
          nmatch=nmatch+1
          obsprec(iday,istn)=in2cm*rdobpr
          obsprec(iday,istn)=AMIN1(maxprcm,obsprec(iday,istn))
          EXIT
        END IF
      END DO
      301     CONTINUE
    END DO
    350 CONTINUE
    351 CONTINUE
    CLOSE(iunit)
    WRITE(6,'(a,i6,a)') ' Read ',nread,' precip obs.'
    WRITE(6,'(i12,a)') nmatch,' obs matched inside grid.'
  END DO
!
!-------------------------------------------------------------------
!
!    Normal end
!
!-------------------------------------------------------------------
!
  RETURN
!
!-------------------------------------------------------------------
!
!    Error destinations
!
!-------------------------------------------------------------------
!
910 CONTINUE
  WRITE(6,'(a,a,/a,i6)')                                                &
  'Error opening list of precip stations,',prcplst,' iostat = ',ios
  STOP

END SUBROUTINE rdnceppr

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!

SUBROUTINE rdextsst(nx,ny,infilename,                                   & 2
                     grid,blats,blons,res,nisst,njsst,sst_flag)
!
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

  INTEGER :: nx,ny
  REAL :: grid(nx,ny), blats, blons, res

  INTEGER :: nisst, njsst
  LOGICAL :: sst_flag

  CHARACTER (LEN=80) dummy
  CHARACTER (LEN=10) infilename

  blats = 0.
  blons = -180.
  res   = 1.

  OPEN (unit=50,file=infilename,form='FORMATTED',                       &
        status='OLD',iostat=ios,err=100)

  100 sst_flag = .true.
!  print *,ios,' ',infilename,nx,ny
  print *,'INFILE = ',infilename
  IF (ios .ne. 0) THEN
    sst_flag = .false.
    RETURN
  END IF
!
  read (50,'(//////,71x,i4,1x,i4)') nisst,njsst
  print *,'Total grid dimensions = ',nisst,njsst
  read (50,'(10x,2I5,11x,2I5)') ix1, ix2, iy1, iy2
  print *,'Subset domain = ',ix1,ix2,iy1,iy2
  read (50,'(///)')
  read (50,110) (kk,k = ix1, ix2 )
  110 format (9x,8(3X,I4,2X),/(8x,8(3X,I4,2X)))
!
!     Loop through rows from the top down.
!
  DO j = iy2, iy1, -1
    read (50,120) jj,(grid(k,j),k=ix1,ix2)
  END DO
  120 format (4x,I3,1X,8F9.4,/(8X,8F9.4))

  RETURN
END SUBROUTINE rdextsst