      SUBROUTINE SEAIC6U(C,EVAPR,GC,GTP,GT,GTA,RES,RUNOFF,SIC,SICP,SNO, 
     1                   QFN,HSEA,BEGO,BWGO,IDAY,GMT,LDAY,DELT,I,J,
     2                   ICEMOD,INCD)
C---------------------------------------------------------------------- 
C 
C     * MAY 20/95 - M.LAZARE.ADD DIAGNOSTIC CALCULATION FOR SNOW EVAP.
C     * JUN 28/94 - G.FLATO. PREVIOUS VERSION SEAIC6O.
C 
C     * THERMODYNAMIC SEA ICE MODEL (AFTER SEMTNER,JPO MAY,1976)
C     * DETERMINES CHANGES TO ICE AND SNOW AMOUNTS OVER PACK ICE. 
C 
C---------------------------------------------------------------------- 
C     * DICTIONARY OF VARIABLES
C
C  ASIC  - initial sea ice amount in (kg m-2)
C  ASNO  - initial snow amount (water equivalent) in (kg m-2)
C  AGT   - initial ground temperature in (K)
C  BEGO  - net energy balance at ice/ocean interface (W m-2)
C  BGT   - snow (ground) temperature (K) after surface melt.
C  BND   - 1.5 * DELT
C  BSIC  - sea ice amount (kg m-2) after sublimation/deposition but
C          before surface melt
C  BSNO  - snow amount (kg m-2) after sublimation/deposition but
C          before surface melt
C  BWGO  - contribution to net fresh water flux into ocean due to ice 
C          growth or melt (kg m-2 s-1)
C  C     - SQRT(2) times the 'average' heat capacity of the upper snow/ice
C          surface layer (J m-2 K-1) - calculated in subroutine GRNBL6O
C  CICE  - SQRT(2) times CPACK (J m-2 K-1)
C  CPACK - heat capacity of 10cm 'upper layer' of pack ice (J m-2 K-1)
C  FSIC  - sea ice amount (kg m-2) before bottom melting or freezing
C  FSNO  - snow amount (kg m-2) before bottom melting or freezing
C  DAY   - current time in fractional Julian days
C  DELE  - EVAPR * DELT - total evaporation over time step (kg m-2)
C  DELM  - equivalent heat available for surface melt (kg m-2)
C  DESIC - change in sea ice amount over time step (kg m-2) from 
C          observations (used when ICEMOD=0)
C  DELT  - model time step (s)
C  DT    - time, in fractional Julian days, between current time
C          and next beginning-of-month day
C  DTS   - time between current time and next beginning-of-month day (s)
C  EVAPR - ground evaporation rate (kg m-2 s-1)
C  FFLXT - fresh water flux into ocean due to surface melt (kg m-2)
C  FFLXB - fresh water flux into ocean due to bottom growth or melt (kg m-2)
C  FO    - heat flux from ocean, taken as max(0.,RES) to eliminate negative 
C          values of RES
C  GC    - ground cover type (1.=pack ice, 0.=open water, -1.=land)
C  GMT   - Greenich meridian time of day (s)
C  GT    - ground temperature (K)
C  GTFSW - freezing point of sea water (K)
C  GTMSI - melting point of sea ice = TFREZ - 0.1 (K)
C  GTP   - "target" ground temperature (K) used in non-interactive ocean.
C  HF    - latent heat of fusion (J kg-1) of snow or sea ice = HS-HV
C  HS    - latent heat of sublimation (J kg-1)
C  HV    - latent heat of vaporization (J kg-1)
C  HSEA  - conductive heat flux up through sea ice (W m-2)
C  ICEMOD- switch =1 if interactive icemodel is being used, =0 if not
C  IDAY  - current Julian day
C  INCD  - date increment switch (0=perpetual, 1=increment)
C  LDAY  - Julian day of next beginning-of-month day
C  QFN   - snow evaporation rate (water equivilent) (kg m-2 s-1)
C  RES   - residual (ground energy balance climatology) (W m-2), used to 
C          mimic oceanic heat transport when slab ocean is used or for
C          'lakes' when 3-D ocean is used.
C  RUNOFF- over ocean, freshwater supplied from liquid precipitation or
C          surface melt (kg m-2)
C  SIC   - sea ice amount (kg m-2)
C  SICMIN- minimum sea ic amount, below which grid cell is considered to 
C          be open water (kg m-2)
C  SICN  - interpolated sea ice amount based on monthly observations used 
C          when ICEMOD=0 (kg m-2)
C  SICP  - "target" sea ice amount (kg m-2) defined at beginning of month 
C  SNO   - snow amount (kg m-2)
C  TFREZ - freezing temperature of fresh water or melting point of snow (K)
C---------------------------------------------------------------------- 
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)

      COMMON /PARAM1/ PI,RVORD,TFREZ,HS,HV,DAYLNT 
      COMMON /PARAM3/CSNO,CPACK,GTFSW,RKHI,SBC,SNOMAX 
      COMMON /PARAM5/ CONI,DENI,SICMIN,DFSIC,CONF 
