      PROGRAM LLOFSIC                                                                                                               
C     PROGRAM LLOFSIC (TS,    ZSNO,    GC,       SIC,       OUTPUT,     )       D2                                                  
C    1           TAPE1=TS, TAPE2=ZSNO, TAPE3=GC, TAPE4=SIC, TAPE6=OUTPUT)                                                           
C     -------------------------------------------------------                   D2                                                  
C                                                                               D2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       D2                                                  
C     MAR 05/99 - M.LAZARE.                                                     D2                                                  
C                 THIS IS SIMILAR TO LLCOSIC EXCEPT:                            D2                                                  
C                 - FOSTER+DAVEY SNOW DEPTH PASSED IN.                          D2                                                  
C                   (NOT USED FOR THE MOMENT, HOWEVER).                         D2                                                  
C                 - RESIDUAL TERM ADDED TO EQUATIONS (14/10 FOR S/N).           D2                                                  
C                 - WORKS FOR OFFSET L-L GRIDS INSTEAD OF "REGULAR"             D2                                                  
C                   GRIDS AS IN LLCOSIC.                                        D2                                                  
C                 - EXTRA ITTERATION ADDED AT END FOR PERMANENTLY ICE-          D2                                                  
C                   COVERED POINTS.                                             D2                                                  
C                 - GROWTH/MELT FORMULAS MODIFIED TO PREDICT MASS AT END        D2                                                  
C                   OF MONTH RATHER THAN MONTHLY MEAN, AND PROPER LATENT        D2                                                  
C                   HEAT OF FUSION CONSTANT USED IN CALCULATION OF "AR",        D2                                                  
C                   RESULTING IN A NEW, HIGHER VALUE AND DIFFERENT              D2                                                  
C                   FUNCTIONAL RELATIONSHIP ON GROWTH SIDE.                     D2                                                  
C                 - BOUND OF +2 ON TEMPERATURE ABOVE T0 IMPOSED ON MELT         D2                                                  
C                   SIDE, TO GIVE MORE REALISTIC MELT RATES OVER SUMMER         D2                                                  
C                   SEASON.                                                     D2                                                  
C     JAN 29/92 - E. CHAN - PREVIOUS VERSION LLCOSIC.                                                                               
C                                                                               D2                                                  
CLLOFSIC - COMPUTES SEA-ICE AMOUNT FROM SNOW DEPTH, GC AND SURFACE TEMP 2  1    D1                                                  
C                                                                               D3                                                  
CAUTHORS - M.LAZARE, F.MAJAESS                                                  D3                                                  
C                                                                               D3                                                  
CPURPOSE - COMPUTES SEA-ICE AMOUNT (SIC) FROM  GROUND COVER (GC), SURFACE       D3                                                  
C          TEMPERATURE (TS) AND SNOW DEPTH (OFFSET GRIDSAS INPUT)               D3                                                  
C          FOR EACH GRID POINT, ONE OF THREE CASES OCCURS:                      D3                                                  
C            1 - NO  SEA-ICE  COVER  ALL  YEAR (I.E. ALL GC VALUES DIFFERENT    D3                                                  
C                FROM ONE), SIC VALUES ARE SET TO ZEROS.                        D3                                                  
C            2 - SEA-ICE COVER FOR PART(S) OF THE  YEAR (I.E GC VALUES EQUAL    D3                                                  
C                TO AND DIFFERENT FROM ONES). IN THIS CASE THE  SIC VALUE OF    D3                                                  
C                THE FIRST MONTH (STARTING FROM JANUARY) WITH GC VALUE EQUAL    D3                                                  
C                TO ONE IS COMPUTED.                                            D3                                                  
C                           _                                                   D3                                                  
C                       T0   IS THE FREEZING TEMPERATURE OF SEA WATER           D3                                                  
C                               (271.2 DEG K),                                  D3                                                  
C                       SIC AND SNO ARE SEA-ICE MASS AND SNOW MASS,             D3                                                  
C                       RESPECTIVELY (KG/M2),                                   D3                                                  
C                       DENI AND DENS ARE SEA-ICE AND SNOW DENSITY,             D3                                                  
C                       RESPECTIVELY (KG/M3),                                   D3                                                  
C                       CONI AND CONS ARE SEA-ICE AND SNOW THERMAL              D3                                                  
C                       CONDUCTIVITY, RESPECTIVELY (W/M-S-K),                   D3                                                  
C                       DIM  IS THE NUMBER OF DAYS IN THE CURRENT MONTH,        D3                                                  
C                       R    IS A PARAMETER EQUAL TO 1062.8                     D3                                                  
C                       R = 2.*RHOI*CONI*86400/HF                               D3                                                  
C                       HF   = LATENT HEAT OF FUSION                            D3                                                  
C                AND    M0   IS THE INITIAL ESTIMATE OF  SEA-ICE  AMOUNT  FOR   D3                                                  
C                            THE MONTH WHOSE SEA-ICE IS TO BE COMPUTED.         D3                                                  
C                            M0 IS SET TO ZERO INITIALLY.                       D3                                                  
C                                                                               D3                                                  
C                THE  SAME COMPUTATIONAL  PROCEDURE IS  REPEATED FOR  ALL THE   D3                                                  
C                FOLLOWING  MONTH(S).                                           D3                                                  
C                WITH GC VALUE(S) EQUAL TO ONE(S)  UNTIL (IF) GC VALUE REVERT   D3                                                  
C                BACK TO A VALUE DIFFERENT FROM ONE.                            D3                                                  
C                THE LAST MONTH COMPUTED SIC VALUE, IN A GIVEN PERIOD OF SEA-   D3                                                  
C                ICE COVERAGE, IS  SET EQUAL TO  HALF  THAT OF  THE  PREVIOUS   D3                                                  
C                MONTH (IF IT IS GREATER THAN OR EQUAL TO 90).                  D3                                                  
C                THE SAME PROCESS IS REPEATED WITH THE NEXT SET, (IF ANY), OF   D3                                                  
C                MONTHS WITH GC VALUE(S) EQUAL TO ONE(S).                       D3                                                  
C            3 - SEA-ICE COVER ALL YEAR (I.E ALL GC VALUES ARE EQUAL TO ONE),   D3                                                  
C                IN THIS CASE  THE MONTH  WITH  LOWEST NEIGHBOURING  POINT(S)   D3                                                  
C                SEA-ICE AMOUNT AVERAGE IS USED TO CHOOSE THE STARTING  MONTH   D3                                                  
C                FROM WHICH  THE COMPUTATION IS TO  START WHERE,  THE INITIAL   D3                                                  
C                GUESS OF SEA-ICE AMOUNT (M0) IS TAKEN:                         D3                                                  
C                   FROM THE  SIC AVERAGE  OF NEIGBOURING POINTS WITH  SIC>0,   D3                                                  
C                OR IF THERE IS  AT LEAST ONE MONTH WITH  ZERO AVERAGE  OR AT   D3                                                  
C                   LEAST ONE  OPEN WATER NEIGHBOURING POINT IT IS SET TO 45,   D3                                                  
C                OR IF NONE OF THE ABOVE CONDITIONS IS MET AND ALL NEIGHBOUR-   D3                                                  
C                   ING POINTS ARE  LAND/UNCOMPUTED  ALL YEAR AROUND  SEA-ICE   D3                                                  
C                   COVERED; IN THIS CASE  M0  IS COMPUTED  FROM  THE AVERAGE   D3                                                  
C                   OF COMPUTED SEA-ICE  AMOUNT FOR POINTS  AROUND  THE  SAME   D3                                                  
C                   LATITUDINAL CIRCLE.                                         D3                                                  
C                                                                               D3                                                  
C                NOTE - M0 >= 45 (KG/M**2).                                     D3                                                  
C                                                                               D3                                                  
C                THE AMOUNT OF  SEA-ICE  (SIC)  FOR  THE FOLLOWING  MONTH  IS   D3                                                  
C                COMPUTED BY APPLYING THE SAME FORMULA AS IN CASE 2 ABOVE.      D3                                                  
C                                                                               D3                                                  
C                THEN, THE COMPUTATION  IS REPEATED FOR  THE FOLLOWING MONTHS   D3                                                  
C                THE SEA-ICE AMOUNT FOR THE LAST MONTH, (THAT PRECEEDING  THE   D3                                                  
C                MONTH FROM WHICH THE SIC COMPUTATION STARTED), IS SET  EQUAL   D3                                                  
C                TO THE AVERAGE OF THE INITIAL AND FINAL M0 VALUES.             D3                                                  
C                NOTE - THE FINAL VALUE FOR M0 IS UPDATED FROM  THE  PREVIOUS   D3                                                  
C                       MONTH SEA-ICE AMOUNT AND M0 VALUES AS DESCRIBED ABOVE.  D3                                                  
C                                                                               D3                                                  
C          NOTE - SURFACE TEMPERATURE VALUES (TS) SHOULD BE IN DEG KELVIN.      D3                                                  
C                 SIC VALUE >= 45 (KG/M**2) FOR  EVERY MONTH  WITH  GC  VALUE   D3                                                  
C                 EQUAL TO ONE.                                                 D3                                                  
C                 GLOBAL SEA-ICE AMOUNT COMPUTATION OF POINTS  SATISFYING THE   D3                                                  
C                 CONDITIONS OF  CASES 1 AND 2 ABOVE  ARE DONE  FIRST  BEFORE   D3                                                  
C                 DOING THOSE FOR POINTS SATISFYING CASE 3 CONDITIONS.          D3                                                  
C                 ALSO, THE STARTING MONTH FOR  CASE 3 POINTS IN THE NORTHERN   D3                                                  
C                 HEMISPHERE IS FORCED TO OCCUR IN  SEPTEMBER  OR AT A  LATER   D3                                                  
C                 MONTH. THIS  "TUNNING"  CONSTRAINT  CAN  BE  ELIMINATED  BY   D3                                                  
C                 COMMENTING  THE  LINE  FOLLOWING  THE  CONTINUE   STATEMENT   D3                                                  
C                 LABELLED 518 IN THE PROGRAM SOURCE CODE.                      D3                                                  
C                                                                               D3                                                  
CINPUT FILES...                                                                 D3                                                  
C                                                                               D3                                                  
C      TS  = CONTAINS MONTHLY AVERAGE SURFACE TEMPERATURES (DEG K)              D3                                                  
C      ZSNO= CONTAINS MONTHLY AVERAGE SNOW DEPTH (CMS) FROM FOSTER/DAVEY        D3                                                  
C      GC  = CONTAINS MONTHLY         GROUND COVER VALUES                       D3                                                  
C            (-1 FOR LAND, 0 FOR SEA, +1 FOR ICE)                               D3                                                  
C                                                                               D3                                                  
COUTPUT FILE...                                                                 D3                                                  
C                                                                               D3                                                  
C      SIC = CONTAINS MONTHLY COMPUTED SEA-ICE AMOUNTS (KG/M**2)                D3                                                  
C----------------------------------------------------------------------------                                                       
C                                                                                                                                   
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      REAL TS(65341,12), GC(65341,12), SIC(65341,12),                                                                               
     1     ZSNO(65341,12), SNO(65341,12),                                                                                           
     2     TGC(12), TTS(12), TSIC(12), TAVG(12)                                                                                     
