MODULE SEAICE 22,1
!@sum  SEAICE contains all the sea ice related subroutines
!@auth Original Development Team
!@ver  1.0
!@cont PREC_SI,SEA_ICE,ADDICE,SIMELT
      USE CONSTANT, only : lhm,rhoi,byrhoi,rhow,shi,shw,byshi,bylhm,sday
     *     ,rhows

      IMPLICIT NONE
      SAVE
!@param LMI number of temperature layers in ice
      INTEGER, PARAMETER :: LMI = 4
!@param XSI fractions of mass layer in each temp. layer
!@param BYXSI recipricol of XSI
      REAL*8, PARAMETER, DIMENSION(LMI) ::
     *     XSI= (/0.5d0, 0.5d0, 0.5d0, 0.5d0/),
     *     BYXSI= (/ 1./XSI(1), 1./XSI(2), 1./XSI(3), 1./XSI(4) /)
!@param Z1I thickness of first layer ice (m)
      REAL*8, PARAMETER :: Z1I = .1d0
!@param ACE1I ice mass first layer (kg/m^2)
      REAL*8, PARAMETER :: ACE1I = Z1I*RHOI
!@param HC1I heat capacity of first layer ice (J/m^2)
      REAL*8, PARAMETER :: HC1I = ACE1I*SHI
!@param Z2OIM min. thickness of 2nd layer ice (m) (if > 0)
      REAL*8, PARAMETER :: Z2OIM = .1d0    ! .4d0
!@param AC2OIM min. ice mass 2nd layer (kg/m^2)   (if > 0)
      REAL*8, PARAMETER :: AC2OIM = Z2OIM*RHOI
!@param ALAMI,ALAMS lambda coefficient for ice/snow J/(m*degC*sec)
      REAL*8, PARAMETER :: ALAMI=2.1762d0, ALAMS=0.35d0
!@param RHOS density of snow (kg/m^3)
      REAL*8, PARAMETER :: RHOS = 300.0
!@var FLEADOC lead fraction for ocean ice (%)
      REAL*8, PARAMETER :: FLEADOC = 0.06d0
!@var FLEADLK lead fraction for lakes (%)
      REAL*8, PARAMETER :: FLEADLK = 0.
!@param BYRLI,BYRLS reciprical of density*lambda
      REAL*8, PARAMETER :: BYRLI = 1./(RHOI*ALAMI),
     *     BYRLS = 1./(RHOS*ALAMS)
!@param MU coefficient of seawater freezing point w.r.t. salinity
      REAL*8, PARAMETER :: MU = 0.054d0
!@param SSI0 default value for sea ice salinity (kg/kg) (=3.2ppt)
      REAL*8, PARAMETER :: SSI0 = 0.0032d0
!@param FSSS fraction of ocean salinity found in new-formed ice
      REAL*8, PARAMETER :: FSSS = 13d0/35d0
!@var qsfix flag is true if salinity of sea ice is constant
      LOGICAL :: QSFIX = .false.
!@var alpha implicity for heat diffusion in sea ice (1=fully implicit)
      REAL*8, PARAMETER :: ALPHA = 1.0
!@dbparam oi_ustar0 default ice-ocean friction velocity (m/s)
      REAL*8 :: oi_ustar0 = 1d-3  ! 5d-3 ! not used if ice dynamics is
!@dbparam silmfac factor controlling lateral melt of ocean ice
      REAL*8 :: silmfac = 1.4d-8 ! = pi*(3d-6)/0.66/1000
!@var silmpow exponent for temperature dependence of lateral melt
      REAL*8 :: silmpow = 1.36d0
!@dbparam snow_ice =1 to allow for snow ice formation (default=1)
      INTEGER :: snow_ice = 1
!@var osurf_tilt controls calc. of ocean surface tilt for ice dyn:
!@+       from geostrophy (=0) or from free surface (=1, default)
      INTEGER :: osurf_tilt = 1
!@var DEBUG flag
      LOGICAL DEBUG

      CONTAINS


      SUBROUTINE PREC_SI(SNOW,MSI2,HSIL,TSIL,SSIL,PRCP,ENRGP,RUN0,SRUN0, 1

     &     WETSNOW)
!@sum  PREC_SI Adds the precipitation to sea/lake ice
!@auth Gary Russell
!@ver  1.0
      IMPLICIT NONE
!@param SNOMAX maximum allowed snowdepth (1m equivalent) (kg/m^2)
!@param dSNdRN rate of conversion of snow to ice as a function of rain
      REAL*8, PARAMETER :: SNOMAX=1d0*RHOS, dSNdRN=0.
!@var ENRGP total energy of precip (C),(J/m^2)
!@var PRCP amount of precip (kg/m^2)
      REAL*8, INTENT(IN) :: ENRGP, PRCP
!@var SNOW snow mass (kg/m^2)
!@var MSI1 first layer ice mass (= SNOW + ACE1I) (kg/m^2)
!@var MSI2 second layer ice mass (kg/m^2)
      REAL*8, INTENT(INOUT) :: SNOW, MSI2
!@var HSIL enthalpy of ice layers (J/m^2)
!@var SSIL salt in ice layers (kg/m^2)
      REAL*8, INTENT(INOUT), DIMENSION(LMI) :: HSIL, SSIL
!@var TSIL temperature of ice layers (C) (DIAGNOSTIC ONLY)
      REAL*8, INTENT(OUT), DIMENSION(LMI) :: TSIL
!@var WETSNOW true if snow is wet (i.e. has been rained on)
      LOGICAL, INTENT(OUT) :: WETSNOW
!@var RUN0 runoff from ice (kg/m^2)
!@var SRUN0 salt in runoff from ice (kg/m^2)
      REAL*8, INTENT(OUT) :: RUN0, SRUN0
      REAL*8 :: BYMSI2, FMSI2, FMSI3, FHSI2, FHSI3,
     *     CMPRS, SNWF, RAIN, FREZI, MSI1, FSSI2, HFREZ,SNOW1,
     *     FSSI3, MELTI, MELTS, SMELTI, DSNOW, SICE, HCMPRS, HSNOW, HICE


      WETSNOW=.FALSE.
      RUN0=0. ; SRUN0=0.

      MSI1=SNOW+ACE1I

      IF (PRCP.gt.0) THEN

C**** Snowfall is calculated from precip energy (0 deg or colder)
      SNWF = MAX(0d0,MIN(PRCP,-ENRGP*BYLHM))
C**** Rain is remaining precip (0 deg or warmer)
      RAIN = PRCP-SNWF
      WETSNOW = RAIN.GT.1d-5*PRCP  ! i.e. a noticeable fraction of prec

C**** New formulation: separate out snow and ice components
C**** Calculate remaining snow and necessary mass flux using
C**** SNOW =SNOW  + SNWF  - CMPRS - MELTS
C**** ACE1I=ACE1I + FREZI + CMPRS - FMSI2
C**** SALTI=SALTI - SMELTI        - FSSI2
      IF (SNOW.gt.0) THEN ! apply fluxes to snow portion first
        SNOW1 = MIN(SNOW,XSI(1)*MSI1)
        HSNOW = HSIL(1)*MIN(SNOW1/(XSI(1)*MSI1),1d0)
        HICE  = HSIL(1)-HSNOW
        SICE  = SSIL(1)

        MELTS=MAX(0d0,(HSNOW+ENRGP)*BYLHM+SNOW1+SNWF)
        IF (MELTS.gt.SNOW1+SNWF) THEN   ! all snow and some ice melts
          MELTS = SNOW1+SNWF
          MELTI = MAX(0d0,(HICE+HSNOW+ENRGP)*BYLHM/(1.-SICE/(XSI(1)*MSI1
     *         -SNOW1))+XSI(1)*MSI1-SNOW1)
          SMELTI = MELTI*SICE/(XSI(1)*MSI1-SNOW1)
          FREZI = 0.
          HICE = HICE+HSNOW+ENRGP
          HSNOW = 0.
          SNOW1 = 0.
          CMPRS = 0.

