/[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.8 - (hide annotations) (download)
Sun Feb 9 18:57:01 2014 UTC (11 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64u
Changes since 1.7: +47 -52 lines
- fix initialisation of hFacZOpenJ (instead of initialising hFacZOpenI 2 times)
- cleaning;

1 jmc 1.8 C $Header: /u/gcmpack/MITgcm/pkg/mom_common/mom_calc_hfacz.F,v 1.7 2012/06/17 14:18:00 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 jmc 1.8 I bi, bj, k,
12     O hFacZ, r_hFacZ,
13     I myThid )
14 adcroft 1.1
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.8 #ifdef ALLOW_DEPTH_CONTROL
25     # ifdef ALLOW_AUTODIFF_TAMC
26     # include "tamc.h"
27     # include "tamc_keys.h"
28     # endif
29     #else /* ALLOW_DEPTH_CONTROL */
30     # ifdef ALLOW_EXCH2
31     # include "W2_EXCH2_SIZE.h"
32     # include "W2_EXCH2_TOPOLOGY.h"
33     # endif
34     #endif /* ALLOW_DEPTH_CONTROL */
35 heimbach 1.2
36 adcroft 1.1 C !INPUT PARAMETERS: ===================================================
37 jmc 1.8 C bi,bj :: tile indices
38     C k :: vertical level
39     C myThid :: my Thread Id number
40     INTEGER bi, bj, k
41 adcroft 1.1 INTEGER myThid
42    
43     C !OUTPUT PARAMETERS: ==================================================
44 jmc 1.8 C hFacZ :: fractional thickness at vorticity points
45     C r_hFacZ :: reciprocal
46 adcroft 1.1 _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47     _RS r_hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48    
49     C !LOCAL VARIABLES: ====================================================
50 jmc 1.8 C i,j :: loop indices
51 jmc 1.3 C hZoption :: forward mode option to select the way hFacZ is computed:
52     C 0 : = minimum of 4 hFacW,hFacS arround (consistent with
53     C definition of partial cell & mask near topography)
54     C 1 : = minimum of 2 average (hFacW)_j,(hFacS)_i
55     C 2 : = average of 4 hFacW,hFacS arround (consistent with
56     C how free surface affects hFacW,hFacS it using r* and
57     C without topography)
58 jmc 1.8 INTEGER i,j
59 heimbach 1.2 #ifdef ALLOW_DEPTH_CONTROL
60 jmc 1.8 _RL hFacZOpen (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61 heimbach 1.2 _RL hFacZOpenI(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62     _RL hFacZOpenJ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63     # ifdef USE_SMOOTH_MIN
64 jmc 1.7 _RS SMOOTHMIN_RS
65     EXTERNAL SMOOTHMIN_RS
66 heimbach 1.2 # endif /* USE_SMOOTH_MIN */
67 jmc 1.8 #else /* ALLOW_DEPTH_CONTROL */
68 jmc 1.4 _RS hFacZOpen
69 jmc 1.3 INTEGER hZoption
70     LOGICAL northWestCorner, northEastCorner,
71     & southWestCorner, southEastCorner
72     INTEGER myFace
73 jmc 1.8 # ifdef ALLOW_EXCH2
74 jmc 1.3 INTEGER myTile
75 jmc 1.8 # endif /* ALLOW_EXCH2 */
76 jmc 1.3 PARAMETER ( hZoption = 0 )
77 heimbach 1.2 #endif /* ALLOW_DEPTH_CONTROL */
78 jmc 1.8 CEOP
79    
80     C-- Calculate open water fraction at vorticity points
81    
82     #ifdef ALLOW_DEPTH_CONTROL
83 adcroft 1.1
84 heimbach 1.2 #ifdef ALLOW_AUTODIFF_TAMC
85     act1 = bi - myBxLo(myThid)
86     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
87     act2 = bj - myByLo(myThid)
88     max2 = myByHi(myThid) - myByLo(myThid) + 1
89     act3 = myThid - 1
90     max3 = nTx*nTy
91     act4 = ikey_dynamics - 1
92     ikey = (act1 + 1) + act2*max1
93     & + act3*max1*max2
94     & + act4*max1*max2*max3
95     kkey = (ikey-1)*Nr + k
96     #endif /* ALLOW_AUTODIFF_TAMC */
97    
98 jmc 1.7 DO j=1-OLy,sNy+OLy
99     DO i=1-OLx,sNx+OLx
100 heimbach 1.2 hFacZ(i,j) =0.
101     r_hFacZ(i,j) =0.
102     hFacZOpen(i,j) =0.
103 jmc 1.8 hFacZOpenI(i,j)=0.
104 heimbach 1.2 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 jmc 1.7 DO j=2-OLy,sNy+OLy
113     DO i=2-OLx,sNx+OLx
114 heimbach 1.2 hFacZOpenJ(i,j)=
115     #ifdef USE_SMOOTH_MIN
116 jmc 1.7 & SMOOTHMIN_RS(_hFacW(i ,j ,k,bi,bj),
117 heimbach 1.2 #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 jmc 1.7 & SMOOTHMIN_RS(_hFacS(i ,j ,k,bi,bj),
125 heimbach 1.2 #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 jmc 1.8 #ifdef ALLOW_AUTODIFF_TAMC
133 heimbach 1.2 CADJ STORE hFacZOpenI(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
134     CADJ STORE hFacZOpenJ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
135 jmc 1.8 #endif /* ALLOW_AUTODIFF_TAMC */
136 jmc 1.7 DO j=2-OLy,sNy+OLy
137     DO i=2-OLx,sNx+OLx
138 heimbach 1.2 hFacZ(i,j) =
139     #ifdef USE_SMOOTH_MIN
140 jmc 1.7 & SMOOTHMIN_RS(hFacZOpenI(i,j),hFacZOpenJ(i,j))
141 heimbach 1.2 #else
142     & MIN(hFacZOpenI(i,j),hFacZOpenJ(i,j))
143     #endif /* USE_SMOOTH_MIN */
144     & *maskW(i,j,k,bi,bj)*maskW(i,j-1,k,bi,bj)
145     & *maskS(i,j,k,bi,bj)*maskS(i-1,j,k,bi,bj)
146     ENDDO
147     ENDDO
148 jmc 1.8 #ifdef ALLOW_AUTODIFF_TAMC
149 heimbach 1.2 CADJ STORE hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
150 jmc 1.8 #endif /* ALLOW_AUTODIFF_TAMC */
151 jmc 1.7 DO j=2-OLy,sNy+OLy
152     DO i=2-OLx,sNx+OLx
153 heimbach 1.2 IF (hFacZ(i,j).EQ.0.) THEN
154 jmc 1.8 r_hFacZ(i,j) = 0. _d 0
155 heimbach 1.2 ELSE
156 jmc 1.8 r_hFacZ(i,j) = 1. _d 0/hFacZ(i,j)
157 heimbach 1.2 ENDIF
158     ENDDO
159     ENDDO
160 jmc 1.8 #ifdef ALLOW_AUTODIFF_TAMC
161 heimbach 1.2 CADJ STORE r_hFacZ(:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
162 jmc 1.8 #endif /* ALLOW_AUTODIFF_TAMC */
163 heimbach 1.2
164     #else /* not ALLOW_DEPTH_CONTROL */
165    
166 jmc 1.3 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
167    
168     C- Initialize hFacZ:
169 jmc 1.7 c DO j=1-OLy,sNy+OLy
170     c DO i=1-OLx,sNx+OLx
171 jmc 1.3 c hFacZ(i,j)=0.
172     c ENDDO
173     c ENDDO
174    
175     C-- 1rst row & column are not computed: fill with zero
176 jmc 1.7 DO i=1-OLx,sNx+OLx
177     hFacZ(i,1-OLy)=0.
178 adcroft 1.1 ENDDO
179 jmc 1.7 DO j=2-OLy,sNy+OLy
180     hFacZ(1-OLx,j)=0.
181 jmc 1.3 ENDDO
182    
183     C-- Calculate open water fraction at vorticity points
184    
185     IF ( hZoption.EQ.2 ) THEN
186 jmc 1.7 DO j=2-OLy,sNy+OLy
187     DO i=2-OLx,sNx+OLx
188 jmc 1.3 c hFacZOpen=
189     c & ( _hFacW(i, j ,k,bi,bj)*rAw(i, j ,bi,bj)
190     c & +_hFacW(i,j-1,k,bi,bj)*rAw(i,j-1,bi,bj) )
191     c & + ( _hFacS( i ,j,k,bi,bj)*rAs( i ,j,bi,bj)
192     c & +_hFacS(i-1,j,k,bi,bj)*rAs(i-1,j,bi,bj) )
193     c hFacZ(i,j) = 0.25 _d 0 * hFacZOpen*recip_rAz(i,j,bi,bj)
194     hFacZOpen=
195     & ( _hFacW(i, j ,k,bi,bj)
196     & +_hFacW(i,j-1,k,bi,bj) )
197     & + ( _hFacS( i ,j,k,bi,bj)
198     & +_hFacS(i-1,j,k,bi,bj) )
199     hFacZ(i,j) = 0.25 _d 0 * hFacZOpen
200     ENDDO
201     ENDDO
202     ELSEIF ( hZoption.EQ.1 ) THEN
203 jmc 1.7 DO j=2-OLy,sNy+OLy
204     DO i=2-OLx,sNx+OLx
205 jmc 1.3 c hFacZOpen=MIN(
206     c & _hFacW(i, j ,k,bi,bj)*rAw(i, j ,bi,bj)
207     c & + _hFacW(i,j-1,k,bi,bj)*rAw(i,j-1,bi,bj)
208     c & , _hFacS( i ,j,k,bi,bj)*rAs( i ,j,bi,bj)
209     c & + _hFacS(i-1,j,k,bi,bj)*rAs(i-1,j,bi,bj)
210     c & )
211     c hFacZ(i,j) = 0.5 _d 0 * hFacZOpen*recip_rAz(i,j,bi,bj)
212     hFacZOpen=MIN(
213     & _hFacW(i, j ,k,bi,bj)
214     & + _hFacW(i,j-1,k,bi,bj)
215     & , _hFacS( i ,j,k,bi,bj)
216     & + _hFacS(i-1,j,k,bi,bj)
217     & )
218     hFacZ(i,j) = 0.5 _d 0 * hFacZOpen
219     ENDDO
220     ENDDO
221     ELSE
222 jmc 1.7 DO j=2-OLy,sNy+OLy
223     DO i=2-OLx,sNx+OLx
224 jmc 1.3 hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
225     & _hFacW(i,j-1,k,bi,bj))
226     hFacZOpen=MIN(_hFacS(i,j,k,bi,bj),hFacZOpen)
227     hFacZOpen=MIN(_hFacS(i-1,j,k,bi,bj),hFacZOpen)
228     hFacZ(i,j)=hFacZOpen
229     ENDDO
230     ENDDO
231     ENDIF
232    
233 jmc 1.8 C-----------------------------------------
234 jmc 1.3 C Special stuff for Cubed Sphere
235     IF ( useCubedSphereExchange .AND. hZoption.GE.1 ) THEN
236    
237     #ifdef ALLOW_EXCH2
238 jmc 1.6 myTile = W2_myTileList(bi,bj)
239 jmc 1.3 myFace = exch2_myFace(myTile)
240     southWestCorner = exch2_isWedge(myTile).EQ.1
241 jmc 1.8 & .AND. exch2_isSedge(myTile).EQ.1
242 jmc 1.3 southEastCorner = exch2_isEedge(myTile).EQ.1
243     & .AND. exch2_isSedge(myTile).EQ.1
244     northEastCorner = exch2_isEedge(myTile).EQ.1
245     & .AND. exch2_isNedge(myTile).EQ.1
246     northWestCorner = exch2_isWedge(myTile).EQ.1
247     & .AND. exch2_isNedge(myTile).EQ.1
248     #else
249     myFace = bi
250     southWestCorner = .TRUE.
251     southEastCorner = .TRUE.
252     northWestCorner = .TRUE.
253     northEastCorner = .TRUE.
254     #endif /* ALLOW_EXCH2 */
255    
256     IF ( southWestCorner ) THEN
257     i=1
258     j=1
259     IF ( hZoption.EQ.1 ) THEN
260     hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
261     & _hFacW(i,j-1,k,bi,bj))
262     hFacZOpen=MIN(_hFacS(i,j,k,bi,bj),hFacZOpen)
263     hFacZ(i,j)=hFacZOpen
264     ELSE
265     IF ( MOD(myFace,2).EQ.1 ) THEN
266     hFacZOpen=
267     & ( _hFacW(i,j-1,k,bi,bj)
268     & +_hFacS( i ,j,k,bi,bj) )
269     & + _hFacW(i, j ,k,bi,bj)
270     ELSE
271     hFacZOpen=
272     & ( _hFacW(i, j ,k,bi,bj)
273     & +_hFacW(i,j-1,k,bi,bj) )
274     & + _hFacS( i ,j,k,bi,bj)
275     ENDIF
276     hFacZ(i,j) = hFacZOpen / 3. _d 0
277     ENDIF
278     ENDIF
279    
280     IF ( southEastCorner ) THEN
281 jmc 1.8 i=sNx+1
282     j=1
283 jmc 1.3 C- to get the same truncation, independent from the face Nb:
284     IF ( hZoption.EQ.1 ) THEN
285     hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
286     & _hFacW(i,j-1,k,bi,bj))
287     hFacZOpen=MIN(_hFacS(i-1,j,k,bi,bj),hFacZOpen)
288     hFacZ(i,j)=hFacZOpen
289     ELSE
290     IF ( myFace.EQ.4 ) THEN
291     hFacZOpen=
292     & ( _hFacS(i-1,j,k,bi,bj)
293     & +_hFacW(i,j-1,k,bi,bj) )
294     & + _hFacW(i, j ,k,bi,bj)
295     ELSEIF ( myFace.EQ.6 ) THEN
296     hFacZOpen=
297     & ( _hFacW(i,j-1,k,bi,bj)
298     & +_hFacW(i, j ,k,bi,bj) )
299     & + _hFacS(i-1,j,k,bi,bj)
300     ELSE
301     hFacZOpen=
302     & ( _hFacW(i, j ,k,bi,bj)
303     & +_hFacS(i-1,j,k,bi,bj) )
304     & + _hFacW(i,j-1,k,bi,bj)
305     ENDIF
306     hFacZ(i,j) = hFacZOpen / 3. _d 0
307     ENDIF
308     ENDIF
309 adcroft 1.1
310 jmc 1.3 IF ( northWestCorner ) THEN
311     i=1
312     j=sNy+1
313     C- to get the same truncation, independent from the face Nb:
314     IF ( hZoption.EQ.1 ) THEN
315     hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
316     & _hFacW(i,j-1,k,bi,bj))
317     hFacZOpen=MIN(_hFacS(i,j,k,bi,bj),hFacZOpen)
318     hFacZ(i,j)=hFacZOpen
319     ELSE
320     IF ( myFace.EQ.1 ) THEN
321     hFacZOpen=
322     & ( _hFacS( i ,j,k,bi,bj)
323     & +_hFacW(i, j ,k,bi,bj) )
324     & + _hFacW(i,j-1,k,bi,bj)
325     ELSEIF ( myFace.EQ.5 ) THEN
326     hFacZOpen=
327     & ( _hFacW(i, j ,k,bi,bj)
328     & +_hFacW(i,j-1,k,bi,bj) )
329     & + _hFacS( i ,j,k,bi,bj)
330     ELSE
331     hFacZOpen=
332     & ( _hFacW(i,j-1,k,bi,bj)
333     & +_hFacS( i ,j,k,bi,bj) )
334     & + _hFacW(i, j ,k,bi,bj)
335     ENDIF
336     hFacZ(i,j) = hFacZOpen / 3. _d 0
337     ENDIF
338     ENDIF
339    
340     IF ( northEastCorner ) THEN
341     i=sNx+1
342     j=sNy+1
343     IF ( hZoption.EQ.1 ) THEN
344     hFacZOpen=MIN(_hFacW(i,j,k,bi,bj),
345     & _hFacW(i,j-1,k,bi,bj))
346     hFacZOpen=MIN(_hFacS(i-1,j,k,bi,bj),hFacZOpen)
347     hFacZ(i,j)=hFacZOpen
348     ELSE
349     IF ( MOD(myFace,2).EQ.1 ) THEN
350     hFacZOpen=
351     & ( _hFacW(i,j-1,k,bi,bj)
352     & +_hFacW(i, j ,k,bi,bj) )
353     & + _hFacS(i-1,j,k,bi,bj)
354     ELSE
355     hFacZOpen=
356     & ( _hFacW(i, j ,k,bi,bj)
357     & +_hFacS(i-1,j,k,bi,bj) )
358     & + _hFacW(i,j-1,k,bi,bj)
359     ENDIF
360     hFacZ(i,j) = hFacZOpen / 3. _d 0
361     ENDIF
362     ENDIF
363    
364     ENDIF
365 jmc 1.8 C Special Cubed Sphere block ends here
366     C-----------------------------------------
367 jmc 1.3
368     C-- Calculate reciprol:
369 jmc 1.7 DO j=1-OLy,sNy+OLy
370     DO i=1-OLx,sNx+OLx
371 adcroft 1.1 IF (hFacZ(i,j).EQ.0.) THEN
372 jmc 1.3 r_hFacZ(i,j) = 0.
373 adcroft 1.1 ELSE
374 jmc 1.3 r_hFacZ(i,j) = 1. _d 0/hFacZ(i,j)
375 adcroft 1.1 ENDIF
376     ENDDO
377     ENDDO
378 jmc 1.3
379     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
380 heimbach 1.2 #endif /* ALLOW_DEPTH_CONTROL */
381 adcroft 1.1
382     RETURN
383     END

  ViewVC Help
Powered by ViewVC 1.1.22