      SUBROUTINE GRNBL6O(GT,SNO,SIC,C,HSEA,EVAPR,
     1                   GC,PRESS,GTA,GTLAST,TS,TFLXS,QFLXS,FSC,
     2                   CSZ,FDL,SNOFAC,SOIL,WL,WF,PVEG,
     3                   DELT,ICLW,ILG,IS,IF)
C.......................................................................
C     * JAN 18/95 - M.LAZARE.
C     * JAN 05/88 - M.LAZARE. PREVIOUS VERSION GRENBL5.
C
C     * PURPOSE:
C     * OVER LAND/SEA-ICE, EVALUATE THE GROUND ENERGY BALANCE, TAKING INTO 
C     * ACCOUNT SOLAR FLUX, ATM. AND TERRESTRIAL FLUXES, SENSIBLE AND 
C     * EVAPORATION HEAT LOSS, HEAT LOST TO UNDERLYING SEA WATER.
C     * METHOD: THERMAL-INERTIA.
C     * SCHEME: FORWARD 1-DELT STEP.
C
C     * NOTE THAT THIS ROUTINE IS A VECTORIZATION OF GRENBL5, BASED ON
C     * CODE ALREADY DEVELOPED FOR OIFPOST (FOR CLASS).
C
C     * DICTIONARY OF VARIABLES:  
C IC  ALFA     AREA RATIO OF LEEDS OR FIRST-YEAR ICE TO TOTAL ICE. 
C O   C       HEAT CAPACITY OF GROUND(J M-2 K-1)
C C   CONF    THERMAL CONDUCTIVITY OF FIRST-YEAR ICE(W M-1 K-1) 
C C   CONI    THERMAL CONDUCTIVITY OF MULTI-YEAR ICE(W M-1 K-1) 
C     CONS    THERMAL CONDUCTIVITY OF SNOW          (W M-1 K-1) 
C C   CPACK   HEAT CAPACITY OF PACK ICE(J M-2 K-1)
C C   CSNO    HEAT CAPACITY OF SNOW(J M-2 K-1)
C C   CVICE   SPECIFIC HEAT OF SEA ICE (J M-3 K-1)
C C   CVMNRL  SPECIFIC HEAT OF MINERAL CONTENT (J M-3 K-1)
C C   CVWTR   SPECIFIC HEAT OF WATER (J M-3 K-1)
C C   DAYLNT  LENGTH OF DAY (S) 
C IC  DELT    TIME STEP OF THE MODEL
C IC  DENI    DENSITY OF ICE (KG M-3) 
C     DENS    DENSITY OF SNOW (KG M-3)
C     DEPTH   EFFECTIVE DEPTH FOR ICE AND SNOW (M)
C IC  DFSIC   DEPTH OF FIRST-YEAR ICE (M) 
C O   EVAPR   EVAPORATION RATE (KG M-2 S-1  I.E.  MM S-1) 
C I   FDL     ATMOSPHERIC NET DOWNWARD LONGWAVE FLUX(W M-2) 
C I   FSG     SOLAR FLUX ABSORBED AT THE GROUND (W/M2)
C I   GC      GROUND COVER (TYPE) 1.=PACK ICE 
C                                 0.=OPEN WATER 
C                                -1.=LAND 
C I/O GT      I, GROUND TEMP. AT TIME T,
C             O, GT AT TIME T+DELT, BUT BEFORE FREEZING/MELTING 
C I   GTA     GROUND TEMPERATURE DIURNALLY AVERAGED TO TIME T.
C O   GTLAST  GT AT TIME T
C     HCAP    VOLUMETRIC HEAT CAPACITY OF SOIL (J M-3 K-1)
C     HCOND   THERMAL CONDUCTIVITY OF SOIL (W M-1 K-1)
C     HLAT    LATENT HEAT LOST TO ATMOSPHERE (W/M2) 
C C   HS      LATENT HEAT OF SUBLIMATION AT 0C (J/KG) 
C O   HSEA    HEAT LOST THROUGH PACK ICE TO UNDERLYING SEA(W/M2)
C     HSENS   SENSIBLE HEAT LOST TO ATM. (W/M2) 
C C   HV      LATENT HEAT OF VAPORIZATION AT 0C(J/KG) 
C     ICLR    COLOUR CODE FOR SOIL
C IC  ICLW    INTERVAL TO CALCULATE LONGWAVE RADIATION
C     ITEX    TEXTURE CODE FOR SOIL 
C C   PI      3.1415926535898 
C I   PRESS   SURFACE PRESSURE (PA) 
C I   PVEG    PRIMARY VEGETATION CLASS
C I   QFLXS   =CD (V*) (QS-QG) EVAPF
C     RCON    RATIO OF HEAT CONDUCTIVITIES (ICE/SNOW) 
C C   RGAS    GAS CONSTANT
C C   RKHI    HEAT CONDUCTIVITY OF PACK ICE(W M-2 K-1)
C C   SBC     STEFAN BOLTZMANN CONSTANT (W M-2 K-4) 
C I   SIC     SEA ICE AMOUNT (KG M-2) 
C C   SICMIN  MINIMUM SEA ICE AMOUNT (KG M-2) 
C     SKRAT   RATIO OF SOIL THERMAL CONDUCTIVITY TO REFERENCE SOIL
C I   SNO     SNOW AMOUNT (KG M-2)
C I   SNOFAC  FRACTIONAL SNOW COVERAGE
C I   SOIL    SOIL CLASS
C I   TFLXS   =CD (V*) (THETAS-THETAG)
C C   TFREZ   FREEZING TEMPERATURE OF WATER 
C I   TS      SCREEN LEVEL TEMPERATURE
C C   GTFSW   FREEZING TEMPERATURE OF SEA ICE 
C I   WF      FRACTION OF THE GROUND WATER HOLDING
C             CAPACITY FILLED BY WATER IN FROZEN FORM 
C I   WL      FRACTION OF THE GROUND WATER HOLDING
C             CAPACITY FILLED BY WATER IN LIQUID FORM 
C     XVOID   VOLUME POROSITY OF SOIL 
C.......................................................................

      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)

