/[MITgcm]/MITgcm/pkg/mom_common/mom_calc_hfacz.F
ViewVC logotype

Diff of /MITgcm/pkg/mom_common/mom_calc_hfacz.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by heimbach, Wed Jun 7 01:55:14 2006 UTC revision 1.3 by jmc, Fri Jun 16 14:44:06 2006 UTC
# Line 18  C Calculates the fractional thickness at Line 18  C Calculates the fractional thickness at
18  C !USES: ===============================================================  C !USES: ===============================================================
19        IMPLICIT NONE        IMPLICIT NONE
20  #include "SIZE.h"  #include "SIZE.h"
21    #include "EEPARAMS.h"
22    #include "PARAMS.h"
23  #include "GRID.h"  #include "GRID.h"
24    #ifdef ALLOW_EXCH2
25    #include "W2_EXCH2_TOPOLOGY.h"
26    #include "W2_EXCH2_PARAMS.h"
27    #endif /* ALLOW_EXCH2 */
28    
29  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
30  # include "EEPARAMS.h"  c # include "EEPARAMS.h"
31  # include "tamc.h"  # include "tamc.h"
32  # include "tamc_keys.h"  # include "tamc_keys.h"
33  #endif  #endif
# Line 41  C  r_hFacZ              :: reciprocal Line 47  C  r_hFacZ              :: reciprocal
47    
48  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
49  C  i,j                  :: loop indices  C  i,j                  :: loop indices
50    C  hZoption   :: forward mode option to select the way hFacZ is computed:
51    C                0 : = minimum of 4 hFacW,hFacS arround (consistent with
52    C                    definition of partial cell & mask near topography)
53    C                1 : = minimum of 2 average (hFacW)_j,(hFacS)_i
54    C                2 : = average of 4 hFacW,hFacS arround (consistent with
55    C                    how free surface affects hFacW,hFacS it using r* and
56    C                    without topography)
57        INTEGER I,J        INTEGER I,J
58  #ifdef ALLOW_DEPTH_CONTROL  #ifdef ALLOW_DEPTH_CONTROL
59        _RL hFacZOpen(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL hFacZOpen(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 52  C  i,j                  :: loop indices Line 65  C  i,j                  :: loop indices
65  # endif /* USE_SMOOTH_MIN */  # endif /* USE_SMOOTH_MIN */
66  #else  #else
67        _RL hFacZOpen        _RL hFacZOpen
68  #endif /* ALLOW_DEPTH_CONTROL */        INTEGER hZoption
69          LOGICAL northWestCorner, northEastCorner,
70         &        southWestCorner, southEastCorner
71          INTEGER myFace
72    #ifdef ALLOW_EXCH2
73          INTEGER myTile
74    #endif /* ALLOW_EXCH2 */
75  CEOP  CEOP
76          PARAMETER ( hZoption = 0 )
77    #endif /* ALLOW_DEPTH_CONTROL */
78    
79  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
80  #ifdef ALLOW_DEPTH_CONTROL  #ifdef ALLOW_DEPTH_CONTROL
# Line 148  CADJ STORE    r_hFacZ(:,:) = comlev1_bib Line 169  CADJ STORE    r_hFacZ(:,:) = comlev1_bib
169    
170  #else /* not ALLOW_DEPTH_CONTROL */  #else /* not ALLOW_DEPTH_CONTROL */
171    
172    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
173    
174    C-    Initialize hFacZ:
175    c     DO j=1-Oly,sNy+Oly
176    c      DO i=1-Olx,sNx+Olx
177    c        hFacZ(i,j)=0.
178    c      ENDDO
179    c     ENDDO
180    
181    C--   1rst row & column are not computed: fill with zero
182        DO i=1-Olx,sNx+Olx        DO i=1-Olx,sNx+Olx
183         hFacZ(i,1-Oly)=0.          hFacZ(i,1-Oly)=0.
        r_hFacZ(i,1-Oly)=0.  
184        ENDDO        ENDDO
   
185        DO j=2-Oly,sNy+Oly        DO j=2-Oly,sNy+Oly
186         hFacZ(1-Olx,j)=0.          hFacZ(1-Olx,j)=0.
187         r_hFacZ(1-Olx,j)=0.        ENDDO
188         DO i=2-Olx,sNx+Olx  
189          hFacZOpen=min(_hFacW(i,j,k,bi,bj),  C--   Calculate open water fraction at vorticity points
190       &                _hFacW(i,j-1,k,bi,bj))  
191          hFacZOpen=min(_hFacS(i,j,k,bi,bj),hFacZOpen)        IF ( hZoption.EQ.2 ) THEN
192          hFacZOpen=min(_hFacS(i-1,j,k,bi,bj),hFacZOpen)          DO j=2-Oly,sNy+Oly
193          hFacZ(i,j)=hFacZOpen           DO i=2-Olx,sNx+Olx
194    c         hFacZOpen=
195    c    &               ( _hFacW(i, j ,k,bi,bj)*rAw(i, j ,bi,bj)
196    c    &                +_hFacW(i,j-1,k,bi,bj)*rAw(i,j-1,bi,bj) )
197    c    &             + ( _hFacS( i ,j,k,bi,bj)*rAs( i ,j,bi,bj)
198    c    &                +_hFacS(i-1,j,k,bi,bj)*rAs(i-1,j,bi,bj) )
199    c         hFacZ(i,j) = 0.25 _d 0 * hFacZOpen*recip_rAz(i,j,bi,bj)
200              hFacZOpen=
201         &               ( _hFacW(i, j ,k,bi,bj)
202         &                +_hFacW(i,j-1,k,bi,bj) )
203         &             + ( _hFacS( i ,j,k,bi,bj)
204         &                +_hFacS(i-1,j,k,bi,bj) )
205              hFacZ(i,j) = 0.25 _d 0 * hFacZOpen
206             ENDDO
207            ENDDO
208          ELSEIF ( hZoption.EQ.1 ) THEN
209            DO j=2-Oly,sNy+Oly
210             DO i=2-Olx,sNx+Olx
211    c         hFacZOpen=MIN(
212    c    &                  _hFacW(i, j ,k,bi,bj)*rAw(i, j ,bi,bj)
213    c    &                + _hFacW(i,j-1,k,bi,bj)*rAw(i,j-1,bi,bj)
214    c    &                , _hFacS( i ,j,k,bi,bj)*rAs( i ,j,bi,bj)
215    c    &                + _hFacS(i-1,j,k,bi,bj)*rAs(i-1,j,bi,bj)
216    c    &                 )
217    c         hFacZ(i,j) = 0.5 _d 0 * hFacZOpen*recip_rAz(i,j,bi,bj)
218              hFacZOpen=MIN(
219         &                  _hFacW(i, j ,k,bi,bj)
220         &                + _hFacW(i,j-1,k,bi,bj)
221         &                , _hFacS( i ,j,k,bi,bj)
222         &                + _hFacS(i-1,j,k,bi,bj)
223         &                 )
224              hFacZ(i,j) = 0.5 _d 0 * hFacZOpen
225             ENDDO
226            ENDDO
227          ELSE
228            DO j=2-Oly,sNy+Oly
229             DO i=2-Olx,sNx+Olx
230              hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
231         &                  _hFacW(i,j-1,k,bi,bj))
232              hFacZOpen=MIN(_hFacS(i,j,k,bi,bj),hFacZOpen)
233              hFacZOpen=MIN(_hFacS(i-1,j,k,bi,bj),hFacZOpen)
234              hFacZ(i,j)=hFacZOpen
235             ENDDO
236            ENDDO
237          ENDIF
238    
239    C---+----1----+----2----+----3----+----4
240    C     Special stuff for Cubed Sphere
241          IF ( useCubedSphereExchange .AND. hZoption.GE.1 ) THEN
242    
243    #ifdef ALLOW_EXCH2
244            myTile = W2_myTileList(bi)
245            myFace = exch2_myFace(myTile)
246            southWestCorner = exch2_isWedge(myTile).EQ.1
247         &               .AND. exch2_isSedge(myTile).EQ.1
248            southEastCorner = exch2_isEedge(myTile).EQ.1
249         &              .AND. exch2_isSedge(myTile).EQ.1
250            northEastCorner = exch2_isEedge(myTile).EQ.1
251         &              .AND. exch2_isNedge(myTile).EQ.1
252            northWestCorner = exch2_isWedge(myTile).EQ.1
253         &              .AND. exch2_isNedge(myTile).EQ.1
254    #else
255            myFace = bi
256            southWestCorner = .TRUE.
257            southEastCorner = .TRUE.
258            northWestCorner = .TRUE.
259            northEastCorner = .TRUE.
260    #endif /* ALLOW_EXCH2 */
261    
262    
263            IF ( southWestCorner ) THEN
264             i=1
265             j=1
266             IF ( hZoption.EQ.1 ) THEN
267              hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
268         &                  _hFacW(i,j-1,k,bi,bj))
269              hFacZOpen=MIN(_hFacS(i,j,k,bi,bj),hFacZOpen)
270              hFacZ(i,j)=hFacZOpen
271             ELSE
272              IF ( MOD(myFace,2).EQ.1 ) THEN
273                hFacZOpen=
274         &               ( _hFacW(i,j-1,k,bi,bj)
275         &                +_hFacS( i ,j,k,bi,bj) )
276         &               + _hFacW(i, j ,k,bi,bj)
277              ELSE
278                hFacZOpen=
279         &               ( _hFacW(i, j ,k,bi,bj)
280         &                +_hFacW(i,j-1,k,bi,bj) )
281         &               + _hFacS( i ,j,k,bi,bj)
282              ENDIF
283              hFacZ(i,j) = hFacZOpen / 3. _d 0
284             ENDIF
285            ENDIF
286    
287            IF ( southEastCorner ) THEN
288             I=sNx+1
289             J=1
290    C-    to get the same truncation, independent from the face Nb:
291             IF ( hZoption.EQ.1 ) THEN
292              hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
293         &                  _hFacW(i,j-1,k,bi,bj))
294              hFacZOpen=MIN(_hFacS(i-1,j,k,bi,bj),hFacZOpen)
295              hFacZ(i,j)=hFacZOpen
296             ELSE
297              IF ( myFace.EQ.4 ) THEN
298                hFacZOpen=
299         &               ( _hFacS(i-1,j,k,bi,bj)
300         &                +_hFacW(i,j-1,k,bi,bj) )
301         &               + _hFacW(i, j ,k,bi,bj)
302              ELSEIF ( myFace.EQ.6 ) THEN
303                hFacZOpen=
304         &               ( _hFacW(i,j-1,k,bi,bj)
305         &                +_hFacW(i, j ,k,bi,bj) )
306         &               + _hFacS(i-1,j,k,bi,bj)
307              ELSE
308                hFacZOpen=
309         &               ( _hFacW(i, j ,k,bi,bj)
310         &                +_hFacS(i-1,j,k,bi,bj) )
311         &               + _hFacW(i,j-1,k,bi,bj)
312              ENDIF
313              hFacZ(i,j) = hFacZOpen / 3. _d 0
314             ENDIF
315            ENDIF
316    
317            IF ( northWestCorner ) THEN
318             i=1
319             j=sNy+1
320    C-    to get the same truncation, independent from the face Nb:
321             IF ( hZoption.EQ.1 ) THEN
322              hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
323         &                  _hFacW(i,j-1,k,bi,bj))
324              hFacZOpen=MIN(_hFacS(i,j,k,bi,bj),hFacZOpen)
325              hFacZ(i,j)=hFacZOpen
326             ELSE
327              IF ( myFace.EQ.1 ) THEN
328                hFacZOpen=
329         &               ( _hFacS( i ,j,k,bi,bj)
330         &                +_hFacW(i, j ,k,bi,bj) )
331         &               + _hFacW(i,j-1,k,bi,bj)
332              ELSEIF ( myFace.EQ.5 ) THEN
333                hFacZOpen=
334         &               ( _hFacW(i, j ,k,bi,bj)
335         &                +_hFacW(i,j-1,k,bi,bj) )
336         &               + _hFacS( i ,j,k,bi,bj)
337              ELSE
338                hFacZOpen=
339         &               ( _hFacW(i,j-1,k,bi,bj)
340         &                +_hFacS( i ,j,k,bi,bj) )
341         &               + _hFacW(i, j ,k,bi,bj)
342              ENDIF
343              hFacZ(i,j) = hFacZOpen / 3. _d 0
344             ENDIF
345            ENDIF
346    
347            IF ( northEastCorner ) THEN
348             i=sNx+1
349             j=sNy+1
350             IF ( hZoption.EQ.1 ) THEN
351              hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
352         &                  _hFacW(i,j-1,k,bi,bj))
353              hFacZOpen=MIN(_hFacS(i-1,j,k,bi,bj),hFacZOpen)
354              hFacZ(i,j)=hFacZOpen
355             ELSE
356              IF ( MOD(myFace,2).EQ.1 ) THEN
357                hFacZOpen=
358         &               ( _hFacW(i,j-1,k,bi,bj)
359         &                +_hFacW(i, j ,k,bi,bj) )
360         &               + _hFacS(i-1,j,k,bi,bj)
361              ELSE
362                hFacZOpen=
363         &               ( _hFacW(i, j ,k,bi,bj)
364         &                +_hFacS(i-1,j,k,bi,bj) )
365         &               + _hFacW(i,j-1,k,bi,bj)
366              ENDIF
367              hFacZ(i,j) = hFacZOpen / 3. _d 0
368             ENDIF
369            ENDIF
370    
371          ENDIF
372    C---+----1----+----2----+----3----+----4
373    
374    C--   Calculate reciprol:
375          DO j=1-Oly,sNy+Oly
376           DO i=1-Olx,sNx+Olx
377          IF (hFacZ(i,j).EQ.0.) THEN          IF (hFacZ(i,j).EQ.0.) THEN
378           r_hFacZ(i,j)=0.           r_hFacZ(i,j) = 0.
379          ELSE          ELSE
380           r_hFacZ(i,j)=1./hFacZ(i,j)           r_hFacZ(i,j) = 1. _d 0/hFacZ(i,j)
381          ENDIF          ENDIF
382         ENDDO         ENDDO
383        ENDDO        ENDDO
384    
385    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
386  #endif /* ALLOW_DEPTH_CONTROL */  #endif /* ALLOW_DEPTH_CONTROL */
387    
388        RETURN        RETURN
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22