!@sum  SEAICE_DRV contains drivers for SEAICE related routines
!@auth Gavin Schmidt
!@ver  1.0
!@cont PRECIP_SI,GROUND_SI


      SUBROUTINE PRECIP_SI 1,11
!@sum  PRECIP_SI driver for applying precipitation to sea ice fraction
!@auth Original Development team
!@ver  1.0
!@calls seaice:prec_si
      USE CONSTANT, only : teeny,rhoi,grav
      USE MODEL_COM, only : im,jm,fland,kocean,itoice,itlkice,focean
     *     ,jday,p,ptop
      USE GEOM, only : imaxj,dxyp,bydxyp
      USE FLUXES, only : runpsi,prec,eprec,srunpsi,gtemp,apress,fwsim

      USE SEAICE, only : prec_si, ace1i, lmi,xsi,debug
      USE SEAICE_COM, only : rsi,msi,snowi,hsi,ssi,flag_dsws,pond_melt

      USE DAGCOM, only : aj,areg,aij,jreg,ij_f0oi,ij_erun2,j_imelt
     *     ,j_smelt
      USE DOMAIN_DECOMP, only : GRID
      USE DOMAIN_DECOMP, only : GET
      IMPLICIT NONE

      REAL*8, DIMENSION(LMI) :: HSIL,TSIL,SSIL
      REAL*8 SNOW,MSI2,PRCP,ENRGP,RUN0,POICE,SRUN0

      INTEGER I,J,JR,ITYPE,N
      LOGICAL WETSNOW
      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)

      DO J=J_0, J_1
      DO I=1,IMAXJ(J)
        JR=JREG(I,J)
      POICE=    RSI(I,J) *(1.-FLAND(I,J))
      RUNPSI(I,J)=0
      SRUNPSI(I,J)=0

      IF (POICE.gt.0) THEN

        IF (FOCEAN(I,J).gt.0) THEN
          ITYPE=ITOICE
        ELSE
          ITYPE=ITLKICE
        END IF
        PRCP=PREC(I,J)
        ENRGP=EPREC(I,J)      ! energy of precip
        SNOW=SNOWI(I,J)
        MSI2=MSI(I,J)
        HSIL(:) = HSI(:,I,J)      ! sea ice temperatures
        SSIL(:) = SSI(:,I,J)      ! sea ice salt


        AIJ(I,J,IJ_F0OI)=AIJ(I,J,IJ_F0OI)+ENRGP*POICE
C**** CALL SUBROUTINE FOR CALCULATION OF PRECIPITATION OVER SEA ICE

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

     *       WETSNOW)

        SNOWI(I,J)  =SNOW
        RUNPSI(I,J) =RUN0
        SRUNPSI(I,J)=SRUN0
        HSI(:,I,J)=HSIL(:)
        SSI(:,I,J)=SSIL(:)

        FLAG_DSWS(I,J)=FLAG_DSWS(I,J).or.WETSNOW
C**** reset flag if there was fresh snow (i.e. prcp but no rain!)
        IF (.not. WETSNOW .and. PRCP.gt.0.) FLAG_DSWS(I,J)=.FALSE.
C**** pond_melt accumulates in melt season only
        if ((J.gt.JM/2 .and. (jday.ge.152 .and. jday.lt.212)) .or.
     *       (J.lt.JM/2 .and. (jday.ge.334 .or. jday.lt.31))) then
          pond_melt(i,j)=pond_melt(i,j)+0.3d0*RUN0
        end if

C**** set gtemp array
        MSI(I,J)=MSI2
        GTEMP(1:2,2,I,J)=TSIL(1:2)

        FWSIM(I,J) = RSI(I,J)*(ACE1I+SNOW+MSI2-SUM(SSIL(1:LMI)))

C**** Accumulate diagnostics for ice fraction
        AJ(J,J_IMELT,ITYPE)=AJ(J,J_IMELT,ITYPE)+RUN0 *POICE
        AJ(J,J_SMELT,ITYPE)=AJ(J,J_SMELT,ITYPE)+SRUN0*POICE
c       AJ(J,J_HMELT,ITYPE)=AJ(J,J_HMELT,ITYPE)+ERUN *POICE  ! ==0
C**** Accumulate regional diagnostics
        AREG(JR,J_IMELT)=AREG(JR,J_IMELT)+RUN0 *POICE*DXYP(J)
        AREG(JR,J_SMELT)=AREG(JR,J_SMELT)+SRUN0*POICE*DXYP(J)
c       AREG(JR,J_HMELT)=AREG(JR,J_HMELT)+ERUN *POICE*DXYP(J) ! ==0

      END IF

C**** Calculate pressure anomaly at surface
      APRESS(I,J) = 100.*(     *     RSI(I,J)*(SNOWI(I,J)+ACE1I+MSI(I,J))*GRAV

      END DO
      END DO
      APRESS(2:IM,1)  = APRESS(1,1)
      APRESS(2:IM,JM) = APRESS(1,JM)
C****
      END SUBROUTINE PRECIP_SI


      SUBROUTINE UNDERICE 1,13
!@sum  underice calculates basal fluxes under sea and lake ice
!@+    saves the resulting fluxes
!@auth Gavin Schmidt
!@ver  1.0
!@calls iceocean,icelake
      USE CONSTANT, only : rhow,rhows,lhm,omega,rhoi,byshi,shw
      USE MODEL_COM, only : im,jm,focean,dtsrc,qcheck,kocean
      USE GEOM, only : sinp,imaxj,dxyp

      USE SEAICE, only : lmi,xsi,icelake,iceocean,ac2oim,alami,alpha
     *     ,tfrez,debug
      USE SEAICE_COM, only : msi,hsi,ssi,rsi

      USE FLUXES, only : fmsi_io,fhsi_io,fssi_io,ui2rho,gtemp,sss,mlhc

      USE LAKES_COM, only : tlake,mwl,flake,gml,mldlk

      USE DOMAIN_DECOMP, only : GRID
      USE DOMAIN_DECOMP, only : GET
      IMPLICIT NONE
      INTEGER I,J,N
      REAL*8 coriol,ustar,Tm,Sm,Si,Ti,dh,mflux,hflux,sflux,fluxlim
     *     ,mlsh   !,mfluxmax

      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)

      DO J=J_0, J_1