C                                                                                                                                   
      REAL AM0F(65341)                                                                                                              
      INTEGER NPERM(65341),M0MINM(65341),NITER(65341)                                                                               
      REAL DIM(12)                                                                                                                  
C                                                                                                                                   
      INTEGER NDAT(65341), ITIM(12), IWTR(12)                                                                                       
      LOGICAL OK                                                                                                                    
      COMMON/ICOM/IBUF(8),IDAT(65341)                                                                                               
C                                                                                                                                   
      DATA MAXX/65341/                                                                                                              
C                                                                                                                                   
C     * INITIALIZE T0 AND R VALUES.                                                                                                 
C                                                                                                                                   
      DATA DIM /31.E0, 28.E0, 31.E0, 30.E0, 31.E0, 30.E0, 31.E0, 31.E0,                                                             
     &          30.E0, 31.E0, 30.E0, 31.E0/                                                                                         
      DATA DENI, CONI, HF/913.E0, 2.25E0, 3.34E5/                                                                                   
      DATA AR, AT0, SNOLIM, SNOMAX/1062.8E0, 271.2E0, 10.E0, 10.E0/                                                                 
C--------------------------------------------------------------------                                                               
      NFF=5                                                                                                                         
      CALL JCLPNT(NFF,1,2,3,4,6)                                                                                                    
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
      REWIND 3                                                                                                                      
      AMAXDF0=0.0E0                                                                                                                 
      AMAXDF1=0.0E0                                                                                                                 
      AMAXDF2=0.0E0                                                                                                                 
      AMAXDF3=0.0E0                                                                                                                 
C                                                                                                                                   
C     * READ IN GC, ZSNO AND TS VALUES FOR THE 12 MONTHS PERIOD.                                                                    
C                                                                                                                                   
      DO 10 M=1,12                                                                                                                  
        CALL GETFLD2(1,TS(1,M),NC4TO8("GRID"),-1,NC4TO8("  TS"),1,                                                                  
     +                                               IBUF,MAXX,OK)                                                                  
        IF(.NOT.OK) CALL                           XIT('LLOFSIC',-1)                                                                
        WRITE(6,6025) IBUF                                                                                                          
        ITIM(M)=IBUF(2)                                                                                                             
        IF(M.EQ.1)THEN                                                                                                              
          ILON=IBUF(5)                                                                                                              
          JLAT=IBUF(6)                                                                                                              
        ENDIF                                                                                                                       
        IF((IBUF(5).NE.ILON).OR.(IBUF(6).NE.JLAT))                                                                                  
     1    CALL                                     XIT('LLOFSIC',-2)                                                                
C                                                                                                                                   
        CALL GETFLD2(2,ZSNO(1,M),NC4TO8("GRID"),ITIM(M),NC4TO8(" SNO"),                                                             
     +                                                  1,IBUF,MAXX,OK)                                                             
        IF(.NOT.OK) CALL                           XIT('LLOFSIC',-3)                                                                
        WRITE(6,6025) IBUF                                                                                                          
        IF((IBUF(5).NE.ILON).OR.(IBUF(6).NE.JLAT))                                                                                  
     1    CALL                                     XIT('LLOFSIC',-4)                                                                
C                                                                                                                                   
        CALL GETFLD2(3,GC(1,M),NC4TO8("GRID"),ITIM(M),NC4TO8("  GC"),                                                               
     +                                                1,IBUF,MAXX,OK)                                                               
        IF(.NOT.OK) CALL                           XIT('LLOFSIC',-5)                                                                
        WRITE(6,6025) IBUF                                                                                                          
        IF((IBUF(5).NE.ILON).OR.(IBUF(6).NE.JLAT))                                                                                  
     1    CALL                                     XIT('LLOFSIC',-6)                                                                
  10  CONTINUE                                                                                                                      
C                                                                                                                                   
      NWDS=ILON*JLAT                                                                                                                
C                                                                                                                                   
      WRITE(6,6040) NWDS                                                                                                            
C                                                                                                                                   
C     * USE NDAT ARRAY AS A WORKING ARRAY HOLDING FOR EACH POINT A                                                                  
C     *      ZERO IF THE SEA-ICE AMOUNT IS COMPUTED,                                                                                
C     * OR A NEGATIVE INTEGER WHOSE ABSOLUTE VALUE SHOULD CORRESPOND TO                                                             
C     *      THE NUMBER OF NEIGHBOURING POINTS WHOSE SEA-ICE AMOUNT IS                                                              
C     *      NOT COMPUTED YET.                                                                                                      
C     * THE VALUES IN NDAT ARE USED IN COMPUTING M0 ESTIMATE FOR POINTS                                                             
C     * WITH ALL YEAR COVERED SEA-ICE CASE.                                                                                         
C     *                                                                                                                             
C     * INITIALIZE NDAT VALUES TO -9.                                                                                               
C                                                                                                                                   
      DO 15 N=1,NWDS                                                                                                                
        NDAT(N)=-9                                                                                                                  
C*********************************************************************                                                              
C       * INITIALIZE NEW FIELDS.                                                                                                    
C                                                                                                                                   
        M0MINM(N)=0                                                                                                                 
        AM0F(N)=0.E0                                                                                                                
C*********************************************************************                                                              
  15  CONTINUE                                                                                                                      
C                                                                                                                                   
C     * INITIALIZE THE SEA-ICE AMOUNT ARRAY "SIC".                                                                                  
C     * CONVERT SNOW DEPTH (CMS) TO SNOW MASS (KG/M2).                                                                              
C                                                                                                                                   
      DO 17 M=1,12                                                                                                                  
        DO 16 N=1,NWDS                                                                                                              
          SIC(N,M)=-5.0E0                                                                                                           
          DSNO=ZSNO(N,M)*0.01E0                                                                                                     
C         SNO(N,M)=(275./0.54)*(EXP(0.54*DSNO)-1.)                                                                                  
          SNO(N,M)=0.E0                                                                                                             
C                                                                                                                                   
C *********************************************                                                                                     
C         NOTE!: NO SNOW IN THIS CASE!!!                                                                                            
C**********************************************                                                                                     
  16    CONTINUE                                                                                                                    
  17  CONTINUE                                                                                                                      
C                                                                                                                                   
      NTTIME=0                                                                                                                      
      JLATF=JLAT/2                                                                                                                  
