/[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.6 - (hide annotations) (download)
Sun Jun 28 01:08:25 2009 UTC (14 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint62, checkpoint63, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61s, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.5: +2 -2 lines
add bj in exch2 arrays and S/R

1 jmc 1.6 C $Header: /u/gcmpack/MITgcm/pkg/mom_common/mom_calc_hfacz.F,v 1.5 2009/05/12 19:56:36 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 jmc 1.3 c # include "EEPARAMS.h"
31 heimbach 1.2 # include "tamc.h"
32     # include "tamc_keys.h"
33     #endif
34    
35 adcroft 1.1 C !INPUT PARAMETERS: ===================================================
36     C bi,bj :: tile indices
37     C k :: vertical level
38     C myThid :: thread number
39     INTEGER bi,bj,k
40     INTEGER myThid
41    
42     C !OUTPUT PARAMETERS: ==================================================
43     C hFacZ :: fractional thickness at vorticity points
44     C r_hFacZ :: reciprocal
45     _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46     _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47    
48     C !LOCAL VARIABLES: ====================================================
49     C i,j :: loop indices
50 jmc 1.3 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 adcroft 1.1 INTEGER I,J
58 heimbach 1.2 #ifdef ALLOW_DEPTH_CONTROL
59     _RL hFacZOpen(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60     _RL hFacZOpenI(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61     _RL hFacZOpenJ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62     # ifdef USE_SMOOTH_MIN
63     _RS smoothMin_R4
64     EXTERNAL smoothMin_R4
65     # endif /* USE_SMOOTH_MIN */
66     #else
67 jmc 1.4 _RS hFacZOpen
68 jmc 1.3 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
76     PARAMETER ( hZoption = 0 )
77 heimbach 1.2 #endif /* ALLOW_DEPTH_CONTROL */
78 adcroft 1.1
79 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
80     #ifdef ALLOW_DEPTH_CONTROL
81     act1 = bi - myBxLo(myThid)
82     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
83     act2 = bj - myByLo(myThid)
84     max2 = myByHi(myThid) - myByLo(myThid) + 1
85     act3 = myThid - 1
86     max3 = nTx*nTy
87     act4 = ikey_dynamics - 1
88     ikey = (act1 + 1) + act2*max1
89     & + act3*max1*max2
90     & + act4*max1*max2*max3
91     kkey = (ikey-1)*Nr + k
92     #endif /* ALLOW_DEPTH_CONTROL */
93     #endif /* ALLOW_AUTODIFF_TAMC */
94    
95 adcroft 1.1 C-- Calculate open water fraction at vorticity points
96    
97 heimbach 1.2 #ifdef ALLOW_DEPTH_CONTROL
98     DO j=1-Oly,sNy+Oly
99     DO i=1-Olx,sNx+Olx
100     hFacZ(i,j) =0.
101     r_hFacZ(i,j) =0.
102     hFacZOpen(i,j) =0.
103     hFacZOpenJ(i,j)=0.
104     hFacZOpenJ(i,j)=0.
105     ENDDO
106     ENDDO
107    
108     #ifdef ALLOW_AUTODIFF_TAMC
109     CADJ STORE hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
110     CADJ STORE r_hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
111     #endif /* ALLOW_AUTODIFF_TAMC */
112     DO j=2-Oly,sNy+Oly
113     DO i=2-Olx,sNx+Olx
114     hFacZOpenJ(i,j)=
115     #ifdef USE_SMOOTH_MIN
116     & smoothMin_R4(_hFacW(i ,j ,k,bi,bj),
117     #else
118     & MIN(_hFacW(i ,j ,k,bi,bj),
119     #endif /* USE_SMOOTH_MIN */
120     & _hFacW(i ,j-1,k,bi,bj))
121     & *maskW(i,j,k,bi,bj)*maskW(i,j-1,k,bi,bj)
122     hFacZOpenI(i,j)=
123     #ifdef USE_SMOOTH_MIN
124     & smoothMin_R4(_hFacS(i ,j ,k,bi,bj),
125     #else
126     & MIN(_hFacS(i ,j ,k,bi,bj),
127     #endif /* USE_SMOOTH_MIN */
128     & _hFacS(i-1,j ,k,bi,bj))
129     & *maskS(i,j,k,bi,bj)*maskS(i-1,j,k,bi,bj)
130     ENDDO
131     ENDDO
132     #ifdef ALLOW_AUTODIFF_TAMC
133     #ifdef ALLOW_DEPTH_CONTROL
134     CADJ STORE hFacZOpenI(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
135     CADJ STORE hFacZOpenJ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
136     #endif /* ALLOW_DEPTH_CONTROL */
137     #endif /* ALLOW_AUTODIFF_TAMC */
138     DO j=2-Oly,sNy+Oly
139     DO i=2-Olx,sNx+Olx
140     hFacZ(i,j) =
141     #ifdef USE_SMOOTH_MIN
142     & smoothMin_R4(hFacZOpenI(i,j),hFacZOpenJ(i,j))
143     #else
144     & MIN(hFacZOpenI(i,j),hFacZOpenJ(i,j))
145     #endif /* USE_SMOOTH_MIN */
146     & *maskW(i,j,k,bi,bj)*maskW(i,j-1,k,bi,bj)
147     & *maskS(i,j,k,bi,bj)*maskS(i-1,j,k,bi,bj)
148     ENDDO
149     ENDDO
150     #ifdef ALLOW_AUTODIFF_TAMC
151     #ifdef ALLOW_DEPTH_CONTROL
152     CADJ STORE hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
153     #endif /* ALLOW_DEPTH_CONTROL */
154     #endif /* ALLOW_AUTODIFF_TAMC */
155     DO j=2-Oly,sNy+Oly
156     DO i=2-Olx,sNx+Olx
157     IF (hFacZ(i,j).EQ.0.) THEN
158     r_hFacZ(i,j)=0.
159     ELSE
160     r_hFacZ(i,j)=1./hFacZ(i,j)
161     ENDIF
162     ENDDO
163     ENDDO
164     #ifdef ALLOW_AUTODIFF_TAMC
165     #ifdef ALLOW_DEPTH_CONTROL
166     CADJ STORE r_hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
167     #endif /* ALLOW_DEPTH_CONTROL */
168     #endif /* ALLOW_AUTODIFF_TAMC */
169    
170     #else /* not ALLOW_DEPTH_CONTROL */
171    
172 jmc 1.3 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 adcroft 1.1 DO i=1-Olx,sNx+Olx
183 jmc 1.3 hFacZ(i,1-Oly)=0.
184 adcroft 1.1 ENDDO
185 jmc 1.3 DO j=2-Oly,sNy+Oly
186     hFacZ(1-Olx,j)=0.
187     ENDDO
188    
189     C-- Calculate open water fraction at vorticity points
190    
191     IF ( hZoption.EQ.2 ) THEN
192     DO j=2-Oly,sNy+Oly
193     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 jmc 1.6 myTile = W2_myTileList(bi,bj)
245 jmc 1.3 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 adcroft 1.1
317 jmc 1.3 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 adcroft 1.1 IF (hFacZ(i,j).EQ.0.) THEN
378 jmc 1.3 r_hFacZ(i,j) = 0.
379 adcroft 1.1 ELSE
380 jmc 1.3 r_hFacZ(i,j) = 1. _d 0/hFacZ(i,j)
381 adcroft 1.1 ENDIF
382     ENDDO
383     ENDDO
384 jmc 1.3
385     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
386 heimbach 1.2 #endif /* ALLOW_DEPTH_CONTROL */
387 adcroft 1.1
388     RETURN
389     END

  ViewVC Help
Powered by ViewVC 1.1.22