C**** Coriolis parameter (on tracer grid)
        coriol = ABS(2.*OMEGA*SINP(J))
        DO I=1,IMAXJ(J)
          IF ((FOCEAN(I,J)+FLAKE(I,J))*RSI(I,J).gt.0) THEN
C**** Set mixed layer conditions
            Tm = GTEMP(1,1,I,J)

            dh = 0.5*(XSI(LMI)*MSI(I,J))/RHOI
c            mfluxmax = (MSI(I,J)-AC2OIM)/dtsrc
            IF (FOCEAN(I,J).gt.0) THEN
C**** Ice lowest layer conditions
              Si = 1d3*SSI(LMI,I,J)/(XSI(LMI)*MSI(I,J))
              Ti = ((HSI(LMI,I,J)-SSI(LMI,I,J)*LHM)/(XSI(LMI)*MSI(I,J))
     *             +LHM)*BYSHI
C**** for the Salinity thermodynamics case (or something similar)
c             Ti = TICE(HSI(LMI,I,J),SSI(LMI,I,J),XSI(LMI)*MSI(I,J))
              IF (KOCEAN.eq.1) THEN
C**** should we calculate ocean rho(Tm,Sm) here?
                Ustar = MAX(5d-4,SQRT(UI2rho(I,J)/RHOWS))
                Sm = SSS(I,J)
                mlsh = MLHC(I,J)
                call iceocean(Ti,Si,Tm,Sm,dh,Ustar,Coriol,dtsrc,mlsh,

     *               mflux,sflux,hflux)
              ELSE ! for fixed SST assume freezing temp at base,implicit
                hflux=alami*(Ti-tfrez(sss(i,j)))/(dh+alpha*dtsrc*alami
     *               *byshi/(XSI(LMI)*MSI(I,J)))
                mflux=0.
                sflux=0.

              END IF
            ELSE   ! for lakes (no salinity so solve directly)
              Ti = (HSI(LMI,I,J)/(XSI(LMI)*MSI(I,J))+LHM)*BYSHI
              mlsh=SHW*MLDLK(I,J)*RHOW
              sflux = 0.
              call icelake(Ti,Tm,dh,dtsrc,mlsh,

     *             mflux,hflux)

C**** Limit lake-to-ice flux if lake is too shallow (< 40cm)
              IF (MWL(I,J).lt.0.4d0*RHOW*FLAKE(I,J)*DXYP(J)) THEN
                FLUXLIM=-GML(I,J)/(DTSRC*FLAKE(I,J)*DXYP(J))
                IF (hflux.lt.FLUXLIM) hflux = FLUXLIM
                if (mflux.lt.0) then
                  mflux = 0.

                end if
                if (qcheck) print*,"Flux limiting",I,J,MWL(I,J)/
     *               (RHOW*FLAKE(I,J)*DXYP(J)),FLUXLIM*DTSRC
              END IF
            END IF
            FMSI_IO(I,J) = mflux*dtsrc   ! positive down
            FHSI_IO(I,J) = hflux*dtsrc
            FSSI_IO(I,J) = sflux*dtsrc

          ELSE
            FMSI_IO(I,J) = 0.
            FHSI_IO(I,J) = 0.
            FSSI_IO(I,J) = 0.

          END IF
        END DO
      END DO

C****
      RETURN
      END SUBROUTINE UNDERICE


      SUBROUTINE MELT_SI 1,13
!@sum  MELT_SI driver for lateral melt of sea ice
!@auth Gary Russell/Gavin Schmidt
!@ver  1.0
!@calls SEAICE:SIMELT
      USE CONSTANT, only : sday
      USE MODEL_COM, only : im,jm,kocean,focean,itoice,itlkice ! ,itime
     *     ,itocean,itlake,dtsrc                               ! ,nday
      USE GEOM, only : dxyp,imaxj
      USE DAGCOM, only : aj,j_imelt,j_hmelt,j_smelt,areg,jreg
      USE SEAICE, only : simelt,tfrez
      USE SEAICE_COM, only : rsi,hsi,msi,lmi,snowi,ssi

      USE LAKES_COM, only : flake
      USE FLUXES, only : sss,melti,emelti,smelti,gtemp,mlhc,fwsim

      USE DOMAIN_DECOMP, only : GRID
      USE DOMAIN_DECOMP, only : GET
      IMPLICIT NONE
      REAL*8, DIMENSION(LMI) :: HSIL,TSIL,SSIL
      REAL*8 MSI2,ROICE,SNOW,ENRGUSED,RUN0,SALT,POCEAN,TFO
     *     ,PWATER,Tm,DT,ENRGMAX

      INTEGER I,J,ITYPE,JR,ITYPEO
      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**** CALCULATE LATERAL MELT ONCE A DAY (ALSO ELIMINATE SMALL AMOUNTS)
C**** We could put this in daily but it then we need an extra routine to
C**** add fluxes to oceans/lakes.
cc      IF (MOD(ITIME,NDAY).eq.0) THEN
cc        DT=SDAY    ! if called more frequently this should change
        DT=DTsrc    ! now do this every physics time step
        DO J=J_0, J_1
        DO I=1,IMAXJ(J)
          PWATER=FOCEAN(I,J)+FLAKE(I,J)
          POCEAN=FOCEAN(I,J)
          RUN0=0. ; ENRGUSED=0. ; SALT=0.