c         ACE1I = ACE1I - MELTI
        ELSE ! some snow remains
          MELTI =0.
          SMELTI=0.
          FREZI=MIN(RAIN,MAX(          HICE  = HICE - LHM*FREZI
          HSNOW = HSNOW + ENRGP + LHM*FREZI
          IF (SNWF.ge.PRCP .and. SNOW+SNWF .GT. SNOMAX) THEN
C**** TOO MUCH SNOW HAS ACCUMULATED, SOME SNOW IS COMPACTED INTO ICE
            CMPRS = SNOW+SNWF-0.9d0*SNOMAX
          ELSE
C**** RAIN and MELT COMPRESSES SNOW INTO ICE
            CMPRS = MIN(dSNdRN*(RAIN+MELTS), SNOW+SNWF-MELTS)
          END IF
          HCMPRS = HSNOW*CMPRS/(SNOW1+SNWF-MELTS)
          HSNOW = HSNOW - HCMPRS
          HICE  = HICE  + HCMPRS

          SNOW1  = SNOW1 + SNWF - MELTS - CMPRS
c         ACE1I = ACE1I + FREZI + CMPRS
        END IF
      ELSE
        HSNOW = MIN(ENRGP,0d0)  ! new snow
        HICE  = HSIL(1)+ENRGP-HSNOW
        SICE  = SSIL(1)
        MELTS = 0.
        CMPRS = 0.
        MELTI = MAX(0d0,HICE*BYLHM/(1.-SICE/(XSI(1)*MSI1))+XSI(2)*MSI1)
        SMELTI= MELTI*SICE/(XSI(2)*MSI1)
        FREZI = MIN(RAIN,MAX(-HICE*BYLHM-XSI(2)*MSI1+SICE,0d0))
        SNOW1 = SNWF
c       ACE1I = ACE1I - MELTI + FREZI

      END IF
      SICE = SICE - SMELTI

C**** Mass fluxes required to keep first layer ice = ACE1I
      FMSI2 = CMPRS+FREZI-MELTI ! either up or down

C***** Calculate consequential heat/salt flux between layers
      IF (FMSI2.gt.0) THEN      ! downward flux to thermal layer 3
        FHSI2 = FMSI2*HSIL(2)/(XSI(2)*MSI1)
        FSSI2 = FMSI2*SSIL(2)/(XSI(2)*MSI1)

      ELSE                      ! upward flux
        FHSI2  = FMSI2*HSIL(3)/(XSI(3)*MSI2)
        FSSI2  = FMSI2*SSIL(3)/(XSI(3)*MSI2)

      END IF

C**** Calculate total snow and ice in upper layers
      IF (ACE1I.gt.XSI(2)*MSI1) THEN ! some ice originally in layer 1
        SNOW = SNOW1
        HICE = HICE + HSIL(2)
        SICE = SICE + SSIL(2)

      ELSE ! some snow orignally in second layer
        SNOW = SNOW1 + (XSI(2)*MSI1-ACE1I)
        SICE = SSIL(2)
        HSNOW = HSNOW + (HSIL(2)-LHM*SICE)*(XSI(2)*MSI1-ACE1I)/(XSI(2)
     $       *MSI1)
        HICE = HICE + (HSIL(2)-LHM*SICE)*ACE1I/(XSI(2)*MSI1) + LHM*SICE

      END IF

C**** reconstitute upper layers
      MSI1 = ACE1I + SNOW
      IF (ACE1I.gt.XSI(2)*MSI1) THEN ! some ice in first layer
        HSIL(1) = (HICE-FHSI2)*(ACE1I-XSI(2)*MSI1)/ACE1I + HSNOW
        SSIL(1) = (SICE-FSSI2)*(ACE1I-XSI(2)*MSI1)/ACE1I ! + SSNOW

      ELSE
        HSIL(1) = HSNOW*XSI(1)*MSI1/(MSI1-ACE1I)
        SSIL(1) = 0.            ! SSNOW*XSI(1)*MSI1/(MSI1-ACE1I)

      END IF
      HSIL(2) = HICE - HSIL(1) - FHSI2 + HSNOW
      SSIL(2) = SICE - SSIL(1) - FSSI2 ! + SSNOW


C**** Calculate mass/heat flux between layers 3 and 4
      FMSI3 = XSI(4)*FMSI2
      IF (FMSI3.gt.0) THEN      ! downward flux to layer 4
        FHSI3 = FMSI3*HSIL(3)/(XSI(3)*MSI2)
        FSSI3 = FMSI3*SSIL(3)/(XSI(3)*MSI2)

      ELSE                      ! upward flux
        FHSI3 = FMSI3*HSIL(4)/(XSI(4)*MSI2)
        FSSI3 = FMSI3*SSIL(4)/(XSI(4)*MSI2)

      END IF

C**** Adjust amounts resulting from fluxes.
      MSI2 = MSI2 + FMSI2

      HSIL(3)=HSIL(3)+(FHSI2-FHSI3)
      HSIL(4)=HSIL(4)+ FHSI3

      SSIL(3)=SSIL(3)+(FSSI2-FSSI3)
      SSIL(4)=SSIL(4)+ FSSI3



C**** Diagnostics for output
      RUN0 = MELTS+MELTI+(RAIN-FREZI) ! runoff mass to ocean
      IF (RUN0.lt.1d-13) RUN0=0. ! clean up roundoff errors
      SRUN0 = SMELTI             ! salt in runoff
c     HFLUX= 0                  ! energy of runoff (currently at 0 deg)


      END IF

 ! ice temp L=1,2,3,4
      TSIL(1:2) = ((HSIL(1:2)-LHM*SSIL(1:2))/(XSI(1:2)*MSI1)+LHM)*BYSHI
      TSIL(3:4) = ((HSIL(3:4)-LHM*SSIL(3:4))/(XSI(3:4)*MSI2)+LHM)*BYSHI

      RETURN
      END SUBROUTINE PREC_SI


      SUBROUTINE SEA_ICE(DTSRCE,SNOW,ROICE,HSIL,SSIL,MSI2,F0DT,F1DT,EVAP 1,1
     &     ,SROX,

     *     FMOC,FHOC,FSOC,RUN,ERUN,SRUN,WETSNOW,MELT12)
!@sum  SEA_ICE applies surface fluxes to ice covered areas
!@auth Gary Russell
!@ver  1.0
      IMPLICIT NONE

      REAL*8, PARAMETER :: dSNdML =0.
!@var DTSRCE source time step (s)
      REAL*8, INTENT(IN) :: DTSRCE
!@var F0DT heat flux on the ice top surface (J/m^2)
      REAL*8, INTENT(IN) :: F0DT
!@var F1DT heat flux between the 1st and 2nd ice layers (J/m^2)
      REAL*8, INTENT(IN) :: F1DT
!@var SROX solar radiation at top and bottom ice surface (J/m^2)
      REAL*8, INTENT(INOUT) :: SROX(2)
!@var EVAP evaporation/dew on the top ice surface (kg/m^2)
      REAL*8, INTENT(IN) :: EVAP
!@var FSRI fraction of solar radiation that goes through the bottom
      REAL*8 :: FSRI(LMI)
!@var FMOC,FHOC,FSOC basal fluxes down of mass,heat,salt (J or kg/m^2)
      REAL*8, INTENT(IN) :: FMOC, FHOC, FSOC
!@var RUN,ERUN,SRUN runoff fluxes down of mass,heat,salt (J or kg/m^2)
      REAL*8, INTENT(OUT) :: RUN,ERUN,SRUN
!@var WETSNOW whether snow is wet or not (Used for albedo calc)
      LOGICAL, INTENT(INOUT) :: WETSNOW
!@var MELT12 save surface melt (Used for albedo calc)
      REAL*8, INTENT(OUT) :: MELT12

!@var ROICE sea ice fraction of open water
      REAL*8, INTENT(IN) :: ROICE
!@var HSIL enthalpy of ice layers (J/m^2)
!@var SSIL salt in ice layers (kg/m^2)
      REAL*8, INTENT(INOUT), DIMENSION(LMI) :: HSIL, SSIL
!@var SNOW snow mass (kg/m^2)
!@var MSI1 first layer ice mass (= SNOW + ACE1I) (kg/m^2)
!@var MSI2 second layer ice mass (kg/m^2)
      REAL*8, INTENT(INOUT) :: SNOW, MSI2
!@var TSIL temperature of ice layers (C)
      REAL*8, DIMENSION(LMI) :: TSIL
      REAL*8 MSI1, MELTS, MELTI, MELT3, MELT4, DSNOW
      REAL*8 DEW, CMPRS, DEWS, DEWI
      REAL*8 SMELTI, SMELT3, SMELT4, FSSI2, FSSI3, SICE
      REAL*8 FMSI2, FMSI3, FHSI2, FHSI3
      REAL*8 HC1, HC2, HC3, HC4, HICE, HSNOW, HCMPRS
      REAL*8 dF1dTI, dF2dTI, dF3dTI, dF4dTI, F1, F2, F3, FO

      FMSI2=0. ; FHSI2=0. ; FHSI3=0. ; FSSI2=0. ; FSSI3=0.
      MELTS=0. ; MELTI=0. ; MELT3=0. ; MELT4=0.
      SMELTI=0 ; SMELT3=0.; SMELT4=0.

C****
      MSI1 = SNOW+ACE1I ! snow and first (physical) layer ice mass
C**** Calculate solar fractions
      IF (SROX(1).gt.0) THEN
        call solar_ice_frac(snow,msi2,wetsnow,fsri,4)
      ELSE
        FSRI = 0.
      END IF
C****
C**** OCEAN ICE, CALCULATE TSIL FROM ENTHALPY
C****
 ! ice temp L=1,2,3,4
      TSIL(1:2) = ((HSIL(1:2)-LHM*SSIL(1:2))/(XSI(1:2)*MSI1)+LHM)*BYSHI
      TSIL(3:4) = ((HSIL(3:4)-LHM*SSIL(3:4))/(XSI(3:4)*MSI2)+LHM)*BYSHI

      HC1 = SHI*XSI(1)*MSI1 ! heat capacity of ice L=1 (J/(degC*m^2))
      HC2 = SHI*XSI(2)*MSI1 ! heat capacity of ice L=2 (J/(degC*m^2))
      HC3 = SHI*XSI(3)*MSI2 ! heat capacity of ice L=3 (J/(degC*m^2))
      HC4 = SHI*XSI(4)*MSI2 ! heat capacity of ice L=4 (J/(degC*m^2))
C**** CALCULATE AND APPLY DIFFUSIVE AND SURFACE ENERGY FLUXES
C**** First layer flux calculated already as part of surface calculation
c      dF1dTI = 2.*DTSRCE/(ACE1I*BYRLI+SNOW*BYRLS)
      dF2dTI = ALAMI*RHOI*DTSRCE/(0.5*XSI(2)*MSI1+0.5*XSI(3)*MSI2)
C****          temperature derivative from F2 diffusive flux
      dF3dTI = ALAMI*RHOI*DTSRCE*2./MSI2
C****          temperature derivative from F3 diffusive flux
      dF4dTI = ALAMI*RHOI*DTSRCE*2.*BYXSI(4)/MSI2
C****          temperature derivative from F4 diffusive flux
C**** EXPLCIIT DIFFUSIVE FLUXES
CEXP  F2 = dF2dTI*(TSIL(2)-TSIL(3))+SROX(1)*FSRI(2)
CEXP  F3 = dF3dTI*(TSIL(3)-TSIL(4))+SROX(1)*FSRI(3)
CEXP  FO = dF4dTI*(TSIL(4)-TFO)    +SROX(1)*FSRI(4)
C**** DIFFUSIVE FLUXES FROM PARTIALLY IMPLICIT METHOD
c      F1DT=(dF1dTI*(HC1*(TSIL(1)-TSIL(2))+ALPHA*F0DT)
c     *    +HC1*SROX(1)*FSRI(1))/(HC1+ALPHA*dF1dTI) ! already calculated
      F2=(dF2dTI*(HC2*(TSIL(2)-TSIL(3))+ALPHA*F1DT)+HC2*SROX(1)*FSRI(2))
     *     /(HC2+ALPHA*dF2dTI)
      F3=(dF3dTI*(HC3*(TSIL(3)-TSIL(4))+ALPHA*F2)+HC3*SROX(1)*FSRI(3))/
     *     (HC3+ALPHA*dF3dTI)
c      FO=(dF4dTI*(HC4*(TSIL(4)-TFO)+ALPHA*F3)+HC4*SROX(1)*FSRI(4))/
c     *     (HC4+ALPHA*dF4dTI)
c      FHOC=(dF4dTI*ALPHA*F3+HC4*(SROX(1)*FSRI(4)+FHOC))/
c     *     (HC4+ALPHA*dF4dTI)
C**** Add solar flux through bottom to basal heat flux already computed
      SROX(2) = SROX(1)*FSRI(4)

      HSIL(1) = HSIL(1)+(F0DT-F1DT)
      HSIL(2) = HSIL(2)+(F1DT-F2)
      HSIL(3) = HSIL(3)+(F2-F3)
      HSIL(4) = HSIL(4)+(F3-SROX(2)-FHOC)

C**** add basal salt flux
      SSIL(4) = SSIL(4)-FSOC


      DEW = -EVAP               ! dew/evap to the surface
      SICE=SSIL(1)+SSIL(2)      ! total upper layer salt (only in ice)

C**** Add DEW directly to ice (which can be in either layer)
C**** DEWI adds to ice, DEWS adds to snow (if applicable)
C**** Calculate melting in first two thermal layers
C**** MELTS,MELTI are the amounts of snow and ice melt
      IF (SNOW*XSI(2).gt.XSI(1)*ACE1I) THEN ! first layer is all snow
        DEWI = MAX(0d0,DEW)    ! >0 i.e. dew to second layer ice
        DEWS = DEW - DEWI      ! <0 remaining evap applied to snow
        HSNOW = HSIL(1) + HSIL(2)*MAX(1.-(ACE1I+DEWI)/(XSI(2)*MSI1+DEWI)
     *       ,0d0)

      ELSE  ! first layer is snow and some ice
        DEWS = -MIN(SNOW,-(DEW-MAX(0d0,DEW))) ! <0 evap upto snow amount
        DEWI = DEW - DEWS       ! either sign, evap/dew applied to ice
        HSNOW = HSIL(1)*MIN((SNOW+DEWS)/(XSI(1)*MSI1+DEW),1d0)

      END IF



C**** calculate melt in snow and ice layers
      HICE = HSIL(1)+HSIL(2)-HSNOW
      MELTS = MAX(0d0,HSNOW*BYLHM+SNOW+DEWS)  ! snow melt amount
      MELTI = MAX(0d0,HICE*BYLHM/(1.-SICE/(ACE1I+DEWI))+ACE1I+DEWI)
      SMELTI= MELTI*SICE/(ACE1I+DEWI)
C**** CMPRS is amount of snow turned to ice during melting
C**** Calculate remaining snow and necessary mass flux using
C**** SNOW =SNOW - MELTS + DEWS - CMPRS
C**** ACE1I=ACE1I- MELTI + DEWI + CMPRS - FMSI2
C**** SICE =SICE- SMELTI                - FSSI2
C****
      IF (SNOW.GT.MELTS-DEWS) THEN ! some snow remains
        CMPRS = MIN(dSNdML*MELTS, SNOW+DEWS-MELTS) ! > 0.
        DSNOW = - (MELTS-DEWS+CMPRS)
        HCMPRS = HSNOW*CMPRS/(SNOW+DEWS-MELTS)

      ELSE ! all snow and some ice melts or evaporates
        CMPRS = 0.
        DSNOW = -SNOW
        HCMPRS=0.

      END IF

C**** Mass fluxes required to keep first layer ice = ACE1I
      FMSI2 = -MELTI+DEWI+CMPRS  ! either up or down

C**** Check for melting in levels 3 and 4
      MELT3 = MAX(0d0,HSIL(3)*BYLHM/(1.-SSIL(3)/(XSI(3)*MSI2))+XSI(3)
     *     *MSI2)
      MELT4 = MAX(0d0,HSIL(4)*BYLHM/(1.-SSIL(4)/(XSI(4)*MSI2))+XSI(4)
     *     *MSI2-FMOC)
      SMELT3 = MELT3*SSIL(3)/(XSI(3)*MSI2)
      SMELT4 = MELT4*SSIL(4)/(XSI(4)*MSI2-FMOC)

C***** Calculate consequential heat flux between layers
      IF (FMSI2.gt.0) THEN ! downward flux to thermal layer 3
        FHSI2 = FMSI2*(HICE+HCMPRS)/(ACE1I+DEWI-MELTI+CMPRS)
        FSSI2 = FMSI2*(SICE-SMELTI)/(ACE1I+DEWI-MELTI+CMPRS)

      ELSE                 ! upward flux
        FHSI2  = FMSI2*HSIL(3)/(XSI(3)*MSI2-MELT3)
        FSSI2  = FMSI2*(SSIL(3)-SMELT3)/(XSI(3)*MSI2-MELT3)

      END IF

C**** reconstitute upper layers
      SNOW = MAX(0d0,SNOW+DSNOW)
      MSI1 = ACE1I + SNOW
      IF (ACE1I.gt.XSI(2)*MSI1) THEN ! some ice in first layer
        HSIL(1) = (HICE+HCMPRS-FHSI2)*(ACE1I-XSI(2)*MSI1)/ACE1I + HSNOW
     *       -HCMPRS
        SSIL(1) = (SICE-SMELTI-FSSI2)*(ACE1I-XSI(2)*MSI1)/ACE1I ! +SSNOW

      ELSE
        HSIL(1) = (HSNOW-HCMPRS)*XSI(1)*MSI1/(MSI1-ACE1I)
        SSIL(1) = 0.            ! SSNOW*XSI(1)*MSI1/(MSI1-ACE1I)

      END IF
      HSIL(2) = HICE          - HSIL(1) - FHSI2 + HSNOW
      SSIL(2) = SICE - SMELTI - SSIL(1) - FSSI2 ! + SSNOW


C**** Calculate mass/heat flux between layers 3 and 4
      FMSI3 = XSI(3)*(MELT4+FMOC) + XSI(4)*(FMSI2-MELT3)
      IF (FMSI3.gt.0) THEN ! downward flux to layer 4
        FHSI3 = FMSI3*HSIL(3)/(XSI(3)*MSI2-MELT3)
        FSSI3 = FMSI3*(SSIL(3)-SMELT3)/(XSI(3)*MSI2-MELT3)

      ELSE                 ! upward flux
        FHSI3 = FMSI3*HSIL(4)/(XSI(4)*MSI2-MELT4-FMOC)
        FSSI3 = FMSI3*(SSIL(4)-SMELT4)/(XSI(4)*MSI2-MELT4-FMOC)

      END IF

C**** Apply fluxes to lower layers
      MSI2=MSI2-(MELT3+MELT4)+FMSI2-FMOC

      HSIL(3)=HSIL(3)+(FHSI2-FHSI3)
      HSIL(4)=HSIL(4)+ FHSI3

      SSIL(3)=SSIL(3)-SMELT3 +(FSSI2-FSSI3)
      SSIL(4)=SSIL(4)-SMELT4 + FSSI3


C**** Calculate additional runoff output fluxes
      RUN = MELTS + MELTI + MELT3 + MELT4 ! mass flux to ocean
      SRUN=        SMELTI +SMELT3 +SMELT4 ! salt flux to ocean
c HMELT currently assumed to be zero since melting is at 0 deg
      ERUN= SROX(2)        ! + HMELT


c**** Decide WETSNOW for albedo calculations and save MELT12 for
c**** possible melt ponds
      MELT12=MELTS+MELTI
      WETSNOW = WETSNOW .or. (MELT12.gt.0)
C****
      RETURN
      END SUBROUTINE SEA_ICE


      SUBROUTINE ADDICE (SNOW,ROICE,HSIL,SSIL,MSI2,TSIL,ENRGFO 1
     &     ,ACEFO,ACEFI,ENRGFI,SALTO,SALTI,

     *     DMIMP,DHIMP,DSIMP,FLEAD,QFIXR)
!@sum  ADDICE adds ice formed in the ocean to ice variables
!@auth Gary Russell/Gavin Schmidt
!@ver  1.0
      IMPLICIT NONE

      REAL*8, PARAMETER, DIMENSION(LMI) :: YSI =
     *     (/XSI(1)* ACE1I/(ACE1I+AC2OIM),XSI(2)*ACE1I/(ACE1I+AC2OIM),
     *       XSI(3)*AC2OIM/(ACE1I+AC2OIM),XSI(4)*AC2OIM/(ACE1I+AC2OIM)/)
c     REAL*8, PARAMETER :: Z2OIX = 5.d0-Z1I,     ! max ice thickness l=2
c    *                     BYZICX=1./(Z1I+Z2OIX)
!@var QFIXR true if RSI and MSI2 are fixed (ie. for fixed SST run)
      LOGICAL, INTENT(IN) :: QFIXR
!@var FLEAD minimum lead fraction for ice (%)
      REAL*8, INTENT(IN) :: FLEAD
!@var ACEFO ice mass formed in ocean for open water fraction (kg/m^2)
!@var ACEFI ice mass formed in ocean for ice covered fraction (kg/m^2)
!@var ENRGFO energy of ice formed in ocean for open water frac (J/m^2)
!@var ENRGFI energy of ice formed in ocean for ice covered frac (J/m^2)
!@var SALTO salt in ice formed in ocean for open water frac (kg/m^2)
!@var SALTI salt in ice formed in ocean for ice covered frac (kg/m^2)
      REAL*8, INTENT(IN) ::  ENRGFI, ENRGFO, ACEFO, ACEFI, SALTO, SALTI
!@var ROICE  ice fraction over open water
!@var SNOW snow amount on ice (kg/m^2)
!@var MSI2 second mass layer ice thickness (kg/m^2)
      REAL*8, INTENT(INOUT) :: ROICE, SNOW, MSI2
!@var HSIL enthalpy of ice layers (J/m^2)
!@var SSIL salt in ice layers (kg/m^2)
      REAL*8, INTENT(INOUT), DIMENSION(LMI) :: HSIL,SSIL
!@var TSIL temperature of ice layers (C)
      REAL*8, INTENT(OUT), DIMENSION(LMI) :: TSIL
!@var DMIMP,DSIMP,DHIMP implicit mass,salt and heat fluxes required
!@+   to maintain minimum ice thickness if ice fraction is fixed
      REAL*8, INTENT(OUT) :: DMIMP,DSIMP,DHIMP

      REAL*8, DIMENSION(LMI) :: FRI
      REAL*8 FMSI4, FHSI3, FHSI4, FSSI3, FSSI4, HSNOW, HICE, SICE
      REAL*8 ROICEN, OPNOCN, DRSI, MSI1

      DMIMP=0. ; DHIMP=0. ; DSIMP=0.


      MSI1=SNOW+ACE1I
      IF (.not. QFIXR) THEN
      IF (ROICE.LE.0. .and. ACEFO.gt.0) THEN
C**** Create new ice in ice-free ocean
        ROICE=ACEFO/(ACE1I+AC2OIM)
        SNOW=0.
C****   TSIL=(ENRGFO/ACEFO + (1-S)*LHM)*BYSHI ! but hsi is primary var.
        HSIL(1:2) =(ENRGFO/ACEFO)*XSI(1:2)*ACE1I
        HSIL(3:4) =(ENRGFO/ACEFO)*XSI(3:4)*AC2OIM
        SSIL(1:2) = (SALTO/ACEFO)*XSI(1:2)*ACE1I
        SSIL(3:4) = (SALTO/ACEFO)*XSI(3:4)*AC2OIM

        MSI1=ACE1I
        MSI2=AC2OIM
      ELSEIF (ROICE.gt.0) THEN
C**** Create new ice in partially ice-covered ocean
      IF (ACEFI.gt.0) THEN
C**** CALCULATE ADVECTIVE HEAT FLUX FROM LAYER 3 TO LAYER 4 OF ICE
CC      FMSI3 = -XSI(3)*ACEFI ! < 0.
CC      FMSI4 = -ACEFI
        FHSI3 = -HSIL(4)*ACEFI*(XSI(3)/XSI(4))/MSI2
        FSSI3 = -SSIL(4)*ACEFI*(XSI(3)/XSI(4))/MSI2

      ELSE
        FHSI3 = 0. ; FSSI3 = 0.

      END IF

      IF (ACEFO .EQ. 0.) THEN
C**** NEW ICE IS ONLY FORMED BELOW OLD SEA ICE
        HSIL(3) = HSIL(3)-FHSI3
        HSIL(4) = HSIL(4)+(FHSI3+ENRGFI)
        SSIL(3) = SSIL(3)-FSSI3
        SSIL(4) = SSIL(4)+(FSSI3+SALTI)

        MSI2 = MSI2+ACEFI       ! new ice mass of physical layer 2
      ELSE
C**** NEW ICE IS FORMED ON OPEN OCEAN AND POSSIBLY BELOW OLD SEA ICE
        DRSI = (1.-ROICE)*ACEFO/(ACE1I+AC2OIM) ! new ice on the open oc.
        ROICEN=ROICE+DRSI
C**** New formulation:
C**** Split all upper mass layer fields into ice and snow components
        IF (ACE1I.gt.XSI(2)*MSI1) THEN ! some ice in first layer
          HSNOW = (HSIL(1)-LHM*SSIL(1))*SNOW/(XSI(1)*MSI1)
        ELSE
          HSNOW = HSIL(1) + (HSIL(2)-LHM*SSIL(2))*(SNOW-XSI(1)*MSI1)
     *         /(XSI(2)*MSI1)
        END IF
        HICE  = HSIL(1)+HSIL(2)-HSNOW
CC      SSNOW = 0.   ! always zero
        SICE  = SSIL(1)+SSIL(2)

C**** distribute snow variables over new ice extent
        SNOW = SNOW*(ROICE/ROICEN)
        HSNOW = HSNOW*(ROICE/ROICEN)
CC      SSNOW = SSNOW*(ROICE/ROICEN)  ! always zero

C**** Add new ice to ice variables
        HICE=((1.-ROICE)*ENRGFO*ACE1I/(ACE1I+AC2OIM)+ROICE*HICE)/ROICEN
        SICE=((1.-ROICE)*SALTO *ACE1I/(ACE1I+AC2OIM)+ROICE*SICE)/ROICEN

C**** reconstitute upper thermal layers
        MSI1 = SNOW + ACE1I     ! mass of layer 1
CC      MSI1 = (DRSI*ACE1I+ROICE*MSI1)/(ROICE+DRSI) ! mass of layer 1
        IF (ACE1I.gt.XSI(2)*MSI1) THEN ! some ice in first layer
          HSIL(1) = HICE*(ACE1I-XSI(2)*MSI1)/ACE1I + HSNOW
          SSIL(1) = SICE*(ACE1I-XSI(2)*MSI1)/ACE1I ! + SSNOW

        ELSE
          HSIL(1) = HSNOW*XSI(1)*MSI1/SNOW
          SSIL(1) = 0.          ! SSNOW*XSI(1)*MSI1/SNOW

        END IF
        HSIL(2) = HICE - HSIL(1) + HSNOW
        SSIL(2) = SICE - SSIL(1) ! + SSNOW

C**** add new ice to old ice (not including snow)
        MSI2 = (DRSI*AC2OIM+ROICE*(MSI2+ACEFI))/ROICEN ! layer 2
C**** COMBINE OPEN OCEAN AND SEA ICE FRACTIONS TO FORM NEW VARIABLES
        HSIL(3)=((1.-ROICE)*ENRGFO*YSI(3)+ROICE*(HSIL(3)-FHSI3))/ROICEN
        HSIL(4)=((1.-ROICE)*ENRGFO*YSI(4)+ROICE*(HSIL(4)+FHSI3+ENRGFI)
     *       )/ROICEN
        SSIL(3)=((1.-ROICE)*SALTO*YSI(3)+ROICE*(SSIL(3)-FSSI3))/ROICEN
        SSIL(4)=((1.-ROICE)*SALTO*YSI(4)+ROICE*(SSIL(4)+FSSI3+SALTI))
     *       /ROICEN

        ROICE = ROICEN          ! new ice concentration
      END IF
C**** COMPRESS THE ICE HORIZONTALLY IF TOO THIN OR LEAD FRAC. TOO SMALL
      OPNOCN=MIN(0.1d0,FLEAD*RHOI/(ROICE*(ACE1I+MSI2)))
      IF ((ROICE*(ACE1I+MSI2)).gt.5.*RHOI) OPNOCN=0.  ! no leads for h>5
      IF (MSI2.LT.AC2OIM .or. ROICE.GT.1.-OPNOCN) THEN
      ROICEN = MIN(ROICE*(ACE1I+MSI2)/(ACE1I+AC2OIM),1.-OPNOCN)
      DRSI = ROICEN-ROICE ! < 0. compressed ice concentration
CC    FMSI3 = XSI(3)*FMSI4 ! < 0. upward ice mass flux into layer 3
      FMSI4 = (MSI1+MSI2)*(DRSI/ROICEN) ! upward ice mass into layer 4
      FHSI3 = HSIL(4)*FMSI4*(XSI(3)/XSI(4))/MSI2 ! upward heat flux
      FHSI4 = (HSIL(1)+HSIL(2)+HSIL(3)+HSIL(4))*(DRSI/ROICEN)
      HSIL(3) = HSIL(3)-FHSI3
      HSIL(4) = HSIL(4)+(FHSI3-FHSI4)
      FSSI3 = SSIL(4)*FMSI4*(XSI(3)/XSI(4))/MSI2 ! upward salt flux
      FSSI4 = (SSIL(1)+SSIL(2)+SSIL(3)+SSIL(4))*(DRSI/ROICEN)
      SSIL(3) = SSIL(3)-FSSI3
      SSIL(4) = SSIL(4)+(FSSI3-FSSI4)

      MSI2 = MSI2-FMSI4 ! new ice mass of second physical layer
CC    SNOW = SNOW   ! snow thickness is conserved
      ROICE = ROICEN
C****
      END IF
C****
      END IF
C**** Clean up ice fraction (if rsi>(1-OPNOCN)-1d-3) => rsi=(1-OPNOCN))
      IF (ROICE.gt.0) THEN
      OPNOCN=MIN(0.1d0,FLEAD*RHOI/(ROICE*(ACE1I+MSI2)))    ! -BYZICX)
      IF ((ROICE*(ACE1I+MSI2)).gt.5.*RHOI) OPNOCN=0.  ! no leads for h>5
      IF (ROICE.gt.(1.-OPNOCN-1d-3)) THEN
        ROICEN = 1.-OPNOCN
        DRSI = ROICEN-ROICE ! +ve
        FMSI4 = (ACE1I+MSI2)*(DRSI/ROICEN) ! >0 ice mass out of layer 4
        FHSI4 = HSIL(4)*FMSI4/(XSI(4)*MSI2)
        FSSI4 = SSIL(4)*FMSI4/(XSI(4)*MSI2)
        FHSI3 = HSIL(3)*FMSI4/MSI2 ! downward heat flux into layer 4
        FSSI3 = SSIL(3)*FMSI4/MSI2 ! downward salt flux into layer 4

        MSI2=MSI2-FMSI4         ! new ice mass of second physical layer
        FRI(1)=ACE1I/(ACE1I+MSI2)
        FRI(3:4)=XSI(3:4)*MSI2/(ACE1I+MSI2)
