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

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

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


Revision 1.7 - (hide annotations) (download)
Sun Jun 17 14:18:00 2012 UTC (12 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64t, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63o, checkpoint64
Changes since 1.6: +28 -29 lines
- rename SMOOTH*_R4,R8 function to the corresponding type (_RS & _RL)

1 jmc 1.7 C $Header: /u/gcmpack/MITgcm/pkg/mom_common/mom_calc_hfacz.F,v 1.6 2009/06/28 01:08:25 jmc Exp $
2 heimbach 1.2 C $Name: $
3 adcroft 1.1
4     #include "MOM_COMMON_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: MOM_CALC_HFACZ
8    
9     C !INTERFACE: ==========================================================
10     SUBROUTINE MOM_CALC_HFACZ(
11     I bi,bj,k,
12     O hFacZ,r_hFacZ,
13     I myThid)
14    
15     C !DESCRIPTION:
16     C Calculates the fractional thickness at vorticity points
17    
18     C !USES: ===============================================================
19     IMPLICIT NONE
20     #include "SIZE.h"
21 jmc 1.3 #include "EEPARAMS.h"
22     #include "PARAMS.h"
23 adcroft 1.1 #include "GRID.h"
24 jmc 1.3 #ifdef ALLOW_EXCH2
25 jmc 1.5 #include "W2_EXCH2_SIZE.h"
26 jmc 1.3 #include "W2_EXCH2_TOPOLOGY.h"
27     #endif /* ALLOW_EXCH2 */
28 adcroft 1.1
29 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
30     # include "tamc.h"
31     # include "tamc_keys.h"
32     #endif
33    
34 adcroft 1.1 C !INPUT PARAMETERS: ===================================================
35     C bi,bj :: tile indices
36     C k :: vertical level
37     C myThid :: thread number
38     INTEGER bi,bj,k
39     INTEGER myThid
40    
41     C !OUTPUT PARAMETERS: ==================================================
42     C hFacZ :: fractional thickness at vorticity points
43     C r_hFacZ :: reciprocal
44     _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
45     _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46    
47     C !LOCAL VARIABLES: ====================================================
48     C i,j :: loop indices
49 jmc 1.3 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 adcroft 1.1 INTEGER I,J
57 heimbach 1.2 #ifdef ALLOW_DEPTH_CONTROL
58     _RL hFacZOpen(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59     _RL hFacZOpenI(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60     _RL hFacZOpenJ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61     # ifdef USE_SMOOTH_MIN
62 jmc 1.7 _RS SMOOTHMIN_RS
63     EXTERNAL SMOOTHMIN_RS
64 heimbach 1.2 # endif /* USE_SMOOTH_MIN */
65     #else
66 jmc 1.4 _RS hFacZOpen
67 jmc 1.3 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
75     PARAMETER ( hZoption = 0 )
76 heimbach 1.2 #endif /* ALLOW_DEPTH_CONTROL */
77 adcroft 1.1
78 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
79     #ifdef ALLOW_DEPTH_CONTROL
80     act1 = bi - myBxLo(myThid)
81     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
82     act2 = bj - myByLo(myThid)
83     max2 = myByHi(myThid) - myByLo(myThid) + 1
84     act3 = myThid - 1
85     max3 = nTx*nTy
86     act4 = ikey_dynamics - 1
87     ikey = (act1 + 1) + act2*max1
88     & + act3*max1*max2
89     & + act4*max1*max2*max3
90     kkey = (ikey-1)*Nr + k
91     #endif /* ALLOW_DEPTH_CONTROL */
92     #endif /* ALLOW_AUTODIFF_TAMC */
93    
94 adcroft 1.1 C-- Calculate open water fraction at vorticity points
95    
96 heimbach 1.2 #ifdef ALLOW_DEPTH_CONTROL
97 jmc 1.7 DO j=1-OLy,sNy+OLy
98     DO i=1-OLx,sNx+OLx
99 heimbach 1.2 hFacZ(i,j) =0.
100     r_hFacZ(i,j) =0.
101     hFacZOpen(i,j) =0.
102     hFacZOpenJ(i,j)=0.
103     hFacZOpenJ(i,j)=0.
104     ENDDO
105     ENDDO
106    
107     #ifdef ALLOW_AUTODIFF_TAMC
108     CADJ STORE hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
109     CADJ STORE r_hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
110     #endif /* ALLOW_AUTODIFF_TAMC */
111 jmc 1.7 DO j=2-OLy,sNy+OLy
112     DO i=2-OLx,sNx+OLx
113 heimbach 1.2 hFacZOpenJ(i,j)=
114     #ifdef USE_SMOOTH_MIN
115 jmc 1.7 & SMOOTHMIN_RS(_hFacW(i ,j ,k,bi,bj),
116 heimbach 1.2 #else
117     & MIN(_hFacW(i ,j ,k,bi,bj),
118     #endif /* USE_SMOOTH_MIN */
119     & _hFacW(i ,j-1,k,bi,bj))
120     & *maskW(i,j,k,bi,bj)*maskW(i,j-1,k,bi,bj)
121     hFacZOpenI(i,j)=
122     #ifdef USE_SMOOTH_MIN
123 jmc 1.7 & SMOOTHMIN_RS(_hFacS(i ,j ,k,bi,bj),
124 heimbach 1.2 #else
125     & MIN(_hFacS(i ,j ,k,bi,bj),
126     #endif /* USE_SMOOTH_MIN */
127     & _hFacS(i-1,j ,k,bi,bj))
128     & *maskS(i,j,k,bi,bj)*maskS(i-1,j,k,bi,bj)
129     ENDDO
130     ENDDO
131     #ifdef ALLOW_AUTODIFF_TAMC
132     #ifdef ALLOW_DEPTH_CONTROL
133     CADJ STORE hFacZOpenI(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
134     CADJ STORE hFacZOpenJ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
135     #endif /* ALLOW_DEPTH_CONTROL */
136     #endif /* ALLOW_AUTODIFF_TAMC */
137 jmc 1.7 DO j=2-OLy,sNy+OLy
138     DO i=2-OLx,sNx+OLx
139 heimbach 1.2 hFacZ(i,j) =
140     #ifdef USE_SMOOTH_MIN
141 jmc 1.7 & SMOOTHMIN_RS(hFacZOpenI(i,j),hFacZOpenJ(i,j))
142 heimbach 1.2 #else
143     & MIN(hFacZOpenI(i,j),hFacZOpenJ(i,j))
144     #endif /* USE_SMOOTH_MIN */
145     & *maskW(i,j,k,bi,bj)*maskW(i,j-1,k,bi,bj)
146     & *maskS(i,j,k,bi,bj)*maskS(i-1,j,k,bi,bj)
147     ENDDO
148     ENDDO
149     #ifdef ALLOW_AUTODIFF_TAMC
150     #ifdef ALLOW_DEPTH_CONTROL
151     CADJ STORE hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
152     #endif /* ALLOW_DEPTH_CONTROL */
153     #endif /* ALLOW_AUTODIFF_TAMC */
154 jmc 1.7 DO j=2-OLy,sNy+OLy
155     DO i=2-OLx,sNx+OLx
156 heimbach 1.2 IF (hFacZ(i,j).EQ.0.) THEN
157     r_hFacZ(i,j)=0.
158     ELSE
159     r_hFacZ(i,j)=1./hFacZ(i,j)
160     ENDIF
161     ENDDO
162     ENDDO
163     #ifdef ALLOW_AUTODIFF_TAMC
164     #ifdef ALLOW_DEPTH_CONTROL
165     CADJ STORE r_hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
166     #endif /* ALLOW_DEPTH_CONTROL */
167     #endif /* ALLOW_AUTODIFF_TAMC */
168    
169     #else /* not ALLOW_DEPTH_CONTROL */
170    
171 jmc 1.3 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
172    
173     C- Initialize hFacZ:
174 jmc 1.7 c DO j=1-OLy,sNy+OLy
175     c DO i=1-OLx,sNx+OLx
176 jmc 1.3 c hFacZ(i,j)=0.
177     c ENDDO
178     c ENDDO
179    
180     C-- 1rst row & column are not computed: fill with zero
181 jmc 1.7 DO i=1-OLx,sNx+OLx
182     hFacZ(i,1-OLy)=0.
183 adcroft 1.1 ENDDO
184 jmc 1.7 DO j=2-OLy,sNy+OLy
185     hFacZ(1-OLx,j)=0.
186 jmc 1.3 ENDDO
187    
188     C-- Calculate open water fraction at vorticity points
189    
190     IF ( hZoption.EQ.2 ) THEN
191 jmc 1.7 DO j=2-OLy,sNy+OLy
192     DO i=2-OLx,sNx+OLx
193 jmc 1.3 c hFacZOpen=
194     c & ( _hFacW(i, j ,k,bi,bj)*rAw(i, j ,bi,bj)
195     c & +_hFacW(i,j-1,k,bi,bj)*rAw(i,j-1,bi,bj) )
196     c & + ( _hFacS( i ,j,k,bi,bj)*rAs( i ,j,bi,bj)
197     c & +_hFacS(i-1,j,k,bi,bj)*rAs(i-1,j,bi,bj) )
198     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 jmc 1.7 DO j=2-OLy,sNy+OLy
209     DO i=2-OLx,sNx+OLx
210 jmc 1.3 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 jmc 1.7 DO j=2-OLy,sNy+OLy
228     DO i=2-OLx,sNx+OLx
229 jmc 1.3 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 jmc 1.6 myTile = W2_myTileList(bi,bj)
244 jmc 1.3 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 adcroft 1.1
316 jmc 1.3 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 jmc 1.7 DO j=1-OLy,sNy+OLy
375     DO i=1-OLx,sNx+OLx
376 adcroft 1.1 IF (hFacZ(i,j).EQ.0.) THEN
377 jmc 1.3 r_hFacZ(i,j) = 0.
378 adcroft 1.1 ELSE
379 jmc 1.3 r_hFacZ(i,j) = 1. _d 0/hFacZ(i,j)
380 adcroft 1.1 ENDIF
381     ENDDO
382     ENDDO
383 jmc 1.3
384     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
385 heimbach 1.2 #endif /* ALLOW_DEPTH_CONTROL */
386 adcroft 1.1
387     RETURN
388     END

  ViewVC Help
Powered by ViewVC 1.1.22