REAL FUNCTION FRH2O(TKELV,SMC,SH2O,SMCMAX,B,PSIS)

      IMPLICIT NONE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC  PURPOSE:  CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT
CC  IF TEMPERATURE IS BELOW 273.15K (T0).  REQUIRES NEWTON-TYPE ITERATION
CC  TO SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF
CC  KOREN ET AL. (1999, JGR, VOL 104(D16), 19569-19585).
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C New version (FEB 2001): much faster and more accurate newton iteration
c achieved by first taking log of eqn cited above -- less than 4
c (typically 1 or 2) iterations achieves convergence.  Also, explicit
c 1-step solution option for special case of parameter Ck=0, which reduces
c the original implicit equation to a simpler explicit form, known as the
c ""Flerchinger Eqn". Improved handling of solution in the limit of
c freezing point temperature T0.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C INPUT:
C
C   TKELV.........Temperature (Kelvin)
C   SMC...........Total soil moisture content (volumetric)
C   SH2O..........Liquid soil moisture content (volumetric)
C   SMCMAX........Saturation soil moisture content (from REDPRM)
C   B.............Soil type "B" parameter (from REDPRM)
C   PSIS..........Saturated soil matric potential (from REDPRM)
C
C OUTPUT:
C   FRH2O.........supercooled liquid water content.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      REAL B
      REAL BLIM
      REAL BX
      REAL CK
      REAL DENOM
      REAL DF
      REAL DH2O
      REAL DICE
      REAL DSWL
      REAL ERROR
      REAL FK
      REAL FRH2O
      REAL GS
      REAL HLICE
      REAL PSIS
      REAL SH2O
      REAL SMC
      REAL SMCMAX
      REAL SWL
      REAL SWLK
      REAL TKELV
      REAL T0

      INTEGER NLOG
      INTEGER KCOUNT

      PARAMETER (CK=8.0)
C      PARAMETER (CK=0.0)
      PARAMETER (BLIM=5.5)
C      PARAMETER (BLIM=7.0)
      PARAMETER (ERROR=0.005)

      PARAMETER (HLICE=3.335E5)
      PARAMETER (GS = 9.81)
      PARAMETER (DICE=920.0)
      PARAMETER (DH2O=1000.0)
      PARAMETER (T0=273.15)

C  ###   LIMITS ON PARAMETER B: B < 5.5  (use parameter BLIM)  ####
C  ###   SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT  ####
C  ###   IS NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES    ####
C##################################################################
C
      BX = B
      IF ( B .GT. BLIM ) BX = BLIM
C------------------------------------------------------------------

C INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
      NLOG=0
      KCOUNT=0

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      IF (TKELV .GT. (T0 - 1.E-3)) THEN

        FRH2O=SMC

      ELSE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       IF (CK .NE. 0.0) THEN

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCC OPTION 1: ITERATED SOLUTION FOR NONZERO CK CCCCCCCCCCC
CCCCCCCCCCCC IN KOREN ET AL, JGR, 1999, EQN 17 CCCCCCCCCCCCCCCCC
C
C INITIAL GUESS FOR SWL (frozen content)
        SWL = SMC-SH2O
C KEEP WITHIN BOUNDS.
        IF(SWL .GT. (SMC-0.02)) THEN
          SWL=SMC-0.02
c         SWL=MAX(0.,SMC-0.02)
        ENDIF
        IF(SWL .LT. 0.)  THEN
          SWL=0.
        ENDIF
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  START OF ITERATIONS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        DO WHILE (NLOG .LT. 10 .AND. KCOUNT .EQ. 0)
         NLOG = NLOG+1
         DF = ALOG(( PSIS*GS/HLICE ) * ( ( 1.+CK*SWL )**2. ) *
     &        ( SMCMAX/(SMC-SWL) )**BX) - ALOG(-(TKELV-T0)/TKELV)
         DENOM = 2. * CK / ( 1.+CK*SWL ) + BX / ( SMC - SWL )
         SWLK = SWL - DF/DENOM
         DSWL=ABS(SWLK-SWL)
         SWL=SWLK
C KEEP WITHIN BOUNDS.
         IF(SWL .GT. (SMC-0.02)) THEN
           SWL=SMC-0.02
c          SWL=MAX(0.,SMC-0.02)
         ENDIF
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
CC WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
         IF ( DSWL .LE. ERROR )  THEN
           KCOUNT=KCOUNT+1
         END IF
        END DO
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  END OF ITERATIONS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C LOWER BOUND FOR SWL (ICE CONTENT)
        IF(SWL .LT. 0.)  THEN
          SWL=0.
        ENDIF
C UPPER BOUND FOR SWL ALREADY APPLIED WITHIN DO-BLOCK
        FRH2O = SMC - SWL
C
CCCCCCCCCCCCCCCCCCCCCCCC END OPTION 1 CCCCCCCCCCCCCCCCCCCCCCCCCCC

       ENDIF

       IF (KCOUNT .EQ. 0) THEN
         Print*,'Flerchinger used. Iterations=',NLOG

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCC OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0 CCCCCCCC
CCCCCCCCCCCCC IN KOREN ET AL., JGR, 1999, EQN 17  CCCCCCCCCCCCCCC
C
        FK=(((HLICE/(GS*(-PSIS)))*((TKELV-T0)/TKELV))**(-1/BX))*SMCMAX
        IF (FK .LT. 0.02) FK = 0.02
        FRH2O = MIN ( FK, SMC )
C
CCCCCCCCCCCCCCCCCCCCCCCCC END OPTION 2 CCCCCCCCCCCCCCCCCCCCCCCCCC

       ENDIF

      ENDIF

      RETURN
      END