C     * OUTPUT FIELDS:
C 
      REAL GT     (ILG), SNO    (ILG), SIC    (ILG), C      (ILG),
     1     HSEA   (ILG), EVAPR  (ILG)
C
C     * INPUT FIELDS:
C
      REAL GC     (ILG), PRESS  (ILG), GTA    (ILG), GTLAST (ILG),
     1     TS     (ILG), TFLXS  (ILG), QFLXS  (ILG), FSC    (ILG),
     2     CSZ    (ILG), FDL    (ILG), SNOFAC (ILG), SOIL   (ILG),
     3     WL     (ILG), WF     (ILG), PVEG   (ILG)
       
      COMMON /PARAMS/WW,TW,A,ASQ,GRAV,RGAS,RGOCP,RGOASQ,CPRES 
      COMMON /PARAM1/PI,RVORD,TFREZ,HS,HV,DAYLNT
      COMMON /PARAM3/CSNO,CPACK,GTFSW,RKHI,SBC,SNOMAX 
      COMMON /PARAM5/CONI,DENI,SICMIN,DFSIC,CONF 
      DATA THIRD/0.33333333333333/
      DATA CVWTR, CVICE, CVMNRL /4.186E+6, 2.05E+6, 2.38E+6/
C
C     * DEFINITION OF THERMAL DIFFUSIVITY FUNCTION (THERMAL CONDUCTIVITY
C     * DIVIDED BY VOLUMETRIC HEAT CAPACITY) FOR A REFERENCE SANDY LOAM 
C     * SOIL. 
C 
      FSK(X) = (2.9E-7*X + 4.0E-9)/(((1.-0.6*X)*X + 0.09)*(0.23+X)) 
C-------------------------------------------------------------------
      SNOLIM=(DENI-275.)/0.54
      DO 100 I=IS,IF
C
C       * INITIALIZE OUTPUT FIELDS.
C
        HSEA(I) =0.
        C(I)    =0.
        EVAPR(I)=0.
        
C
C       * CALCULATIONS EXCLUSIVELY FOR SEA-ICE.
C 
        IF(GC(I).GE.0.5)                                       THEN
C 
C         * HEAT CAPACITY OF PACK ICE (NORMALIZED BY SQRT(2)).
C 
          C(I)=(SNOFAC(I)*CSNO+(1.-SNOFAC(I))*CPACK)*SQRT(2.) 
C 
C         * HEAT FLUX THROUGH PACK ICE FROM UNDERLYING SEA WATER (HSEA). 
C         * SEE SUBROUTINE OIFPREP FOR SNOW DENSITY JUSTIFICATION COMMENTS. 
C         * CONVERT ANY SNOW MASS HAVING DENSITY GREATER THAN DENI TO 
C         * SEA-ICE (BASED ON USUAL EXPONENTIAL FORMULA WITH A=0.54 AND 
C         * SURFACE DENSITY OF 275 KG/M3).
C 
          IF(SNO(I).GT.SNOLIM) THEN 
            SIC(I)=SIC(I)+(SNO(I)-SNOLIM)
            SNO(I)=SNOLIM
          ENDIF
          IF(SNO(I).GT.SNOMAX) THEN 
            DENS=0.54*SNO(I)/LOG(1.+0.54*SNO(I)/275.)
          ELSE 
            DENS=100. 
          ENDIF
          CONS=2.805E-6*DENS**2
          RCONMF=CONI/CONF 
          RCON=CONI/CONS 
          DSNO=SNO(I)/DENS
          DSIC=SIC(I)/DENI
          FACT=DSIC+RCON*DSNO
          ALFA=MIN((45./SIC(I))**1.25,1.) 
          DEPTH=RCONMF*DFSIC*FACT/(ALFA*FACT+(1.-ALFA)*RCONMF*DFSIC) 
          HSEA(I)=CONI*(GTFSW-GTA(I))/DEPTH
      ENDIF