C**** New formulation:
C**** Split all upper mass layer fields into ice and snow components
        IF (ACE1I.gt.XSI(2)*MSI1) THEN ! some ice in first layer
          HSNOW = (HSIL(1)-LHM*SSIL(1))*SNOW/(XSI(1)*MSI1)
        ELSE
          HSNOW = HSIL(1) + (HSIL(2)-LHM*SSIL(2))*(SNOW-XSI(1)*MSI1)
     *         /(XSI(2)*MSI1)
        END IF
        HICE  = HSIL(1)+HSIL(2)-HSNOW
CC      SSNOW = 0.   ! always zero
        SICE  = SSIL(1)+SSIL(2)

C**** distribute snow variables over new ice extent
        SNOW = SNOW*(ROICE/ROICEN)
        HSNOW = HSNOW*(ROICE/ROICEN)
CC      SSNOW = SSNOW*(ROICE/ROICEN)  ! always zero

C**** Add new ice to ice variables
        HICE = (ROICE/ROICEN)*(FHSI4*FRI(1)+HICE)
        SICE = (ROICE/ROICEN)*(FSSI4*FRI(1)+SICE)

C**** reconstitute upper thermal layers
        MSI1 = SNOW + ACE1I     ! mass of layer 1
CC      MSI1 = (DRSI*ACE1I+ROICE*MSI1)/(ROICE+DRSI) ! mass of layer 1
        IF (ACE1I.gt.XSI(2)*MSI1) THEN ! some ice in first layer
          FRI(1)=(ACE1I-XSI(2)*MSI1)/ACE1I
          HSIL(1) = HICE*FRI(1) + HSNOW
          SSIL(1) = SICE*FRI(1) ! + SSNOW

        ELSE
          HSIL(1) = HSNOW*XSI(1)*MSI1/SNOW
          SSIL(1) = 0.          ! SSNOW*XSI(1)*MSI1/SNOW

        END IF
        HSIL(2) = HICE - HSIL(1) + HSNOW
        SSIL(2) = SICE - SSIL(1) ! + SSNOW

        HSIL(3) = (ROICE/ROICEN)*(HSIL(3)+FHSI4*FRI(3)-FHSI3)
        HSIL(4) = (ROICE/ROICEN)*(HSIL(4)+FHSI4*FRI(4)+FHSI3-FHSI4)
        SSIL(3) = (ROICE/ROICEN)*(SSIL(3)+FSSI4*FRI(3)-FSSI3)
        SSIL(4) = (ROICE/ROICEN)*(SSIL(4)+FSSI4*FRI(4)+FSSI3-FSSI4)

        ROICE = ROICEN          ! new ice concentration
      END IF
      END IF
      ELSE
