pkg/seaice/seaice_init.F

C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init.F,v 1.18 2004/03/02 08:45:46 dimitri Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"

CStartOfInterface
      SUBROUTINE SEAICE_INIT( myThid )
C     /==========================================================\
C     | SUBROUTINE SEAICE_INIT                                   |
C     | o Initialization of sea ice model.                       |
C     |==========================================================|
C     \==========================================================/
      IMPLICIT NONE
 
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SEAICE.h"
#include "SEAICE_GRID.h"
#include "SEAICE_DIAGS.h"
#include "SEAICE_PARAMS.h"
#ifdef ALLOW_EXCH2
#include "W2_EXCH2_TOPOLOGY.h"
#include "W2_EXCH2_PARAMS.h"
#endif /* ALLOW_EXCH2 */

C     === Routine arguments ===
C     myThid - Thread no. that called this routine.
      INTEGER myThid
CEndOfInterface
 
#ifdef ALLOW_SEAICE
C     === Local variables ===
C     i,j,k,bi,bj - Loop counters

      INTEGER i, j, k, bi, bj
      _RS  mask_uice
      INTEGER myIter, myTile

#ifdef ALLOW_TIMEAVE
C     Initialize averages to zero
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        CALL TIMEAVE_RESET(FUtave   ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(FVtave   ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(EmPmRtave,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(QNETtave ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(QSWtave  ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(UICEtave ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(VICEtave ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(HEFFtave ,1,bi,bj,myThid)
        CALL TIMEAVE_RESET(AREAtave ,1,bi,bj,myThid)
        DO k=1,Nr
         SEAICE_TimeAve(k,bi,bj)=ZERO
        ENDDO
       ENDDO
      ENDDO
#endif /* ALLOW_TIMEAVE */

cph(
cph   make sure TAF sees proper initialisation
cph   to avoid partial recomputation issues
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
c
        DO K=1,3
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           HEFF(I,J,k,bi,bj)=SEAICE_initialHEFF
           AREA(I,J,k,bi,bj)=ZERO
           UICE(I,J,k,bi,bj)=ZERO
           VICE(I,J,k,bi,bj)=ZERO
          ENDDO
         ENDDO
        ENDDO
c
        DO J=1-OLy,sNy+OLy
         DO I=1-OLx,sNx+OLx
          HSNOW(I,J,bi,bj)=0.2 _d 0
          ZETA(I,J,bi,bj)=ZERO
         ENDDO
        ENDDO
c
       ENDDO
      ENDDO
cph)

C--   Initialize grid info
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO J=1-OLy,sNy+OLy
         DO I=1-OLx,sNx+OLx
          CSTICE(i,j,bi,bj) =cos(yC(I,J,bi,bj)*deg2rad)
          CSUICE(i,j,bi,bj) =cos(yG(I,J,bi,bj)*deg2rad)
          SINEICE(i,j,bi,bj)=sin(yC(I,J,bi,bj)*deg2rad)
          TNGTICE(i,j,bi,bj)=SINEICE(i,j,bi,bj)/CSTICE(i,j,bi,bj)
          SINEICE(i,j,bi,bj)=sin(yG(I,J,bi,bj)*deg2rad)
          TNGICE(i,j,bi,bj) =SINEICE(i,j,bi,bj)/CSUICE(i,j,bi,bj)
          DXTICE(i,j,bi,bj)=dxF(i,j,bi,bj)/CSTICE(i,j,bi,bj)
          DXUICE(i,j,bi,bj)=dxV(i,j,bi,bj)/CSUICE(i,j,bi,bj)
          DYTICE(i,j,bi,bj)=dyF(i,j,bi,bj)
          DYUICE(i,j,bi,bj)=dyU(i,j,bi,bj)
         ENDDO
        ENDDO
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          HEFFM(i,j,bi,bj)=ONE
          IF (_hFacC(i,j,1,bi,bj).eq.0.) HEFFM(i,j,bi,bj)=ZERO
         ENDDO
        ENDDO
        DO J=2-OLy,sNy+OLy
         DO I=2-OLx,sNx+OLx
          UVM(i,j,bi,bj)=ZERO
          mask_uice=HEFFM(I,J,  bi,bj)+HEFFM(I-1,J-1,bi,bj)
     &             +HEFFM(I,J-1,bi,bj)+HEFFM(I-1,J,  bi,bj)
          IF(mask_uice.GT.3.5) UVM(I,J,bi,bj)=ONE
         ENDDO
        ENDDO

#ifdef ALLOW_EXCH2
C--   Special stuff for cubed sphere: assume grid is rectangular and
C     set UV mask to zero except for Arctic and Antarctic cube faces.
        IF (useCubedSphereExchange) THEN
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           CSTICE(i,j,bi,bj) = ONE
           CSUICE(i,j,bi,bj) = ONE
           TNGTICE(i,j,bi,bj)= ZERO
           TNGICE(i,j,bi,bj) = ZERO
           DXTICE(i,j,bi,bj) = dxF(i,j,bi,bj)
           DXUICE(i,j,bi,bj) = dxV(i,j,bi,bj)
          ENDDO
         ENDDO
         myTile = W2_myTileList(bi)
         IF ( exch2_myFace(myTile) .EQ. 1 .OR.
     &        exch2_myFace(myTile) .EQ. 2 .OR.
     &        exch2_myFace(myTile) .EQ. 4 .OR.
     &        exch2_myFace(myTile) .EQ. 5 ) THEN
          DO J=1-OLy,sNy+OLy
           DO I=1-OLx,sNx+OLx
            UVM(i,j,bi,bj)=ZERO
           ENDDO
          ENDDO
         ELSEIF ( exch2_isWedge(myTile) .EQ. 1 ) THEN
          I=1
          DO J=1-OLy,sNy+OLy
           UVM(i,j,bi,bj)=ZERO
          ENDDO
         ELSEIF ( exch2_isSedge(myTile) .EQ. 1 ) THEN
          J=1
          DO I=1-OLx,sNx+OLx
           UVM(i,j,bi,bj)=ZERO
          ENDDO
         ENDIF
        ENDIF
#endif /* ALLOW_EXCH2 */

        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          TICE(I,J,bi,bj)=273.0 _d 0
#ifdef SEAICE_MULTILEVEL
          DO k=1,MULTDIM
           TICES(I,J,k,bi,bj)=273.0 _d 0
          ENDDO
#endif
          UICEC(I,J,bi,bj)=ZERO
          VICEC(I,J,bi,bj)=ZERO
          AMASS(I,J,bi,bj)=1000.0 _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO

C--   Update overlap regions
      _EXCH_XY_R8(UVM, myThid)

      myIter=0
      CALL PLOT_FIELD_XYRL( CSTICE  , 'Current CSTICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( CSUICE  , 'Current CSUICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( TNGTICE  , 'Current TNGTICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( TNGICE  , 'Current TNGICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( SINEICE  , 'Current SINEICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( DXTICE  , 'Current DXTICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( DXUICE  , 'Current DXUICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( DYTICE  , 'Current DYTICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( DYUICE  , 'Current DYUICE  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( HEFFM  , 'Current HEFFM  ' ,
     &                        myIter, myThid )
      CALL PLOT_FIELD_XYRL( UVM  , 'Current UVM  ' ,
     &                        myIter, myThid )

C--   Set model variables to initial/restart conditions
      IF ( nIter0 .NE. 0 ) THEN

         CALL SEAICE_READ_PICKUP ( myThid )

      ELSE

         DO bj=myByLo(myThid),myByHi(myThid)
            DO bi=myBxLo(myThid),myBxHi(myThid)
               DO j=1-OLy,sNy+OLy
                  DO i=1-OLx,sNx+OLx
                     HSNOW(I,J,bi,bj)=0.2 _d 0
                     YNEG(I,J,bi,bj)=ZERO
                     TMIX(I,J,bi,bj)=TICE(I,J,bi,bj)
                     DO k=1,3
                        HEFF(I,J,k,bi,bj)=SEAICE_initialHEFF
                        AREA(I,J,k,bi,bj)=HEFFM(i,j,bi,bj)
                        UICE(I,J,k,bi,bj)=ZERO
                        VICE(I,J,k,bi,bj)=ZERO
                     ENDDO
                  ENDDO
               ENDDO
            ENDDO
         ENDDO

C--   Read initial sea-ice thickness from file if available.
         IF ( HeffFile .NE. ' ' ) THEN
            _BEGIN_MASTER( myThid )
            CALL READ_FLD_XY_RL( HeffFile, ' ', ZETA, 0, myThid )
            _END_MASTER(myThid)
            _EXCH_XY_R8(ZETA,myThid)
            DO bj=myByLo(myThid),myByHi(myThid)
               DO bi=myBxLo(myThid),myBxHi(myThid)
                  DO j=1-OLy,sNy+OLy
                     DO i=1-OLx,sNx+OLx
                        DO k=1,3
                         HEFF(I,J,k,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
                         IF ( ZETA(i,j,bi,bj).EQ.ZERO )
     &                        AREA(I,J,k,bi,bj) = ZERO
                        ENDDO
                     ENDDO
                  ENDDO
               ENDDO
            ENDDO
         ENDIF
         
      ENDIF

C---  Complete initialization
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          ZETA(I,J,bi,bj)=HEFF(I,J,1,bi,bj)*(1.0 _d 11)
          ETA(I,J,bi,bj)=ZETA(I,J,bi,bj)/4.0 _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#endif /* ALLOW_SEAICE */

      RETURN
      END
1	adcroft	1.19	C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init.F,v 1.18 2004/03/02 08:45:46 dimitri Exp $
2	edhill	1.10	C $Name: $
3	heimbach	1.2
4			#include "SEAICE_OPTIONS.h"
5
6			CStartOfInterface
7			SUBROUTINE SEAICE_INIT( myThid )
8			C /==========================================================\
9			C \| SUBROUTINE SEAICE_INIT \|
10			C \| o Initialization of sea ice model. \|
11			C \|==========================================================\|
12			C \==========================================================/
13			IMPLICIT NONE
14
15			C === Global variables ===
16			#include "SIZE.h"
17			#include "EEPARAMS.h"
18			#include "PARAMS.h"
19			#include "GRID.h"
20			#include "SEAICE.h"
21			#include "SEAICE_GRID.h"
22			#include "SEAICE_DIAGS.h"
23			#include "SEAICE_PARAMS.h"
24	adcroft	1.19	#ifdef ALLOW_EXCH2
25	dimitri	1.15	#include "W2_EXCH2_TOPOLOGY.h"
26			#include "W2_EXCH2_PARAMS.h"
27	adcroft	1.19	#endif /* ALLOW_EXCH2 */
28	heimbach	1.2
29			C === Routine arguments ===
30			C myThid - Thread no. that called this routine.
31			INTEGER myThid
32			CEndOfInterface
33
34			#ifdef ALLOW_SEAICE
35			C === Local variables ===
36			C i,j,k,bi,bj - Loop counters
37
38			INTEGER i, j, k, bi, bj
39			_RS mask_uice
40	dimitri	1.15	INTEGER myIter, myTile
41	heimbach	1.2
42			#ifdef ALLOW_TIMEAVE
43			C Initialize averages to zero
44			DO bj = myByLo(myThid), myByHi(myThid)
45			DO bi = myBxLo(myThid), myBxHi(myThid)
46			CALL TIMEAVE_RESET(FUtave ,1,bi,bj,myThid)
47			CALL TIMEAVE_RESET(FVtave ,1,bi,bj,myThid)
48			CALL TIMEAVE_RESET(EmPmRtave,1,bi,bj,myThid)
49			CALL TIMEAVE_RESET(QNETtave ,1,bi,bj,myThid)
50			CALL TIMEAVE_RESET(QSWtave ,1,bi,bj,myThid)
51			CALL TIMEAVE_RESET(UICEtave ,1,bi,bj,myThid)
52			CALL TIMEAVE_RESET(VICEtave ,1,bi,bj,myThid)
53			CALL TIMEAVE_RESET(HEFFtave ,1,bi,bj,myThid)
54			CALL TIMEAVE_RESET(AREAtave ,1,bi,bj,myThid)
55			DO k=1,Nr
56	dimitri	1.3	SEAICE_TimeAve(k,bi,bj)=ZERO
57	heimbach	1.2	ENDDO
58			ENDDO
59			ENDDO
60			#endif /* ALLOW_TIMEAVE */
61
62	heimbach	1.9	cph(
63			cph make sure TAF sees proper initialisation
64			cph to avoid partial recomputation issues
65			DO bj=myByLo(myThid),myByHi(myThid)
66			DO bi=myBxLo(myThid),myBxHi(myThid)
67			c
68			DO K=1,3
69	dimitri	1.17	DO J=1-OLy,sNy+OLy
70			DO I=1-OLx,sNx+OLx
71	heimbach	1.9	HEFF(I,J,k,bi,bj)=SEAICE_initialHEFF
72			AREA(I,J,k,bi,bj)=ZERO
73			UICE(I,J,k,bi,bj)=ZERO
74			VICE(I,J,k,bi,bj)=ZERO
75			ENDDO
76			ENDDO
77			ENDDO
78			c
79	dimitri	1.17	DO J=1-OLy,sNy+OLy
80			DO I=1-OLx,sNx+OLx
81	heimbach	1.9	HSNOW(I,J,bi,bj)=0.2 _d 0
82			ZETA(I,J,bi,bj)=ZERO
83			ENDDO
84			ENDDO
85			c
86			ENDDO
87			ENDDO
88			cph)
89
90	dimitri	1.15	C-- Initialize grid info
91	heimbach	1.2	DO bj=myByLo(myThid),myByHi(myThid)
92			DO bi=myBxLo(myThid),myBxHi(myThid)
93	dimitri	1.17	DO J=1-OLy,sNy+OLy
94			DO I=1-OLx,sNx+OLx
95	dimitri	1.7	CSTICE(i,j,bi,bj) =cos(yC(I,J,bi,bj)*deg2rad)
96			CSUICE(i,j,bi,bj) =cos(yG(I,J,bi,bj)*deg2rad)
97	dimitri	1.8	SINEICE(i,j,bi,bj)=sin(yC(I,J,bi,bj)*deg2rad)
98			TNGTICE(i,j,bi,bj)=SINEICE(i,j,bi,bj)/CSTICE(i,j,bi,bj)
99	dimitri	1.7	SINEICE(i,j,bi,bj)=sin(yG(I,J,bi,bj)*deg2rad)
100			TNGICE(i,j,bi,bj) =SINEICE(i,j,bi,bj)/CSUICE(i,j,bi,bj)
101	heimbach	1.2	DXTICE(i,j,bi,bj)=dxF(i,j,bi,bj)/CSTICE(i,j,bi,bj)
102	dimitri	1.7	DXUICE(i,j,bi,bj)=dxV(i,j,bi,bj)/CSUICE(i,j,bi,bj)
103	heimbach	1.2	DYTICE(i,j,bi,bj)=dyF(i,j,bi,bj)
104	dimitri	1.7	DYUICE(i,j,bi,bj)=dyU(i,j,bi,bj)
105	heimbach	1.2	ENDDO
106			ENDDO
107			DO j=1-OLy,sNy+OLy
108			DO i=1-OLx,sNx+OLx
109	dimitri	1.3	HEFFM(i,j,bi,bj)=ONE
110			IF (_hFacC(i,j,1,bi,bj).eq.0.) HEFFM(i,j,bi,bj)=ZERO
111	heimbach	1.2	ENDDO
112			ENDDO
113	dimitri	1.17	DO J=2-OLy,sNy+OLy
114			DO I=2-OLx,sNx+OLx
115	dimitri	1.3	UVM(i,j,bi,bj)=ZERO
116	dimitri	1.7	mask_uice=HEFFM(I,J, bi,bj)+HEFFM(I-1,J-1,bi,bj)
117			& +HEFFM(I,J-1,bi,bj)+HEFFM(I-1,J, bi,bj)
118	dimitri	1.3	IF(mask_uice.GT.3.5) UVM(I,J,bi,bj)=ONE
119	heimbach	1.2	ENDDO
120			ENDDO
121	dimitri	1.15
122	adcroft	1.19	#ifdef ALLOW_EXCH2
123	dimitri	1.15	C-- Special stuff for cubed sphere: assume grid is rectangular and
124			C set UV mask to zero except for Arctic and Antarctic cube faces.
125			IF (useCubedSphereExchange) THEN
126	dimitri	1.17	DO J=1-OLy,sNy+OLy
127			DO I=1-OLx,sNx+OLx
128	dimitri	1.15	CSTICE(i,j,bi,bj) = ONE
129			CSUICE(i,j,bi,bj) = ONE
130			TNGTICE(i,j,bi,bj)= ZERO
131			TNGICE(i,j,bi,bj) = ZERO
132			DXTICE(i,j,bi,bj) = dxF(i,j,bi,bj)
133			DXUICE(i,j,bi,bj) = dxV(i,j,bi,bj)
134			ENDDO
135			ENDDO
136			myTile = W2_myTileList(bi)
137			IF ( exch2_myFace(myTile) .EQ. 1 .OR.
138			& exch2_myFace(myTile) .EQ. 2 .OR.
139			& exch2_myFace(myTile) .EQ. 4 .OR.
140			& exch2_myFace(myTile) .EQ. 5 ) THEN
141	dimitri	1.17	DO J=1-OLy,sNy+OLy
142			DO I=1-OLx,sNx+OLx
143	dimitri	1.15	UVM(i,j,bi,bj)=ZERO
144			ENDDO
145			ENDDO
146			ELSEIF ( exch2_isWedge(myTile) .EQ. 1 ) THEN
147			I=1
148	dimitri	1.17	DO J=1-OLy,sNy+OLy
149	dimitri	1.15	UVM(i,j,bi,bj)=ZERO
150			ENDDO
151			ELSEIF ( exch2_isSedge(myTile) .EQ. 1 ) THEN
152			J=1
153	dimitri	1.17	DO I=1-OLx,sNx+OLx
154	dimitri	1.15	UVM(i,j,bi,bj)=ZERO
155			ENDDO
156			ENDIF
157			ENDIF
158	adcroft	1.19	#endif /* ALLOW_EXCH2 */
159	heimbach	1.9
160	heimbach	1.2	DO j=1-OLy,sNy+OLy
161			DO i=1-OLx,sNx+OLx
162			TICE(I,J,bi,bj)=273.0 _d 0
163	dimitri	1.8	#ifdef SEAICE_MULTILEVEL
164	heimbach	1.13	DO k=1,MULTDIM
165	dimitri	1.8	TICES(I,J,k,bi,bj)=273.0 _d 0
166			ENDDO
167			#endif
168	dimitri	1.3	UICEC(I,J,bi,bj)=ZERO
169			VICEC(I,J,bi,bj)=ZERO
170			AMASS(I,J,bi,bj)=1000.0 _d 0
171	heimbach	1.2	ENDDO
172			ENDDO
173			ENDDO
174			ENDDO
175
176			C-- Update overlap regions
177			_EXCH_XY_R8(UVM, myThid)
178	dimitri	1.8
179			myIter=0
180			CALL PLOT_FIELD_XYRL( CSTICE , 'Current CSTICE ' ,
181			& myIter, myThid )
182			CALL PLOT_FIELD_XYRL( CSUICE , 'Current CSUICE ' ,
183			& myIter, myThid )
184			CALL PLOT_FIELD_XYRL( TNGTICE , 'Current TNGTICE ' ,
185			& myIter, myThid )
186			CALL PLOT_FIELD_XYRL( TNGICE , 'Current TNGICE ' ,
187			& myIter, myThid )
188			CALL PLOT_FIELD_XYRL( SINEICE , 'Current SINEICE ' ,
189			& myIter, myThid )
190			CALL PLOT_FIELD_XYRL( DXTICE , 'Current DXTICE ' ,
191			& myIter, myThid )
192			CALL PLOT_FIELD_XYRL( DXUICE , 'Current DXUICE ' ,
193			& myIter, myThid )
194			CALL PLOT_FIELD_XYRL( DYTICE , 'Current DYTICE ' ,
195			& myIter, myThid )
196			CALL PLOT_FIELD_XYRL( DYUICE , 'Current DYUICE ' ,
197			& myIter, myThid )
198			CALL PLOT_FIELD_XYRL( HEFFM , 'Current HEFFM ' ,
199			& myIter, myThid )
200			CALL PLOT_FIELD_XYRL( UVM , 'Current UVM ' ,
201			& myIter, myThid )
202	heimbach	1.2
203			C-- Set model variables to initial/restart conditions
204	dimitri	1.5	IF ( nIter0 .NE. 0 ) THEN
205	dimitri	1.18
206	dimitri	1.5	CALL SEAICE_READ_PICKUP ( myThid )
207	dimitri	1.18
208	dimitri	1.5	ELSE
209	dimitri	1.18
210	heimbach	1.2	DO bj=myByLo(myThid),myByHi(myThid)
211			DO bi=myBxLo(myThid),myBxHi(myThid)
212			DO j=1-OLy,sNy+OLy
213			DO i=1-OLx,sNx+OLx
214			HSNOW(I,J,bi,bj)=0.2 _d 0
215	dimitri	1.3	YNEG(I,J,bi,bj)=ZERO
216	heimbach	1.2	TMIX(I,J,bi,bj)=TICE(I,J,bi,bj)
217			DO k=1,3
218	dimitri	1.6	HEFF(I,J,k,bi,bj)=SEAICE_initialHEFF
219	heimbach	1.2	AREA(I,J,k,bi,bj)=HEFFM(i,j,bi,bj)
220	dimitri	1.3	UICE(I,J,k,bi,bj)=ZERO
221			VICE(I,J,k,bi,bj)=ZERO
222	heimbach	1.2	ENDDO
223			ENDDO
224			ENDDO
225			ENDDO
226			ENDDO
227	dimitri	1.5
228			C-- Read initial sea-ice thickness from file if available.
229	dimitri	1.18	IF ( HeffFile .NE. ' ' ) THEN
230			_BEGIN_MASTER( myThid )
231			CALL READ_FLD_XY_RL( HeffFile, ' ', ZETA, 0, myThid )
232			_END_MASTER(myThid)
233			_EXCH_XY_R8(ZETA,myThid)
234			DO bj=myByLo(myThid),myByHi(myThid)
235			DO bi=myBxLo(myThid),myBxHi(myThid)
236			DO j=1-OLy,sNy+OLy
237			DO i=1-OLx,sNx+OLx
238			DO k=1,3
239			HEFF(I,J,k,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
240			IF ( ZETA(i,j,bi,bj).EQ.ZERO )
241			& AREA(I,J,k,bi,bj) = ZERO
242			ENDDO
243	dimitri	1.5	ENDDO
244			ENDDO
245			ENDDO
246			ENDDO
247	dimitri	1.18	ENDIF
248
249	heimbach	1.2	ENDIF
250
251			C--- Complete initialization
252			DO bj=myByLo(myThid),myByHi(myThid)
253			DO bi=myBxLo(myThid),myBxHi(myThid)
254			DO j=1-OLy,sNy+OLy
255			DO i=1-OLx,sNx+OLx
256			ZETA(I,J,bi,bj)=HEFF(I,J,1,bi,bj)*(1.0 _d 11)
257	dimitri	1.3	ETA(I,J,bi,bj)=ZETA(I,J,bi,bj)/4.0 _d 0
258	heimbach	1.2	ENDDO
259			ENDDO
260			ENDDO
261			ENDDO
262
263	dimitri	1.5	#endif /* ALLOW_SEAICE */
264	heimbach	1.2
265			RETURN
266			END