/[MITgcm]/MITgcm/pkg/obcs/obcs_init_fixed.F
ViewVC logotype

Contents of /MITgcm/pkg/obcs/obcs_init_fixed.F

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


Revision 1.19 - (show annotations) (download)
Tue Nov 25 01:08:20 2014 UTC (9 years, 6 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65l, checkpoint65m, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65o, HEAD
Changes since 1.18: +8 -5 lines
add call to new S/R OBCS_SET_CONNECT

1 C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_init_fixed.F,v 1.18 2012/09/18 20:09:17 jmc Exp $
2 C $Name: $
3
4 #include "OBCS_OPTIONS.h"
5
6 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
7 CBOP
8 C !ROUTINE: OBCS_INIT_FIXED
9
10 C !INTERFACE:
11 SUBROUTINE OBCS_INIT_FIXED( myThid )
12
13 C !DESCRIPTION:
14 C *==========================================================*
15 C | SUBROUTINE OBCS_INIT_FIXED
16 C | o Initialise OBCs fixed arrays
17 C *==========================================================*
18
19 C !USES:
20 IMPLICIT NONE
21
22 C === Global variables ===
23 #include "SIZE.h"
24 #include "EEPARAMS.h"
25 #include "PARAMS.h"
26 #include "GRID.h"
27 #include "OBCS_PARAMS.h"
28 #include "OBCS_GRID.h"
29
30 C !INPUT/OUTPUT PARAMETERS:
31 C myThid :: my Thread Id. number
32 INTEGER myThid
33 CEOP
34
35 #ifdef ALLOW_OBCS
36 C !LOCAL VARIABLES:
37 C msgBuf :: Informational/error message buffer
38 C OB_ApplX :: number of grid points (in X dir) overwritten by obcs_apply
39 C OB_ApplY :: number of grid points (in Y dir) overwritten by obcs_apply
40 C bi,bj :: tile indices
41 C i, j :: Loop counters
42 CHARACTER*(MAX_LEN_MBUF) msgBuf, errMsg
43 INTEGER OB_ApplX
44 INTEGER OB_ApplY
45 INTEGER bi, bj
46 INTEGER i, j
47 INTEGER im, jm
48 INTEGER iB, jB
49 LOGICAL flag
50 INTEGER ioUnit
51
52 #ifdef ALLOW_DEBUG
53 IF (debugMode) CALL DEBUG_ENTER('OBCS_INIT_FIXED',myThid)
54 #endif
55
56 C== Set Interior mask at Cell Center:
57
58 DO bj = myByLo(myThid), myByHi(myThid)
59 DO bi = myBxLo(myThid), myBxHi(myThid)
60 DO j=1-OLy,sNy+OLy
61 DO i=1-OLx,sNx+OLx
62 OBCS_insideMask(i,j,bi,bj) = 1.
63 ENDDO
64 ENDDO
65 ENDDO
66 ENDDO
67
68 IF ( insideOBmaskFile.EQ.' ' ) THEN
69 C-- If no maskFile is provided, set Interior mask from OB list of indices
70
71 DO bj = myByLo(myThid), myByHi(myThid)
72 DO bi = myBxLo(myThid), myBxHi(myThid)
73 DO j=1,sNy
74 C- Eastern boundary
75 IF ( OB_Ie(j,bi,bj).NE.OB_indexNone ) THEN
76 flag = .TRUE.
77 DO i=OB_Ie(j,bi,bj),sNx
78 flag = flag .AND.
79 & kSurfC(i,j,bi,bj).LE.Nr .AND. i.NE.OB_Iw(j,bi,bj)
80 IF ( flag ) OBCS_insideMask(i,j,bi,bj) = 0.
81 ENDDO
82 ENDIF
83 C- Western boundary
84 IF ( OB_Iw(j,bi,bj).NE.OB_indexNone ) THEN
85 flag = .TRUE.
86 DO i=OB_Iw(j,bi,bj),1,-1
87 flag = flag .AND.
88 & kSurfC(i,j,bi,bj).LE.Nr .AND. i.NE.OB_Ie(j,bi,bj)
89 IF ( flag ) OBCS_insideMask(i,j,bi,bj) = 0.
90 ENDDO
91 ENDIF
92 ENDDO
93 DO i=1,sNx
94 C- Northern boundary
95 IF ( OB_Jn(i,bi,bj).NE.OB_indexNone ) THEN
96 flag = .TRUE.
97 DO j=OB_Jn(i,bi,bj),sNy
98 flag = flag .AND.
99 & kSurfC(i,j,bi,bj).LE.Nr .AND. j.NE.OB_Js(i,bi,bj)
100 IF ( flag ) OBCS_insideMask(i,j,bi,bj) = 0.
101 ENDDO
102 ENDIF
103 C- Southern boundary
104 IF ( OB_Js(i,bi,bj).NE.OB_indexNone ) THEN
105 flag = .TRUE.
106 DO j=OB_Js(i,bi,bj),1,-1
107 flag = flag .AND.
108 & kSurfC(i,j,bi,bj).LE.Nr .AND. j.NE.OB_Jn(i,bi,bj)
109 IF ( flag ) OBCS_insideMask(i,j,bi,bj) = 0.
110 ENDDO
111 ENDIF
112 ENDDO
113
114 C-- end bi,bj loops
115 ENDDO
116 ENDDO
117
118 ELSE
119 C-- Read in Interior mask from file :
120
121 CALL READ_FLD_XY_RS( insideOBmaskFile, ' ', OBCS_insideMask,
122 & 0, myThid )
123
124 DO bj = myByLo(myThid), myByHi(myThid)
125 DO bi = myBxLo(myThid), myBxHi(myThid)
126 DO j=1,sNy
127 DO i=1,sNx
128 IF ( OBCS_insideMask(i,j,bi,bj).NE.0. )
129 & OBCS_insideMask(i,j,bi,bj) = 1.
130 ENDDO
131 ENDDO
132 ENDDO
133 ENDDO
134
135 C-- end computing/reading Interior mask
136 ENDIF
137
138 C-- Fill in the overlap:
139 _EXCH_XY_RS( OBCS_insideMask, myThid )
140
141 C== Set interior mask at U & V location (grid-cell Wester & Southern edges)
142 C leave OB edges inside (maskIn=1) (e.g., Eastern OB: maskInW(OB_Ie)=1 )
143 C so that velocity normal-component at OB is still in Interior region.
144 DO bj = myByLo(myThid), myByHi(myThid)
145 DO bi = myBxLo(myThid), myBxHi(myThid)
146 DO j=2-OLy,sNy+OLy
147 DO i=2-OLx,sNx+OLx
148 maskInW(i,j,bi,bj) = maskInW(i,j,bi,bj)
149 & *MAX( OBCS_insideMask(i-1,j,bi,bj),
150 & OBCS_insideMask(i,j,bi,bj) )
151 maskInS(i,j,bi,bj) = maskInS(i,j,bi,bj)
152 & *MAX( OBCS_insideMask(i,j-1,bi,bj),
153 & OBCS_insideMask(i,j,bi,bj) )
154 ENDDO
155 ENDDO
156 ENDDO
157 ENDDO
158
159 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
160
161 C-- After exchange, set tiled index arrays OB_Jn/Js/Ie/Iw in overlap region
162 C issue a warning if not consistent (similar to OBCS_CHECK but for overlap)
163
164 c IF ( .TRUE. ) THEN
165 IF ( OBCS_indexStatus .LT. 2 ) THEN
166 ioUnit = standardMessageUnit
167 WRITE(msgBuf,'(2A)')
168 & 'OBCS_INIT_FIXED: Setting OB indices in Overlap'
169 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
170 WRITE(errMsg,'(2A)') 'S/R OBCS_INIT_FIXED: ',
171 & 'Inside Mask and OB locations disagree :'
172 flag = .TRUE.
173 DO bj = myByLo(myThid), myByHi(myThid)
174 DO bi = myBxLo(myThid), myBxHi(myThid)
175
176 C- Eastern boundary
177 DO j=1-OLy,sNy+OLy
178 DO i=1,sNx+1
179 IF ( OBCS_insideMask(i,j,bi,bj).LT.
180 & OBCS_insideMask(i-1,j,bi,bj)
181 & .AND. ( j.LT.1 .OR. j.GT.sNy )
182 & .AND. kSurfW(i,j,bi,bj).LE.Nr ) THEN
183 IF ( OB_Ie(j,bi,bj).EQ.OB_indexNone ) THEN
184 OB_Ie(j,bi,bj) = i
185 WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
186 & ' Sets OBE(j,bi,bj=',j,',',bi,',',bj,')=', OB_Ie(j,bi,bj)
187 CALL PRINT_MESSAGE( msgBuf,ioUnit,SQUEEZE_RIGHT,myThid )
188 ELSEIF ( OB_Ie(j,bi,bj).NE.i ) THEN
189 IF ( flag ) CALL PRINT_ERROR( errMsg, myThid )
190 flag = .FALSE.
191 WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
192 & ' OBE(j,bi,bj=',j,',',bi,',',bj,')=', OB_Ie(j,bi,bj),
193 & ' but from insideMask expects I=', i
194 CALL PRINT_ERROR( msgBuf, myThid )
195 ENDIF
196 ENDIF
197 ENDDO
198 ENDDO
199 C- Western boundary
200 DO j=1-OLy,sNy+OLy
201 DO i=0,sNx
202 IF ( OBCS_insideMask(i,j,bi,bj).LT.
203 & OBCS_insideMask(i+1,j,bi,bj)
204 & .AND. ( j.LT.1 .OR. j.GT.sNy )
205 & .AND. kSurfW(i+1,j,bi,bj).LE.Nr ) THEN
206 IF ( OB_Iw(j,bi,bj).EQ.OB_indexNone ) THEN
207 OB_Iw(j,bi,bj) = i
208 WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
209 & ' Sets OBW(j,bi,bj=',j,',',bi,',',bj,')=', OB_Iw(j,bi,bj)
210 CALL PRINT_MESSAGE( msgBuf,ioUnit,SQUEEZE_RIGHT,myThid )
211 ELSEIF ( OB_Iw(j,bi,bj).NE.i ) THEN
212 IF ( flag ) CALL PRINT_ERROR( errMsg, myThid )
213 flag = .FALSE.
214 WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
215 & ' OBW(j,bi,bj=',j,',',bi,',',bj,')=', OB_Iw(j,bi,bj),
216 & ' but from insideMask expects I=', i
217 CALL PRINT_ERROR( msgBuf, myThid )
218 ENDIF
219 ENDIF
220 ENDDO
221 ENDDO
222 C- Northern boundary
223 DO j=1,sNy+1
224 DO i=1-OLx,sNx+OLx
225 IF ( OBCS_insideMask(i,j,bi,bj).LT.
226 & OBCS_insideMask(i,j-1,bi,bj)
227 & .AND. ( i.LT.1 .OR. i.GT.sNx )
228 & .AND. kSurfS(i,j,bi,bj).LE.Nr ) THEN
229 IF ( OB_Jn(i,bi,bj).EQ.OB_indexNone ) THEN
230 OB_Jn(i,bi,bj) = j
231 WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
232 & ' Sets OBN(i,bi,bj=',i,',',bi,',',bj,')=', OB_Jn(i,bi,bj)
233 CALL PRINT_MESSAGE( msgBuf,ioUnit,SQUEEZE_RIGHT,myThid )
234 ELSEIF ( OB_Jn(i,bi,bj).NE.j ) THEN
235 IF ( flag ) CALL PRINT_ERROR( errMsg, myThid )
236 flag = .FALSE.
237 WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
238 & ' OBN(i,bi,bj=',i,',',bi,',',bj,')=', OB_Jn(i,bi,bj),
239 & ' but from insideMask expects J=', j
240 CALL PRINT_ERROR( msgBuf, myThid )
241 ENDIF
242 ENDIF
243 ENDDO
244 ENDDO
245 C- Southern boundary
246 DO j=0,sNy
247 DO i=1-OLx,sNx+OLx
248 IF ( OBCS_insideMask(i,j,bi,bj).LT.
249 & OBCS_insideMask(i,j+1,bi,bj)
250 & .AND. ( i.LT.1 .OR. i.GT.sNx )
251 & .AND. kSurfS(i,j+1,bi,bj).LE.Nr ) THEN
252 IF ( OB_Js(i,bi,bj).EQ.OB_indexNone ) THEN
253 OB_Js(i,bi,bj) = j
254 WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
255 & ' Sets OBS(i,bi,bj=',i,',',bi,',',bj,')=', OB_Js(i,bi,bj)
256 CALL PRINT_MESSAGE( msgBuf,ioUnit,SQUEEZE_RIGHT,myThid )
257 ELSEIF ( OB_Js(i,bi,bj).NE.j ) THEN
258 IF ( flag ) CALL PRINT_ERROR( errMsg, myThid )
259 flag = .FALSE.
260 WRITE(msgBuf,'(A,I5,2(A,I3),2(A,I5))')
261 & ' OBS(i,bi,bj=',i,',',bi,',',bj,')=', OB_Js(i,bi,bj),
262 & ' but from insideMask expects J=', j
263 CALL PRINT_ERROR( msgBuf, myThid )
264 ENDIF
265 ENDIF
266 ENDDO
267 ENDDO
268
269 ENDDO
270 ENDDO
271 WRITE(msgBuf,'(2A)')
272 & 'OBCS_INIT_FIXED: Setting OB indices in Overlap <= done'
273 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
274 ENDIF
275
276 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
277
278 C== After EXCH: reset Interior mask to zero beyond OB: this is necessary
279 C when EXCH are not disabled (e.g. with EXCH1) between tile Edges
280 C that are closed by OB.
281 C Do it over OLx,OLy grid points beyond OB, in agreement with OBCS code
282 C (apply_tracer) which overwrites tracer over the same width.
283 OB_ApplX = OLx
284 OB_ApplY = OLy
285
286 DO bj = myByLo(myThid), myByHi(myThid)
287 DO bi = myBxLo(myThid), myBxHi(myThid)
288
289 C-- Set Interior mask to zero beyond OB
290 C- Eastern boundary
291 DO j=1-OLy,sNy+OLy
292 jm = MAX( j-1, 1-OLy )
293 iB = OB_Ie(j,bi,bj)
294 IF ( iB.NE.OB_indexNone ) THEN
295 DO i=iB,iB+OB_ApplX-1
296 OBCS_insideMask(i,j,bi,bj) = 0.
297 ENDDO
298 DO i=iB+1,iB+OB_ApplX-1
299 maskInW(i,j,bi,bj) = 0.
300 ENDDO
301 IF ( OB_Ie(jm,bi,bj).NE.OB_indexNone ) THEN
302 iB = MAX( iB, OB_Ie(jm,bi,bj) )
303 DO i=iB,iB+OB_ApplX-1
304 maskInS(i,j,bi,bj) = 0.
305 ENDDO
306 ENDIF
307 ENDIF
308 ENDDO
309 C- Western boundary
310 DO j=1-OLy,sNy+OLy
311 jm = MAX( j-1, 1-OLy )
312 iB = OB_Iw(j,bi,bj)
313 IF ( iB.NE.OB_indexNone ) THEN
314 DO i=1-OB_ApplX+iB,iB
315 OBCS_insideMask(i,j,bi,bj) = 0.
316 ENDDO
317 DO i=2-OB_ApplX+iB,iB
318 maskInW(i,j,bi,bj) = 0.
319 ENDDO
320 IF ( OB_Iw(jm,bi,bj).NE.OB_indexNone ) THEN
321 iB = MIN( iB, OB_Iw(jm,bi,bj) )
322 DO i=1-OB_ApplX+iB,iB
323 maskInS(i,j,bi,bj) = 0.
324 ENDDO
325 ENDIF
326 ENDIF
327 ENDDO
328 C- Northern boundary
329 DO i=1-OLx,sNx+OLx
330 im = MAX( i-1, 1-OLx )
331 jB = OB_Jn(i,bi,bj)
332 IF ( jB.NE.OB_indexNone ) THEN
333 DO j=jB,jB+OB_ApplY-1
334 OBCS_insideMask(i,j,bi,bj) = 0.
335 ENDDO
336 DO j=jB+1,jB+OB_ApplY-1
337 maskInS(i,j,bi,bj) = 0.
338 ENDDO
339 IF ( OB_Jn(im,bi,bj).NE.OB_indexNone ) THEN
340 jB = MAX( jB, OB_Jn(im,bi,bj) )
341 DO j=jB,jB+OB_ApplY-1
342 maskInW(i,j,bi,bj) = 0.
343 ENDDO
344 ENDIF
345 ENDIF
346 ENDDO
347 C- Southern boundary
348 DO i=1-OLx,sNx+OLx
349 im = MAX( i-1, 1-OLx )
350 jB = OB_Js(i,bi,bj)
351 IF ( jB.NE.OB_indexNone ) THEN
352 DO j=1-OB_ApplY+jB,jB
353 OBCS_insideMask(i,j,bi,bj) = 0.
354 ENDDO
355 DO j=2-OB_ApplY+jB,jB
356 maskInS(i,j,bi,bj) = 0.
357 ENDDO
358 IF ( OB_Js(im,bi,bj).NE.OB_indexNone ) THEN
359 jB = MIN( jB, OB_Js(im,bi,bj) )
360 DO j=1-OB_ApplY+jB,jB
361 maskInW(i,j,bi,bj) = 0.
362 ENDDO
363 ENDIF
364 ENDIF
365 ENDDO
366
367 C-- Apply mask to maskInC :
368 DO j=1-OLy,sNy+OLy
369 DO i=1-OLx,sNx+OLx
370 maskInC(i,j,bi,bj) = maskInC(i,j,bi,bj)
371 & *OBCS_insideMask(i,j,bi,bj)
372 ENDDO
373 ENDDO
374
375 C-- end bi,bj loops
376 ENDDO
377 ENDDO
378
379 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
380
381 C-- Set OB active tiles:
382 DO bj = myByLo(myThid), myByHi(myThid)
383 DO bi = myBxLo(myThid), myBxHi(myThid)
384 tileHasOBE(bi,bj) = .FALSE.
385 tileHasOBW(bi,bj) = .FALSE.
386 tileHasOBN(bi,bj) = .FALSE.
387 tileHasOBS(bi,bj) = .FALSE.
388 DO j=1-OLy,sNy+OLy
389 tileHasOBE(bi,bj) = tileHasOBE(bi,bj) .OR.
390 & ( OB_Ie(j,bi,bj).NE.OB_indexNone )
391 tileHasOBW(bi,bj) = tileHasOBW(bi,bj) .OR.
392 & ( OB_Iw(j,bi,bj).NE.OB_indexNone )
393 ENDDO
394 DO i=1-OLx,sNx+OLx
395 tileHasOBN(bi,bj) = tileHasOBN(bi,bj) .OR.
396 & ( OB_Jn(i,bi,bj).NE.OB_indexNone )
397 tileHasOBS(bi,bj) = tileHasOBS(bi,bj) .OR.
398 & ( OB_Js(i,bi,bj).NE.OB_indexNone )
399 ENDDO
400 ENDDO
401 ENDDO
402
403 C-- Set domain connected-piece Id for OB grid points:
404 CALL OBCS_SET_CONNECT( myThid )
405
406 #ifdef ALLOW_DEBUG
407 IF (debugMode) CALL DEBUG_LEAVE('OBCS_INIT_FIXED',myThid)
408 #endif
409
410 #endif /* ALLOW_OBCS */
411 RETURN
412 END

  ViewVC Help
Powered by ViewVC 1.1.22