ifenty/seaiceAdjointCode/seaice_init_varia.F

C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.6 2007/04/27 15:50:14 jmc Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"

CStartOfInterface
      SUBROUTINE SEAICE_INIT_VARIA( myThid )
C     /==========================================================\
C     | SUBROUTINE SEAICE_INIT_VARIA                             |
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"
CML#include "SEAICE_GRID.h"
#include "SEAICE_DIAGS.h"
#include "SEAICE_PARAMS.h"
#include "FFIELDS.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
 
C     === Local variables ===
C     i,j,k,bi,bj - Loop counters

      INTEGER i, j, k, bi, bj
      _RS  mask_uice
      INTEGER myIter, myTile
cif(  Helper variable for determining the fraction of sw radiation
cif(  penetrating the model's shallowest layer
      _RL swfracba(two)
      _RL FACTORM,dummyTime
      INTEGER IMAX
      
      IMAX = 2
      FACTORM = -1.0
      dummyTime = 1.0

      swfracba(1) = abs(rF(1)) 
      swfracba(2) = abs(rF(2))
      CALL SWFRAC(
     I       IMAX,FACTORM,
     I       dummyTime,myThid,
     U       swfracba)

      SWFRACB = swfracba(2)

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)=ZERO
           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)=ZERO
          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
          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=1-OLy+1,sNy+OLy
         DO I=1-OLx+1,sNx+OLx
#ifndef SEAICE_CGRID
          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
#else
          seaiceMaskU(I,J,bi,bj)=   0.0 _d 0
          seaiceMaskV(I,J,bi,bj)=   0.0 _d 0
          mask_uice=HEFFM(I,J,bi,bj)+HEFFM(I-1,J  ,bi,bj)
          IF(mask_uice.GT.1.5) seaiceMaskU(I,J,bi,bj)=ONE
          mask_uice=HEFFM(I,J,bi,bj)+HEFFM(I  ,J-1,bi,bj)
          IF(mask_uice.GT.1.5) seaiceMaskV(I,J,bi,bj)=ONE
#endif /* not SEAICE_CGRID */
         ENDDO
        ENDDO

#ifdef ALLOW_EXCH2
#ifdef SEAICE_CGRID
#else
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
         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
#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_MULTICATEGORY
          DO k=1,MULTDIM
           TICES(I,J,k,bi,bj)=273.0 _d 0
          ENDDO
#endif /* SEAICE_MULTICATEGORY */
          UICEC      (I,J,bi,bj)=ZERO
          VICEC      (I,J,bi,bj)=ZERO
          DWATN      (I,J,bi,bj)=ZERO
#ifndef SEAICE_CGRID
          DAIRN      (I,J,bi,bj)=ZERO
          AMASS      (I,J,bi,bj)=1000.0 _d 0
#else
          seaiceMassC(I,J,bi,bj)=1000.0 _d 0
          seaiceMassU(I,J,bi,bj)=1000.0 _d 0
          seaiceMassV(I,J,bi,bj)=1000.0 _d 0
#endif
          GWATX      (I,J,bi,bj)=ZERO
          GWATY      (I,J,bi,bj)=ZERO
         ENDDO
        ENDDO

C--   Choose a proxy level for geostrophic velocity,
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
#ifdef SEAICE_TEST_ICE_STRESS_1
          KGEO(I,J,bi,bj) = 1
