CLEAN/code/shelfice_update_masks_JJ.F

C $Header: /u/gcmpack/MITgcm_contrib/shelfice_remeshing/CLEAN/code/shelfice_update_masks_JJ.F,v 1.1 2015/12/07 17:08:44 dgoldberg Exp $
C $Name:  $

#include "SHELFICE_OPTIONS.h"
#ifdef ALLOW_CTRL
# include "CTRL_OPTIONS.h"
#endif

CBOP
C     !ROUTINE: SHELFICE_UPDATE_MASKS
C     !INTERFACE:
      SUBROUTINE SHELFICE_UPDATE_MASKS_JJ(
     I     rF, recip_drF,
     U     hFacC,
     I     myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE SHELFICE_UPDATE_MASKS
C     | o modify topography factor hFacC according to ice shelf
C     |   topography
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "SURFACE.h"
#ifdef ALLOW_SHELFICE
# include "SHELFICE.h"
#endif /* ALLOW_SHELFICE */

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     rF        :: R-coordinate of face of cell (units of r).
C     recip_drF :: Recipricol of cell face separation along Z axis ( units of r ).
C     hFacC     :: Fraction of cell in vertical which is open (see GRID.h)
C     myThid    :: Number of this instance of SHELFICE_UPDATE_MASKS
      _RS rF        (1:Nr+1)
      _RS recip_drF (1:Nr)
      _RS hFacC     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:Nr,nSx,nSy)
      
      INTEGER myThid

#ifdef ALLOW_SHELFICE
C     !LOCAL VARIABLES:
C     == Local variables ==
C     bi,bj   :: tile indices
C     I,J,K   :: Loop counters
      INTEGER bi, bj
      INTEGER I, J, K
      _RL hFacCtmp
      _RL hFacMnSz