C                                                                                                                                   
  20  CONTINUE                                                                                                                      
      NTTIME=NTTIME+1                                                                                                               
      NPTSC=0                                                                                                                       
C                                                                                                                                   
      WRITE(6,6045) NTTIME, NPTSC                                                                                                   
C                                                                                                                                   
C     * ALTERNATE PROCESSING BETWEEN NORTHERN AND SOUTHERN HEMISPHERE                                                               
C     * STARTING WITH THE SOUTHERN ONE. GRID POINTS ARE PROCESSED FROM                                                              
C     * THE EQUATOR TOWARDS THE POLE.                                                                                               
C                                                                                                                                   
      DO 730 KK=1,2                                                                                                                 
       IF(KK.EQ.1)THEN                                                                                                              
        LAT1=JLATF                                                                                                                  
        LAT2=1                                                                                                                      
        LAT3=-1                                                                                                                     
       ELSE                                                                                                                         
        LAT1=JLATF+1                                                                                                                
        LAT2=JLAT                                                                                                                   
        LAT3=1                                                                                                                      
       ENDIF                                                                                                                        
C                                                                                                                                   
       DO 720 J=LAT1,LAT2,LAT3                                                                                                      
       IF((J.LT.1).OR.(J.GT.JLAT)) CALL            XIT('LLOFSIC',-7)                                                                
       NTIME=0                                                                                                                      
  25   CONTINUE                                                                                                                     
       ISTRT=ILON                                                                                                                   
       IFNSH=1                                                                                                                      
       NTIME=NTIME+1                                                                                                                
  26   NSICNC=0                                                                                                                     
       NSICPC=0                                                                                                                     
       NSICNW=0                                                                                                                     
       IINC=ISTRT                                                                                                                   
       ISTRT=IFNSH                                                                                                                  
       IFNSH=IINC                                                                                                                   
       IF(ISTRT.LE.IFNSH)THEN                                                                                                       
         IINC=1                                                                                                                     
       ELSE                                                                                                                         
         IINC=-1                                                                                                                    
       ENDIF                                                                                                                        
C                                                                                                                                   
C     * WITHIN A GIVEN LATITUDE, GRID POINTS SEA-ICE COMPUTATION IS PROCESSED                                                       
C     * ON THE GREENWICH MERIDIAN POINT, (ZERO DEGREE LONGITUDE), THEN                                                              
C     * ALTERNATELY ON POINTS LOCATED TO THE EAST AND WEST OF IT.                                                                   
C                                                                                                                                   
       DO 710 I=ISTRT,IFNSH,IINC                                                                                                    
C                                                                                                                                   
C       * FOR EACH GRID POINT...                                                                                                    
C                                                                                                                                   
        IF(I.EQ.1)THEN                                                                                                              
          IM1=ILON                                                                                                                  
        ELSE                                                                                                                        
          IM1=I-1                                                                                                                   
        ENDIF                                                                                                                       
        IF(I.EQ.ILON)THEN                                                                                                           
          IP1=1                                                                                                                     
        ELSE                                                                                                                        
          IP1=I+1                                                                                                                   
        ENDIF                                                                                                                       
        N=(J-1)*ILON+I                                                                                                              
C                                                                                                                                   
        IF(N.GT.NWDS)THEN                                                                                                           
          PRINT *,' I,J,N,NWDS,IM1,IP1,LAT1,LAT2,LAT3,ISTRT,IFNSH,',                                                                
     1     'IINC= ',I,J,N,NWDS,IM1,IP1,LAT1,LAT2,LAT3,ISTRT,IFNSH,IINC                                                              
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * SKIP POINT(S) WHOSE SEA-ICE AMOUNT ALREADY COMPUTED.                                                                      
C                                                                                                                                   
        IF (NDAT(N).GE.0) GO TO 710                                                                                                 
C                                                                                                                                   
C       * DO COMPUTATIONS VIA TEMPORARY ARRAYS TGC, TTS AND TSIC.                                                                   
C                                                                                                                                   
        DO 100 M=1,12                                                                                                               
          TGC(M) = GC(N,M)                                                                                                          
          TTS(M) = TS(N,M)                                                                                                          
 100    CONTINUE                                                                                                                    
C                                                                                                                                   
C       * CHECK GROUND COVER VALUES IN TGC.                                                                                         
C                                                                                                                                   
        ISUM=0                                                                                                                      
        DO 110 M=1,12                                                                                                               
          IF(TGC(M).NE.1.0E0)THEN                                                                                                   
            TSIC(M)=0.0E0                                                                                                           
          ELSE                                                                                                                      
            TSIC(M)=-1.0E0                                                                                                          
            ISUM=ISUM+1                                                                                                             
          ENDIF                                                                                                                     
 110    CONTINUE                                                                                                                    
C                                                                                                                                   
C       * NO ICE-COVER CASE                                                                                                         
C                                                                                                                                   
        IF(ISUM.EQ.0) GO TO 600                                                                                                     
C                                                                                                                                   
C       * ICE-COVER ALL YEAR (SKIP COMPUTATION IF THIS IS THE FIRST TIME                                                            
C       *                     THROUGH LOOP 710 I.E. COMPUTE SIC FIRST FOR                                                           
C       *                     THOSE POINTS WITH NO ICE COVER OR ICE-COVERED                                                         
C       *                     FOR PART(S) OF THE YEAR)                                                                              
C                                                                                                                                   
        IF(ISUM.EQ.12)THEN                                                                                                          
          NPERM(N)=1                                                                                                                
          IF((NTIME.EQ.1).OR.(NTTIME.LE.1))THEN                                                                                     
            GO TO 710                                                                                                               
          ELSE                                                                                                                      
            GO TO 500                                                                                                               
          ENDIF                                                                                                                     
C*********************************************************************                                                              
C         * SET CONDITIONS FOR PERMANENT ICE-COVERAGE.                                                                              
C                                                                                                                                   
        ELSE                                                                                                                        
          NPERM(N)=0                                                                                                                
C*********************************************************************                                                              
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * ICE-COVER DURING PART(S) OF THE YEAR CASE                                                                                 
C                                                                                                                                   
 200    MPTR=0                                                                                                                      
        DO 210 M=1,12                                                                                                               
          MP1=M+1                                                                                                                   
          IF(MP1.GT.12)MP1=MP1-12                                                                                                   
          IF((TSIC(M).EQ.0.0E0).AND.(TSIC(MP1).LT.0.0E0))THEN                                                                       
            MPTR=M                                                                                                                  
            GO TO 220                                                                                                               
          ENDIF                                                                                                                     
 210    CONTINUE                                                                                                                    
        DO 215 M=1,12                                                                                                               
          MP1=M+1                                                                                                                   
          IF(MP1.GT.12)MP1=MP1-12                                                                                                   
          IF((TSIC(M).GT.0.E0).AND.(TSIC(MP1).EQ.0.E0))THEN                                                                         
            MM1=M-1                                                                                                                 
            IF(MM1.LE.0)MM1=MM1+12                                                                                                  
C                                                                                                                                   
C           * ADJUST THE LAST COMPUTED SEA-ICE AMOUNT.                                                                              
C                                                                                                                                   
            IF(TSIC(MM1).GE.90.E0)THEN                                                                                              
              TSIC(M)=TSIC(MM1)/2.E0                                                                                                
            ELSE                                                                                                                    
C             IF((TSIC(MM1).EQ.0.0).AND.(TSIC(M).GE.90.))THEN                                                                       
C               TSIC(M)=0.5*TSIC(M)                                                                                                 
              IF((TSIC(MM1).EQ.0.0E0).AND.(TSIC(M).GE.45.E0))THEN                                                                   
                TSIC(M)=TSIC(M)                                                                                                     
              ELSE                                                                                                                  
                TSIC(M)=45.0E0                                                                                                      
              ENDIF                                                                                                                 
            ENDIF                                                                                                                   
          ENDIF                                                                                                                     
 215    CONTINUE                                                                                                                    
        GO TO 600                                                                                                                   
 220    CONTINUE                                                                                                                    
        AM0=0.0E0                                                                                                                   
        IF(N.GT.NWDS/2) THEN                                                                                                        
          RES=10.E0                                                                                                                 
        ELSE                                                                                                                        
          RES=14.E0                                                                                                                 
        ENDIF                                                                                                                       
        IF(MPTR.LE.0) CALL                         XIT('LLOFSIC',-8)                                                                
        DO 230 M=MPTR,(MPTR+11)                                                                                                     
          MP1=M+1                                                                                                                   
          IF(MP1.GT.12)MP1=MP1-12                                                                                                   
          IF(TSIC(MP1).GE.0.0E0)GO TO 200                                                                                           
          IF(SNO(N,MP1).GT.SNOMAX) THEN                                                                                             
            DENS=0.54E0*SNO(N,MP1)/LOG(1.E0+0.54E0*SNO(N,MP1)/275.E0)                                                               
          ELSE                                                                                                                      
            DENS=275.E0                                                                                                             
          ENDIF                                                                                                                     
          CONS=2.805E-6*DENS**2                                                                                                     
          RESF=RES*86400.E0*DIM(MP1)/(2.E0*HF)                                                                                      