C**** Call simelt if (lake and v. small ice) or (q-flux ocean, some ice)
C**** now include lat melt for lakes and any RSI < 1
          IF (     *         .0)) .or. (KOCEAN.eq.1.and.POCEAN*RSI(I,J).gt.0) ) THEN
            JR=JREG(I,J)
            IF (POCEAN.gt.0) THEN
              ITYPE =ITOICE
              ITYPEO=ITOCEAN
              TFO = tfrez(sss(i,j))
            ELSE
              ITYPE =ITLKICE
              ITYPEO=ITLAKE
              TFO = 0.
            END IF
            Tm=GTEMP(1,1,I,J)
            ROICE=RSI(I,J)
            MSI2=MSI(I,J)
            SNOW=SNOWI(I,J)     ! snow mass
            HSIL(:)= HSI(:,I,J) ! sea ice enthalpy
            SSIL(:)= SSI(:,I,J) ! sea ice salt
            ENRGMAX= MAX(Tm-TFO,0d0)*MLHC(I,J) ! max energy for melt

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

     *           ,ENRGMAX,ENRGUSED,RUN0,SALT)

C**** accumulate diagnostics
           AJ(J,J_HMELT,ITYPE)=AJ(J,J_HMELT,ITYPE)-ENRGUSED*ROICE*PWATER
           AJ(J,J_SMELT,ITYPE)=AJ(J,J_SMELT,ITYPE)+    SALT*ROICE*PWATER
           AJ(J,J_IMELT,ITYPE)=AJ(J,J_IMELT,ITYPE)+    RUN0*ROICE*PWATER
           AJ(J,J_HMELT,ITYPEO)=AJ(J,J_HMELT,ITYPEO)-ENRGUSED*(1.-ROICE)
     *          *PWATER
           AJ(J,J_SMELT,ITYPEO)=AJ(J,J_SMELT,ITYPEO)+    SALT*(1.-ROICE)
     *          *PWATER
           AJ(J,J_IMELT,ITYPEO)=AJ(J,J_IMELT,ITYPEO)+    RUN0*(1.-ROICE)
     *          *PWATER
           AREG(JR,J_HMELT)=AREG(JR,J_HMELT)-ENRGUSED*PWATER*DXYP(J)
           AREG(JR,J_SMELT)=AREG(JR,J_SMELT)+    SALT*PWATER*DXYP(J)
           AREG(JR,J_IMELT)=AREG(JR,J_IMELT)+    RUN0*PWATER*DXYP(J)
C**** Update prognostic sea ice variables
            RSI(I,J)=ROICE
            MSI(I,J)=MSI2
            SNOWI(I,J)=SNOW
            HSI(:,I,J)=HSIL(:)
            SSI(:,I,J)=SSIL(:)

          END IF
C**** Save fluxes (in kg, J etc.), positive into ocean
          MELTI(I,J) = RUN0*PWATER*DXYP(J)
          EMELTI(I,J)=-ENRGUSED*PWATER*DXYP(J)
          SMELTI(I,J)= SALT*PWATER*DXYP(J)

        END DO
        END DO
cc      ELSE
cc        MELTI=0. ; EMELTI=0. ; SMELTI=0.
cc#ifdef TRACERS_WATER
cc        TRMELTI=0.
cc#endif
cc      END IF
C****
      RETURN
      END SUBROUTINE MELT_SI


      SUBROUTINE GROUND_SI 1,14