C**** Ensure that MSI2 does not get too small for fixed-SST case.
        IF (ROICE.gt.0. and. MSI2.lt.AC2OIM) then
          DMIMP=AC2OIM-MSI2
          DHIMP=SUM(HSIL(3:4)*XSI(3:4)*DMIMP)
          DSIMP=SUM(SSIL(3:4)*XSI(3:4)*DMIMP)
          HSIL(3:4)=HSIL(3:4)*AC2OIM/MSI2
          SSIL(3:4)=SSIL(3:4)*AC2OIM/MSI2

          MSI2=AC2OIM
        END IF
      END IF
C**** Calculate temperatures for diagnostics and radiation
      TSIL(1:2)=((HSIL(1:2)-SSIL(1:2)*LHM)/(XSI(1:2)*MSI1)+LHM)*BYSHI
      TSIL(3:4)=((HSIL(3:4)-SSIL(3:4)*LHM)/(XSI(3:4)*MSI2)+LHM)*BYSHI

      RETURN
      END SUBROUTINE ADDICE


      SUBROUTINE SIMELT(DT,ROICE,SNOW,MSI2,HSIL,SSIL,POCEAN,Tm,TFO,TSIL, 1

     &     ENRGMAX,ENRGUSED,RUN0,SALT)
!@sum  SIMELT melts sea ice laterally and if it is too small
!@+    Note: all amounts are with respect to the ocean/lake fraction
!@auth Original Development Team
!@ver  1.1
      IMPLICIT NONE