C---------------------------------------------------------------------- 
      ASIC=SIC
      ASNO=SNO
      AGT=GT
      GTMSI=TFREZ-0.1 
      HF=HS-HV
      BEGO=-1.*HSEA 
C     RUNOFF=0.      (COMMENTED OUT SINCE INITIALIZED IN VGRFAC5/6).
      DELE=0. 
      DELM=0. 
      CICE=CPACK*SQRT(2.) 
C 
C     * SUBLIMATION/DEPOSITION. 
C     * NOTE: ALL OF EVAPORATIVE FLUX IS APPLIED TO SNOW WHEN SNOW IS
C     * PRESENT.
C 
      DELE=EVAPR*DELT 
      IF(SNO.GT.DELE) THEN
         SNO=SNO-DELE 
         QFN=EVAPR
      ELSE
         SIC=SIC-(DELE-SNO) 
         QFN=SNO/DELT
         SNO=0. 
         IF(ICEMOD.NE.0) THEN 
            SIC=MAX(SIC,0.) 
            IF(SIC.LT.SICMIN) GO TO 100 
         ENDIF
      ENDIF 
C
C     * SAVE ICE AND SNOW MASS BEFORE SURFACE MELT
C
      BSIC=SIC
      BSNO=SNO
C 
C     * CALCULATE MELTING AND RUNOFF IF GT IS GREATER THAN FREEZING 
C     * AND ADJUST GT FOR LATENT HEAT LOSS. 
C 
      IF(GT.GT.TFREZ .AND. SNO.GT.0.0) THEN 
C 
C        * WHEN THERE IS SNOW THE SURFACE TEMPERATURE MUST BE ABOVE 
C        * THE MELTING POINT OF SNOW (TFREZ) FOR MELTING TO OCCUR.
C        * IF ALL SNOW IS MELTED, THE RESULTING GT IS CALCULATED. 
C 
         DELM=(GT-TFREZ)*C/HF 
         IF(SNO.GE.DELM) THEN 
            SNO=SNO-DELM
            RUNOFF=RUNOFF+DELM
            GT=TFREZ
         ELSE 
            RUNOFF=RUNOFF+SNO
            GT=GT-HF*SNO/C
            SNO=0.0 
         ENDIF
      ENDIF 
      BGT=GT
C 
      IF(GT.GT.GTMSI .AND. SNO.EQ.0.0) THEN 
C 
C        * WHEN THERE ISN'T ANY SNOW THE SURFACE TEMPERATURE NEEDS ONLY 
C        * BE ABOVE THE MELTING POINT OF ICE (GTMSI) FOR MELTING
C        * TO OCCUR.
C 
         DELM=(GT-GTMSI)*CICE/HF
         IF(DELM.LT.SIC) THEN 
            SIC=SIC-DELM
            GT=GTMSI
         ELSE 
            GT=GT-HF*SIC/CICE 
            SIC=0.
         ENDIF
         IF(ICEMOD.NE.0) THEN 
            RUNOFF=RUNOFF+MIN(DELM,BSIC)
            IF(SIC.LT.SICMIN) GO TO 100 
         ENDIF
      ENDIF 
C
C     * CALCULATE FRESH WATER FLUX DUE TO SURFACE MELT - SAVE ICE
C     * AND SNOW AMOUNT BEFORE GROWTH/MELT
C
      FFLXT=((BSIC-SIC)+(BSNO-SNO))/DELT
      FSIC=SIC 
      FSNO=SNO 