C                                                                                                                                   
          X=(DENI*CONI)*SNO(N,MP1)/(DENS*CONS)                                                                                      
          IF(AT0.GE.TTS(MP1))THEN                                                                                                   
            DT=AT0-TTS(MP1)                                                                                                         
            DTFAC=AR*DIM(MP1)*DT                                                                                                    
            AM1=SQRT((X+AM0-RESF)**2 + DTFAC) - (X+RESF)                                                                            
          ELSE                                                                                                                      
            DT=MIN(TTS(MP1)-AT0,2.E0)                                                                                               
            XMELT=SQRT(AR*DIM(MP1)*DT)                                                                                              
            IF(XMELT.GT.SNO(N,MP1)) THEN                                                                                            
              AM1=AM0+SNO(N,MP1)-XMELT                                                                                              
            ELSE                                                                                                                    
              AM1=AM0                                                                                                               
            ENDIF                                                                                                                   
          ENDIF                                                                                                                     
          TSIC(MP1)=0.5E0*(AM0+AM1)                                                                                                 
          AM0=AM1                                                                                                                   
          IF(TSIC(MP1).LT.45.0E0)TSIC(MP1)=45.0E0                                                                                   
          IF(AM0.LT.45.0E0)AM0=45.0E0                                                                                               
 230    CONTINUE                                                                                                                    
        GO TO 200                                                                                                                   
C                                                                                                                                   
C       * ICE-COVERED ALL YEAR CASE COMPUTATIONS                                                                                    
C                                                                                                                                   
 500    CONTINUE                                                                                                                    
C       PRINT *,'  ICE COVERED ALL YEAR CASE,J,I= ',J,', ',I                                                                        
C                                                                                                                                   
C       * DEFER SIC COMPUTATION FOR THE CURRENT POINT IF THERE IS UNCOMPUTED                                                        
C       * NEIGHBOURING POINT(S) WITH MORE NEIGHBOURS WHOSE SIC ALREADY                                                              
C       * COMPUTED.                                                                                                                 
C                                                                                                                                   
        IF(J.GT.1)THEN                                                                                                              
         JB=(J-2)*ILON                                                                                                              
         IF((NDAT(JB+IM1).LT.0).AND.(NDAT(JB+IM1).GT.NDAT(N)))GO TO 710                                                             
         IF((NDAT(JB+ I ).LT.0).AND.(NDAT(JB+ I ).GT.NDAT(N)))GO TO 710                                                             
         IF((NDAT(JB+IP1).LT.0).AND.(NDAT(JB+IP1).GT.NDAT(N)))GO TO 710                                                             
        ENDIF                                                                                                                       
C                                                                                                                                   
        JB=(J-1)*ILON                                                                                                               
        IF((NDAT(JB+IM1).LT.0).AND.(NDAT(JB+IM1).GT.NDAT(N)))GO TO 710                                                              
        IF((NDAT(JB+IP1).LT.0).AND.(NDAT(JB+IP1).GT.NDAT(N)))GO TO 710                                                              
C                                                                                                                                   
        IF(J.LT.JLAT)THEN                                                                                                           
         JB=J*ILON                                                                                                                  
         IF((NDAT(JB+IM1).LT.0).AND.(NDAT(JB+IM1).GT.NDAT(N)))GO TO 710                                                             
         IF((NDAT(JB+ I ).LT.0).AND.(NDAT(JB+ I ).GT.NDAT(N)))GO TO 710                                                             
         IF((NDAT(JB+IP1).LT.0).AND.(NDAT(JB+IP1).GT.NDAT(N)))GO TO 710                                                             
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * COMPUTE MONTHLY AVERAGES OF NEIGHBOURING POINTS WITH SIC>0.                                                               
C                                                                                                                                   
        ISKP=0                                                                                                                      
 503    ICM=0                                                                                                                       
        MXLAND=0                                                                                                                    
        MXWTR=0                                                                                                                     
        DO 505 M=1,12                                                                                                               
          NPTS=0                                                                                                                    
          ASUM=-1.0E0                                                                                                               
          ITLAND=0                                                                                                                  
          IWTR(M)=0                                                                                                                 
C                                                                                                                                   
          JB=(J-1)*ILON                                                                                                             
          IF(GC(JB+IM1,M).LT.-0.1E0)THEN                                                                                            
            ITLAND=ITLAND+1                                                                                                         
          ELSE                                                                                                                      
            IF(ABS(GC(JB+IM1,M)).LT.0.1E0) IWTR(M)=IWTR(M)+1                                                                        
          ENDIF                                                                                                                     
          IF((NDAT(JB+IM1).GE.0).AND.(SIC(JB+IM1,M).GE.0.0E0))THEN                                                                  
            IF(ASUM.LT.0.0E0) ASUM=0.0E0                                                                                            
            IF(SIC(JB+IM1,M).GT.0.0E0)THEN                                                                                          
              NPTS=NPTS+1                                                                                                           
              ASUM=ASUM+SIC(JB+IM1,M)                                                                                               
            ENDIF                                                                                                                   
          ENDIF                                                                                                                     
          IF(GC(JB+IP1,M).LT.-0.1E0)THEN                                                                                            
            ITLAND=ITLAND+1                                                                                                         
          ELSE                                                                                                                      
            IF(ABS(GC(JB+IP1,M)).LT.0.1E0) IWTR(M)=IWTR(M)+1                                                                        
          ENDIF                                                                                                                     
          IF((NDAT(JB+IP1).GE.0).AND.(SIC(JB+IP1,M).GE.0.0E0))THEN                                                                  
            IF(ASUM.LT.0.0E0) ASUM=0.0E0                                                                                            
            IF(SIC(JB+IP1,M).GT.0.0E0)THEN                                                                                          
              NPTS=NPTS+1                                                                                                           
              ASUM=ASUM+SIC(JB+IP1,M)                                                                                               
            ENDIF                                                                                                                   
          ENDIF                                                                                                                     
C                                                                                                                                   
          IF(J.GT.1)THEN                                                                                                            
            JB=(J-2)*ILON                                                                                                           
            IF(GC(JB+IM1,M).LT.-0.1E0)THEN                                                                                          
              ITLAND=ITLAND+1                                                                                                       
            ELSE                                                                                                                    
              IF(ABS(GC(JB+IM1,M)).LT.0.1E0) IWTR(M)=IWTR(M)+1                                                                      
            ENDIF                                                                                                                   
            IF(GC(JB+ I ,M).LT.-0.1E0)THEN                                                                                          
              ITLAND=ITLAND+1                                                                                                       
            ELSE                                                                                                                    
              IF(ABS(GC(JB+ I ,M)).LT.0.1E0) IWTR(M)=IWTR(M)+1                                                                      
            ENDIF                                                                                                                   
            IF(GC(JB+IP1,M).LT.-0.1E0)THEN                                                                                          
              ITLAND=ITLAND+1                                                                                                       
            ELSE                                                                                                                    
              IF(ABS(GC(JB+IP1,M)).LT.0.1E0) IWTR(M)=IWTR(M)+1                                                                      
            ENDIF                                                                                                                   
          ENDIF                                                                                                                     
C                                                                                                                                   
          IF(J.LT.JLAT)THEN                                                                                                         
            JB=J*ILON                                                                                                               
            IF(GC(JB+IM1,M).LT.-0.1E0)THEN                                                                                          
              ITLAND=ITLAND+1                                                                                                       
            ELSE                                                                                                                    
              IF(ABS(GC(JB+IM1,M)).LT.0.1E0) IWTR(M)=IWTR(M)+1                                                                      
            ENDIF                                                                                                                   
            IF(GC(JB+ I ,M).LT.-0.1E0)THEN                                                                                          
              ITLAND=ITLAND+1                                                                                                       
            ELSE                                                                                                                    
              IF(ABS(GC(JB+ I ,M)).LT.0.1E0) IWTR(M)=IWTR(M)+1                                                                      
            ENDIF                                                                                                                   
            IF(GC(JB+IP1,M).LT.-0.1E0)THEN                                                                                          
              ITLAND=ITLAND+1                                                                                                       
            ELSE                                                                                                                    
              IF(ABS(GC(JB+IP1,M)).LT.0.1E0) IWTR(M)=IWTR(M)+1                                                                      
            ENDIF                                                                                                                   
          ENDIF                                                                                                                     