#else /* SEAICE_TEST_ICE_STRESS_1 */
          IF (klowc(i,j,bi,bj) .LT. 2) THEN
           KGEO(I,J,bi,bj) = 1
          ELSE
           KGEO(I,J,bi,bj) = 2
           DO WHILE ( abs(rC(KGEO(I,J,bi,bj))) .LT. 50.0 .AND.
     &          KGEO(I,J,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
              KGEO(I,J,bi,bj) = KGEO(I,J,bi,bj) + 1
           ENDDO
          ENDIF
#endif /* SEAICE_TEST_ICE_STRESS_1 */
         ENDDO
        ENDDO

       ENDDO
      ENDDO

C--   Update overlap regions
#ifdef SEAICE_CGRID
      CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid)
#else
      _EXCH_XY_R8(UVM, myThid)
#endif

C--   Now lets look at all these beasts
      IF ( debugLevel .GE. debLevB ) THEN
         myIter=0
         CALL PLOT_FIELD_XYRL( HEFFM   , 'Current HEFFM   ' ,
     &        myIter, myThid )
#ifdef SEAICE_CGRID
         CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU',
     &        myIter, myThid )
         CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV',
     &        myIter, myThid )
#else
         CALL PLOT_FIELD_XYRL( UVM     , 'Current UVM     ' ,
     &        myIter, myThid )
#endif
      ENDIF

#if (defined (SEAICE_CGRID) && defined (SEAICE_ALLOW_EVP))
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           stressDivergenceX(I,J,bi,bj) = 0. _d 0
           stressDivergenceY(I,J,bi,bj) = 0. _d 0
           seaice_sigma1 (I,J,bi,bj) = 0. _d 0
           seaice_sigma2 (I,J,bi,bj) = 0. _d 0
           seaice_sigma12(I,J,bi,bj) = 0. _d 0
          ENDDO
         ENDDO
        ENDDO
       ENDDO
#endif /* SEAICE_ALLOW_EVP and SEAICE_CGRID */

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*HEFFM(i,j,bi,bj)
           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*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)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDIF

       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
            IF(HEFF(I,J,k,bi,bj).GT.ZERO) THEN
             AREA(I,J,k,bi,bj)=ONE
            ELSE
             HSNOW(I,J,bi,bj)=ZERO
            ENDIF
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO
         
      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
        IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           sIceLoad(i,j,bi,bj) = HEFF(I,J,1,bi,bj)*SEAICE_rhoIce
     &                         + HSNOW(I,J,bi,bj)* 330. _d 0
          
          ENDDO
         ENDDO
        ENDIF
       ENDDO
      ENDDO

      RETURN
      END