C- Update etaN
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO J = 1-OLy,sNy+OLy
          DO I = 1-OLx,sNx+OLx
           IF ( R_shelfice(I,J,bi,bj) .LT. 0.0) THEN
            IF (etah(I,J,bi,bj) .GT. SHELFICESplitThreshold ) THEN
             K  = MAX(1,kTopC(I,J,bi,bj))
             etaN(I,J,bi,bj) = etaN(I,J,bi,bj) - 1/recip_drF(K)
             etaH(I,J,bi,bj) = etaH(I,J,bi,bj) - 1/recip_drF(K)
            R_shelfIce(I,J,bi,bj) = R_shelfIce(I,J,bi,bj)+1/recip_drF(K)
             uVel(I,J,K-1,bi,bj)=uVel(I,J,K,bi,bj)
             uVel(I+1,J,K-1,bi,bj)=uVel(I+1,J,K,bi,bj)
             vVel(I,J,K-1,bi,bj)=vVel(I,J,K,bi,bj)
             vVel(I,J+1,K-1,bi,bj)=vVel(I,J+1,K,bi,bj)

             gvnm1(I,J,K-1,bi,bj)=gvnm1(I,J,K,bi,bj)
             gvnm1(I,J+1,K-1,bi,bj)=gvnm1(I,J+1,K,bi,bj)
             gunm1(I,J,K-1,bi,bj)=gunm1(I,J,K,bi,bj)
             gunm1(I+1,J,K-1,bi,bj)=gunm1(I,J,K,bi,bj)

             salt(I,J,K-1,bi,bj)=salt(I,J,K,bi,bj)
             theta(I,J,K-1,bi,bj)=theta(I,J,K,bi,bj)
             hfacC(I,J,K,bi,bj)=1.0

           ENDIF
             IF (R_shelfice(i,j,bi,bj) .NE. R_grounding(i,j,bi,bj))THEN
             IF (etah(I,J,bi,bj) .LT. SHELFICEMergeThreshold ) THEN
             K = MAX(1,kTopC(I,J,bi,bj))

             salt(I,J,K+1,bi,bj)=((salt(I,J,K,bi,bj)*(1/recip_drF(K)+
     &       etaN(I,J,bi,bj)))+(salt(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/(
     &       1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))

             theta(I,J,K+1,bi,bj)=((theta(I,J,K,bi,bj)*(1/recip_drF(K)+
     &      etaN(I,J,bi,bj)))+(theta(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/(
     &       1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))

             vVel(I,J,K+1,bi,bj)=((vVel(I,J,K,bi,bj)*(1/recip_drF(K)+
     &       etaN(I,J,bi,bj)))+(vVel(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/(
     &       1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))

            vVel(I,J+1,K+1,bi,bj)=((vVel(I,J+1,K,bi,bj)*(1/recip_drF(K)+
     &      etaN(I,J,bi,bj)))+(vVel(I,J+1,K+1,bi,bj)*1/recip_drF(K+1)))/
     &       (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))

             uVel(I,J,K+1,bi,bj)=((uVel(I,J,K,bi,bj)*(1/recip_drF(K)+
     &       etaN(I,J,bi,bj)))+(uVel(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/(
     &      1/recip_ drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))

            uVel(I+1,J,K+1,bi,bj)=((uVel(I+1,J,K,bi,bj)*(1/recip_drF(K)+
     &      etaN(I,J,bi,bj)))+(uVel(I+1,J,K+1,bi,bj)*1/recip_drF(K+1)))/
     &       (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))

             etaN(I,J,bi,bj) = etaN(I,J,bi,bj) +1/recip_drF(K)
             etaH(I,J,bi,bj) = etaH(I,J,bi,bj) +1/recip_drF(K)
           R_shelfice(I,J,bi,bj) = R_shelfice(I,J,bi,bj) -1/recip_drF(K)

         gunm1(I+1,J,K+1,bi,bj)=((gunm1(I+1,J,K,bi,bj)*(1/recip_drF(K)+
     &     etaN(I,J,bi,bj)))+(gunm1(I+1,J,K+1,bi,bj)*1/recip_drF(K+1)))/
     &       (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))

         gunm1(I,J,K+1,bi,bj)=((gunm1(I,J,K,bi,bj)*(1/recip_drF(K)+
     &     etaN(I,J,bi,bj)))+(gunm1(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/
     &       (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))

      gvnm1(I,J+1,K+1,bi,bj)=((gvnm1(I,J+1,K,bi,bj)*(1/recip_drF(K)+
     &     etaN(I,J,bi,bj)))+(gvnm1(I,J+1,K+1,bi,bj)*1/recip_drF(K+1)))/
     &       (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))

      gvnm1(I,J,K+1,bi,bj)=((gvnm1(I,J,K,bi,bj)*(1/recip_drF(K)+
     &     etaN(I,J,bi,bj)))+(gvnm1(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/
     &       (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))

             hfacC(I,J,K,bi,bj)=1.0
             ENDIF 
            ENDIF
           ENDIF
          ENDDO
         ENDDO
        ENDDO
       ENDDO

     
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO J = 1-OLy,sNy+OLy
          DO I = 1-OLx,sNx+OLx
            etaH(I,J,bi,bj)=etaN(I,J,bi,bj)
            etaHnm1(I,J,bi,bj)=etaH(I,J,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO


       DO   bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO J = 1-OLy,sNy+OLy
          DO I = 1-OLx,sNx+OLx
           K = MAX(1,kTopC(I,J,bi,bj))
            hfac_surfc(I,J,bi,bj)= ((etaH(I,J,bi,bJ) +(1/recip_drF(K)))
     &     *recip_drF(K))
          ENDDO
         ENDDO
        ENDDO
       ENDDO

     
C- fill in the overlap (+ BARRIER):
      _EXCH_XY_RS(R_shelfIce, myThid )

C--   Calculate lopping factor hFacC : Remove part outside of the domain
C     taking into account the Reference (=at rest) Surface Position Ro_shelfIce
      DO bj=myByLo(myThid), myByHi(myThid)
       DO bi=myBxLo(myThid), myBxHi(myThid)

C--   compute contributions of shelf ice to looping factors
        DO K=1, Nr
         hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
C      o Non-dimensional distance between grid boundary and model surface
           hFacCtmp = (rF(k)-R_shelfIce(I,J,bi,bj))*recip_drF(K)
C      o Reduce the previous fraction : substract the outside part.
           hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
C      o set to zero if empty Column :
           hFacCtmp = max( hFacCtmp, 0. _d 0)
C      o Impose minimum fraction and/or size (dimensional)
           IF (hFacCtmp.LT.hFacMnSz) THEN
            IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
             hFacC(I,J,K,bi,bj)=0.
            ELSE
             hFacC(I,J,K,bi,bj)=hFacMnSz
            ENDIF
           ELSE
             hFacC(I,J,K,bi,bj)=hFacCtmp
           ENDIF
          ENDDO
         ENDDO
        ENDDO

#ifdef ALLOW_SHIFWFLX_CONTROL
C     maskSHI is a hack to play along with the general ctrl-package
C     infrastructure, where only the k=1 layer of a 3D mask is used
C     for 2D fields. We cannot use maskInC instead, because routines
C     like ctrl_get_gen and ctrl_set_unpack_xy require 3D masks.
        DO K=1,Nr
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           maskSHI(I,J,K,bi,bj) = 0. _d 0
          ENDDO
         ENDDO
        ENDDO
        DO K=1,Nr
         DO J=1-OLy,sNy+OLy
          DO I=1-OLx,sNx+OLx
           IF ( ABS(R_shelfice(I,J,bi,bj)) .GT. 0. _d 0
     &          .AND. hFacC(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
            maskSHI(I,J,K,bi,bj) = 1. _d 0
            maskSHI(I,J,1,bi,bj) = 1. _d 0
           ENDIF
          ENDDO
         ENDDO
        ENDDO
#endif /* ALLOW_SHIFWFLX_CONTROL */

C - end bi,bj loops.
       ENDDO
      ENDDO
#endif /* ALLOW_SHELFICE */
      RETURN
      END
1	dgoldberg	1.2	C $Header: /u/gcmpack/MITgcm_contrib/shelfice_remeshing/CLEAN/code/shelfice_update_masks_JJ.F,v 1.1 2015/12/07 17:08:44 dgoldberg Exp $
2	dgoldberg	1.1	C $Name: $
3
4			#include "SHELFICE_OPTIONS.h"
5			#ifdef ALLOW_CTRL
6			# include "CTRL_OPTIONS.h"
7			#endif
8
9			CBOP
10			C !ROUTINE: SHELFICE_UPDATE_MASKS
11			C !INTERFACE:
12			SUBROUTINE SHELFICE_UPDATE_MASKS_JJ(
13			I rF, recip_drF,
14			U hFacC,
15			I myThid )
16			C !DESCRIPTION: \bv
17			C ==========================================================
18			C \| SUBROUTINE SHELFICE_UPDATE_MASKS
19			C \| o modify topography factor hFacC according to ice shelf
20			C \| topography
21			C ==========================================================
22			C \ev
23
24			C !USES:
25			IMPLICIT NONE
26			C === Global variables ===
27			#include "SIZE.h"
28			#include "EEPARAMS.h"
29			#include "PARAMS.h"
30			#include "DYNVARS.h"
31			#include "SURFACE.h"
32			#ifdef ALLOW_SHELFICE
33			# include "SHELFICE.h"
34			#endif /* ALLOW_SHELFICE */
35
36			C !INPUT/OUTPUT PARAMETERS:
37			C == Routine arguments ==
38			C rF :: R-coordinate of face of cell (units of r).
39			C recip_drF :: Recipricol of cell face separation along Z axis ( units of r ).
40			C hFacC :: Fraction of cell in vertical which is open (see GRID.h)
41			C myThid :: Number of this instance of SHELFICE_UPDATE_MASKS
42			_RS rF (1:Nr+1)
43			_RS recip_drF (1:Nr)
44			_RS hFacC (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:Nr,nSx,nSy)
45
46			INTEGER myThid
47
48			#ifdef ALLOW_SHELFICE
49			C !LOCAL VARIABLES:
50			C == Local variables ==
51			C bi,bj :: tile indices
52			C I,J,K :: Loop counters
53			INTEGER bi, bj
54			INTEGER I, J, K
55			_RL hFacCtmp
56			_RL hFacMnSz
57
58			C- Update etaN
59			DO bj = myByLo(myThid), myByHi(myThid)
60			DO bi = myBxLo(myThid), myBxHi(myThid)
61			DO J = 1-OLy,sNy+OLy
62			DO I = 1-OLx,sNx+OLx
63			IF ( R_shelfice(I,J,bi,bj) .LT. 0.0) THEN
64			IF (etah(I,J,bi,bj) .GT. SHELFICESplitThreshold ) THEN
65			K = MAX(1,kTopC(I,J,bi,bj))
66			etaN(I,J,bi,bj) = etaN(I,J,bi,bj) - 1/recip_drF(K)
67			etaH(I,J,bi,bj) = etaH(I,J,bi,bj) - 1/recip_drF(K)
68			R_shelfIce(I,J,bi,bj) = R_shelfIce(I,J,bi,bj)+1/recip_drF(K)
69			uVel(I,J,K-1,bi,bj)=uVel(I,J,K,bi,bj)
70			uVel(I+1,J,K-1,bi,bj)=uVel(I+1,J,K,bi,bj)
71			vVel(I,J,K-1,bi,bj)=vVel(I,J,K,bi,bj)
72			vVel(I,J+1,K-1,bi,bj)=vVel(I,J+1,K,bi,bj)
73
74			gvnm1(I,J,K-1,bi,bj)=gvnm1(I,J,K,bi,bj)
75			gvnm1(I,J+1,K-1,bi,bj)=gvnm1(I,J+1,K,bi,bj)
76			gunm1(I,J,K-1,bi,bj)=gunm1(I,J,K,bi,bj)
77			gunm1(I+1,J,K-1,bi,bj)=gunm1(I,J,K,bi,bj)
78
79			salt(I,J,K-1,bi,bj)=salt(I,J,K,bi,bj)
80			theta(I,J,K-1,bi,bj)=theta(I,J,K,bi,bj)
81			hfacC(I,J,K,bi,bj)=1.0
82
83			ENDIF
84			IF (R_shelfice(i,j,bi,bj) .NE. R_grounding(i,j,bi,bj))THEN
85			IF (etah(I,J,bi,bj) .LT. SHELFICEMergeThreshold ) THEN
86			K = MAX(1,kTopC(I,J,bi,bj))
87
88			salt(I,J,K+1,bi,bj)=((salt(I,J,K,bi,bj)*(1/recip_drF(K)+
89			& etaN(I,J,bi,bj)))+(salt(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/(
90			& 1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
91
92			theta(I,J,K+1,bi,bj)=((theta(I,J,K,bi,bj)*(1/recip_drF(K)+
93			& etaN(I,J,bi,bj)))+(theta(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/(
94			& 1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
95
96			vVel(I,J,K+1,bi,bj)=((vVel(I,J,K,bi,bj)*(1/recip_drF(K)+
97			& etaN(I,J,bi,bj)))+(vVel(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/(
98			& 1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
99
100			vVel(I,J+1,K+1,bi,bj)=((vVel(I,J+1,K,bi,bj)*(1/recip_drF(K)+
101			& etaN(I,J,bi,bj)))+(vVel(I,J+1,K+1,bi,bj)*1/recip_drF(K+1)))/
102			& (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
103
104			uVel(I,J,K+1,bi,bj)=((uVel(I,J,K,bi,bj)*(1/recip_drF(K)+
105			& etaN(I,J,bi,bj)))+(uVel(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/(
106			& 1/recip_ drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
107
108			uVel(I+1,J,K+1,bi,bj)=((uVel(I+1,J,K,bi,bj)*(1/recip_drF(K)+
109			& etaN(I,J,bi,bj)))+(uVel(I+1,J,K+1,bi,bj)*1/recip_drF(K+1)))/
110			& (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
111
112			etaN(I,J,bi,bj) = etaN(I,J,bi,bj) +1/recip_drF(K)
113			etaH(I,J,bi,bj) = etaH(I,J,bi,bj) +1/recip_drF(K)
114			R_shelfice(I,J,bi,bj) = R_shelfice(I,J,bi,bj) -1/recip_drF(K)
115
116			gunm1(I+1,J,K+1,bi,bj)=((gunm1(I+1,J,K,bi,bj)*(1/recip_drF(K)+
117			& etaN(I,J,bi,bj)))+(gunm1(I+1,J,K+1,bi,bj)*1/recip_drF(K+1)))/
118			& (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
119
120			gunm1(I,J,K+1,bi,bj)=((gunm1(I,J,K,bi,bj)*(1/recip_drF(K)+
121			& etaN(I,J,bi,bj)))+(gunm1(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/
122			& (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
123
124			gvnm1(I,J+1,K+1,bi,bj)=((gvnm1(I,J+1,K,bi,bj)*(1/recip_drF(K)+
125			& etaN(I,J,bi,bj)))+(gvnm1(I,J+1,K+1,bi,bj)*1/recip_drF(K+1)))/
126			& (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
127
128			gvnm1(I,J,K+1,bi,bj)=((gvnm1(I,J,K,bi,bj)*(1/recip_drF(K)+
129			& etaN(I,J,bi,bj)))+(gvnm1(I,J,K+1,bi,bj)*1/recip_drF(K+1)))/
130			& (1/recip_drF(K)+1/recip_drF(K+1)+etaN(I,J,bi,bj))
131
132			hfacC(I,J,K,bi,bj)=1.0
133			ENDIF
134			ENDIF
135			ENDIF
136			ENDDO
137			ENDDO
138			ENDDO
139			ENDDO
140
141
142			DO bj = myByLo(myThid), myByHi(myThid)
143			DO bi = myBxLo(myThid), myBxHi(myThid)
144			DO J = 1-OLy,sNy+OLy
145			DO I = 1-OLx,sNx+OLx
146			etaH(I,J,bi,bj)=etaN(I,J,bi,bj)
147			etaHnm1(I,J,bi,bj)=etaH(I,J,bi,bj)
148			ENDDO
149			ENDDO
150			ENDDO
151			ENDDO
152
153
154
155	dgoldberg	1.2	DO bj = myByLo(myThid), myByHi(myThid)
156	dgoldberg	1.1	DO bi = myBxLo(myThid), myBxHi(myThid)
157			DO J = 1-OLy,sNy+OLy
158			DO I = 1-OLx,sNx+OLx
159	dgoldberg	1.2	K = MAX(1,kTopC(I,J,bi,bj))
160	dgoldberg	1.1	hfac_surfc(I,J,bi,bj)= ((etaH(I,J,bi,bJ) +(1/recip_drF(K)))
161			& *recip_drF(K))
162			ENDDO
163			ENDDO
164			ENDDO
165			ENDDO
166
167
168
169
170
171			C- fill in the overlap (+ BARRIER):
172			_EXCH_XY_RS(R_shelfIce, myThid )
173
174			C-- Calculate lopping factor hFacC : Remove part outside of the domain
175			C taking into account the Reference (=at rest) Surface Position Ro_shelfIce
176			DO bj=myByLo(myThid), myByHi(myThid)
177			DO bi=myBxLo(myThid), myBxHi(myThid)
178
179			C-- compute contributions of shelf ice to looping factors
180			DO K=1, Nr
181			hFacMnSz=max( hFacMin, min(hFacMinDr*recip_drF(k),1. _d 0) )
182			DO J=1-OLy,sNy+OLy
183			DO I=1-OLx,sNx+OLx
184			C o Non-dimensional distance between grid boundary and model surface
185			hFacCtmp = (rF(k)-R_shelfIce(I,J,bi,bj))*recip_drF(K)
186			C o Reduce the previous fraction : substract the outside part.
187			hFacCtmp = hFacC(I,J,K,bi,bj) - max( hFacCtmp, 0. _d 0)
188			C o set to zero if empty Column :
189			hFacCtmp = max( hFacCtmp, 0. _d 0)
190			C o Impose minimum fraction and/or size (dimensional)
191			IF (hFacCtmp.LT.hFacMnSz) THEN
192			IF (hFacCtmp.LT.hFacMnSz*0.5) THEN
193			hFacC(I,J,K,bi,bj)=0.
194			ELSE
195			hFacC(I,J,K,bi,bj)=hFacMnSz
196			ENDIF
197			ELSE
198			hFacC(I,J,K,bi,bj)=hFacCtmp
199			ENDIF
200			ENDDO
201			ENDDO
202			ENDDO
203
204			#ifdef ALLOW_SHIFWFLX_CONTROL
205			C maskSHI is a hack to play along with the general ctrl-package
206			C infrastructure, where only the k=1 layer of a 3D mask is used
207			C for 2D fields. We cannot use maskInC instead, because routines
208			C like ctrl_get_gen and ctrl_set_unpack_xy require 3D masks.
209			DO K=1,Nr
210			DO J=1-OLy,sNy+OLy
211			DO I=1-OLx,sNx+OLx
212			maskSHI(I,J,K,bi,bj) = 0. _d 0
213			ENDDO
214			ENDDO
215			ENDDO
216			DO K=1,Nr
217			DO J=1-OLy,sNy+OLy
218			DO I=1-OLx,sNx+OLx
219			IF ( ABS(R_shelfice(I,J,bi,bj)) .GT. 0. _d 0
220			& .AND. hFacC(I,J,K,bi,bj) .NE. 0. _d 0 ) THEN
221			maskSHI(I,J,K,bi,bj) = 1. _d 0
222			maskSHI(I,J,1,bi,bj) = 1. _d 0
223			ENDIF
224			ENDDO
225			ENDDO
226			ENDDO
227			#endif /* ALLOW_SHIFWFLX_CONTROL */
228
229			C - end bi,bj loops.
230			ENDDO
231			ENDDO
232			#endif /* ALLOW_SHELFICE */
233			RETURN
234			END