C                                                                                                                                   
C         * IF IT IS MORE DESIRABLE TO COMPUTE THE MONTHLY AVERAGES, (IF                                                            
C         * POSSIBLE), FROM SEA-ICE COVERED ADJACENT POINTS ON THE SAME                                                             
C         * LATITUDE, THE USER MUST ACTIVATE THE FIVE FORTRAN STATEMENT                                                             
C         * LINES OF CODE WHICH FOLLOW AND ARE COMMENTED.                                                                           
C                                                                                                                                   
C         IF((ASUM.GT.0.0).AND.(ISKP.LT.6))THEN                                                                                     
C           ISKP=6                                                                                                                  
C           GO TO 503                                                                                                               
C         ENDIF                                                                                                                     
C                                                                                                                                   
C         * RESTRICT MONTHLY AVERAGES COMPUTATION TO POINTS ON THE SAME                                                             
C         * LATITUDE WITH SIC>0; (IF POSSIBLE).                                                                                     
C                                                                                                                                   
C         IF(ISKP.GT.0) GO TO 504                                                                                                   
C                                                                                                                                   
          IF(J.GT.1)THEN                                                                                                            
            JB=(J-2)*ILON                                                                                                           
            IF((NDAT(JB+IM1).GE.0).AND.(SIC(JB+IM1,M).GE.0.0E0))THEN                                                                
              IF(ASUM.LT.0.0E0) ASUM=0.0E0                                                                                          
              IF(SIC(JB+IM1,M).GT.0.0E0)THEN                                                                                        
                NPTS=NPTS+1                                                                                                         
                ASUM=ASUM+SIC(JB+IM1,M)                                                                                             
              ENDIF                                                                                                                 
            ENDIF                                                                                                                   
            IF((NDAT(JB+ I ).GE.0).AND.(SIC(JB+ I ,M).GE.0.0E0))THEN                                                                
              IF(ASUM.LT.0.0E0) ASUM=0.0E0                                                                                          
              IF(SIC(JB+ I ,M).GT.0.0E0)THEN                                                                                        
                NPTS=NPTS+1                                                                                                         
                ASUM=ASUM+SIC(JB+ I ,M)                                                                                             
              ENDIF                                                                                                                 
            ENDIF                                                                                                                   
            IF((NDAT(JB+IP1).GE.0).AND.(SIC(JB+IP1,M).GE.0.0E0))THEN                                                                
              IF(ASUM.LT.0.0E0) ASUM=0.0E0                                                                                          
              IF(SIC(JB+IP1,M).GT.0.0E0)THEN                                                                                        
                NPTS=NPTS+1                                                                                                         
                ASUM=ASUM+SIC(JB+IP1,M)                                                                                             
              ENDIF                                                                                                                 
            ENDIF                                                                                                                   
          ENDIF                                                                                                                     
C                                                                                                                                   
          IF(J.LT.JLAT)THEN                                                                                                         
            JB=J*ILON                                                                                                               
            IF((NDAT(JB+IM1).GE.0).AND.(SIC(JB+IM1,M).GE.0.0E0))THEN                                                                
              IF(ASUM.LT.0.0E0) ASUM=0.0E0                                                                                          
              IF(SIC(JB+IM1,M).GT.0.0E0)THEN                                                                                        
                NPTS=NPTS+1                                                                                                         
                ASUM=ASUM+SIC(JB+IM1,M)                                                                                             
              ENDIF                                                                                                                 
            ENDIF                                                                                                                   
            IF((NDAT(JB+ I ).GE.0).AND.(SIC(JB+ I ,M).GE.0.0E0))THEN                                                                
              IF(ASUM.LT.0.0E0) ASUM=0.0E0                                                                                          
              IF(SIC(JB+ I ,M).GT.0.0E0)THEN                                                                                        
                NPTS=NPTS+1                                                                                                         
                ASUM=ASUM+SIC(JB+ I ,M)                                                                                             
              ENDIF                                                                                                                 
            ENDIF                                                                                                                   
            IF((NDAT(JB+IP1).GE.0).AND.(SIC(JB+IP1,M).GE.0.0E0))THEN                                                                
              IF(ASUM.LT.0.0E0) ASUM=0.0E0                                                                                          
              IF(SIC(JB+IP1,M).GT.0.0E0)THEN                                                                                        
                NPTS=NPTS+1                                                                                                         
                ASUM=ASUM+SIC(JB+IP1,M)                                                                                             
              ENDIF                                                                                                                 
            ENDIF                                                                                                                   
          ENDIF                                                                                                                     
C                                                                                                                                   
 504      IF(NPTS.GE.1)THEN                                                                                                         
            TAVG(M)=ASUM/NPTS                                                                                                       
            ICM=ICM+1                                                                                                               
          ELSE                                                                                                                      
            TAVG(M)=ASUM                                                                                                            
          ENDIF                                                                                                                     
          IF(ITLAND.GT.MXLAND) MXLAND=ITLAND                                                                                        
          IF(IWTR(M).GT.MXWTR) MXWTR=IWTR(M)                                                                                        
 505    CONTINUE                                                                                                                    
C                                                                                                                                   
C       * SET INITIAL M0 BASED ON ICM VALUE.                                                                                        
C                                                                                                                                   
        IF(ICM.EQ.0)THEN                                                                                                            
C                                                                                                                                   
          IF((NTIME.GT.ILON).OR.((MXWTR.EQ.0).AND.                                                                                  
     1                           ((MXLAND.EQ.0).OR.                                                                                 
     2                            (NTIME.LE.(ILON-5))))) THEN                                                                       
C                                                                                                                                   
C         * DEFER COMPUTATION, IF NTIME.LE.ILON AND ALL NEIGHBOURING POINTS                                                         
C         * ARE SEA-ICE COVERED ALL YEAR AROUND AND NONE IS COMPUTED,                                                               
C         * OTHERWISE ABORT.                                                                                                        
C                                                                                                                                   
            IF(NTIME.LE.ILON)THEN                                                                                                   
              NSICNC=NSICNC+1                                                                                                       
              GO TO 710                                                                                                             
            ELSE                                                                                                                    
              PRINT *,' ALL YEAR SEA-ICE COVERED AND NONE IS COMPUTED',                                                             
     1                ' J,I,NTIME= ',J,', ',I,', NTIME= ',NTIME                                                                     
              CALL                                 XIT('LLOFSIC',-9)                                                                
            ENDIF                                                                                                                   
C                                                                                                                                   
          ELSE                                                                                                                      
C                                                                                                                                   
C           * COMPUTE POINTS WITH MORE NEIGHBOURING OPEN WATER POINTS FIRST.                                                        
C                                                                                                                                   
            IF((MXWTR.GT.0).AND.(((9-MXWTR)+1).GE.NTIME))THEN                                                                       
              NSICNW=NSICNW+1                                                                                                       
              GO TO 710                                                                                                             
            ENDIF                                                                                                                   
C                                                                                                                                   
C           * FIND MAXIMUM TS AND THE MONTH WHERE IT IS OCCURING (STARTING                                                          
C           * FROM JANUARY)                                                                                                         
C                                                                                                                                   
            MPTR=1                                                                                                                  
            AMAXTS=TTS(1)                                                                                                           
            DO 508 M=2,12                                                                                                           
              IF(TTS(M).GE.AMAXTS)THEN                                                                                              
                MPTR=M                                                                                                              
                AMAXTS=TTS(M)                                                                                                       
              ENDIF                                                                                                                 
 508        CONTINUE                                                                                                                
C                                                                                                                                   
            MPTRS=MPTR                                                                                                              
            DO 509 M=MPTRS,(MPTRS+11)                                                                                               
              MM=M                                                                                                                  
              IF(MM.GT.12) MM=MM-12                                                                                                 
              MMP1=MM+1                                                                                                             
              IF(MMP1.GT.12) MMP1=MMP1-12                                                                                           
              IF((TTS(MM).EQ.AMAXTS).AND.(TTS(MMP1).LT.AMAXTS))THEN                                                                 
                MPTR=MM                                                                                                             
                GO TO 510                                                                                                           
              ENDIF                                                                                                                 
 509        CONTINUE                                                                                                                
C                                                                                                                                   
 510        IF(MXWTR.GT.0)THEN                                                                                                      
C                                                                                                                                   
C             * THE CASE OF AT LEAST ONE NEIGHBOURING POINT WITH SEA WATER                                                          
C             * COVERAGE DURING THE YEAR.                                                                                           
C             * FIND THE MONTH WITH MAXIMUM SEA-WATER COVERAGE STARTING                                                             
C             * FROM THE MONTH WITH HIGHEST TS.                                                                                     
C                                                                                                                                   
              MPTRW=MPTR+1                                                                                                          
              IF(MPTRW.GT.12) MPTRW=MPTRW-12                                                                                        
              DO 511 M=(MPTR+1),(MPTR+12)                                                                                           
                MM=M                                                                                                                
                IF(MM.GT.12) MM=MM-12                                                                                               
                MMP1=MM+1                                                                                                           
                IF(MMP1.GT.12) MMP1=MMP1-12                                                                                         
                IF((IWTR(MM).EQ.MXWTR).AND.(IWTR(MMP1).LT.MXWTR))THEN                                                               
                  MPTRW=MM                                                                                                          
                  GO TO 501                                                                                                         
                ENDIF                                                                                                               
 511          CONTINUE                                                                                                              