!@sum  GROUND_SI driver for applying surface + base fluxes to sea ice
!@auth Gary Russell/Gavin Schmidt
!@ver  1.0
!@calls SEAICE:SEA_ICE
      USE CONSTANT, only : grav,rhows,rhow
      USE MODEL_COM, only : im,jm,dtsrc,fland,kocean,focean
     *     ,itoice,itlkice,jday,p,ptop
      USE GEOM, only : imaxj,dxyp
      USE FLUXES, only : e0,e1,evapor,runosi,erunosi,srunosi,solar
     *     ,fmsi_io,fhsi_io,fssi_io,apress,gtemp,sss

      USE SEAICE, only : sea_ice,ssidec,lmi,xsi,ace1i,qsfix,debug
     *     ,snowice, snow_ice, rhos
      USE SEAICE_COM, only : rsi,msi,snowi,hsi,ssi,pond_melt,flag_dsws

      USE LAKES_COM, only : mwl,gml,flake
      USE DAGCOM, only : aj,areg,aij,jreg,ij_erun2,ij_rsoi,ij_msi2
     *     ,j_imelt,j_hmelt,j_smelt,j_rsnow,ij_rsit,ij_rsnw,ij_snow
     *     ,ij_mltp
      USE DOMAIN_DECOMP, only : GRID
      USE DOMAIN_DECOMP, only : GET
      IMPLICIT NONE

      REAL*8, DIMENSION(LMI) :: HSIL,SSIL
      REAL*8 SNOW,ROICE,MSI2,F0DT,F1DT,EVAP,SROX(2)
     *     ,FMOC,FHOC,FSOC,POICE,PWATER,SCOVI
      REAL*8 MFLUX,HFLUX,SFLUX,RUN,ERUN,SRUN,MELT12
      REAL*8 MSNWIC,HSNWIC,SSNWIC,SM,TM
      INTEGER I,J,JR,ITYPE
      LOGICAL WETSNOW

      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)

      debug=.false.
      DO J=J_0, J_1
      DO I=1,IMAXJ(J)
      PWATER=FOCEAN(I,J)+FLAKE(I,J)   ! 1.-FLAND(I,J)
      ROICE=RSI(I,J)
      POICE=ROICE*PWATER
      JR=JREG(I,J)
      SOLAR(3,I,J)=0
      RUNOSI(I,J)=0
      ERUNOSI(I,J)=0
      SRUNOSI(I,J)=0

      IF (POICE.gt.0) THEN

        F0DT=E0(I,J,2) ! heat flux to the top ice surface (J/m^2)
        F1DT=E1(I,J,2) ! heat flux between 1st and 2nd ice layer (J/m^2)
        EVAP=EVAPOR(I,J,2) ! evaporation/dew at the ice surface (kg/m^2)
        SROX(1)=SOLAR(2,I,J) ! solar radiation absrbd by sea ice (J/m^2)
        FMOC=fmsi_io(i,j)  ! mass flux at base (kg/m^2)
        FSOC=fssi_io(i,j)  ! salt flux at base (kg/m^2)
        FHOC=fhsi_io(i,j)  ! heat flux at base (J/M^2)
        SNOW= SNOWI(I,J)  ! snow mass (kg/m^2)
        MSI2= MSI(I,J)
        HSIL(:) = HSI(:,I,J)  ! sea ice enthalpy
        SSIL(:) = SSI(:,I,J)  ! sea ice salt
        WETSNOW=FLAG_DSWS(I,J)  ! wetness of snow
        Tm=GTEMP(1,1,I,J)    ! ocean mixed layer temperature (C)

        IF (FOCEAN(I,J).gt.0) THEN
          ITYPE=ITOICE
          Sm=SSS(I,J)           ! ocean mixed layer salinity (psu)
        ELSE
          ITYPE=ITLKICE
          Sm=0.                 ! lakes always fresh
        END IF

        AIJ(I,J,IJ_RSOI) =AIJ(I,J,IJ_RSOI) +POICE
        AIJ(I,J,IJ_MSI2) =AIJ(I,J,IJ_MSI2) +MSI2*POICE

        CALL SEA_ICE(DTSRC,SNOW,ROICE,HSIL,SSIL,MSI2,F0DT,F1DT,EVAP,SROX

     *       ,FMOC,FHOC,FSOC,RUN,ERUN,SRUN,WETSNOW,MELT12)

C**** Decay sea ice salinity
        if (.not. qsfix .and. FOCEAN(I,J).gt.0) then
          CALL SSIDEC(SNOW,MSI2,HSIL,SSIL,DTsrc,

     *         MFLUX,HFLUX,SFLUX)
        else
          MFLUX=0. ; SFLUX=0. ; HFLUX=0.

        end if

C**** Calculate snow-ice possibility
        if (snow_ice .eq. 1 .and. FOCEAN(I,J).gt.0) then
          call snowice(Tm,Sm,SNOW,MSI2,HSIL,SSIL,qsfix,

     *         MSNWIC,HSNWIC,SSNWIC)
        else
          MSNWIC=0. ; SSNWIC=0. ; HSNWIC=0.

        end if

C**** RESAVE PROGNOSTIC QUANTITIES
        SNOWI(I,J)=SNOW
        HSI(:,I,J)=HSIL(:)
        SSI(:,I,J)=SSIL(:)
        MSI(I,J) = MSI2

        FLAG_DSWS(I,J)=WETSNOW
C**** pond_melt accumulates in melt season only
        if ((J.gt.JM/2 .and. (jday.ge.152 .and. jday.lt.212)) .or.
     *       (J.lt.JM/2 .and. (jday.ge.334 .or. jday.lt.31))) then
          pond_melt(i,j)=pond_melt(i,j)+0.3d0*MELT12
          pond_melt(i,j)=MIN(pond_melt(i,j),0.5*(MSI2+SNOW+ACE1I))
        end if

C**** Net fluxes to ocean
        RUNOSI(I,J) = FMOC + RUN  + MFLUX + MSNWIC
        ERUNOSI(I,J)= FHOC + ERUN + HFLUX + HSNWIC
        SRUNOSI(I,J)= FSOC + SRUN + SFLUX + SSNWIC
        SOLAR(3,I,J)= SROX(2)


C**** ACCUMULATE DIAGNOSTICS
          SCOVI=0.
C**** snow cover diagnostic now matches that seen by the radiation
          IF (SNOW.GT.0) SCOVI=MIN(1d0,SNOW/(RHOS*0.1d0))*POICE

          AIJ(I,J,IJ_RSNW)=AIJ(I,J,IJ_RSNW)+SCOVI
          AIJ(I,J,IJ_SNOW)=AIJ(I,J,IJ_SNOW)+SNOW*POICE
          AIJ(I,J,IJ_RSIT)=AIJ(I,J,IJ_RSIT)+POICE
          AIJ(I,J,IJ_MLTP)=AIJ(I,J,IJ_MLTP)+pond_melt(i,j)*POICE

          AJ(J,J_RSNOW,ITYPE)=AJ(J,J_RSNOW,ITYPE)+SCOVI
          AJ(J,J_IMELT,ITYPE)=AJ(J,J_IMELT,ITYPE)+(FMOC+RUN+MFLUX
     *         +MSNWIC)*POICE
          AJ(J,J_HMELT,ITYPE)=AJ(J,J_HMELT,ITYPE)+(FHOC+ERUN+HFLUX
     *         +HSNWIC)*POICE
          AJ(J,J_SMELT,ITYPE)=AJ(J,J_SMELT,ITYPE)+(FSOC+SRUN+SFLUX
     *         +SSNWIC)*POICE

          AREG(JR,J_RSNOW)=AREG(JR,J_RSNOW)+SCOVI*DXYP(J)
          AREG(JR,J_IMELT)=AREG(JR,J_IMELT)+(FMOC+RUN+MFLUX+MSNWIC)
     *         *POICE*DXYP(J)
          AREG(JR,J_HMELT)=AREG(JR,J_HMELT)+(FHOC+ERUN+HFLUX+HSNWIC)
     *         *POICE*DXYP(J)
          AREG(JR,J_SMELT)=AREG(JR,J_SMELT)+(FSOC+SRUN+SFLUX+SSNWIC)
     *         *POICE*DXYP(J)

      END IF