!@var POCEAN ocean fraction (zero if lake)
      REAL*8, INTENT(IN) :: POCEAN
!@var TFO freezing temperature of water (C)
      REAL*8, INTENT(IN) :: TFO
!@var ROICE,SNOW,MSI2 ice variables (%,kg/m^2,kg/m^2)
      REAL*8, INTENT(INOUT) :: ROICE, SNOW, MSI2
!@var Tm mixed layer ocean temperature (C)
      REAL*8, INTENT(IN) :: Tm
!@var DT time step (s)
      REAL*8, INTENT(IN) :: DT
!@var HSIL ice enthalpy  (J/m^2)
!@var SSIL ice salt (kg/m^2)
      REAL*8, INTENT(INOUT), DIMENSION(LMI) :: HSIL, SSIL
!@var TSIL ice temperature (C)
      REAL*8, INTENT(OUT), DIMENSION(LMI) :: TSIL
!@var ENRGMAX max energy available for melting ice (J/m^2)
      REAL*8, INTENT(IN) :: ENRGMAX
!@var ENRGUSED energy used to melt ice (J/m^2)
      REAL*8, INTENT(OUT) :: ENRGUSED
!@var RUN0 amount of sea ice melt (kg/m^2)
!@var SALT amount of salt in sea ice melt (kg/m^2)
      REAL*8, INTENT(OUT) :: RUN0, SALT
!@var DRSI change in RSI fraction
!@var DTEMP temperature diff. used in lateral melt calculation.
      REAL*8 :: DRSI,DTEMP


C**** Estimate DRSI
      DRSI=0.
      IF (ROICE.lt.1d-3) THEN
        DRSI=ROICE
      ELSE      ! IF (POCEAN.gt.0) THEN ! now for lakes too
C**** Estimate lateral melt using parameterisation from Makyut/Steele
C**** (via C. Bitz): Rside=dt*pi/(floesize*eta)*(3d-6)*(delT)^(1.36)
        dtemp=MAX(Tm-TFO,0d0)
        DRSI=DT*SILMFAC*dtemp**SILMPOW
        IF (ROICE-DRSI.lt.1d-3) DRSI=ROICE
        IF (ENRGMAX+DRSI*SUM(HSIL).lt.0) DRSI=-ENRGMAX/SUM(HSIL)
      END IF
C**** Remove DRSI amount of ice
      ENRGUSED=-DRSI*SUM(HSIL) !  [J/m^2]
      RUN0=DRSI*(SNOW + ACE1I + MSI2)
      SALT=DRSI*SUM(SSIL)

      ROICE=ROICE-DRSI
      IF (ROICE.lt.1d-10) THEN
        ROICE=0.                ! deal with possible round off err
C**** set defaults if no ice is left
        SNOW=0.
        MSI2=AC2OIM
        IF (POCEAN.gt.0) THEN
          SSIL(1:2)=SSI0*XSI(1:2)*ACE1I
          SSIL(3:4)=SSI0*XSI(3:4)*AC2OIM
        ELSE
          SSIL(:) = 0.
        END IF
        HSIL(1:2)=(SHI*TFO-LHM)*XSI(1:2)*ACE1I+LHM*SSIL(1:2)
        HSIL(3:4)=(SHI*TFO-LHM)*XSI(3:4)*AC2OIM+LHM*SSIL(3:4)
        TSIL=TFO

      END IF
C****
      RETURN
      END SUBROUTINE SIMELT


      SUBROUTINE SSIDEC(SNOW,MSI2,HSIL,SSIL,DT, 1

     &     MFLUX,HFLUX,SFLUX)
!@sum  SSIDEC decays salinity in sea ice
!@auth Jiping Liu
!@ver  1.0
      IMPLICIT NONE

!@var HSIL enthalpy of ice layers (J/m^2)
!@var SSIL salt in ice layers (kg/m^2)
      REAL*8, INTENT(INOUT), DIMENSION(LMI) :: SSIL,HSIL
!@var MSI2 second mass layer ice thickness (kg/m^2)
!@var MSI1 first mass layer ice thickness (kg/m^2)
!@var DT source time step (s)
      REAL*8, INTENT(INOUT) :: MSI2
      REAL*8, INTENT(IN) :: DT, SNOW
!@var MFLUX,SFLUX,HFLUX mass, salt and heat flux arising from
!@+   sea salinity decay
      REAL*8, INTENT(OUT) :: MFLUX,HFLUX,SFLUX


      REAL*8 DSSI(3:LMI),DMSI(3:LMI),DHSI(3:LMI),TSIL(3:LMI),DS12
     *     ,DM12,DH12
      REAL*8 FMSI1,FMSI2,FMSI3,FSSI1,FSSI2,FSSI3,FHSI1,FHSI2,FHSI3
      REAL*8 HSNOW,HICE,SICE,TICE,MSI1
      INTEGER L
!@var dtssi decay time scale for sea ice salinity (days)
!@var bydtssi decay constant for sea ice salinity (1/s)
      REAL*8, parameter :: dtssi=30.d0, bydtssi=1./(dtssi*sday)

      DSSI = 0. ; DMSI = 0. ; DHSI = 0.
      DS12 = 0. ; DM12 = 0. ; DH12 = 0.

C**** salinity in sea ice decays if it is above threshold level ssi0

C**** check first layer (default ice and snow)
C**** New formulation:
C**** Split all upper mass layer fields into ice and snow components
      MSI1=SNOW+ACE1I
      IF (ACE1I.gt.XSI(2)*MSI1) THEN ! some ice in first layer
        HSNOW = (HSIL(1)-LHM*SSIL(1))*SNOW/(XSI(1)*MSI1)
      ELSE
        HSNOW = HSIL(1) + (HSIL(2)-LHM*SSIL(2))*(XSI(2)*MSI1-ACE1I)
     *       /(XSI(2)*MSI1)
      END IF
      HICE  = HSIL(1)+HSIL(2)-HSNOW
CC    SSNOW = 0.   ! always zero
      SICE  = SSIL(1)+SSIL(2)
      TICE = ((HICE-LHM*SICE)/ACE1I+LHM)*BYSHI


C**** calculate removal of excess salinity
      IF (SICE.GT.ssi0*ACE1I) THEN
        DS12 = (SICE-ACE1I*ssi0)*DT*BYDTSSI
        DM12 = DS12
        DH12 = -DM12*TICE*SHI

        SICE=SICE-DS12
        HICE=HICE-DH12
      END IF

C**** check remaining layers
      DO L=3,LMI
        IF (SSIL(L).GT.ssi0*XSI(L)*MSI2) THEN
          DSSI(L) = (SSIL(L)-ssi0*XSI(L)*MSI2)*DT*BYDTSSI
          DMSI(L) = DSSI(L)
          TSIL(L) = ((HSIL(L)-LHM*SSIL(L))/(XSI(L)*MSI2)+LHM)*BYSHI
          DHSI(L) = -DMSI(L)*TSIL(L)*SHI

          SSIL(L) = SSIL(L) - DSSI(L)
          HSIL(L) = HSIL(L) - DHSI(L)
        END IF
      END DO

C**** Calculate fluxes required to rebalance layers
C**** Mass/heat/salt moves from layer 3 to 2
      FMSI2 = -DM12
      FHSI2 = FMSI2*HSIL(3)/(XSI(3)*MSI2-DMSI(3))
      FSSI2 = FMSI2*SSIL(3)/(XSI(3)*MSI2-DMSI(3))


C**** reconstitute upper layers
      IF (DM12.gt.0) THEN
        IF (ACE1I.gt.XSI(2)*MSI1) THEN ! some ice in first layer
          HSIL(1) = (HICE-FHSI2)*(ACE1I-XSI(2)*MSI1)/ACE1I + HSNOW
          SSIL(1) = (SICE-FSSI2)*(ACE1I-XSI(2)*MSI1)/ACE1I ! + SSNOW

        ELSE
          HSIL(1) = HSNOW*XSI(1)*MSI1/(MSI1-ACE1I)
          SSIL(1) = 0.          ! SSNOW*XSI(1)*MSI1/(MSI1-ACE1I)

        END IF
        HSIL(2) = HICE - HSIL(1) - FHSI2 + HSNOW
        SSIL(2) = SICE - SSIL(1) - FSSI2 ! + SSNOW

      END IF

C**** Mass/heat moves between layers 3 to 4
      FMSI3 = XSI(3)*DMSI(4)+XSI(4)*(FMSI2-DMSI(3))
      IF (FMSI3.gt.0) THEN      ! downward flux to layer 4
        FHSI3 = FMSI3*HSIL(3)/(XSI(3)*MSI2-DMSI(3))
        FSSI3 = FMSI3*SSIL(3)/(XSI(3)*MSI2-DMSI(3))

      ELSE                      ! upward flux
        FHSI3 = FMSI3*HSIL(4)/(XSI(4)*MSI2-DMSI(4))
        FSSI3 = FMSI3*SSIL(4)/(XSI(4)*MSI2-DMSI(4))

      END IF

C**** Apply the second mass layer fluxes
      HSIL(3)=HSIL(3)+(FHSI2-FHSI3)
      HSIL(4)=HSIL(4)+ FHSI3
      SSIL(3)=SSIL(3)+(FSSI2-FSSI3)
      SSIL(4)=SSIL(4)+ FSSI3