C                                                                                                                                   
C             * THE CASE OF AT LEAST ONE NEIGHBOURING POINT WITH SEA WATER                                                          
C             * COVERAGE DURING THE YEAR AND AT LEAST ONE NEIGHBOURING                                                              
C             * COMPUTED SIC POINT.                                                                                                 
C             * FIND THE MONTH WITH MINIMUM MONTHLY SIC COVERAGE STARTING                                                           
C             * FROM THE MONTH WITH HIGHEST TS.                                                                                     
C                                                                                                                                   
 501          MPTRJ=0                                                                                                               
              IF(ICM.GT.0)THEN                                                                                                      
                DO 506 M=(MPTR+1),(MPTR+12)                                                                                         
                  MM=M                                                                                                              
                  IF(MM.GT.12) MM=MM-12                                                                                             
                  IF(TAVG(MM).LT.0.0E0) GO TO 506                                                                                   
                  MMP1=MM+1                                                                                                         
                  IF(MMP1.GT.12) MMP1=MMP1-12                                                                                       
                  IF(TAVG(MMP1).GT.TAVG(MM))THEN                                                                                    
                    MPTRJ=MM                                                                                                        
                    GO TO 512                                                                                                       
                  ENDIF                                                                                                             
 506            CONTINUE                                                                                                            
              ENDIF                                                                                                                 
C                                                                                                                                   
 512          CONTINUE                                                                                                              
              IF((ICM.GT.0).AND.(MPTRJ.GT.0))THEN                                                                                   
                MPTR=MPTRJ-1                                                                                                        
              ELSE                                                                                                                  
                MPTR=MPTRW-1                                                                                                        
              ENDIF                                                                                                                 
              IF(MPTR.LT.0) CALL                   XIT('LLOFSIC',-10)                                                               
              IF(MPTR.EQ.0) MPTR=12                                                                                                 
              AM00=45.0E0                                                                                                           
              GO TO 519                                                                                                             
C                                                                                                                                   
            ELSE                                                                                                                    
C                                                                                                                                   
C           * THE CASE OF ALL NEIGHBOURS ARE LAND , OR LAND AND UNCOMPUTED                                                          
C           * ALL YEAR SEA-ICE COVERED POINTS.                                                                                      
C                                                                                                                                   
              IF((((NTIME.LE.50).AND.(MXLAND.LT.8)).OR.                                                                             
     1                                (NTIME.LT.30)).OR.                                                                            
     2                                  (NPTSC.GT.0))THEN                                                                           
C                                                                                                                                   
C             * DEFER COMPUTATION IF THERE IS AT LEAST ONE UNCOMPUTED ALL                                                           
C             * YEAR SEA-ICE COVERED NEIGHBOURING POINT.                                                                            
C                                                                                                                                   
                NSICPC=NSICPC+1                                                                                                     
                GO TO 710                                                                                                           
                                                                                                                                    
              ELSE                                                                                                                  
                                                                                                                                    
C             * COMPUTE AVERAGE BASED ON SIC POINTS AROUND THE SAME                                                                 
C             * LATITUDE CIRCLE.                                                                                                    
C                                                                                                                                   
                IICM=0                                                                                                              
                MPTRS=MPTR                                                                                                          
                MPTR=0                                                                                                              
                AM00=-1.0E0                                                                                                         
                DO 514 M=1,12                                                                                                       
                  ASUM=0.0E0                                                                                                        
                  NPTS=0                                                                                                            
                  DO 513 II=1,ILON                                                                                                  
                    NN=(J-1)*ILON+II                                                                                                
                    IF((NDAT(NN).GE.0).AND.(SIC(NN,M).GT.0.0E0))THEN                                                                
                      NPTS=NPTS+1                                                                                                   
                      ASUM=ASUM+SIC(NN,M)                                                                                           
                    ENDIF                                                                                                           
 513              CONTINUE                                                                                                          
                  IF(NPTS.GE.1)THEN                                                                                                 
                    TAVG(M)=ASUM/NPTS                                                                                               
                    IICM=IICM+1                                                                                                     
                    IF((AM00.LT.0.E0).OR.(AM00.GE.TAVG(M)))THEN                                                                     
                      MPTR=M                                                                                                        
                      AM00=TAVG(M)                                                                                                  
                    ENDIF                                                                                                           
                  ELSE                                                                                                              
                    TAVG(M)=-1.0E0                                                                                                  
                  ENDIF                                                                                                             
 514            CONTINUE                                                                                                            
                IF(IICM.EQ.0)THEN                                                                                                   
                  IF(NTIME.LE.40)THEN                                                                                               
                    NSICPC=NSICPC+1                                                                                                 
                    GO TO 710                                                                                                       
                  ELSE                                                                                                              
                    IF((MPTRS.LE.0).OR.(MPTRS.GT.12))                                                                               
     1                CALL                         XIT('LLOFSIC',-11)                                                               
                    MPTR=MPTRS+2                                                                                                    
                    IF(MPTR.GT.12)MPTR=MPTR-12                                                                                      
C                   AM00=(((4565.-45.)*(273.-TTS(MPTR)))/40.)+45.0                                                                  
                    AM00=(((4565.E0-45.E0)*(273.E0-AMAXTS))/40.E0)+                                                                 
     &                45.0E0                                                                                                        
                    IF(AM00.LT.45.0E0) AM00=45.0E0                                                                                  
                    PRINT *,' FORMULA USED, J,I,AM00= ',J,', ',I,                                                                   
     1                                                  ', ',AM00                                                                   
C                   CALL                           XIT('LLOFSIC',-12)                                                               
                  ENDIF                                                                                                             
                ELSE                                                                                                                
                  IF(AM00.LT.45.0E0) AM00=45.0E0                                                                                    
                  PRINT *,'  MINN LAT CIRCLE AVRG IS USED J,I,AM00= ',                                                              
     1                    J,', ',I,', ',AM00                                                                                        
                  IF((MPTR.LE.0).OR.(MPTR.GT.12))                                                                                   
     1              CALL                           XIT('LLOFSIC',-13)                                                               
                  MPTR=MPTR-1                                                                                                       
                  IF(MPTR.LE.0)MPTR=MPTR+12                                                                                         
                  GO TO 519                                                                                                         
                ENDIF                                                                                                               
              ENDIF                                                                                                                 
            ENDIF                                                                                                                   
C                                                                                                                                   
          ENDIF                                                                                                                     
C                                                                                                                                   
        ELSE                                                                                                                        
                                                                                                                                    
C                                                                                                                                   
C       * CHECK THE MONTH WITH MINIMUM AVERAGE OF COMPUTED SIC>0.                                                                   
C                                                                                                                                   
          MPTR=0                                                                                                                    
          AM00=-1.0E0                                                                                                               
          IZAVG=0                                                                                                                   
          DO 515 M=1,12                                                                                                             
            IF(TAVG(M).LE.0.0E0)THEN                                                                                                
              IF(TAVG(M).EQ.0.0E0) IZAVG=IZAVG+1                                                                                    
              GO TO 515                                                                                                             
            ENDIF                                                                                                                   
            IF((AM00.LT.0.E0).OR.(AM00.GE.TAVG(M)))THEN                                                                             
              MPTR=M                                                                                                                
              AM00=TAVG(M)                                                                                                          
            ENDIF                                                                                                                   
 515      CONTINUE                                                                                                                  
C                                                                                                                                   
          IF(MPTR.EQ.0) CALL                       XIT('LLOFSIC',-14)                                                               
C                                                                                                                                   
C         * ENSURE COMPUTED SIC MINIMUM OCCUR IN THE SAME MONTH AS THE                                                              
C         * LAST MONTH WITH OPEN WATER NEIGHBOURING POINTS BEFORE THE                                                               
C         * SEA-ICE START FORMING.                                                                                                  
C                                                                                                                                   
          IF(IZAVG.GT.0)THEN                                                                                                        
            MPTRS=MPTR-1                                                                                                            
            IF(MPTRS.LE.0)MPTRS=MPTRS+12                                                                                            
            DO 516 M=MPTRS,(MPTRS-10),-1                                                                                            
              MM=M                                                                                                                  
              IF(MM.LE.0)MM=MM+12                                                                                                   
              IF(TAVG(MM).EQ.0.0E0)THEN                                                                                             
                MPTR=MM-1                                                                                                           
                IF(MPTR.LE.0)MPTR=MPTR+12                                                                                           
                AM00=45.0E0                                                                                                         
                GO TO 518                                                                                                           
              ENDIF                                                                                                                 
 516        CONTINUE                                                                                                                
          ELSE                                                                                                                      