C 
      IF(ICEMOD.EQ.0) THEN
C 
C        * "DMELT" IS THE ICE-THICKNESS CHANGE FROM TOP.
C 
         DMELT=SIC-ASIC 

C        * CALCULATE CHANGE IN SEA-ICE AMOUNT BASED ON OBSERVATIONS.
C 
         IF(INCD.NE.0) THEN 
            DAY=FLOAT(IDAY)+GMT/86400.
            DT=FLOAT(LDAY)-DAY
            IF(DT.LT.0.) DT=DT+365. 
            DTS=DT*86400. 
            SICN=ASIC+(SICP-ASIC)*DELT/DTS
            BND=1.5*DELT
            IF(SICP.LT.SICMIN.AND.DTS.LE.BND) GT=GTFSW
         ELSE 
            SICN=ASIC 
         ENDIF
         DELSIC=SICN-ASIC 
C 
C        * DIAGNOSTIC EQUATION FOR ENERGY BALANCE AT OCEAN INTERFACE, 
C        * I.E. "BEGO". 
C        * RESET SIC TO OBSERVED VALUE SICN AND GT TO VALUE PRIOR TO
C        * SEA-ICE MELTING. 
C 
         BEGO = BEGO + HF*DELSIC/DELT - HF*DMELT/DELT 
         SIC=MAX(SICN,0.) 
         IF(SNO.EQ.0.0.AND.DELSIC.LT.0.) THEN 
            IF(SIC.GT.0.) THEN
               GT=GTMSI 
            ELSE
               GT=BGT-HF*ABS(DELSIC)/CICE 
            ENDIF 
         ENDIF
      ELSE
C 
C        * CALCULATE MELTING/FREEZING AT BOTTOM. IF SEA ICE MELTS TO
C        * SICMIN AND THERE IS SOME SNOW ON TOP OF IT, THEN CONVERT 
C        * ENOUGH SNOW INTO SEA ICE TO RETURN SIC TO SICMIN.
C        * LIMIT RESIDUAL TO NON-NEGATIVE VALUES SINCE OTHERWISE WOULD
C        * BE INDICATIVE OF HEAT ADVECTION AWAY FROM "COLD" POINT.
C 
         FO=MAX(0.,RES) 
         SIC=SIC+DELT*(HSEA-FO)/HF
         BEGO = BEGO + HF*(SIC-FSIC)/DELT 
C 
         IF((SNO+SIC).LT.SICMIN) THEN 
            SIC=SNO+SIC 
            SNO=0.
         ENDIF
         IF(SIC.LT.SICMIN.AND.(SNO+SIC).GE.SICMIN) THEN 
            SNO=SNO-(SICMIN-SIC)
            SIC=SICMIN
         ENDIF
      ENDIF 
C
C     * CALCULATE FRESH WATER FLUX AT BOTTOM DUE TO GROWTH OR MELT
C
      FFLXB=((FSIC-SIC)+(FSNO-SNO))/DELT
C
C     * CALCULATE CONTRIBUTION TO NET FRESH WATER FLUX INTO OCEAN
C
      BWGO=(FFLXT+FFLXB)
C 
  100 IF(ICEMOD.NE.0.AND.SIC.LT.SICMIN) THEN
C 
C        * SEA ICE AT THIS GRID POINT HAS BEEN MELTED TO LESS THAN
C        * SICMIN AND THE GROUND COVER MUST BE ADJUSTED.
C 
         GC=0.0 
         SNO=0.0
         GT=GTFSW 
         GTP=GTFSW
         IF(ASIC.GE.SICMIN) THEN
           WRITE(6,6000) I,J
           WRITE(6,6010) DELE,DELM,AGT,ASIC,ASNO
           WRITE(6,6020) HSEA,RES,GT,SIC,SNO
         ENDIF
      ENDIF 
C---------------------------------------------------------------------
 6000 FORMAT('0SEA ICE DISAPPEARS AT I=',I5,', J=',I5)
 6010 FORMAT(' EVAP,MELT,AGT,ASIC,ASNO:',5F9.3) 
 6020 FORMAT(' HSEA, RES, GT, SIC, SNO:',5F9.3,/) 
      RETURN
      END 