C
C     * CALCULATIONS EXCLUSIVELY FOR GLACIER LAND.
C 
      IF(GC(I).LE.-0.5 .AND. PVEG(I).EQ.1.)                      THEN
        C(I)=(SNOFAC(I)*CSNO+(1.-SNOFAC(I))*CPACK)/SQRT(2.)
      ENDIF
C
C     * CALCULATIONS EXCLUSIVELY FOR NON-GLACIER LAND.
C 
      IF(GC(I).LE.-0.5 .AND. PVEG(I).NE.1.)                      THEN
C 
C 
C       * OBTAIN TEXTURE CODE OF SOIL, RANGING FROM 0 (FINE) TO 9 (COARSE). 
C 
        IC   = INT(SOIL(I))/100
        ICLR = IC/100 
        ITEX = (IC-100*ICLR)/10 
C 
C       * SKRAT IS THE RATIO OF SOIL(I) THERMAL CONDUCTIVITY TO THAT OF A
C       * REFERENCE SANDY LOAM SOIL(I), AS A FUCTION OF SOIL(I) TYPE. 
C 
        SKRAT = 0.8+(REAL(ITEX))*(1.5-0.8)/9.
C 
C       * GET VOLUME POROSITY (VOLUME OCCUPIED BY WATER, AIR OR ICE) AS 
C       * FUNCTION OF SOIL TYPE.
C 
        XVOID = 0.3+0.03*(ITEX+1) 
C 
C       * CALCULATE VOLUMETRIC HEAT CAPACITY OF SOIL (CONTRIBUTION OF AIR 
C       * TO TOTAL IGNORED DUE ITS MUCH SMALLER HEAT CAPACITY). 
C 
        HCAP = XVOID*(WL(I)*CVWTR+WF(I)*CVICE) + (1.-XVOID)*CVMNRL
C 
C       * NOW CALCULATE THERMAL CONDUCTIVITY FOR THE PARTICULAR SOIL
C       * TYPE OF THE GRID POINT. 
C 
        HCOND = HCAP * FSK((WL(I)+WF(I))*XVOID) * SKRAT 
C 
C       * CALCULATE "EFFECTIVE" HEAT CAPACITY WITH CORRECTION FOR SNOW
C       * AND NORMALIZATION.
C 
        C(I)=SQRT(HCOND*HCAP*DAYLNT/(2.*PI)) 
        C(I)=C(I)*(1.-SNOFAC(I))+CSNO*SNOFAC(I) 
        C(I)=(MAX(12.0E4,C(I)))/SQRT(2.) 
      ENDIF
C
C     * GENERAL CALCULATIONS FOR LAND OR SEA-ICE.
C 
        IF(GC(I).NE.0.)                                       THEN
C 
C         * SENSIBLE HEAT. 
C 
          RAUS=PRESS(I)/(RGAS*TS(I)) 
          HSENS=-CPRES*RAUS*TFLXS(I) 
C
C         * LATENT HEAT.
C         * THE FOLLOWING PICE OF CODE WAS EXTRACTED FROM THE
C         * FUNCTION "HTOCP" AND IN-LINED, SO THAT THIS LOOP WOULD
C         * VECTORIZE.
C
          IF(GT(I).GT.TFREZ)    THEN
            HL=HV*(1.-SNOFAC(I)) +HS*SNOFAC(I) 
          ELSE
            HL=HS
          ENDIF
          EVAPR(I)=-RAUS*QFLXS(I) 
          HLAT=HL*EVAPR(I) 
C
C         * SHORTWAVE FLUX INTO THE GROUND.
C
          FSG=MAX(FSC(I)*CSZ(I), 0.)
C 
C         * LONGWAVE FLUX FROM THE GROUND.
C 
          FLG=SBC* GT(I)**4
          IF(ICLW.EQ.0) FLG=0.
C 
C         * TIME STEP THE GROUND TEMPERATURE. 
C
          GTLAST(I)=GT(I)
          GT(I)=GT(I)+(FSG+FDL(I)-FLG-HSENS-HLAT+HSEA(I))*DELT/C(I)
        ENDIF
C
C       * FORCE-RESTORE OVER LAND.
C 
        IF(GC(I).LE.-0.5)                                       THEN
          GT(I)=GT(I)+2.*PI*DELT*(GTA(I)-GT(I))/DAYLNT
        ENDIF
  100 CONTINUE
C---------------------------------------------------------------------
      RETURN
      END 
