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

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

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


Revision 1.31 - (show annotations) (download)
Fri Mar 4 04:53:42 2011 UTC (13 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62w
Changes since 1.30: +3 -1 lines
add insideOBmaskFile (to specify interior OB mask) in 1rst namelist.

1 C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_readparms.F,v 1.30 2011/02/28 15:22:09 jmc Exp $
2 C $Name: $
3
4 #include "OBCS_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: OBCS_READPARMS
8 C !INTERFACE:
9 SUBROUTINE OBCS_READPARMS( myThid )
10
11 C !DESCRIPTION: \bv
12 C *==========================================================*
13 C | SUBROUTINE OBCS_READPARMS
14 C | o Routine to initialize OBCS variables and constants.
15 C *==========================================================*
16 C \ev
17
18 C !USES:
19 IMPLICIT NONE
20 C === Global variables ===
21 #include "SIZE.h"
22 #include "EEPARAMS.h"
23 #include "PARAMS.h"
24 #include "OBCS.h"
25 #ifdef ALLOW_ORLANSKI
26 #include "ORLANSKI.h"
27 #endif
28 #ifdef ALLOW_PTRACERS
29 #include "PTRACERS_SIZE.h"
30 #include "OBCS_PTRACERS.h"
31 #endif /* ALLOW_PTRACERS */
32 #ifdef ALLOW_EXCH2
33 #include "W2_EXCH2_SIZE.h"
34 #include "W2_EXCH2_TOPOLOGY.h"
35 #include "W2_EXCH2_PARAMS.h"
36 #endif /* ALLOW_EXCH2 */
37
38 C !INPUT/OUTPUT PARAMETERS:
39 C === Routine arguments ===
40 INTEGER myThid
41
42 #ifdef ALLOW_OBCS
43
44 C !LOCAL VARIABLES:
45 C === Local variables ===
46 C msgBuf :: Informational/error message buffer
47 C iUnit :: Work variable for IO unit number
48 CHARACTER*(MAX_LEN_MBUF) msgBuf
49 INTEGER iUnit
50 INTEGER I, J
51 INTEGER bi, bj, iG, jG, iGm, jGm
52 #ifdef ALLOW_PTRACERS
53 INTEGER iTracer
54 #endif
55 #ifdef ALLOW_EXCH2
56 INTEGER tN
57 #endif /* ALLOW_EXCH2 */
58
59 C These are input arrays (of integers) that contain the *absolute*
60 C computational index of an open-boundary (OB) point.
61 C A zero (0) element means there is no corresponding OB in that column/row.
62 C The computational coordinate refers to "tracer" cells.
63 C For a northern/southern OB, the OB V point is to the south/north.
64 C For an eastern/western OB, the OB U point is to the west/east.
65 C eg.
66 C OB_Jnorth(3)=34 means that:
67 C T( 3 ,34) is a an OB point
68 C U(3:4,34) is a an OB point
69 C V( 4 ,34) is a an OB point
70 C while
71 C OB_Jsouth(3)=1 means that:
72 C T( 3 ,1) is a an OB point
73 C U(3:4,1) is a an OB point
74 C V( 4 ,2) is a an OB point
75 C
76 C For convenience, negative values for Jnorth/Ieast refer to
77 C points relative to the Northern/Eastern edges of the model
78 C eg. OB_Jnorth(3)=-1 means that the point (3,Ny) is a northern O-B.
79 C
80 C With exch2, the global domain used for specifying the boundary (and
81 C boundary value files) is different for N,S and E,W boundaries:
82 C - for N,S, the facets are stacked in x (like W2_mapIO=-1)
83 C - for E,W, the facets are stacked in y, so that E,W boundaries in
84 C different facets cannot have the same I
85 C
86 C OB_Jnorth(W2_maxXStackNx) :: global index array of northern open-boundary point
87 C OB_Jsouth(W2_maxXStackNx) :: global index array of southern open-boundary point
88 C OB_Ieast(W2_maxYStackNy) :: global index array of eastern open-boundary point
89 C OB_Iwest(W2_maxYStackNy) :: global index array of western open-boundary point
90
91 COMMON/OBCS_GLOBAL/ OB_Jnorth, OB_Jsouth, OB_Ieast, OB_Iwest
92 #ifdef ALLOW_EXCH2
93 INTEGER OB_Jnorth(W2_maxXStackNx)
94 INTEGER OB_Jsouth(W2_maxXStackNx)
95 INTEGER OB_Ieast(W2_maxYStackNy)
96 INTEGER OB_Iwest(W2_maxYStackNy)
97 #else
98 INTEGER OB_Jnorth(Nx)
99 INTEGER OB_Jsouth(Nx)
100 INTEGER OB_Ieast(Ny)
101 INTEGER OB_Iwest(Ny)
102 #endif
103
104 C With exch2, we use different global domains for specifying
105 C N,S resp. E,W boundaries (and for reading in the corresponding data):
106 C
107 C OBNS_Nx :: width of global domain for OB_Jnorth, OB_Jsouth
108 C OBNS_Ny :: height of global domain for OB_Jnorth, OB_Jsouth
109 C OBEW_Nx :: width of global domain for OB_Ieast, OB_Iwest
110 C OBEW_Ny :: height of global domain for OB_Ieast, OB_Iwest
111
112 INTEGER OBNS_Nx, OBNS_Ny
113 INTEGER OBEW_Nx, OBEW_Ny
114
115 #ifdef ALLOW_EXCH2
116 C buf :: used to exchange OB_Jnorth, ...
117 _RS buf(1-Olx:sNx+Olx,1-Oly:sNy+Oly,nSx,nSy)
118 #endif
119 CEOP
120
121 NAMELIST /OBCS_PARM01/
122 & insideOBmaskFile,
123 & OB_Jnorth,OB_Jsouth,OB_Ieast,OB_Iwest,
124 & useOrlanskiNorth,useOrlanskiSouth,
125 & useOrlanskiEast,useOrlanskiWest,
126 & useStevensNorth,useStevensSouth,
127 & useStevensEast,useStevensWest,
128 & OBNuFile,OBNvFile,OBNtFile,OBNsFile,OBNaFile,OBNhFile,
129 & OBSuFile,OBSvFile,OBStFile,OBSsFile,OBSaFile,OBShFile,
130 & OBEuFile,OBEvFile,OBEtFile,OBEsFile,OBEaFile,OBEhFile,
131 & OBWuFile,OBWvFile,OBWtFile,OBWsFile,OBWaFile,OBWhFile,
132 & OBNslFile,OBSslFile,OBEslFile,OBWslFile,
133 & OBNsnFile,OBSsnFile,OBEsnFile,OBWsnFile,
134 & OBNuiceFile,OBSuiceFile,OBEuiceFile,OBWuiceFile,
135 & OBNviceFile,OBSviceFile,OBEviceFile,OBWviceFile,
136 & OBNetaFile, OBSetaFile, OBEetaFile, OBWetaFile,
137 & OBNwFile, OBSwFile, OBEwFile, OBWwFile,
138 #ifdef ALLOW_PTRACERS
139 & OBNptrFile,OBSptrFile,OBEptrFile,OBWptrFile,
140 #endif
141 & useOBCSsponge, useOBCSbalance, useOBCSprescribe,
142 & OBCS_balanceFacN, OBCS_balanceFacS,
143 & OBCS_balanceFacE, OBCS_balanceFacW,
144 & useOBCSYearlyFields, OBCSfixTopo,
145 & OBCS_monitorFreq, OBCS_monSelect, OBCSprintDiags
146
147 #ifdef ALLOW_ORLANSKI
148 NAMELIST /OBCS_PARM02/
149 & CMAX, cvelTimeScale, CFIX, useFixedCEast, useFixedCWest
150 #endif
151
152 #ifdef ALLOW_OBCS_SPONGE
153 NAMELIST /OBCS_PARM03/
154 & Urelaxobcsinner,Urelaxobcsbound,
155 & Vrelaxobcsinner,Vrelaxobcsbound,
156 & spongeThickness
157 #endif
158 #ifdef ALLOW_OBCS_STEVENS
159 NAMELIST /OBCS_PARM04/
160 & TrelaxStevens,SrelaxStevens,
161 & useStevensPhaseVel,useStevensAdvection
162 #endif /* ALLOW_OBCS_STEVENS */
163
164 _BEGIN_MASTER(myThid)
165
166 #ifdef ALLOW_EXCH2
167 OBNS_Nx = exch2_xStack_Nx
168 OBNS_Ny = exch2_xStack_Ny
169 OBEW_Nx = exch2_yStack_Nx
170 OBEW_Ny = exch2_yStack_Ny
171 #else
172 OBNS_Nx = Nx
173 OBNS_Ny = Ny
174 OBEW_Nx = Nx
175 OBEW_Ny = Ny
176 #endif
177
178 C-- Default flags and values for OBCS
179 insideOBmaskFile = ' '
180 DO I=1,OBNS_Nx
181 OB_Jnorth(I)=0
182 OB_Jsouth(I)=0
183 ENDDO
184 DO J=1,OBEW_Ny
185 OB_Ieast(J)=0
186 OB_Iwest(J)=0
187 ENDDO
188 useOrlanskiNorth =.FALSE.
189 useOrlanskiSouth =.FALSE.
190 useOrlanskiEast =.FALSE.
191 useOrlanskiWest =.FALSE.
192 useStevensNorth =.FALSE.
193 useStevensSouth =.FALSE.
194 useStevensEast =.FALSE.
195 useStevensWest =.FALSE.
196 useStevensPhaseVel =.TRUE.
197 useStevensAdvection=.TRUE.
198 useOBCSsponge =.FALSE.
199 useOBCSbalance =.FALSE.
200 OBCS_balanceFacN = 1. _d 0
201 OBCS_balanceFacS = 1. _d 0
202 OBCS_balanceFacE = 1. _d 0
203 OBCS_balanceFacW = 1. _d 0
204 useOBCSprescribe =.FALSE.
205 useOBCSYearlyFields=.FALSE.
206 OBCSfixTopo =.TRUE.
207 OBCS_monitorFreq = monitorFreq
208 OBCS_monSelect = 0
209 OBCSprintDiags = debugLevel.GE.debLevB
210
211 OBNuFile = ' '
212 OBNvFile = ' '
213 OBNtFile = ' '
214 OBNsFile = ' '
215 OBNaFile = ' '
216 OBNslFile = ' '
217 OBNsnFile = ' '
218 OBNuiceFile = ' '
219 OBNviceFile = ' '
220 OBNhFile = ' '
221 OBSuFile = ' '
222 OBSvFile = ' '
223 OBStFile = ' '
224 OBSsFile = ' '
225 OBSaFile = ' '
226 OBShFile = ' '
227 OBSslFile = ' '
228 OBSsnFile = ' '
229 OBSuiceFile = ' '
230 OBSviceFile = ' '
231 OBEuFile = ' '
232 OBEvFile = ' '
233 OBEtFile = ' '
234 OBEsFile = ' '
235 OBEaFile = ' '
236 OBEhFile = ' '
237 OBEslFile = ' '
238 OBEsnFile = ' '
239 OBEuiceFile = ' '
240 OBEviceFile = ' '
241 OBWuFile = ' '
242 OBWvFile = ' '
243 OBWtFile = ' '
244 OBWsFile = ' '
245 OBWaFile = ' '
246 OBWhFile = ' '
247 OBWslFile = ' '
248 OBWsnFile = ' '
249 OBWuiceFile = ' '
250 OBWviceFile = ' '
251 OBNetaFile = ' '
252 OBSetaFile = ' '
253 OBEetaFile = ' '
254 OBWetaFile = ' '
255 OBNwFile = ' '
256 OBSwFile = ' '
257 OBEwFile = ' '
258 OBWwFile = ' '
259 #ifdef ALLOW_PTRACERS
260 DO iTracer = 1, PTRACERS_num
261 OBNptrFile(iTracer) = ' '
262 OBSptrFile(iTracer) = ' '
263 OBEptrFile(iTracer) = ' '
264 OBWptrFile(iTracer) = ' '
265 ENDDO
266 #endif
267
268 C Open and read the data.obcs file
269 WRITE(msgBuf,'(A)') ' OBCS_READPARMS: opening data.obcs'
270 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
271 & SQUEEZE_RIGHT , myThid )
272 CALL OPEN_COPY_DATA_FILE(
273 I 'data.obcs', 'OBCS_READPARMS',
274 O iUnit,
275 I myThid )
276
277 C-- Read parameters from open data file
278 READ(UNIT=iUnit,NML=OBCS_PARM01)
279
280 #ifdef ALLOW_ORLANSKI
281 C Default Orlanski radiation parameters
282 CMAX = 0.45 _d 0 /* maximum allowable phase speed-CFL for AB-II */
283 cvelTimeScale = 2000.0 _d 0 /* Averaging period for phase speed in sec. */
284 CFIX = 0.8 _d 0 /* Fixed boundary phase speed in m/s */
285 useFixedCEast=.FALSE.
286 useFixedCWest=.FALSE.
287 IF (useOrlanskiNorth.OR.
288 & useOrlanskiSouth.OR.
289 & useOrlanskiEast.OR.
290 & useOrlanskiWest)
291 & READ(UNIT=iUnit,NML=OBCS_PARM02)
292 #endif
293
294 #ifdef ALLOW_OBCS_SPONGE
295 C Default sponge layer parameters:
296 C sponge layer is turned off by default
297 spongeThickness = 0
298 Urelaxobcsinner = 0. _d 0
299 Urelaxobcsbound = 0. _d 0
300 Vrelaxobcsinner = 0. _d 0
301 Vrelaxobcsbound = 0. _d 0
302 CML this was the previous default in units of days
303 CML spongeThickness = 2
304 CML Urelaxobcsinner = 5. _d 0
305 CML Urelaxobcsbound = 1. _d 0
306 CML Vrelaxobcsinner = 5. _d 0
307 CML Vrelaxobcsbound = 1. _d 0
308 IF (useOBCSsponge)
309 & READ(UNIT=iUnit,NML=OBCS_PARM03)
310 #endif
311 #ifdef ALLOW_OBCS_STEVENS
312 TrelaxStevens = 0. _d 0
313 SrelaxStevens = 0. _d 0
314 IF ( useStevensNorth .OR. useStevensSouth
315 & .OR. useStevensEast .OR. useStevensWest )
316 & READ(UNIT=iUnit,NML=OBCS_PARM04)
317 #endif
318
319 WRITE(msgBuf,'(A)') ' OBCS_READPARMS: finished reading data.obcs'
320 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
321 & SQUEEZE_RIGHT , myThid )
322
323 C-- Close the open data file
324 CLOSE(iUnit)
325
326 C- Account for periodicity if negative indices were supplied
327 DO J=1,OBEW_Ny
328 IF (OB_Ieast(J).LT.0) OB_Ieast(J)=OB_Ieast(J)+OBEW_Nx+1
329 ENDDO
330 DO I=1,OBNS_Nx
331 IF (OB_Jnorth(I).LT.0) OB_Jnorth(I)=OB_Jnorth(I)+OBNS_Ny+1
332 ENDDO
333 IF ( debugLevel.GE.debLevA ) THEN
334 c write(*,*) 'OB Jn =',OB_Jnorth
335 c write(*,*) 'OB Js =',OB_Jsouth
336 c write(*,*) 'OB Ie =',OB_Ieast
337 c write(*,*) 'OB Iw =',OB_Iwest
338 WRITE(msgBuf,'(A)') ' Northern OB global indices : OB_Jnorth ='
339 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
340 & SQUEEZE_RIGHT, myThid )
341 CALL PRINT_LIST_I( OB_Jnorth, 1, OBNS_Nx, INDEX_I,
342 & .FALSE., .TRUE., standardMessageUnit )
343 WRITE(msgBuf,'(A)') ' Southern OB global indices : OB_Jsouth ='
344 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
345 & SQUEEZE_RIGHT, myThid )
346 CALL PRINT_LIST_I( OB_Jsouth, 1, OBNS_Nx, INDEX_I,
347 & .FALSE., .TRUE., standardMessageUnit )
348 WRITE(msgBuf,'(A)') ' Eastern OB global indices : OB_Ieast ='
349 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
350 & SQUEEZE_RIGHT, myThid )
351 CALL PRINT_LIST_I( OB_Ieast, 1, OBEW_Ny, INDEX_J,
352 & .FALSE., .TRUE., standardMessageUnit )
353 WRITE(msgBuf,'(A)') ' Western OB global indices : OB_Iwest ='
354 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
355 & SQUEEZE_RIGHT, myThid )
356 CALL PRINT_LIST_I( OB_Iwest, 1, OBEW_Ny, INDEX_J,
357 & .FALSE., .TRUE., standardMessageUnit )
358 WRITE(msgBuf,'(A)') ' '
359 CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
360 & SQUEEZE_RIGHT, myThid )
361 ENDIF
362
363 _END_MASTER(myThid)
364 C-- Everyone else must wait for the parameters to be loaded
365 _BARRIER
366
367 C-- Calculate the tiled index arrays OB_Jn/Js/Ie/Iw here from the
368 C global arrays OB_Jnorth/Jsouth/Ieast/Iwest.
369 C Note: This part of the code has been moved from obcs_init_fixed to
370 C routine routine because the OB_Jn/Js/Ie/Iw index arrays are
371 C required by ini_depth which is called before obcs_init_fixed
372 DO bj = myByLo(myThid), myByHi(myThid)
373 DO bi = myBxLo(myThid), myBxHi(myThid)
374
375 DO I=1-Olx,sNx+Olx
376 OB_Jn(I,bi,bj)=0
377 OB_Js(I,bi,bj)=0
378 ENDDO
379
380 DO J=1-Oly,sNy+Oly
381 OB_Ie(J,bi,bj)=0
382 OB_Iw(J,bi,bj)=0
383 ENDDO
384
385 #ifdef ALLOW_EXCH2
386 C We apply OBCS only inside tile and exchange overlaps later
387 tN = W2_myTileList(bi,bj)
388 C 1. N/S boundaries
389 DO J=1,sNy
390 C convert from local y index J to global y index jG
391 c for N/S boundaries, we use faces stacked in x direction
392 jG = exch2_tyXStackLo(tN)+J-1
393 C loop over local x index I
394 DO I=1,sNx
395 iG = exch2_txXStackLo(tN)+I-1
396 IF (jG.EQ.OB_Jnorth(iG)) OB_Jn(I,bi,bj)=J
397 IF (jG.EQ.OB_Jsouth(iG)) OB_Js(I,bi,bj)=J
398 ENDDO
399 ENDDO
400 C 2. E/W boundaries
401 DO J=1,sNy
402 C convert from local y index J to global y index jG
403 c for E/W boundaries, we use faces stacked in y direction
404 jG = exch2_tyYStackLo(tN)+J-1
405 C loop over local x index I
406 DO I=1,sNx
407 iG = exch2_txYStackLo(tN)+I-1
408 IF (iG.EQ.OB_Ieast(jG)) OB_Ie(J,bi,bj)=I
409 IF (iG.EQ.OB_Iwest(jG)) OB_Iw(J,bi,bj)=I
410 ENDDO
411 ENDDO
412
413 #else /* ALLOW_EXCH2 */
414
415 DO J=1-Oly,sNy+Oly
416 C convert from local y index J to global y index jG
417 jG = myYGlobalLo-1+(bj-1)*sNy+J
418 C use periodicity to deal with out of range points caused by the overlaps.
419 C they will be excluded by the mask in any case, but this saves array
420 C out-of-bounds errors here.
421 jGm = 1+mod( jG-1+Ny , Ny )
422 C loop over local x index I
423 DO I=1,sNx
424 iG = myXGlobalLo-1+(bi-1)*sNx+I
425 iGm = 1+mod( iG-1+Nx , Nx )
426 C OB_Ieast(jGm) allows for the eastern boundary to be at variable x locations
427 IF (iG.EQ.OB_Ieast(jGm)) OB_Ie(J,bi,bj)=I
428 IF (iG.EQ.OB_Iwest(jGm)) OB_Iw(J,bi,bj)=I
429 ENDDO
430 ENDDO
431 DO J=1,sNy
432 jG = myYGlobalLo-1+(bj-1)*sNy+J
433 jGm = 1+mod( jG-1+Ny , Ny )
434 DO I=1-Olx,sNx+Olx
435 iG = myXGlobalLo-1+(bi-1)*sNx+I
436 iGm = 1+mod( iG-1+Nx , Nx )
437 C OB_Jnorth(iGm) allows for the northern boundary to be at variable y locations
438 IF (jG.EQ.OB_Jnorth(iGm)) OB_Jn(I,bi,bj)=J
439 IF (jG.EQ.OB_Jsouth(iGm)) OB_Js(I,bi,bj)=J
440 ENDDO
441 ENDDO
442 #endif /* ALLOW_EXCH2 */
443
444 C bi,bj-loops
445 ENDDO
446 ENDDO
447
448 #ifdef ALLOW_EXCH2
449 C exchange with neighbors
450 DO bj = myByLo(myThid), myByHi(myThid)
451 DO bi = myBxLo(myThid), myBxHi(myThid)
452 DO J=1,sNy
453 buf(sNx,J,bi,bj) = OB_Ie(J,bi,bj)
454 buf( 1,J,bi,bj) = OB_Iw(J,bi,bj)
455 ENDDO
456 ENDDO
457 ENDDO
458 CALL EXCH_3D_RS( buf, 1, myThid )
459 DO bj = myByLo(myThid), myByHi(myThid)
460 DO bi = myBxLo(myThid), myBxHi(myThid)
461 DO J=1-Oly,sNy+Oly
462 OB_Ie(J,bi,bj) = buf(sNx,J,bi,bj)
463 OB_Iw(J,bi,bj) = buf( 1,J,bi,bj)
464 ENDDO
465 ENDDO
466 ENDDO
467
468 DO bj = myByLo(myThid), myByHi(myThid)
469 DO bi = myBxLo(myThid), myBxHi(myThid)
470 DO I=1,sNx
471 buf(I,sNy,bi,bj) = OB_Jn(I,bi,bj)
472 buf(I, 1,bi,bj) = OB_Js(I,bi,bj)
473 ENDDO
474 ENDDO
475 ENDDO
476 CALL EXCH_3D_RS( buf, 1, myThid )
477 DO bj = myByLo(myThid), myByHi(myThid)
478 DO bi = myBxLo(myThid), myBxHi(myThid)
479 DO I=1-Olx,sNx+Olx
480 OB_Jn(I,bi,bj) = buf(I,sNy,bi,bj)
481 OB_Js(I,bi,bj) = buf(I, 1,bi,bj)
482 ENDDO
483 ENDDO
484 ENDDO
485 #endif /* ALLOW_EXCH2 */
486
487 #endif /* ALLOW_OBCS */
488 RETURN
489 END

  ViewVC Help
Powered by ViewVC 1.1.22