c     MSI1 = MSI1 - (DMSI(1)+DMSI(2)) - FMSI2 ! stays fixed
      MSI2 = MSI2 - (DMSI(3)+DMSI(4)) + FMSI2

C**** output fluxes and diagnostics
      MFLUX = DM12+SUM(DMSI(3:LMI))         ! mass flux to ocean
      HFLUX = DH12+SUM(DHSI(3:LMI))         ! energy flux to ocean
      SFLUX = DS12+SUM(DSSI(3:LMI))         ! salt flux to ocean

C****
      RETURN
      END SUBROUTINE SSIDEC


      subroutine iceocean(Ti,Si,Tm,Sm,dh,ustar,Coriol,dtsrc,mlsh, 1,1

     &     mflux,sflux,hflux)
!@sum  iceocean calculates fluxes at base of sea ice
!@auth Gavin Schmidt
!@ver  1.0
!@usage At the ice-ocean interface the following equations hold:
!@+                                              Tb = -mu Sb         (1)
!@+ -lam_i(Ti,Si)(Ti-Tb)/dh + rho_m shw g_T (Tb-Tm) = -m Lh(Tib,Sib)
!@+                                                  + m shw (Tib-Tb)(2)
!@+                           rho_m     g_S (Sb-Sm) =  m (Sib-Sb )   (3)
!@+  with  Sib=Si, Tib=Ti m>0,  or Sib = fsss Sb, Tib=Tb m<0        (4)
!@+ This routine iteratively solves for m, Tb, Sb and Sib, Tib using
!@+ Newton's method. The form of the latent heat and conductivity
!@+ can optionally include salinity effects.
!@+ The actual fluxes (positive down) are then:
!@+  mflux = m   ; sflux = 1d-3 m Sib (kg/m^2/s)
!@+  hflux = lam_i(Ti,Si) (Ti-Tb)/dh -m Lh(Tib,Sib)+m shw Tib  (J/m^2/s)
!@+
      implicit none
!@var alamdS salinity coefficient for conductivity (from common?)
      real*8, parameter :: alamdS=0.117d0

!@var G_mole_T,G_mole_S molecular diffusion terms for heat/salt
C****  G_mole_T = 12.5 * 13.8d0**(2d0/3d0) - 6. = 65.9d0
C****  G_mole_S = 12.5 * 2432d0**(2d0/3d0) - 6. = 2255d0
      real*8, parameter :: G_mole_T = 65.9d0 , G_mole_S = 2255d0
!@var Si,Sm salinity in lowest ice layer and mixed layer (psu)
!@var Ti,Tm temperatures in ice and mixed layer (C)
!@var dh distance from center of bottom ice layer to base of ice (m)
      real*8, intent(in) :: Ti,Si,Tm,Sm,dh
!@var coriol Coriolis parameter (used in turbulent flux calc) (1/s)
!@var ustar friction velocity at ice-ocean interface (m/s)
      real*8, intent(in) :: ustar,Coriol
!@var dtsrc source time step (s)
!@var mlsh mixed layer specific heat capactity (J/m^2 C) (UNUSED)
      real*8, intent(in) :: dtsrc,mlsh
!@var mflux,sflux,hflux mass, salt and heat fluxes at base of ice
      real*8, intent(out) :: mflux,sflux,hflux


!@var g_T,g_S turbulent exchange velocities (m/s)
!@var G_turb turbulent diffusion term
!@var Sb,Sb0 final and initial estimate for salinity at interface (psu)
      real*8 :: G_turb,g_T,g_S,Sb,Sb0
!@var dSbdSb differential of Sb with respect to initial Sb0
      real*8 :: dSbdSb
!@var m melt rate (negative implies freezing) (kg/m^2/s)
!@var Tb temperature at interface (C)
!@var Tib,Sib actual temperature/salinity in melting/freezing ice (psu)
      real*8 :: m,Tb,Sib  ! ,Tib
      real*8 :: lh,left2,df3dm,dmdTb,dmdSi,alamdh,f0,df,rsg
      integer, parameter :: niter=5 !@param niter number of iterations
      integer :: i

C**** Calculate turbulent exchange velocities g_T, g_S

C**** G_turb = (1/k) ln (u* E n^2 / f h) + 1 / (2 E n) - (1/k)
C****        = 2.5 ln ( 5300*(u*)^2/coriol ) + 9.62 - 2.5 (assuming n=1)
C**** Note: n should theoretically depend on the buoyancy flux and
C****       therfore this should be included in the iteration below.
      G_turb = 2.5d0 * log ( 5300.*ustar*ustar/coriol ) + 7.12d0
C****  g  = u* / ( G_turb + G_mole )
      g_T = ustar / ( G_turb + G_mole_T )
      g_S = ustar / ( G_turb + G_mole_S )

C**** set conductivity term
C**** Diffusive flux is implicit in ice !  + ml temperature (not used)
      rsg=rhows*shw*g_T    ! /(1.+alpha*dtsrc*rhows*shw*g_T/mlsh)
c no salinity effects
      alamdh = alami/(dh+alpha*dtsrc*alami*byshi/(2.*dh*rhoi))
c S thermo
c      alamdh = (alami+alamdS*Si)/(dh+alpha*dtsrc*
c     *          (alami+alamdS*Si)*byshi/(2.*rhoi*dh))

C**** solve for boundary values (uses Newton-Rapheson with bounds)
C**** Estimate initial boundary salinity
      Sb0 = 0.75d0*Sm+0.25d0*Si
      do i=1,niter
C**** we use the fill tfrez calc here, but assume for simplicity that
C**** the gradient w.r.t. T is still = -mu below.
        Tb = tfrez(Sb0)   ! -mu*Sb0  ! freezing point at salinity Sb0
C**** calculate left hand side of equation 2
        left2 = -alamdh*(Ti-Tb) + rsg*(Tb-Tm)
C**** depending on whether there is freezing or melting, the sea ice
C**** salinity in the melt/freeze fraction is set to be the original
C**** ice salinity or is a function of the boundary salinity.
C**** Similarly for temperature Tib
        if (left2.gt.0) then    ! freezing
          if (qsfix) then   ! keep salinity in ice constant
            Sib = ssi0*1d3
          else              ! it is a function of boundary value
            Sib = Sb0*fsss
          end if
c no salinity effects
c         lh = lhm + Tb*(shw-shi)
c salinity effects only mass
          lh = lhm*(1.-Sib*1d-3) + Tb*(shw-shi)
c S thermo
c         lh = lhm*(1.+mu*Sib/Tb) + (Tb+mu*Sib)*(shw-shi)
          m = -left2/lh

c no salinity effects
c          dmdTb = (left2*(shw-shi)-lh*(alamdh + rsg))/(lh*lh)
c          dmdSi = 0.
c salinity effects only mass
          dmdTb = (left2*(shw-shi)-lh*(alamdh + rsg))/(lh*lh)
          dmdSi = -left2*lhm*1d-3/(lh*lh)
c S thermo
c         dmdTb = (left2*(-lhm*mu*Sib/Tb**2+shw-shi)-lh*(alamdh + rsg
c    *       ))/(lh*lh)
c         dmdSi = (left2*mu*(lhm/Tb+shw-shi))/(lh*lh)

          if (qsfix) then   ! keep salinity in ice constant
            Sb = (rhows*g_S*Sm+m*Sib)/(rhows*g_S+m)
            df3dm= rhows*g_S*(Sib-Sm)/(rhows*g_S+m)**2
            dSbdSb= -mu*dmdTb*df3dm
          else              ! it is a function of boundary value
            Sb = rhows*g_S*Sm/(rhows*g_S+m*(1.-fsss))
            df3dm= -rhows*g_S*Sm*(1.-fsss)/(rhows*g_S+m*(1.-fsss))**2
            dSbdSb= (fsss*dmdSi-mu*dmdTb)*df3dm
          end if
        else                    ! melting
          Sib = Si
c no salinity effects
c         lh = lhm + Tb*shw - Ti*shi
c salinity effects only mass
          lh = lhm*(1.-Sib*1d-3) + Tb*shw - Ti*shi
c S thermo
c         lh = lhm*(1.+mu*Sib/Ti) + (Ti+mu*Sib)*(shw-shi) - shw*(Ti-Tb)
          m = -left2/lh

          Sb = (m*Sib+rhows*g_S*Sm)/(rhows*g_S+m)
          df3dm= (Sib*(rhows*g_S+m)-(m*Sib+rhows*g_S*Sm))/(rhows*g_S+m)
     *         **2
          dmdTb = -(alamdh + rsg)/lh + left2*shw/lh**2
          dSbdSb= -mu*dmdTb*df3dm
        end if

        f0=Sb-Sb0
        df=dSbdSb-1.+1d-20
        Sb = min(max(Sib,Sb0 - f0/df),40d0)
        Sb0=Sb
      end do
C**** define fluxes (positive down)
C**** Cap mass flux at at 90% of bottom layer
      if (m.gt.0.9d0*2.*dh*rhoi/dtsrc) then
        m=0.9d0*2.*dh*rhoi/dtsrc
      end if

      mflux = m                       ! (kg/m^2 s)
      sflux = 1d-3*m*Sib              ! (kg/m^2 s)
      hflux = alamdh*(Ti-Tb) - m*lh +m*shw*Tb ! (J/m^2 s)

C****
      return
      end subroutine iceocean


      subroutine icelake(Ti,Tm,dh,dtsrc,mlsh, 1

     &     mflux,hflux)
!@sum  icelake calculates fluxes at base of lake ice (no salinity)
!@auth Gavin Schmidt
!@ver  1.0
!@usage At the ice-lake interface the following equations hold:
!@+   interface is at freezing point (Tb=0.)   (1)
!@+   -lam_i Ti/dh - rho_m shw g_T Tm = -m Lh(Tib) + m shw Tib   (2)
!@+     with  Tib=Ti m>0,  or Tib=0. m<0        (4)
      implicit none
!@var mflux,hflux mass and heat fluxes at base of ice
      real*8, intent(out) :: mflux,hflux
!@var dtsrc source time step (s)
!@var Ti,Tm temperatures in ice and upper lake layer (C)
!@var dh distance from center of bottom ice layer to base of ice (m)
!@var mlsh mixed layer specific heat capactity (J/m^2 C) (UNUSED)
      real*8, intent(in) :: Ti,Tm,dh,dtsrc,mlsh
