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