subroutine sfch2o(sm,stc,smc,ioneta,sh2o)

      implicit none
      include 'param_o.h'
      include 'const.h'

      integer::i,j,n,k
      real,dimension(0:im+1,0:jm+1,nm)::sm
      integer,dimension(0:im+1,0:jm+1,nm)::ioneta
      real,dimension(0:im+1,0:jm+1,nm,4)::stc,smc,sh2o
      real::bx,fk,frh2o
      integer,PARAMETER:: NSOTYP=9
      real,parameter::blim=5.5,hlice=3.335E5,grav=9.80616,t0=273.15
      real,dimension(NSOTYP)::beta2,psis,smcmax
      DATA BETA2 /4.26,8.72,11.55,4.74,10.73,8.17,6.77,5.25,4.26/
      DATA PSIS /0.04,0.62,0.47,0.14,0.10,0.26,0.14,0.36,0.04/
      DATA SMCMAX /0.421,0.464,0.468,0.434,0.406,  &
                   0.465,0.404,0.439,0.421/

       print *,'ioneta,21,21,6,',ioneta(21,21,6)
! Add sh2o stuff right here
        do n=1,nm
        DO J=1,JM
        DO I=1,IM
        DO K=1,4
! ----------------------------------------------------------------------
! cold start:  determine liquid soil water content (SH2O)
! SH2O <= SMC for T < 273.149K (-0.001C)
!
!new
!
        if (SM(I,J,n) .eq. 1) then
        STC(I,J,n,K)=amax1(273.15,STC(I,J,n,K))
        endif
!
!new
!
            IF (STC(I,J,n,K) .LT. 273.149) THEN
! ----------------------------------------------------------------------
! first guess following explicit solution for Flerchinger Eqn from Koren
! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
!	print *,i,j,n,k,stc(i,j,n,k),ioneta(i,j,n)
              BX = BETA2(IONETA(I,J,n))
              IF ( BETA2(IONETA(I,J,n)) .GT. BLIM ) BX = BLIM
              FK = (((HLICE/(GRAV*(-PSIS(IONETA(I,J,n)))))*  &
                   ((STC(I,J,n,K)-T0)/STC(I,J,n,K)))**  &
                   (-1/BX))*SMCMAX(IONETA(I,J,n))
              IF (FK .LT. 0.02) FK = 0.02
              SH2O(I,J,n,K) = MIN ( FK, SMC(I,J,n,K) )
! ----------------------------------------------------------------------
! now use iterative solution for liquid soil water content using
! FUNCTION FRH2O (from the Eta "NOAH" land-surface model) with the
! initial guess for SH2O from above explicit first guess.
              SH2O(I,J,n,K)=FRH2O(STC(I,J,n,K),SMC(I,J,n,K),SH2O(I,J,n,K),  &
                          SMCMAX(IONETA(I,J,n)),BETA2(IONETA(I,J,n)),  &
                          PSIS(IONETA(I,J,n)))
            ELSE
! ----------------------------------------------------------------------
! SH2O = SMC for T => 273.149K (-0.001C)
              SH2O(I,J,n,K)=SMC(I,J,n,K)
!
            ENDIF
        enddo
        enddo
        enddo
        enddo

	print *,'sh2o,',sh2o(17,31,5,1)

        end subroutine sfch2o