C**** Assume constant g_T = 1.3d-7, g_S = 0.025 * g_T m/s
!@var rsg = rhow * shw * g_T turbulent energy flux (J/m^2 s)
!@var rgS = rhow * g_S turbulent tracer flux (kg/m^2 s)
      real*8, parameter ::  rsg = rhow*shw*1.3d-7, rgS=rhow*3.2d-9
      real*8 left2, lh, m, alamdh


C**** Diffusive flux is implicit in ice temperature
      alamdh = alami/(dh+alpha*dtsrc*byshi*alami/(2.*dh*rhoi))
C**** calculate left hand side of equation 2
      left2 = -alamdh*Ti - rsg*Tm    !/(1.+alpha*dtsrc*rsg/mlsh)
      if (left2.gt.0) then      ! freezing
        lh = lhm
      else                      ! melting
        lh = lhm -  Ti*shi
      end if
C**** define fluxes (positive down)
      m = -left2/lh
C**** Cap mass flux at 90% of bottom layer
      if (m.gt.0.9d0*2.*dh*rhoi/dtsrc) then
        m=0.9d0*2.*dh*rhoi/dtsrc
      end if
      mflux = m                 ! (kg/m^2 s)
      hflux = alamdh*Ti - m*lh  ! (J/m^2 s)

C****
      return
      end subroutine icelake


      subroutine solar_ice_frac(snow,msi2,wetsnow,fsri,lmax) 2
!@sum solar_ice_frac calculates the fraction of solar radiation
!@+   transmitted through the ice, taking albedo into account, but
!@+   assuming constant proportions of visible and near-ir.
!@+   Based on Ebert et al, 1995.
!@auth Gavin Schmidt
      implicit none
!@var kiextvis,kiextnir1 ice extinction coeff. for vis and nir (1/m)
      real*8, parameter :: kiextvis=1.5d0, kiextnir1=18.d0
!@var snow,msi2 snow and second layer ice amounts (kg/m^2)
      real*8, intent(in) :: snow,msi2
!@var wetsnow =true if snow is wet
      logical, intent(in) :: wetsnow
!@var lmax number of ice layers to do calculation for
      integer, intent(in) :: lmax
!@var fsri fraction of solar energy reaching bottom of layer
      real*8, dimension(lmax), intent(out) :: fsri
      real*8, dimension(lmax) :: fsrvis, fsrnir1, hice
!@var ksextvis,ksextnir1 snow extinction coeff. for vis and nir (1/m)
!@var fracvis,fracnir fraction of absorbed solar in vis and nir
      real*8 fracvis,fracnir1,ksextvis,ksextnir1
      real*8 hsnow,hsnow1,hice12
      integer l

      hsnow = snow/rhos
      hice12= ace1i*byrhoi

C**** old code
c!@param KIEXT extinction coefficient for light in sea ice (1/m)
c      REAL*8, PARAMETER :: KIEXT = 1.5d0
c!@param KSEXT extinction coefficient for light in snow (1/m)
c      REAL*8, PARAMETER :: KSEXT = 15d0
c        IF (ACE1I*XSI(1).gt.SNOW*XSI(2)) THEN
cC**** first thermal layer is snow and ice
c          HICE1 = (ACE1I-XSI(2)*(SNOW+ACE1I))*BYRHOI
c          FSRI(1) = EXP(-KSEXT*HSNOW-KIEXT*HICE1)
c        ELSE ! all snow
c          HSNOW1 = (ACE1I+SNOW)*XSI(1)/RHOS
c          FSRI(1) = EXP(-KSEXT*HSNOW1)
c        END IF
c        FSRI(2) = EXP(-KSEXT*HSNOW-KIEXT*HICE12)
c        FSRI(3) = FSRI(2)*EXP(-KIEXT*MSI2*XSI(3)*BYRHOI)
c        FSRI(4) = FSRI(2)*EXP(-KIEXT*MSI2*BYRHOI)

c**** Set fractions of visible and near-ir bands based on weighting the
c**** solar input by the approximate co-albedo (ignoring melt ponds etc)
c**** VIS:  250 - 690 nm,       49.3% of incoming solar at ground
c**** NIR1: 690 - 1190 nm       34.9%       "       "   "    "
c**** NIR2/3 (> 1190 nm assumed not to be transmitted)
c****
      if (hsnow.gt.0.02) then    ! same cutoff as for albedo
c**** Extinction coefficients for snow depend on wetness
c**** Since albedo does as well, absorbed fractions also vary
        if (wetsnow) then
          fracvis =0.20d0 ;  fracnir1=0.33d0
          ksextvis=10.7d0 ; ksextnir1=118.d0
        else
          fracvis =0.06d0 ;  fracnir1=0.31d0
          ksextvis=19.6d0 ; ksextnir1=196.d0
        end if
      else
        fracvis=0.24d0 ; fracnir1=0.43d0
        ksextvis=10.7d0 ; ksextnir1=118.d0
      end if

C**** calculate sucessive fractions assuming each band is attenuated
C**** independently.
      if (ace1i*xsi(1).gt.snow*xsi(2)) then
c**** First thermal layer is snow and ice
        hice(1)    = (ace1i-xsi(2)*(snow+ace1i))*byrhoi
        fsrvis (1) = exp(-ksextvis *hsnow - kiextvis *hice(1))
        fsrnir1(1) = exp(-ksextnir1*hsnow - kiextnir1*hice(1))
        fsri(1)    = fracvis*fsrvis(1) + fracnir1*fsrnir1(1)
      else                      ! all snow
        hsnow1     = (ace1i+snow)*xsi(1)/rhos
        fsrvis (1) = exp(-ksextvis *hsnow1)
        fsrnir1(1) = exp(-ksextnir1*hsnow1)
        fsri(1)    = fracvis*fsrvis(1) + fracnir1*fsrnir1(1)
      end if
      fsrvis (2) = exp(-ksextvis *hsnow - kiextvis *hice12)
      fsrnir1(2) = exp(-ksextnir1*hsnow - kiextnir1*hice12)
      fsri(2)    = fracvis*fsrvis(2) + fracnir1*fsrnir1(2)

      do l=3,lmax   ! only done if lmax>2
        hice(l)    = xsi(l)*msi2*byrhoi
        fsrvis (l) = exp(-kiextvis *hice(l)) * fsrvis (l-1)
        fsrnir1(l) = exp(-kiextnir1*hice(l)) * fsrnir1(l-1)
        fsri(l)    = fracvis*fsrvis(l) + fracnir1*fsrnir1(l)
      end do
c****
      return
      end subroutine solar_ice_frac


      function tfrez(sss,press) 6
!@sum tfrez calculates freezing temperature of sea water
!@auth Gavin Schmidt
      implicit none
!@var sss sea surface salinity (psu)
      real*8, intent(in) :: sss
!@var press gauge pressure (default=0) (Pa)
      real*8, intent(in), optional :: press
!@var tfrez approx. freezing point of sea water (C)
      real*8 tfrez,pr
      real*8, parameter :: dtdp = -7.5d-8
      real*8 :: a01 = -.0575d0, a02 = -2.154996d-4, a03 =1.710523d-3

      pr=0.
      if (present(press)) then
        pr=press
      end if
C**** linear approximation
c      tfrez = -mu*sss + dtdp*pr
C**** UNESCO formula (1983)
      tfrez = (a01 + a02*sss)*sss + a03*sss*sqrt(sss)+ dtdp*pr

      return
      end function tfrez


      subroutine snowice(Tm,Sm,SNOW,MSI2,HSIL,SSIL,qsfix, 1,1

     *         MSNWIC,HSNWIC,SSNWIC)
      IMPLICIT NONE
!@var Tm ocean mixed layer temperature (C)
!@var Sm ocean mixed layer salinity (psu)
      REAL*8, INTENT(IN) :: Tm, Sm
!@var MSI1, MSI2 snow and ice mass in two mass layers (kg/m^2)
      REAL*8, INTENT(INOUT) :: SNOW,MSI2
!@var HSIL, SSIL energy and salt profiles in ice (J, kg)/m^2
      REAL*8, DIMENSION(LMI), INTENT(INOUT) :: HSIL,SSIL
!@var MSNWIC, HSNWIC, SSNWIC mass, energy and salt flux to ocean (<0) 
      REAL*8, INTENT(OUT) :: MSNWIC, HSNWIC, SSNWIC
!@var qsfix true if sea ice salinity is fixed
      LOGICAL, INTENT(IN) :: qsfix

      REAL*8 DSNOW,Z0,MAXM,MAXME,Eoc,Esnow,Ei,Erat,Si,MSI1
      REAL*8 SICE,HICE,HSNOW,FMSI2,FMSI3,FHSI2,FSSI2,FHSI3,FSSI3
      INTEGER L

C**** test for snow ice possibility      
      IF (RHOI*SNOW.gt.(ACE1I+MSI2)*(RHOWS-RHOI)) THEN
        MSI1=SNOW+ACE1I
        Z0=(MSI1+MSI2)/RHOWS-(ACE1I+MSI2)/RHOI
        MAXM=Z0*RHOWS*(RHOI-RHOS)/(RHOWS+RHOS-RHOI)

C**** heat, salt, tracer amounts in snow, first layer ice
        IF (ACE1I.gt.XSI(2)*MSI1) THEN ! some ice in first layer
          HSNOW = (HSIL(1)-LHM*SSIL(1))*(MSI1-ACE1I)/(XSI(1)*MSI1)
        ELSE
          HSNOW = HSIL(1) + (HSIL(2)-LHM*SSIL(2))*(XSI(2)*MSI1-ACE1I)
     *         /(XSI(2)*MSI1)
        END IF
        HICE  = HSIL(1)+HSIL(2)-HSNOW
        SICE  = SSIL(1)+SSIL(2)


C**** Check energy constraint on mass flux
        if (qsfix) then
          Si=1d3*ssi0
        else
          Si=fsss*Sm
        end if
c        Ei=EICE(TFREZ(Sm),Si) 
        Ei=-LHM*(1.-0.001*Si)+ Tfrez(Sm)*shi
        Esnow=HSNOW/SNOW
        Eoc=Tm*SHW
        ERAT=(Esnow-Ei)/(Ei-Eoc)

        MAXME=Z0*RHOI*ERAT/(1.+ERAT*(RHOWS-RHOI)/RHOWS)
        MAXM=MAX(0d0,MIN(MAXME,MAXM))   ! mass of sea water to be added

        DSNOW=Z0*RHOI-MAXM*(RHOWS-RHOI)/RHOWS  ! loss of snow

        FMSI2 = MAXM+DSNOW      ! mass of ice pushed down to layer 3