C**** set total atmopsheric pressure anomaly in case needed by ocean
      APRESS(I,J) = 100.*(     *     *(SNOWI(I,J)+ACE1I+MSI(I,J))*GRAV

      END DO
      END DO
      APRESS(2:IM,1)  = APRESS(1,1)
      APRESS(2:IM,JM) = APRESS(1,JM)
C****
      END SUBROUTINE GROUND_SI


      SUBROUTINE FORM_SI 1,13
!@sum  FORM_SI driver for adding new sea ice
!@auth Original Development team
!@ver  1.0
!@calls seaice:addice
      USE CONSTANT, only : lhm,byshi
      USE MODEL_COM, only : im,jm,focean,kocean,fland
     *     ,itocean,itoice,itlake,itlkice,itime
      USE GEOM, only : imaxj,dxyp

      USE DAGCOM, only : aj,areg,aij,jreg,j_rsi,j_ace1,j_ace2,j_snow
     *     ,j_smelt,j_imelt,j_hmelt,ij_tsi,ij_ssi1,ij_ssi2,j_implh
     *     ,j_implm,ij_smfx
      USE SEAICE, only : ace1i,addice,lmi,fleadoc,fleadlk,xsi,debug
      USE SEAICE_COM, only : rsi,msi,snowi,hsi,ssi

      USE FLUXES, only : dmsi,dhsi,dssi,gtemp,fwsim

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

      REAL*8, DIMENSION(LMI) :: HSIL,TSIL,SSIL
      REAL*8 SNOW,ROICE,MSI2,ENRGFO,ACEFO,ACEFI,ENRGFI,SALTO,SALTI
     *     ,POICE,PWATER,FLEAD,POCEAN,DMIMP,DHIMP,DSIMP

      LOGICAL QFIXR
      INTEGER I,J,JR,ITYPE,ITYPEO,N
      integer :: J_0, J_1
C****
C**** Extract useful local domain parmeters from "grid"
C****
      CALL GET(grid, J_STRT = J_0, J_STOP = J_1)

      debug=.false.
      DO J=J_0, J_1
      DO I=1,IMAXJ(J)
      PWATER=FOCEAN(I,J)+FLAKE(I,J)
      ROICE=RSI(I,J)
      POICE=ROICE*PWATER
      POCEAN=(1.-ROICE)*PWATER
      JR=JREG(I,J)
      IF (PWATER.gt.0) THEN

        SNOW= SNOWI(I,J)      ! snow mass (kg/m^2)
        MSI2= MSI(I,J)
        HSIL(:) = HSI(:,I,J)      ! sea ice enthalpy
        SSIL(:) = SSI(:,I,J)      ! sea ice salt


        IF (FOCEAN(I,J).gt.0) THEN
          FLEAD=FLEADOC
          ITYPE=ITOICE
          ITYPEO=ITOCEAN
          IF (KOCEAN.eq.1) THEN
            QFIXR=.FALSE.
          ELSE
            QFIXR=.TRUE.
          END IF
        ELSE
          FLEAD=FLEADLK
          ITYPE=ITLKICE
          ITYPEO=ITLAKE
          QFIXR=.FALSE.
        END IF

        ACEFO=DMSI(1,I,J)
        ACEFI=DMSI(2,I,J)
        ENRGFO=DHSI(1,I,J)
        ENRGFI=DHSI(2,I,J)
        SALTO=DSSI(1,I,J)
        SALTI=DSSI(2,I,J)

C**** ice formation diagnostics on the atmospheric grid
C**** open ocean diagnostics
        AJ(J,J_SMELT,ITYPEO)=AJ(J,J_SMELT,ITYPEO)-SALTO *POCEAN
        AJ(J,J_HMELT,ITYPEO)=AJ(J,J_HMELT,ITYPEO)-ENRGFO*POCEAN
        AJ(J,J_IMELT,ITYPEO)=AJ(J,J_IMELT,ITYPEO)-ACEFO *POCEAN
C**** Ice-covered ocean diagnostics
        AJ(J,J_SMELT,ITYPE)=AJ(J,J_SMELT,ITYPE)-SALTI *POICE
        AJ(J,J_HMELT,ITYPE)=AJ(J,J_HMELT,ITYPE)-ENRGFI*POICE
        AJ(J,J_IMELT,ITYPE)=AJ(J,J_IMELT,ITYPE)-ACEFI *POICE
C**** regional diagnostics
        AREG(JR,J_IMELT)=AREG(JR,J_IMELT)-
     *       (ACEFO *POCEAN+ACEFI *POICE)*DXYP(J)
        AREG(JR,J_HMELT)=AREG(JR,J_HMELT)-
     *       (ENRGFO*POCEAN+ENRGFI*POICE)*DXYP(J)
        AREG(JR,J_SMELT)=AREG(JR,J_SMELT)-
     *       (SALTO *POCEAN+SALTI *POICE)*DXYP(J)

c        debug=i.eq.20.and.j.eq.37
        CALL ADDICE (SNOW,ROICE,HSIL,SSIL,MSI2,TSIL,ENRGFO,ACEFO,ACEFI,
     *       ENRGFI,SALTO,SALTI,

     *       DMIMP,DHIMP,DSIMP,FLEAD,QFIXR)

C**** RESAVE PROGNOSTIC QUANTITIES
        SNOWI(I,J) = SNOW
        MSI(I,J)=MSI2
        HSI(:,I,J) = HSIL(:)
        SSI(:,I,J) = SSIL(:)

        IF (.not. QFIXR) THEN
          RSI(I,J)=ROICE
        ELSE
C**** save implicit mass-flux diagnostics
          AIJ(I,J,IJ_SMFX)=AIJ(I,J,IJ_SMFX)+ROICE*DMIMP
          AJ(J,J_IMPLM,ITOICE)=AJ(J,J_IMPLM,ITOICE)-DMIMP*POICE
          AJ(J,J_IMPLH,ITOICE)=AJ(J,J_IMPLH,ITOICE)-DHIMP*POICE
        END IF
C**** set gtemp array
        GTEMP(1:2,2,I,J)=TSIL(1:2)

        FWSIM(I,J) = RSI(I,J)*(ACE1I+SNOW+MSI2-SUM(SSIL(1:LMI)))

C**** ACCUMULATE DIAGNOSTICS
        AJ(J,J_RSI, ITYPE)=AJ(J,J_RSI, ITYPE)+        POICE
        AJ(J,J_ACE1,ITYPE)=AJ(J,J_ACE1,ITYPE)+ACE1I  *POICE
        AJ(J,J_ACE2,ITYPE)=AJ(J,J_ACE2,ITYPE)+MSI2   *POICE
        AJ(J,J_SNOW,ITYPE)=AJ(J,J_SNOW,ITYPE)+SNOW   *POICE
        IF (JR.ne.24) THEN
          AREG(JR,J_RSI) =AREG(JR,J_RSI) +        POICE*DXYP(J)
          AREG(JR,J_SNOW)=AREG(JR,J_SNOW)+SNOW   *POICE*DXYP(J)
          AREG(JR,J_ACE1)=AREG(JR,J_ACE1)+ACE1I  *POICE*DXYP(J)
          AREG(JR,J_ACE2)=AREG(JR,J_ACE2)+MSI2   *POICE*DXYP(J)
        END IF
        AIJ(I,J,IJ_TSI)=AIJ(I,J,IJ_TSI)+
     *       POICE*(XSI(3)*TSIL(3)+XSI(4)*TSIL(4))
        AIJ(I,J,IJ_SSI1)=AIJ(I,J,IJ_SSI1)+POICE*(SSIL(1)+SSIL(2))/ACE1I
        AIJ(I,J,IJ_SSI2)=AIJ(I,J,IJ_SSI2)+POICE*(SSIL(3)+SSIL(4))
     *       /MSI(I,J)


      if (TSIL(1).lt.-100.) then
         write(6,*) "Seaice: T < -100. i,j,TSI = ",i,j,TSIL(1:LMI)
         call stop_model("Seaice too cold after ADDICE",255)
      end if

      END IF
      END DO
      END DO
C**** replicate ice values at the north pole
      DO I=2,IM
        RSI(I,JM)=RSI(1,JM)
        MSI(I,JM)=MSI(1,JM)
        HSI(:,I,JM)=HSI(:,1,JM)
        SSI(:,I,JM)=SSI(:,1,JM)
        SNOWI(I,JM)=SNOWI(1,JM)
        GTEMP(1:2,2,I,JM)=GTEMP(1:2,2,1,JM)

        FWSIM(I,JM) = FWSIM(1,JM)
      END DO
C****
      END SUBROUTINE FORM_SI


      SUBROUTINE vflx_OCEAN 1,7
!@sum  vflx_OCEAN saves quantities for OHT calculations
!@auth Original Development Team
!@ver  1.0
      USE MODEL_COM, only : im,jm,focean
      USE DAGCOM, only : oa
      USE SEAICE_COM, only : hsi,snowi
      USE FLUXES, only : fwsim
      USE DOMAIN_DECOMP, only : GRID
      USE DOMAIN_DECOMP, only : GET
      IMPLICIT NONE
      INTEGER I,J
      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****
C****       DATA SAVED IN ORDER TO CALCULATE OCEAN TRANSPORTS
C****
C****       1  SNOWOI (INSTANTANEOUS AT NOON GMT)
C****       2  FWSIM  (INSTANTANEOUS AT NOON GMT)
C****       3  HSIT   (INSTANTANEOUS AT NOON GMT)
C****
      DO J=J_0, J_1
        DO I=1,IM
          IF (FOCEAN(I,J).gt.0) THEN
            OA(I,J,1)=SNOWI(I,J)
            OA(I,J,2)=FWSIM(I,J)
            OA(I,J,3)=SUM(HSI(:,I,J))
          END IF
        END DO
      END DO

      RETURN
C****
      END SUBROUTINE vflx_OCEAN


      SUBROUTINE init_ice(iniOCEAN) 1,15
!@sum  init_ice initialises ice arrays
!@auth Original Development Team
!@ver  1.0
      USE CONSTANT, only : byshi,lhm,shi,rhows
      USE MODEL_COM, only : im,jm,kocean,focean,flake0
      USE SEAICE, only : xsi,ace1i,ac2oim,ssi0,tfrez,oi_ustar0,silmfac
     *     ,lmi,snow_ice
      USE SEAICE_COM, only : rsi,msi,hsi,snowi,ssi,pond_melt,flag_dsws

      USE FLUXES, only : gtemp,ui2rho,fwsim,msicnv

      USE DAGCOM, only : npts,icon_OMSI,icon_OHSI,icon_OSSI,icon_LMSI
     *     ,icon_LHSI,conpt0
      USE PARAM
      USE DOMAIN_DECOMP, only : GRID
      USE DOMAIN_DECOMP, only : GET
      IMPLICIT NONE
      LOGICAL :: QCON(NPTS), T=.TRUE. , F=.FALSE. , iniOCEAN
      CHARACTER CONPT(NPTS)*10
      INTEGER I,J
      REAL*8 MSI1,TFO
      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**** set up a default ice-ocean stress field. This can be changed by
C**** adjusting oi_ustar0 in the parameter list. If ice dynamics
C**** is used, this is overwritten.
      call sync_param("oi_ustar0",oi_ustar0)
      UI2rho = rhows*(oi_ustar0)**2

C**** Adjust degree of lateral melt by changing silmfac
C**** Default is 2.5d-8, but could be changed by a factor of 2.
      call sync_param("silmfac",silmfac)

C**** Decide whether snow_ice formation is allowed
      call sync_param("snow_ice",snow_ice)

C**** clean up ice fraction/sea ice salinity possibly incorrect in I.C.
      DO J=J_0, J_1
      DO I=1,IM
        IF (FOCEAN(I,J)+FLAKE0(I,J).eq.0 .and. RSI(i,j).gt.0) RSI(I,J)=0
        IF (RSI(I,J).gt.0 .and. FLAKE0(I,J).gt.0) SSI(:,I,J)=0.
      END DO
      END DO
      IF (KOCEAN.EQ.0.and.iniOCEAN) THEN
C****   set defaults for no ice case
        DO J=J_0, J_1
        DO I=1,IM
          IF (RSI(I,J).le.0) THEN
            MSI1        =ACE1I
            MSI(I,J)    =AC2OIM
            SNOWI(I,J)  =0.
            IF (FOCEAN(I,J).gt.0) THEN
              SSI(1:2,I,J)=SSI0*XSI(1:2)*ACE1I
              SSI(3:4,I,J)=SSI0*XSI(3:4)*AC2OIM
              TFO = -1.87d0  ! reasonable value, doesn't really matter
            ELSE
              SSI(:,I,J)  = 0.
              TFO = 0.
            END IF
            HSI(1:2,I,J)=(SHI*TFO-LHM)*XSI(1:2)*ACE1I+LHM*SSI(1:2,I,J)
            HSI(3:4,I,J)=(SHI*TFO-LHM)*XSI(3:4)*AC2OIM+LHM*SSI(3:4,I,J)

            pond_melt(i,j) = 0.
            flag_dsws(i,j) = .FALSE.
          END IF
        END DO
        END DO
      END IF
C**** set GTEMP etc. array for ice
      DO J=J_0, J_1
      DO I=1,IM
        MSI1=SNOWI(I,J)+ACE1I
        GTEMP(1:2,2,I,J)=((HSI(1:2,I,J)-SSI(1:2,I,J)*LHM)/
     *       (XSI(1:2)*MSI1)+LHM)*BYSHI

        FWSIM(I,J) = RSI(I,J)*(MSI1+MSI(I,J)-SUM(SSI(1:LMI,I,J)))
        MSICNV(I,J)=0.   ! always initialise to zero
      END DO
      END DO

C**** Set conservation diagnostics for ice mass, energy, salt
      CONPT=CONPT0
      CONPT(3)="LAT. MELT" ; CONPT(4)="PRECIP"
      CONPT(5)="THERMO"    ; CONPT(6)="ADVECT"
      CONPT(8)="OCN FORM"
      QCON=(/ F, F, T, T, T, T, F, T, T, F, F/)
      CALL SET_CON(QCON,CONPT,"OICE MAS","(KG/M^2)        ",
     *     "(10**-9 KG/SM^2)",1d0,1d9,icon_OMSI)
      QCON=(/ F, F, T, T, T, T, F, T, T, F, F/)
      CALL SET_CON(QCON,CONPT,"OICE ENR","(10**6 J/M^2)   ",
     *     "(10**-3 W/M^2)  ",1d-6,1d3,icon_OHSI)
      QCON=(/ F, F, T, T, T, T, F, T, T, F, F/)
      CALL SET_CON(QCON,CONPT,"OICE SLT","(10**-3 KG/M^2) ",
     *     "(10**-12KG/SM^2)",1d3,1d12,icon_OSSI)
      CONPT(8)="LK FORM"
      QCON=(/ F, F, T, T, T, F, F, T, T, F, F/)
      CALL SET_CON(QCON,CONPT,"LKICE MS","(KG/M^2)        ",
     *     "(10**-9 KG/SM^2)",1d0,1d9,icon_LMSI)
      QCON=(/ F, F, T, T, T, F, F, T, T, F, F/)
      CALL SET_CON(QCON,CONPT,"LKICE EN","(10**6 J/M^2)   ",
     *     "(10**-3 W/M^2)  ",1d-6,1d3,icon_LHSI)
C****
      END SUBROUTINE init_ice


      SUBROUTINE conserv_OMSI(ICE),4
!@sum  conserv_MSI calculates total amount of snow and ice over ocean
!@auth Gavin Schmidt
!@ver  1.0
      USE MODEL_COM, only : im,jm,fim,focean
      USE GEOM, only : imaxj
      USE SEAICE, only : ace1i
      USE SEAICE_COM, only : rsi,msi,snowi
      IMPLICIT NONE
!@var ICE total ocean snow and ice mass (kg/m^2)
      REAL*8, DIMENSION(JM) :: ICE
      INTEGER I,J

      DO J=1,JM
        ICE(J)=0
        DO I=1,IMAXJ(J)
          ICE(J)=ICE(J)+RSI(I,J)*(MSI(I,J)+ACE1I+SNOWI(I,J))
     *         *FOCEAN(I,J)
        END DO
      END DO
      ICE(1) =FIM*ICE(1)
      ICE(JM)=FIM*ICE(JM)
      RETURN
C****
      END SUBROUTINE conserv_OMSI


      SUBROUTINE conserv_OHSI(EICE),3
!@sum  conserv_HSI calculates total ice energy over ocean
!@auth Gavin Schmidt
!@ver  1.0
      USE MODEL_COM, only : im,jm,fim,focean
      USE GEOM, only : imaxj
      USE SEAICE_COM, only : rsi,hsi
      IMPLICIT NONE
!@var EICE total ocean snow and ice energy (J/m^2)
      REAL*8, DIMENSION(JM) :: EICE
      INTEGER I,J

      DO J=1,JM
        EICE(J)=0
        DO I=1,IMAXJ(J)
          EICE(J)=EICE(J)+RSI(I,J)*FOCEAN(I,J)*SUM(HSI(:,I,J))
        END DO
      END DO
      EICE(1) =FIM*EICE(1)
      EICE(JM)=FIM*EICE(JM)
      RETURN
C****
      END SUBROUTINE conserv_OHSI


      SUBROUTINE conserv_OSSI(SALT),3
!@sum  conserv_SSI calculates total amount of salt in ocean ice
!@auth Gavin Schmidt
!@ver  1.0
      USE MODEL_COM, only : im,jm,fim,focean
      USE GEOM, only : imaxj
      USE SEAICE_COM, only : rsi,ssi,lmi
      IMPLICIT NONE
!@var SALT total salt in ocean ice (kg/m^2)
      REAL*8, DIMENSION(JM) :: SALT
      INTEGER I,J,L

      DO J=1,JM
        SALT(J)=0
        DO I=1,IMAXJ(J)
          IF (FOCEAN(I,J).gt.0) THEN
            SALT(J)=SALT(J)+FOCEAN(I,J)*RSI(I,J)*SUM(SSI(:,I,J))
          END IF
        END DO
      END DO
      SALT(1) =FIM*SALT(1)
      SALT(JM)=FIM*SALT(JM)
      RETURN
C****
      END SUBROUTINE conserv_OSSI


      SUBROUTINE conserv_LMSI(ICE),5
!@sum  conserv_LMSI calculates total amount of snow and ice over lakes
!@auth Gavin Schmidt
!@ver  1.0
      USE MODEL_COM, only : im,jm,fim
      USE GEOM, only : imaxj
      USE SEAICE, only : ace1i
      USE SEAICE_COM, only : rsi,msi,snowi
      USE LAKES_COM, only : flake
      IMPLICIT NONE
!@var ICE total lake snow and ice mass (kg/m^2)
      REAL*8, DIMENSION(JM) :: ICE
      INTEGER I,J

      DO J=1,JM
        ICE(J)=0
        DO I=1,IMAXJ(J)
          ICE(J)=ICE(J)+RSI(I,J)*(MSI(I,J)+ACE1I+SNOWI(I,J))
     *         *FLAKE(I,J)
        END DO
      END DO
      ICE(1) =FIM*ICE(1)
      ICE(JM)=FIM*ICE(JM)
      RETURN
C****
      END SUBROUTINE conserv_LMSI


      SUBROUTINE conserv_LHSI(EICE),4
!@sum  conserv_LHSI calculates total ice energy over lakes
!@auth Gavin Schmidt
!@ver  1.0
      USE MODEL_COM, only : im,jm,fim
      USE GEOM, only : imaxj
      USE SEAICE_COM, only : rsi,hsi
      USE LAKES_COM, only : flake
      IMPLICIT NONE
!@var EICE total lake snow and ice energy (J/m^2)
      REAL*8, DIMENSION(JM) :: EICE
      INTEGER I,J

      DO J=1,JM
        EICE(J)=0
        DO I=1,IMAXJ(J)
          EICE(J)=EICE(J)+RSI(I,J)*FLAKE(I,J)*SUM(HSI(:,I,J))
        END DO
      END DO
      EICE(1) =FIM*EICE(1)
      EICE(JM)=FIM*EICE(JM)
      RETURN
C****
      END SUBROUTINE conserv_LHSI


      SUBROUTINE daily_ice 1,6
!@sum daily_ice performs ice processes that are needed everyday
!@auth Gavin Schmidt
      USE MODEL_COM, only : jm,jday
      USE GEOM, only : imaxj
      USE SEAICE_COM, only : flag_dsws,pond_melt
      USE DOMAIN_DECOMP, only : GRID
      USE DOMAIN_DECOMP, only : GET
      IMPLICIT NONE
      INTEGER I,J
      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**** Every day adjust pond-melt
      DO J=J_0, J_1
      DO I=1,IMAXJ(J)

C**** pond_melt decreases linearly in shoulder season
        if (J.gt.JM/2 .and. (jday.ge.212 .and. jday.lt.244)) then
          pond_melt(i,j)=pond_melt(i,j)*(1.-1./(244.-jday))
        elseif (J.lt.JM/2 .and. (jday.ge.31 .and. jday.lt.60)) then
          pond_melt(i,j)=pond_melt(i,j)*(1.-1./(60.-jday))
C**** allow pond_melt to accumulate in melt season only
C**** otherwise pond_melt is zero
        elseif (.not.((J.gt.JM/2 .and. (jday.ge.152 .and. jday.lt.212))
     *         .or. (J.lt.JM/2 .and. (jday.ge.334 .or. jday.lt.31))) )
     *         then
          pond_melt(i,j)=0.
        end if
      END DO
      END DO

      RETURN
      END SUBROUTINE daily_ice