C                                                                                                                                   
C           * DEFER COMPUTATION UNTIL ALL POINTS WITH OPEN WATER NEIGHBOURING                                                       
C           * POINTS ARE COMPUTED.                                                                                                  
C                                                                                                                                   
            IF((NTIME.LE.2).OR.((NTIME.LE.10).AND.(MXWTR.EQ.0)))                                                                    
     1                                                      GO TO 710                                                               
            MPTRP1=MPTR+1                                                                                                           
            IF(MPTRP1.GT.12)MPTRP1=MPTRP1-12                                                                                        
            IF(TAVG(MPTRP1).EQ.TAVG(MPTR))THEN                                                                                      
              PRINT *,' NOT LAST MINN. MNTHLY SIC AVRG,J,I= ',J,', ',I                                                              
              PRINT *,' TAVG(1-12)= ',(TAVG(M),M=1,12)                                                                              
            ELSE                                                                                                                    
              MPTR=MPTR-1                                                                                                           
              IF(MPTR.LE.0)MPTR=MPTR+12                                                                                             
            ENDIF                                                                                                                   
            IF(MXWTR.GT.0) AM00=45.0E0                                                                                              
C                                                                                                                                   
          ENDIF                                                                                                                     
C                                                                                                                                   
 518      CONTINUE                                                                                                                  
          IF((IZAVG.GT.0).AND.(KK.EQ.2).AND.(MPTR.LT.8)) MPTR=8                                                                     
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * COMPUTE GRID SEA-ICE AMOUNT IN "TSIC" ARRAY.                                                                              
C                                                                                                                                   
 519    CONTINUE                                                                                                                    
        IF((MPTR.LE.0).OR.(MPTR.GT.12)) CALL       XIT('LLOFSIC',-15)                                                               
        IF(AM00.LT.45.0E0)THEN                                                                                                      
          PRINT *,' M0 (NGBRS)!!, J,I,AM00= ',J,', ',I,', ',AM00                                                                    
          AM00=45.0E0                                                                                                               
        ENDIF                                                                                                                       
        AM0=AM00                                                                                                                    
        IF(N.GT.NWDS/2) THEN                                                                                                        
          RES=10.E0                                                                                                                 
        ELSE                                                                                                                        
          RES=14.E0                                                                                                                 
        ENDIF                                                                                                                       
        DO 520 M=MPTR,(MPTR+10)                                                                                                     
          MP1=M+1                                                                                                                   
          IF(MP1.GT.12)MP1=MP1-12                                                                                                   
C                                                                                                                                   
          RESF=RES*86400.E0*DIM(MP1)/(2.E0*HF)                                                                                      
          IF(SNO(N,MP1).GT.SNOMAX) THEN                                                                                             
            DENS=0.54E0*SNO(N,MP1)/LOG(1.E0+0.54E0*SNO(N,MP1)/275.E0)                                                               
          ELSE                                                                                                                      
            DENS=275.E0                                                                                                             
          ENDIF                                                                                                                     
          CONS=2.805E-6*DENS**2                                                                                                     
          X=(DENI*CONI)*SNO(N,MP1)/(DENS*CONS)                                                                                      
          IF(AT0.GE.TTS(MP1))THEN                                                                                                   
            DT=AT0-TTS(MP1)                                                                                                         
            DTFAC=AR*DIM(MP1)*DT                                                                                                    
            AM1=SQRT((X+AM0-RESF)**2 + DTFAC) - (X+RESF)                                                                            
          ELSE                                                                                                                      
            DT=MIN(TTS(MP1)-AT0,2.E0)                                                                                               
            XMELT=SQRT(AR*DIM(MP1)*DT)                                                                                              
            IF(XMELT.GT.SNO(N,MP1)) THEN                                                                                            
              AM1=AM0+SNO(N,MP1)-XMELT                                                                                              
            ELSE                                                                                                                    
              AM1=AM0                                                                                                               
            ENDIF                                                                                                                   
          ENDIF                                                                                                                     
          TSIC(MP1)=0.5E0*(AM0+AM1)                                                                                                 
          AM0=AM1                                                                                                                   
          IF(TSIC(MP1).LT.45.0E0)TSIC(MP1)=45.0E0                                                                                   
          IF(AM0.LT.45.0E0) AM0=45.0E0                                                                                              
 520    CONTINUE                                                                                                                    
C*********************************************************************                                                              
C       * SET VALUES FOR EXTRA YEAR INTEGRATION AT END FOR                                                                          
C       * PERMANENTLY-ICE COVERED POINTS.                                                                                           
C                                                                                                                                   
        M0MINM(N)=MPTR                                                                                                              
        AM0F(N)=AM0                                                                                                                 
C*********************************************************************                                                              
        MPTRM1=MPTR-1                                                                                                               
        IF(MPTRM1.LE.0)MPTRM1=MPTRM1+12                                                                                             
        MPTRP1=MPTR+1                                                                                                               
        IF(MPTRP1.GT.12)MPTRP1=MPTRP1-12                                                                                            
        TSIC(MPTR)=(AM0+AM00)/2.0E0                                                                                                 
        IF(TSIC(MPTR).LT.45.0E0) TSIC(MPTR)=45.0E0                                                                                  
        TSICVAL=(TSIC(MPTRM1)+TSIC(MPTRP1))/2.0E0                                                                                   
        IF(TSICVAL.LT.45.0E0) TSICVAL=45.0E0                                                                                        
C                                                                                                                                   
C       TTSICV=TSIC(MPTR)                                                                                                           
C       TSIC(MPTR)=TSICVAL                                                                                                          
C       TSICVAL=TTSICV                                                                                                              
C                                                                                                                                   
        IF(ABS(TSIC(MPTRM1)-TSIC(MPTRP1)).LT.1.E-50) GO TO 600                                                                      
C                                                                                                                                   
C       * CHECK FOR LOWEST ABSOLUTE ERROR                                                                                           
C                                                                                                                                   
        AMAXDIF=ABS(AM00-TSIC(MPTR))/ABS(TSIC(MPTRM1)-TSIC(MPTRP1))                                                                 
        IF(AMAXDIF.GT.AMAXDF0) AMAXDF0=AMAXDIF                                                                                      
        AMAXDIF=ABS(TSIC(MPTRM1)-TSIC(MPTRP1))                                                                                      
        IF(AMAXDIF.GT.AMAXDF1) AMAXDF1=AMAXDIF                                                                                      
        AMAXDIF=ABS(AM00-TSICVAL)/ABS(TSIC(MPTRM1)-TSIC(MPTRP1))                                                                    
        IF(AMAXDIF.GT.AMAXDF2) AMAXDF2=AMAXDIF                                                                                      
        AMAXDIF=ABS(TSIC(MPTR)-TSICVAL)/ABS(TSIC(MPTRM1)-TSIC(MPTRP1))                                                              
        IF(AMAXDIF.GT.AMAXDF3) AMAXDF3=AMAXDIF                                                                                      
C                                                                                                                                   
C       * TRANSFER COMPUTED SEA-ICE AMOUNT FROM TSIC TO SIC                                                                         
C                                                                                                                                   
 600    CONTINUE                                                                                                                    
C                                                                                                                                   
        DO 610 M=1,12                                                                                                               
          SIC(N,M)=TSIC(M)                                                                                                          
 610    CONTINUE                                                                                                                    
        NDAT(N)=0                                                                                                                   
C                                                                                                                                   
C       * INCREASE NDAT VALUE FOR NEIGHBOURING POINTS (NDAT.LT.-1) BY 1.                                                            
C                                                                                                                                   
        IF(J.GT.1)THEN                                                                                                              
          JB=(J-2)*ILON                                                                                                             
          IF(NDAT(JB+IM1).LT.-1) NDAT(JB+IM1)=NDAT(JB+IM1)+1                                                                        
          IF(NDAT(JB+ I ).LT.-1) NDAT(JB+ I )=NDAT(JB+ I )+1                                                                        
          IF(NDAT(JB+IP1).LT.-1) NDAT(JB+IP1)=NDAT(JB+IP1)+1                                                                        
        ENDIF                                                                                                                       
C                                                                                                                                   
        JB=(J-1)*ILON                                                                                                               
        IF(NDAT(JB+IM1).LT.-1) NDAT(JB+IM1)=NDAT(JB+IM1)+1                                                                          
        IF(NDAT(JB+IP1).LT.-1) NDAT(JB+IP1)=NDAT(JB+IP1)+1                                                                          