C**** Add new ice 
        SICE=SICE+MAXM*0.001d0*Si
        HICE=HICE+MAXM*Eoc+DSNOW*Esnow

        FSSI2 = FMSI2*SICE/(ACE1I+FMSI2)
        FHSI2 = FMSI2*HICE/(ACE1I+FMSI2)
C**** Update ice and snow concentrations

        SICE=SICE-FSSI2
        HICE=HICE-FHSI2
        HSNOW=HSNOW*(SNOW-DSNOW)/SNOW   ! snow temp remains const
        SNOW=SNOW-DSNOW
        MSI1=ACE1I+SNOW

C**** Reconstitute upper layers
        IF (ACE1I.gt.XSI(2)*MSI1) THEN ! some ice in first layer
          HSIL(1) = HICE*(ACE1I-XSI(2)*MSI1)/ACE1I + HSNOW
          SSIL(1) = SICE*(ACE1I-XSI(2)*MSI1)/ACE1I ! + SSNOW

        ELSE
          HSIL(1) = HSNOW*XSI(1)*MSI1/(MSI1-ACE1I)
          SSIL(1) = 0.          ! SSNOW*XSI(1)*MSI1/(MSI1-ACE1I)

        END IF
        HSIL(2) = HICE - HSIL(1) + HSNOW
        SSIL(2) = SICE - SSIL(1) ! + SSNOW

C**** lower level fluxes
        FMSI3 = FMSI2*XSI(4)    ! mass of ice pushed down to layer 4
        FSSI3 = FMSI3*(SSIL(3)+FSSI2)/(XSI(3)*MSI2+FMSI2)
        FHSI3 = FMSI3*(HSIL(3)+FHSI2)/(XSI(3)*MSI2+FMSI2)


C**** Update lower layers
        MSI2=MSI2+FMSI2
        SSIL(3)=SSIL(3)+FSSI2-FSSI3
        SSIL(4)=SSIL(4)      +FSSI3
        HSIL(3)=HSIL(3)+FHSI2-FHSI3
        HSIL(4)=HSIL(4)      +FHSI3


C**** output flux (positive down)
        MSNWIC = -MAXM
        HSNWIC = MSNWIC*Eoc
        SSNWIC = 0.001d0*Si*MSNWIC

      ELSE
        MSNWIC = 0. ; HSNWIC = 0. ; SSNWIC = 0.

      END IF
      RETURN
      end subroutine snowice

      END MODULE SEAICE


      MODULE SEAICE_COM 38,3
!@sum  SEAICE_COM contains the model arrays for seaice
!@auth Gavin Schmidt
!@ver  1.0
      USE MODEL_COM, only : im,jm
      USE DOMAIN_DECOMP, only : grid

      USE SEAICE, only : lmi

      IMPLICIT NONE
!@var RSI fraction of open water area covered in ice
      REAL*8, ALLOCATABLE, DIMENSION(:,:) :: RSI
!@var SNOWI snow amount on sea ice (kg/m^2)
      REAL*8, ALLOCATABLE, DIMENSION(:,:) :: SNOWI
!@var MSI mass of ice second layer (layer 1=const) (kg/m^2)
C**** Note that MSI includes the mass of salt in sea ice
      REAL*8, ALLOCATABLE, DIMENSION(:,:) :: MSI
!@var HSI enthalpy of each ice layer (J/m^2)
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: HSI
!@var SSI sea ice salt content (kg/m^2)
      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: SSI
!@var pond_melt amount of melt pond mass (kg/m^2)
C**** Note this is a virtual meltpond and is only used for
C**** albedo calculations
      REAL*8, ALLOCATABLE, DIMENSION(:,:) :: pond_melt
!@var flag_dsws true if snow on ice is wet (ie. rain or surface melt)
      LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: flag_dsws



      END MODULE SEAICE_COM


      SUBROUTINE ALLOC_SEAICE_COM(grid) 1,6
!@sum  To allocate arrays whose sizes now need to be determined at 
!@+    run-time
!@auth Rodger Abel
!@ver  1.0
      USE DOMAIN_DECOMP, ONLY : DYN_GRID
      USE DOMAIN_DECOMP, ONLY : GET
      USE MODEL_COM, ONLY : IM, JM

      USE SEAICE, only : LMI

      USE SEAICE_COM, ONLY : RSI, SNOWI, MSI, HSI, SSI, pond_melt, 
     *     flag_dsws

      IMPLICIT NONE
      TYPE (DYN_GRID), INTENT(IN) :: grid

      INTEGER :: J_1H, J_0H
      INTEGER :: IER

      CALL GET(grid, J_STRT_HALO=J_0H, J_STOP_HALO=J_1H)

      ALLOCATE( RSI(IM, J_0H:J_1H),
     *     SNOWI(IM, J_0H:J_1H),
     *     MSI(IM, J_0H:J_1H),
     *     pond_melt(IM, J_0H:J_1H),
     *     flag_dsws(IM, J_0H:J_1H),
     *     STAT=IER)

      ALLOCATE( HSI(LMI, IM, J_0H:J_1H),
     *     SSI(LMI, IM, J_0H:J_1H),
     *     STAT=IER)



      RETURN
      END SUBROUTINE ALLOC_SEAICE_COM


      SUBROUTINE io_seaice(kunit,iaction,ioerr) 2,2
!@sum  io_seaice reads and writes seaice variables to file
!@auth Gavin Schmidt
!@ver  1.0
      USE MODEL_COM, only : ioread,iowrite,lhead,irsfic,irsficno,irerun
      USE SEAICE_COM
      IMPLICIT NONE

      INTEGER kunit   !@var kunit unit number of read/write
      INTEGER iaction !@var iaction flag for reading or writing to file
!@var IOERR 1 (or -1) if there is (or is not) an error in i/o
      INTEGER, INTENT(INOUT) :: IOERR
!@var HEADER Character string label for individual records
      CHARACTER*80 :: HEADER, MODULE_HEADER = "SICE02"


      write (MODULE_HEADER(lhead+1:80),'(a14,i1,a20,i1,a23)')
     * 'R8 F(im,jm),H(',lmi,',im,jm),snw,msi,ssi(',lmi,
     *     '),pond_melt,L flag_dsws'

      SELECT CASE (IACTION)
      CASE (:IOWRITE)            ! output to standard restart file
         WRITE (kunit,err=10) MODULE_HEADER,RSI,HSI,SNOWI,MSI,SSI
     *     ,POND_MELT,FLAG_DSWS

      CASE (IOREAD:)            ! input from restart file
         READ (kunit,err=10) HEADER,RSI,HSI,SNOWI,MSI,SSI
     *     ,POND_MELT,FLAG_DSWS
        IF (HEADER(1:LHEAD).NE.MODULE_HEADER(1:LHEAD)) THEN
          PRINT*,"Discrepancy in module version ",HEADER,MODULE_HEADER
          GO TO 10
        END IF

      END SELECT

      RETURN
 10   IOERR=1
      RETURN
      END SUBROUTINE io_seaice


      SUBROUTINE CHECKI(SUBR) 1,16
!@sum  CHECKI Checks whether Ice values are reasonable
!@auth Original Development Team
!@ver  1.0
      USE CONSTANT, only : lhm,shi
      USE MODEL_COM
      USE GEOM, only : imaxj

      USE SEAICE, only : lmi,xsi,ace1i
      USE SEAICE_COM, only : rsi,msi,hsi,snowi,ssi

      USE LAKES_COM, only : flake
      USE FLUXES
      USE DOMAIN_DECOMP, only : GRID
      USE DOMAIN_DECOMP, only : GET
      IMPLICIT NONE

!@var SUBR identifies where CHECK was called from
      CHARACTER*6, INTENT(IN) :: SUBR
!@var QCHECKI true if errors found in seaice
      LOGICAL QCHECKI
      INTEGER I,J,L
      REAL*8 TICE


      integer :: J_0, J_1
C**** 
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT = J_0, J_STOP = J_1)

C**** Check for NaN/INF in ice data
      CALL CHECK3(RSI,IM,JM,1,SUBR,'rsi')
      CALL CHECK3(MSI,IM,JM,1,SUBR,'msi')
      CALL CHECK3(HSI,LMI,IM,JM,SUBR,'hsi')
      CALL CHECK3(SSI,LMI,IM,JM,SUBR,'ssi')
      CALL CHECK3(SNOWI,IM,JM,1,SUBR,'sni')

      QCHECKI = .FALSE.
C**** Check for reasonable values for ice variables
      DO J=J_0, J_1
        DO I=1,IMAXJ(J)
          IF (RSI(I,J).lt.0 .or. RSI(I,j).gt.1 .or. MSI(I,J).lt.0) THEN
            WRITE(6,*) 'After ',SUBR,': I,J,RSI,MSI=',I,J,RSI(I,J)
     *           ,MSI(I,J)
            QCHECKI = .TRUE.
          END IF
          IF (          DO L=1,LMI
            IF (L.le.2) TICE = ((HSI(L,I,J)-SSI(L,I,J)*LHM)/(XSI(L)
     *           *(ACE1I+SNOWI(I,J)))+LHM)/SHI
            IF (L.gt.2) TICE = ((HSI(L,I,J)-SSI(L,I,J)*LHM)/(XSI(L)
     *           *MSI(I,J))+LHM)/SHI
            IF (HSI(L,I,J).gt.0.or.TICE.gt.1d-10.or.TICE.lt.-80.) THEN
              WRITE(6,'(3a,3i3,6e12.4/1X,6e12.4)')
     *             'After ',SUBR,': I,J,L,TSI=',I,J,L,TICE,RSI(I,J)
              WRITE(6,*) HSI(:,I,J),MSI(I,J),SNOWI(I,J),SSI(:,I,J)
              IF (TICE.gt.1d-3.or.TICE.lt.-100.) QCHECKI = .TRUE.
            END IF
            IF (SSI(L,I,J).lt.0) THEN
              WRITE(6,*) 'After ',SUBR,': I,J,L,SSI=',I,J,L,SSI(:,I
     *             ,J),MSI(I,J),SNOWI(I,J),RSI(I,J)
              QCHECKI = .TRUE.
            END IF
          END DO
          IF (SNOWI(I,J).lt.0) THEN
            WRITE(6,*) 'After ',SUBR,': I,J,SNOWI=',I,J,SNOWI(I,J)
            QCHECKI = .TRUE.
          END IF
          IF (MSI(I,J).gt.10000) THEN
            WRITE(6,*) 'After ',SUBR,': I,J,MSI=',I,J,MSI(I,J),RSI(I,J)
c            QCHECKI = .TRUE.
          END IF
          END IF
        END DO
      END DO



      IF (QCHECKI)
     &     call stop_model("CHECKI: Ice variables out of bounds",255)

      END SUBROUTINE CHECKI