1	gforget	1.1	C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_init_varia.F,v 1.6 2007/04/27 15:50:14 jmc Exp $
2			C $Name: $
3
4			#include "SEAICE_OPTIONS.h"
5
6			CStartOfInterface
7			SUBROUTINE SEAICE_INIT_VARIA( myThid )
8			C /==========================================================\
9			C \| SUBROUTINE SEAICE_INIT_VARIA \|
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			CML#include "SEAICE_GRID.h"
22			#include "SEAICE_DIAGS.h"
23			#include "SEAICE_PARAMS.h"
24			#include "FFIELDS.h"
25			#ifdef ALLOW_EXCH2
26			#include "W2_EXCH2_TOPOLOGY.h"
27			#include "W2_EXCH2_PARAMS.h"
28			#endif /* ALLOW_EXCH2 */
29
30			C === Routine arguments ===
31			C myThid - Thread no. that called this routine.
32			INTEGER myThid
33			CEndOfInterface
34
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			INTEGER myIter, myTile
41			cif( Helper variable for determining the fraction of sw radiation
42			cif( penetrating the model's shallowest layer
43			_RL swfracba(two)
44			_RL FACTORM,dummyTime
45			INTEGER IMAX
46
47			IMAX = 2
48			FACTORM = -1.0
49			dummyTime = 1.0
50
51			swfracba(1) = abs(rF(1))
52			swfracba(2) = abs(rF(2))
53			CALL SWFRAC(
54			I IMAX,FACTORM,
55			I dummyTime,myThid,
56			U swfracba)
57
58			SWFRACB = swfracba(2)
59
60			cph(
61			cph make sure TAF sees proper initialisation
62			cph to avoid partial recomputation issues
63			DO bj=myByLo(myThid),myByHi(myThid)
64			DO bi=myBxLo(myThid),myBxHi(myThid)
65			c
66			DO K=1,3
67			DO J=1-OLy,sNy+OLy
68			DO I=1-OLx,sNx+OLx
69			HEFF(I,J,k,bi,bj)=ZERO
70			AREA(I,J,k,bi,bj)=ZERO
71			UICE(I,J,k,bi,bj)=ZERO
72			VICE(I,J,k,bi,bj)=ZERO
73			ENDDO
74			ENDDO
75			ENDDO
76			c
77			DO J=1-OLy,sNy+OLy
78			DO I=1-OLx,sNx+OLx
79			HSNOW(I,J,bi,bj)=ZERO
80			ZETA(I,J,bi,bj)=ZERO
81			ENDDO
82			ENDDO
83			c
84			ENDDO
85			ENDDO
86			cph)
87
88			C-- Initialize grid info
89			DO bj=myByLo(myThid),myByHi(myThid)
90			DO bi=myBxLo(myThid),myBxHi(myThid)
91			DO j=1-OLy,sNy+OLy
92			DO i=1-OLx,sNx+OLx
93			HEFFM(i,j,bi,bj)=ONE
94			IF (_hFacC(i,j,1,bi,bj).eq.0.) HEFFM(i,j,bi,bj)=ZERO
95			ENDDO
96			ENDDO
97			DO J=1-OLy+1,sNy+OLy
98			DO I=1-OLx+1,sNx+OLx
99			#ifndef SEAICE_CGRID
100			UVM(i,j,bi,bj)=ZERO
101			mask_uice=HEFFM(I,J, bi,bj)+HEFFM(I-1,J-1,bi,bj)
102			& +HEFFM(I,J-1,bi,bj)+HEFFM(I-1,J, bi,bj)
103			IF(mask_uice.GT.3.5) UVM(I,J,bi,bj)=ONE
104			#else
105			seaiceMaskU(I,J,bi,bj)= 0.0 _d 0
106			seaiceMaskV(I,J,bi,bj)= 0.0 _d 0
107			mask_uice=HEFFM(I,J,bi,bj)+HEFFM(I-1,J ,bi,bj)
108			IF(mask_uice.GT.1.5) seaiceMaskU(I,J,bi,bj)=ONE
109			mask_uice=HEFFM(I,J,bi,bj)+HEFFM(I ,J-1,bi,bj)
110			IF(mask_uice.GT.1.5) seaiceMaskV(I,J,bi,bj)=ONE
111			#endif /* not SEAICE_CGRID */
112			ENDDO
113			ENDDO
114
115			#ifdef ALLOW_EXCH2
116			#ifdef SEAICE_CGRID
117			#else
118			C-- Special stuff for cubed sphere: assume grid is rectangular and
119			C set UV mask to zero except for Arctic and Antarctic cube faces.
120			IF (useCubedSphereExchange) THEN
121			myTile = W2_myTileList(bi)
122			IF ( exch2_myFace(myTile) .EQ. 1 .OR.
123			& exch2_myFace(myTile) .EQ. 2 .OR.
124			& exch2_myFace(myTile) .EQ. 4 .OR.
125			& exch2_myFace(myTile) .EQ. 5 ) THEN
126			DO J=1-OLy,sNy+OLy
127			DO I=1-OLx,sNx+OLx
128			UVM(i,j,bi,bj)=ZERO
129			ENDDO
130			ENDDO
131			ELSEIF ( exch2_isWedge(myTile) .EQ. 1 ) THEN
132			I=1
133			DO J=1-OLy,sNy+OLy
134			UVM(i,j,bi,bj)=ZERO
135			ENDDO
136			ELSEIF ( exch2_isSedge(myTile) .EQ. 1 ) THEN
137			J=1
138			DO I=1-OLx,sNx+OLx
139			UVM(i,j,bi,bj)=ZERO
140			ENDDO
141			ENDIF
142			ENDIF
143			#endif
144			#endif /* ALLOW_EXCH2 */
145
146			DO j=1-OLy,sNy+OLy
147			DO i=1-OLx,sNx+OLx
148			TICE(I,J,bi,bj)=273.0 _d 0
149			#ifdef SEAICE_MULTICATEGORY
150			DO k=1,MULTDIM
151			TICES(I,J,k,bi,bj)=273.0 _d 0
152			ENDDO
153			#endif /* SEAICE_MULTICATEGORY */
154			UICEC (I,J,bi,bj)=ZERO
155			VICEC (I,J,bi,bj)=ZERO
156			DWATN (I,J,bi,bj)=ZERO
157			#ifndef SEAICE_CGRID
158			DAIRN (I,J,bi,bj)=ZERO
159			AMASS (I,J,bi,bj)=1000.0 _d 0
160			#else
161			seaiceMassC(I,J,bi,bj)=1000.0 _d 0
162			seaiceMassU(I,J,bi,bj)=1000.0 _d 0
163			seaiceMassV(I,J,bi,bj)=1000.0 _d 0
164			#endif
165			GWATX (I,J,bi,bj)=ZERO
166			GWATY (I,J,bi,bj)=ZERO
167			ENDDO
168			ENDDO
169
170			C-- Choose a proxy level for geostrophic velocity,
171			DO j=1-OLy,sNy+OLy
172			DO i=1-OLx,sNx+OLx
173			#ifdef SEAICE_TEST_ICE_STRESS_1
174			KGEO(I,J,bi,bj) = 1
175			#else /* SEAICE_TEST_ICE_STRESS_1 */
176			IF (klowc(i,j,bi,bj) .LT. 2) THEN
177			KGEO(I,J,bi,bj) = 1
178			ELSE
179			KGEO(I,J,bi,bj) = 2
180			DO WHILE ( abs(rC(KGEO(I,J,bi,bj))) .LT. 50.0 .AND.
181			& KGEO(I,J,bi,bj) .LT. (klowc(i,j,bi,bj)-1) )
182			KGEO(I,J,bi,bj) = KGEO(I,J,bi,bj) + 1
183			ENDDO
184			ENDIF
185			#endif /* SEAICE_TEST_ICE_STRESS_1 */
186			ENDDO
187			ENDDO
188
189			ENDDO
190			ENDDO
191
192			C-- Update overlap regions
193			#ifdef SEAICE_CGRID
194			CALL EXCH_UV_XY_RL(seaiceMaskU,seaiceMaskV,.FALSE.,myThid)
195			#else
196			_EXCH_XY_R8(UVM, myThid)
197			#endif
198
199			C-- Now lets look at all these beasts
200			IF ( debugLevel .GE. debLevB ) THEN
201			myIter=0
202			CALL PLOT_FIELD_XYRL( HEFFM , 'Current HEFFM ' ,
203			& myIter, myThid )
204			#ifdef SEAICE_CGRID
205			CALL PLOT_FIELD_XYRL( seaiceMaskU, 'Current seaiceMaskU',
206			& myIter, myThid )
207			CALL PLOT_FIELD_XYRL( seaiceMaskV, 'Current seaiceMaskV',
208			& myIter, myThid )
209			#else
210			CALL PLOT_FIELD_XYRL( UVM , 'Current UVM ' ,
211			& myIter, myThid )
212			#endif
213			ENDIF
214
215			#if (defined (SEAICE_CGRID) && defined (SEAICE_ALLOW_EVP))
216			DO bj=myByLo(myThid),myByHi(myThid)
217			DO bi=myBxLo(myThid),myBxHi(myThid)
218			DO j=1-OLy,sNy+OLy
219			DO i=1-OLx,sNx+OLx
220			stressDivergenceX(I,J,bi,bj) = 0. _d 0
221			stressDivergenceY(I,J,bi,bj) = 0. _d 0
222			seaice_sigma1 (I,J,bi,bj) = 0. _d 0
223			seaice_sigma2 (I,J,bi,bj) = 0. _d 0
224			seaice_sigma12(I,J,bi,bj) = 0. _d 0
225			ENDDO
226			ENDDO
227			ENDDO
228			ENDDO
229			#endif /* SEAICE_ALLOW_EVP and SEAICE_CGRID */
230
231			C-- Set model variables to initial/restart conditions
232			IF ( nIter0 .NE. 0 ) THEN
233
234			CALL SEAICE_READ_PICKUP ( myThid )
235
236			ELSE
237
238			DO bj=myByLo(myThid),myByHi(myThid)
239			DO bi=myBxLo(myThid),myBxHi(myThid)
240			DO j=1-OLy,sNy+OLy
241			DO i=1-OLx,sNx+OLx
242			HSNOW(I,J,bi,bj)=0.2*HEFFM(i,j,bi,bj)
243			YNEG(I,J,bi,bj)=ZERO
244			TMIX(I,J,bi,bj)=TICE(I,J,bi,bj)
245			DO k=1,3
246			HEFF(I,J,k,bi,bj)=SEAICE_initialHEFF*HEFFM(i,j,bi,bj)
247			UICE(I,J,k,bi,bj)=ZERO
248			VICE(I,J,k,bi,bj)=ZERO
249			ENDDO
250			ENDDO
251			ENDDO
252			ENDDO
253			ENDDO
254
255			C-- Read initial sea-ice thickness from file if available.
256			IF ( HeffFile .NE. ' ' ) THEN
257			_BEGIN_MASTER( myThid )
258			CALL READ_FLD_XY_RL( HeffFile, ' ', ZETA, 0, myThid )
259			_END_MASTER(myThid)
260			_EXCH_XY_R8(ZETA,myThid)
261			DO bj=myByLo(myThid),myByHi(myThid)
262			DO bi=myBxLo(myThid),myBxHi(myThid)
263			DO j=1-OLy,sNy+OLy
264			DO i=1-OLx,sNx+OLx
265			DO k=1,3
266			HEFF(I,J,k,bi,bj) = MAX(ZETA(i,j,bi,bj),ZERO)
267			ENDDO
268			ENDDO
269			ENDDO
270			ENDDO
271			ENDDO
272			ENDIF
273
274			DO bj=myByLo(myThid),myByHi(myThid)
275			DO bi=myBxLo(myThid),myBxHi(myThid)
276			DO j=1-OLy,sNy+OLy
277			DO i=1-OLx,sNx+OLx
278			DO k=1,3
279			IF(HEFF(I,J,k,bi,bj).GT.ZERO) THEN
280			AREA(I,J,k,bi,bj)=ONE
281			ELSE
282			HSNOW(I,J,bi,bj)=ZERO
283			ENDIF
284			ENDDO
285			ENDDO
286			ENDDO
287			ENDDO
288			ENDDO
289
290			ENDIF
291
292			C--- Complete initialization
293			DO bj=myByLo(myThid),myByHi(myThid)
294			DO bi=myBxLo(myThid),myBxHi(myThid)
295			DO j=1-OLy,sNy+OLy
296			DO i=1-OLx,sNx+OLx
297			ZETA(I,J,bi,bj)=HEFF(I,J,1,bi,bj)*(1.0 _d 11)
298			ETA(I,J,bi,bj)=ZETA(I,J,bi,bj)/4.0 _d 0
299			ENDDO
300			ENDDO
301			IF ( useRealFreshWaterFlux .AND. .NOT.useThSIce ) THEN
302			DO j=1-OLy,sNy+OLy
303			DO i=1-OLx,sNx+OLx
304			sIceLoad(i,j,bi,bj) = HEFF(I,J,1,bi,bj)*SEAICE_rhoIce
305			& + HSNOW(I,J,bi,bj)* 330. _d 0
306
307			ENDDO
308			ENDDO
309			ENDIF
310			ENDDO
311			ENDDO
312
313			RETURN
314			END