C                                                                                                                                   
        IF(J.LT.JLAT)THEN                                                                                                           
          JB=J*ILON                                                                                                                 
          IF(NDAT(JB+IM1).LT.-1) NDAT(JB+IM1)=NDAT(JB+IM1)+1                                                                        
          IF(NDAT(JB+ I ).LT.-1) NDAT(JB+ I )=NDAT(JB+ I )+1                                                                        
          IF(NDAT(JB+IP1).LT.-1) NDAT(JB+IP1)=NDAT(JB+IP1)+1                                                                        
        ENDIF                                                                                                                       
C                                                                                                                                   
        NPTSC=NPTSC+1                                                                                                               
        IF((NTTIME.GT.1).AND.(NTIME.GT.1)) GO TO 26                                                                                 
C                                                                                                                                   
 710   CONTINUE                                                                                                                     
C                                                                                                                                   
C      * GO BACK AND COMPUTE THE REMAINING POINT(S) ON THE SAME LAT IF ANY.                                                         
C                                                                                                                                   
       IF((NTTIME.GT.1).AND.(NTIME.LT.ILON))THEN                                                                                    
        JB=(J-1)*ILON                                                                                                               
        DO 715 I=1,ILON                                                                                                             
          IF(NDAT(JB+I).LT.0) GO TO 25                                                                                              
 715    CONTINUE                                                                                                                    
        WRITE(6,6035)J,NTIME                                                                                                        
       ENDIF                                                                                                                        
 720   CONTINUE                                                                                                                     
 730  CONTINUE                                                                                                                      
C                                                                                                                                   
C     IF(NSICNC.GT.0) PRINT *,' ALL YEAR SEA-ICE COVERED PTS AND NONE',                                                             
C    1                        ' IS COMPUTED, NSICNC= ',NSICNC                                                                       
C     IF(NSICPC.GT.0) PRINT *,' UNCOMPUTED PARTLY SIC (LAND), ',                                                                    
C    1                        ' NSICPC= ',NSICPC                                                                                    
C     IF(NSICNW.GT.0) PRINT *,' DEFERRED POINTS WITH O.W. NGBS',                                                                    
C    1                        ', NSICNW= ',NSICNW                                                                                   
C                                                                                                                                   
      WRITE(6,6045)NTTIME,NPTSC                                                                                                     
C                                                                                                                                   
C     * GO BACK AND COMPUTE THE REMAINING POINT(S) IF ANY.                                                                          
C                                                                                                                                   
      DO 750 N=1,NWDS                                                                                                               
        IF(NDAT(N).LT.0) GO TO 20                                                                                                   
 750  CONTINUE                                                                                                                      
C************************************************************                                                                       
C     * DO EXTRA YEAR INTEGRATION FOR PERMANENTLY ICE-COVERED                                                                       
C     * POINTS TO ELIMINATE PROBLEMS WITH INITIAL "AM0" SETTING.                                                                    
C                                                                                                                                   
      PRINT *, '0IN EXTRA YEAR INTEGRATION FOR PERMANENT ICE:'                                                                      
C                                                                                                                                   
      DO 775 N=1,NWDS                                                                                                               
        IF(NPERM(N).NE.0)                          THEN                                                                             
          AM0=AM0F(N)                                                                                                               
          MPTR=M0MINM(N)                                                                                                            
          NITER(N)=1                                                                                                                
  769     CONTINUE                                                                                                                  
          AM0O=AM0                                                                                                                  
          NITER(N)=NITER(N)+1                                                                                                       
          IF(N.GT.NWDS/2) THEN                                                                                                      
            RES=10.E0                                                                                                               
          ELSE                                                                                                                      
            RES=14.E0                                                                                                               
          ENDIF                                                                                                                     
          DO 770 M=MPTR,MPTR+10                                                                                                     
            MP1=M+1                                                                                                                 
            IF(MP1.GT.12)MP1=MP1-12                                                                                                 
            TTS(MP1)=TS(N,MP1)                                                                                                      
            RESF=RES*86400.E0*DIM(MP1)/(2.E0*HF)                                                                                    
C                                                                                                                                   
            IF(SNO(N,MP1).GT.SNOMAX) THEN                                                                                           
              DENS=0.54E0*SNO(N,MP1)/LOG(1.E0+0.54E0*SNO(N,MP1)/275.E0)                                                             
            ELSE                                                                                                                    
              DENS=275.E0                                                                                                           
            ENDIF                                                                                                                   
            CONS=2.805E-6*DENS**2                                                                                                   
            X=(DENI*CONI)*SNO(N,MP1)/(DENS*CONS)                                                                                    
            IF(AT0.GE.TTS(MP1))THEN                                                                                                 
              DT=AT0-TTS(MP1)                                                                                                       
              DTFAC=AR*DIM(MP1)*DT                                                                                                  
              AM1=SQRT((X+AM0-RESF)**2 + DTFAC) - (X+RESF)                                                                          
            ELSE                                                                                                                    
              DT=MIN(TTS(MP1)-AT0,2.E0)                                                                                             
              XMELT=SQRT(AR*DIM(MP1)*DT)                                                                                            
              IF(XMELT.GT.SNO(N,MP1)) THEN                                                                                          
                AM1=AM0+SNO(N,MP1)-XMELT                                                                                            
              ELSE                                                                                                                  
                AM1=AM0                                                                                                             
              ENDIF                                                                                                                 
            ENDIF                                                                                                                   
            TSIC(MP1)=0.5E0*(AM0+AM1)                                                                                               
            AM0=AM1                                                                                                                 
            IF(TSIC(MP1).LT.45.0E0) TSIC(MP1)=45.0E0                                                                                
            IF(AM0.LT.45.0E0) AM0=45.0E0                                                                                            
  770     CONTINUE                                                                                                                  
          TSIC(MPTR)=0.5E0*(AM0+AM0O)                                                                                               
C                                                                                                                                   
          IF(AM0-AM0O.LE.50.E0 .OR. NITER(N).GT.50) THEN                                                                            
C                                                                                                                                   
C           * PRINT OUT DEBUG INFO AT END.                                                                                          
C                                                                                                                                   
            PRINT *, '0N,NITER,MPTR,AM0,OLDAM0: ',N,NITER(N),MPTR,AM0,                                                              
     1                 AM0O                                                                                                         
            PRINT *, ' TSIC: ',(TSIC(M),M=1,12)                                                                                     
            PRINT *, ' TTS:  ',(TTS(M),M=1,12)                                                                                      
C           PRINT *, ' SNO:  ',(SNO(N,M),M=1,12)                                                                                    
                                                                                                                                    
          ENDIF                                                                                                                     
          IF(AM0-AM0O.GT.50.E0)     THEN                                                                                            
            PRINT *, '0N,NITER,AM0O,AM0,DIFF = ',N,NITER(N),                                                                        
     1                  AM0O,AM0,AM0-AM0O                                                                                           
            IF(NITER(N).LE.50)          GO TO 769                                                                                   
          ENDIF                                                                                                                     
          IF(NITER(N).GT.50)   CALL                XIT('LLOFSIC',-16)                                                               
C                                                                                                                                   
          DO 772 M=1,12                                                                                                             
            SIC(N,M)=TSIC(M)                                                                                                        
  772     CONTINUE                                                                                                                  
        ENDIF                                                                                                                       
  775 CONTINUE                                                                                                                      
C************************************************************                                                                       
C                                                                                                                                   
C     * WRITE-OUT THE COMPUTED SEA-ICE AMOUNT TO FILE SIC                                                                           
C                                                                                                                                   
      CALL SETLAB(IBUF,NC4TO8("GRID"),0,NC4TO8(" SIC"),1,ILON,JLAT,0,1)                                                             
      DO 800 M=1,12                                                                                                                 
        IBUF(2)=ITIM(M)                                                                                                             
        CALL PUTFLD2(4,SIC(1,M),IBUF,MAXX)                                                                                          
        WRITE(6,6025)IBUF                                                                                                           
 800  CONTINUE                                                                                                                      
      WRITE(6,6050)NTTIME,AMAXDF0,AMAXDF1,AMAXDF2,AMAXDF3                                                                           
      CALL                                         XIT('LLOFSIC',0)                                                                 
C--------------------------------------------------------------------                                                               
 6025 FORMAT(' ',A4,I10,1X,A4,I10,4I6)                                                                                              
 6035 FORMAT('  J= ',I7,', NTIME= ',I7)                                                                                             
 6040 FORMAT('  NWDS= ',I10)                                                                                                        
 6045 FORMAT('  NTTIME= ',I5,', NPTSC= ',I7)                                                                                        
 6050 FORMAT('  NTTIME = ',I5,/,                                                                                                    
     1       '  AMAXDF0= ',E17.6,/,                                                                                                 
     2       '  AMAXDF1= ',E17.6,/,                                                                                                 
     3       '  AMAXDF2= ',E17.6,/,                                                                                                 
     4       '  AMAXDF3= ',E17.6)                                                                                                   
      END                                                                                                                           
