model/src/ini_cartesian_grid.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_cartesian_grid.F,v 1.21 2009/01/27 15:35:27 jmc Exp $
C $Name:  $

c#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: INI_CARTESIAN_GRID
C     !INTERFACE:
      SUBROUTINE INI_CARTESIAN_GRID( myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_CARTESIAN_GRID
C     | o Initialise model coordinate system
C     *==========================================================*
C     | The grid arrays, initialised here, are used throughout
C     | the code in evaluating gradients, integrals and spatial
C     | avarages. This routine
C     | is called separately by each thread and initialises only
C     | the region of the domain it is "responsible" for.
C     | Notes:
C     | Two examples are included. One illustrates the
C     | initialisation of a cartesian grid (this routine).
C     | The other shows the
C     | inialisation of a spherical polar grid. Other orthonormal
C     | grids can be fitted into this design. In this case
C     | custom metric terms also need adding to account for the
C     | projections of velocity vectors onto these grids.
C     | The structure used here also makes it possible to
C     | implement less regular grid mappings. In particular
C     | o Schemes which leave out blocks of the domain that are
C     |   all land could be supported.
C     | o Multi-level schemes such as icosohedral or cubic
C     |   grid projections onto a sphere can also be fitted
C     |   within the strategy we use.
C     |   Both of the above also require modifying the support
C     |   routines that map computational blocks to simulation
C     |   domain blocks.
C     | Under the cartesian grid mode primitive distances in X
C     | and Y are in metres. Disktance in Z are in m or Pa
C     | depending on the vertical gridding mode.
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
c#ifdef ALLOW_EXCH2
c#include "W2_EXCH2_SIZE.h"
c#include "W2_EXCH2_TOPOLOGY.h"
c#include "W2_EXCH2_PARAMS.h"
c#endif /* ALLOW_EXCH2 */

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid  ::  Number of this instance of INI_CARTESIAN_GRID
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER iG, jG, bi, bj, i, j
      _RL xG0, yG0
C     "Long" real for temporary coordinate calculation
C      NOTICE the extended range of indices!!
      _RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
      _RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
C     These functions return the "global" index with valid values beyond
C     halo regions
      INTEGER iGl,jGl
      iGl(i,bi) = 1+MOD(myXGlobalLo-1+(bi-1)*sNx+i+Olx*Nx-1,Nx)
      jGl(j,bj) = 1+MOD(myYGlobalLo-1+(bj-1)*sNy+j+Oly*Ny-1,Ny)
c#ifdef ALLOW_EXCH2
c      INTEGER tN
c#endif /* ALLOW_EXCH2 */
CEOP

C     For each tile ...
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)

C--     "Global" index (place holder)
        jG = myYGlobalLo + (bj-1)*sNy
        iG = myXGlobalLo + (bi-1)*sNx
c#ifdef ALLOW_EXCH2
c        IF ( W2_useE2ioLayOut ) THEN
cC- note: does not work for non-uniform delX or delY
c          tN = W2_myTileList(bi,bj)
c          iG = exch2_txGlobalo(tN)
c          jG = exch2_tyGlobalo(tN)
c        ENDIF
c#endif /* ALLOW_EXCH2 */

C--   First find coordinate of tile corner (meaning outer corner of halo)
        xG0 = xgOrigin
C       Find the X-coordinate of the outer grid-line of the "real" tile
        DO i=1, iG-1
         xG0 = xG0 + delX(i)
        ENDDO
C       Back-step to the outer grid-line of the "halo" region
        DO i=1, Olx
         xG0 = xG0 - delX( 1+MOD(Olx*Nx-1+iG-i,Nx) )
        ENDDO
C       Find the Y-coordinate of the outer grid-line of the "real" tile
        yG0 = ygOrigin
        DO j=1, jG-1
         yG0 = yG0 + delY(j)
        ENDDO
C       Back-step to the outer grid-line of the "halo" region
        DO j=1, Oly
         yG0 = yG0 - delY( 1+MOD(Oly*Ny-1+jG-j,Ny) )
        ENDDO

C--     Calculate coordinates of cell corners for N+1 grid-lines
        DO j=1-Oly,sNy+Oly +1
         xGloc(1-Olx,j) = xG0
         DO i=1-Olx,sNx+Olx
c         xGloc(i+1,j) = xGloc(i,j) + delX(1+mod(Nx-1+iG-1+i,Nx))
          xGloc(i+1,j) = xGloc(i,j) + delX( iGl(i,bi) )
         ENDDO
        ENDDO
        DO i=1-Olx,sNx+Olx +1
         yGloc(i,1-Oly) = yG0
         DO j=1-Oly,sNy+Oly
c         yGloc(i,j+1) = yGloc(i,j) + delY(1+mod(Ny-1+jG-1+j,Ny))
          yGloc(i,j+1) = yGloc(i,j) + delY( jGl(j,bj) )
         ENDDO
        ENDDO

C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          xG(i,j,bi,bj) = xGloc(i,j)
          yG(i,j,bi,bj) = yGloc(i,j)
         ENDDO
        ENDDO

C--     Calculate [xC,yC], coordinates of cell centers
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
C         by averaging
          xC(i,j,bi,bj) = 0.25 _d 0*(
     &     xGloc(i,j)+xGloc(i+1,j)+xGloc(i,j+1)+xGloc(i+1,j+1) )
          yC(i,j,bi,bj) = 0.25 _d 0*(
     &     yGloc(i,j)+yGloc(i+1,j)+yGloc(i,j+1)+yGloc(i+1,j+1) )
         ENDDO
        ENDDO

C--     Calculate [dxF,dyF], lengths between cell faces (through center)
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          dxF(i,j,bi,bj) = delX( iGl(i,bi) )
          dyF(i,j,bi,bj) = delY( jGl(j,bj) )
         ENDDO
        ENDDO

C--     Calculate [dxG,dyG], lengths along cell boundaries
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          dxG(i,j,bi,bj) = delX( iGl(i,bi) )
          dyG(i,j,bi,bj) = delY( jGl(j,bj) )
         ENDDO
        ENDDO

C--     The following arrays are not defined in some parts of the halo
C       region. We set them to zero here for safety. If they are ever
C       referred to, especially in the denominator then it is a mistake!
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          dxC(i,j,bi,bj) = 0.
          dyC(i,j,bi,bj) = 0.
          dxV(i,j,bi,bj) = 0.
          dyU(i,j,bi,bj) = 0.
          rAw(i,j,bi,bj) = 0.
          rAs(i,j,bi,bj) = 0.
         ENDDO
        ENDDO

C--     Calculate [dxC], zonal length between cell centers
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx+1,sNx+Olx ! NOTE range
          dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
         ENDDO
        ENDDO

C--     Calculate [dyC], meridional length between cell centers
        DO j=1-Oly+1,sNy+Oly ! NOTE range
         DO i=1-Olx,sNx+Olx
          dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
         ENDDO
        ENDDO

C--     Calculate [dxV,dyU], length between velocity points (through corners)
        DO j=1-Oly+1,sNy+Oly ! NOTE range
         DO i=1-Olx+1,sNx+Olx ! NOTE range
C         by averaging (method I)
          dxV(i,j,bi,bj) = 0.5 _d 0*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
          dyU(i,j,bi,bj) = 0.5 _d 0*(dyG(i,j,bi,bj)+dyG(i,j-1,bi,bj))
C         by averaging (method II)
c         dxV(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
c         dyU(i,j,bi,bj) = 0.5*(dyC(i,j,bi,bj)+dyC(i-1,j,bi,bj))
         ENDDO
        ENDDO

C--     Calculate vertical face area
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          rA (i,j,bi,bj) = dxF(i,j,bi,bj)*dyF(i,j,bi,bj)
          rAw(i,j,bi,bj) = dxC(i,j,bi,bj)*dyG(i,j,bi,bj)
          rAs(i,j,bi,bj) = dxG(i,j,bi,bj)*dyC(i,j,bi,bj)
          rAz(i,j,bi,bj) = dxV(i,j,bi,bj)*dyU(i,j,bi,bj)
C--     Set trigonometric terms & grid orientation:
          tanPhiAtU(i,j,bi,bj) = 0.
          tanPhiAtV(i,j,bi,bj) = 0.
          angleCosC(i,j,bi,bj) = 1.
          angleSinC(i,j,bi,bj) = 0.
         ENDDO
        ENDDO

C--     Cosine(lat) scaling
        DO j=1-OLy,sNy+OLy
         cosFacU(j,bi,bj)=1.
         cosFacV(j,bi,bj)=1.
         sqcosFacU(j,bi,bj)=1.
         sqcosFacV(j,bi,bj)=1.
        ENDDO

C--   end bi,bj loops
       ENDDO
      ENDDO

C--   Set default (=whole domain) for where relaxation to climatology applies
      _BEGIN_MASTER(myThid)
      IF ( latBandClimRelax.EQ.UNSET_RL ) THEN
        latBandClimRelax = 0.
        DO j=1,Ny
          latBandClimRelax = latBandClimRelax + delY(j)
        ENDDO
        latBandClimRelax = latBandClimRelax*3. _d 0
      ENDIF
      _END_MASTER(myThid)

      RETURN
      END
1	jmc	1.22	C $Header: /u/gcmpack/MITgcm/model/src/ini_cartesian_grid.F,v 1.21 2009/01/27 15:35:27 jmc Exp $
2	adcroft	1.16	C $Name: $
3	cnh	1.1
4	jmc	1.22	c#include "PACKAGES_CONFIG.h"
5	cnh	1.10	#include "CPP_OPTIONS.h"
6	cnh	1.1
7	cnh	1.17	CBOP
8			C !ROUTINE: INI_CARTESIAN_GRID
9			C !INTERFACE:
10	cnh	1.1	SUBROUTINE INI_CARTESIAN_GRID( myThid )
11	jmc	1.22
12	cnh	1.17	C !DESCRIPTION: \bv
13			C ==========================================================
14	jmc	1.21	C \| SUBROUTINE INI_CARTESIAN_GRID
15			C \| o Initialise model coordinate system
16	cnh	1.17	C ==========================================================
17	jmc	1.21	C \| The grid arrays, initialised here, are used throughout
18			C \| the code in evaluating gradients, integrals and spatial
19			C \| avarages. This routine
20			C \| is called separately by each thread and initialises only
21			C \| the region of the domain it is "responsible" for.
22			C \| Notes:
23			C \| Two examples are included. One illustrates the
24			C \| initialisation of a cartesian grid (this routine).
25			C \| The other shows the
26			C \| inialisation of a spherical polar grid. Other orthonormal
27			C \| grids can be fitted into this design. In this case
28			C \| custom metric terms also need adding to account for the
29			C \| projections of velocity vectors onto these grids.
30			C \| The structure used here also makes it possible to
31			C \| implement less regular grid mappings. In particular
32			C \| o Schemes which leave out blocks of the domain that are
33			C \| all land could be supported.
34			C \| o Multi-level schemes such as icosohedral or cubic
35			C \| grid projections onto a sphere can also be fitted
36			C \| within the strategy we use.
37			C \| Both of the above also require modifying the support
38			C \| routines that map computational blocks to simulation
39			C \| domain blocks.
40			C \| Under the cartesian grid mode primitive distances in X
41			C \| and Y are in metres. Disktance in Z are in m or Pa
42			C \| depending on the vertical gridding mode.
43	cnh	1.17	C ==========================================================
44			C \ev
45
46			C !USES:
47	adcroft	1.12	IMPLICIT NONE
48	cnh	1.1	C === Global variables ===
49			#include "SIZE.h"
50			#include "EEPARAMS.h"
51			#include "PARAMS.h"
52			#include "GRID.h"
53	jmc	1.22	c#ifdef ALLOW_EXCH2
54			c#include "W2_EXCH2_SIZE.h"
55			c#include "W2_EXCH2_TOPOLOGY.h"
56			c#include "W2_EXCH2_PARAMS.h"
57			c#endif /* ALLOW_EXCH2 */
58	cnh	1.1
59	cnh	1.17	C !INPUT/OUTPUT PARAMETERS:
60	cnh	1.1	C == Routine arguments ==
61	jmc	1.21	C myThid :: Number of this instance of INI_CARTESIAN_GRID
62	cnh	1.1	INTEGER myThid
63
64	cnh	1.17	C !LOCAL VARIABLES:
65	cnh	1.1	C == Local variables ==
66	jmc	1.22	INTEGER iG, jG, bi, bj, i, j
67	adcroft	1.16	_RL xG0, yG0
68	cnh	1.17	C "Long" real for temporary coordinate calculation
69			C NOTICE the extended range of indices!!
70	adcroft	1.16	_RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
71			_RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
72	cnh	1.17	C These functions return the "global" index with valid values beyond
73			C halo regions
74	adcroft	1.16	INTEGER iGl,jGl
75	jmc	1.22	iGl(i,bi) = 1+MOD(myXGlobalLo-1+(bi-1)sNx+i+OlxNx-1,Nx)
76			jGl(j,bj) = 1+MOD(myYGlobalLo-1+(bj-1)sNy+j+OlyNy-1,Ny)
77			c#ifdef ALLOW_EXCH2
78			c INTEGER tN
79			c#endif /* ALLOW_EXCH2 */
80	cnh	1.17	CEOP
81	adcroft	1.16
82			C For each tile ...
83	cnh	1.1	DO bj = myByLo(myThid), myByHi(myThid)
84			DO bi = myBxLo(myThid), myBxHi(myThid)
85	adcroft	1.16
86			C-- "Global" index (place holder)
87			jG = myYGlobalLo + (bj-1)*sNy
88	cnh	1.1	iG = myXGlobalLo + (bi-1)*sNx
89	jmc	1.22	c#ifdef ALLOW_EXCH2
90			c IF ( W2_useE2ioLayOut ) THEN
91			cC- note: does not work for non-uniform delX or delY
92			c tN = W2_myTileList(bi,bj)
93			c iG = exch2_txGlobalo(tN)
94			c jG = exch2_tyGlobalo(tN)
95			c ENDIF
96			c#endif /* ALLOW_EXCH2 */
97	cnh	1.1
98	adcroft	1.16	C-- First find coordinate of tile corner (meaning outer corner of halo)
99	jmc	1.21	xG0 = xgOrigin
100	adcroft	1.16	C Find the X-coordinate of the outer grid-line of the "real" tile
101			DO i=1, iG-1
102			xG0 = xG0 + delX(i)
103			ENDDO
104			C Back-step to the outer grid-line of the "halo" region
105			DO i=1, Olx
106	jmc	1.22	xG0 = xG0 - delX( 1+MOD(Olx*Nx-1+iG-i,Nx) )
107	adcroft	1.16	ENDDO
108			C Find the Y-coordinate of the outer grid-line of the "real" tile
109	jmc	1.21	yG0 = ygOrigin
110	adcroft	1.16	DO j=1, jG-1
111			yG0 = yG0 + delY(j)
112			ENDDO
113			C Back-step to the outer grid-line of the "halo" region
114			DO j=1, Oly
115	jmc	1.22	yG0 = yG0 - delY( 1+MOD(Oly*Ny-1+jG-j,Ny) )
116	adcroft	1.16	ENDDO
117
118			C-- Calculate coordinates of cell corners for N+1 grid-lines
119	jmc	1.22	DO j=1-Oly,sNy+Oly +1
120			xGloc(1-Olx,j) = xG0
121			DO i=1-Olx,sNx+Olx
122			c xGloc(i+1,j) = xGloc(i,j) + delX(1+mod(Nx-1+iG-1+i,Nx))
123			xGloc(i+1,j) = xGloc(i,j) + delX( iGl(i,bi) )
124	adcroft	1.16	ENDDO
125			ENDDO
126	jmc	1.22	DO i=1-Olx,sNx+Olx +1
127			yGloc(i,1-Oly) = yG0
128			DO j=1-Oly,sNy+Oly
129			c yGloc(i,j+1) = yGloc(i,j) + delY(1+mod(Ny-1+jG-1+j,Ny))
130			yGloc(i,j+1) = yGloc(i,j) + delY( jGl(j,bj) )
131	adcroft	1.16	ENDDO
132			ENDDO
133
134			C-- Make a permanent copy of [xGloc,yGloc] in [xG,yG]
135	jmc	1.22	DO j=1-Oly,sNy+Oly
136			DO i=1-Olx,sNx+Olx
137			xG(i,j,bi,bj) = xGloc(i,j)
138			yG(i,j,bi,bj) = yGloc(i,j)
139	adcroft	1.16	ENDDO
140			ENDDO
141
142			C-- Calculate [xC,yC], coordinates of cell centers
143	jmc	1.22	DO j=1-Oly,sNy+Oly
144			DO i=1-Olx,sNx+Olx
145	adcroft	1.16	C by averaging
146	jmc	1.22	xC(i,j,bi,bj) = 0.25 _d 0*(
147			& xGloc(i,j)+xGloc(i+1,j)+xGloc(i,j+1)+xGloc(i+1,j+1) )
148			yC(i,j,bi,bj) = 0.25 _d 0*(
149			& yGloc(i,j)+yGloc(i+1,j)+yGloc(i,j+1)+yGloc(i+1,j+1) )
150	adcroft	1.16	ENDDO
151			ENDDO
152
153			C-- Calculate [dxF,dyF], lengths between cell faces (through center)
154	jmc	1.22	DO j=1-Oly,sNy+Oly
155			DO i=1-Olx,sNx+Olx
156			dxF(i,j,bi,bj) = delX( iGl(i,bi) )
157			dyF(i,j,bi,bj) = delY( jGl(j,bj) )
158	adcroft	1.16	ENDDO
159			ENDDO
160
161			C-- Calculate [dxG,dyG], lengths along cell boundaries
162	jmc	1.22	DO j=1-Oly,sNy+Oly
163			DO i=1-Olx,sNx+Olx
164			dxG(i,j,bi,bj) = delX( iGl(i,bi) )
165			dyG(i,j,bi,bj) = delY( jGl(j,bj) )
166	adcroft	1.16	ENDDO
167			ENDDO
168
169			C-- The following arrays are not defined in some parts of the halo
170			C region. We set them to zero here for safety. If they are ever
171			C referred to, especially in the denominator then it is a mistake!
172	jmc	1.22	DO j=1-Oly,sNy+Oly
173			DO i=1-Olx,sNx+Olx
174			dxC(i,j,bi,bj) = 0.
175			dyC(i,j,bi,bj) = 0.
176			dxV(i,j,bi,bj) = 0.
177			dyU(i,j,bi,bj) = 0.
178			rAw(i,j,bi,bj) = 0.
179			rAs(i,j,bi,bj) = 0.
180	adcroft	1.16	ENDDO
181			ENDDO
182
183			C-- Calculate [dxC], zonal length between cell centers
184	jmc	1.22	DO j=1-Oly,sNy+Oly
185			DO i=1-Olx+1,sNx+Olx ! NOTE range
186			dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
187	cnh	1.1	ENDDO
188			ENDDO
189	adcroft	1.16
190			C-- Calculate [dyC], meridional length between cell centers
191	jmc	1.22	DO j=1-Oly+1,sNy+Oly ! NOTE range
192			DO i=1-Olx,sNx+Olx
193			dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
194	cnh	1.1	ENDDO
195			ENDDO
196	adcroft	1.16
197			C-- Calculate [dxV,dyU], length between velocity points (through corners)
198	jmc	1.22	DO j=1-Oly+1,sNy+Oly ! NOTE range
199			DO i=1-Olx+1,sNx+Olx ! NOTE range
200	adcroft	1.16	C by averaging (method I)
201	jmc	1.22	dxV(i,j,bi,bj) = 0.5 _d 0*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
202			dyU(i,j,bi,bj) = 0.5 _d 0*(dyG(i,j,bi,bj)+dyG(i,j-1,bi,bj))
203	adcroft	1.16	C by averaging (method II)
204	jmc	1.22	c dxV(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
205			c dyU(i,j,bi,bj) = 0.5*(dyC(i,j,bi,bj)+dyC(i-1,j,bi,bj))
206	cnh	1.1	ENDDO
207			ENDDO
208	adcroft	1.16
209	jmc	1.21	C-- Calculate vertical face area
210	jmc	1.22	DO j=1-Oly,sNy+Oly
211			DO i=1-Olx,sNx+Olx
212			rA (i,j,bi,bj) = dxF(i,j,bi,bj)*dyF(i,j,bi,bj)
213			rAw(i,j,bi,bj) = dxC(i,j,bi,bj)*dyG(i,j,bi,bj)
214			rAs(i,j,bi,bj) = dxG(i,j,bi,bj)*dyC(i,j,bi,bj)
215			rAz(i,j,bi,bj) = dxV(i,j,bi,bj)*dyU(i,j,bi,bj)
216	jmc	1.18	C-- Set trigonometric terms & grid orientation:
217	jmc	1.22	tanPhiAtU(i,j,bi,bj) = 0.
218			tanPhiAtV(i,j,bi,bj) = 0.
219			angleCosC(i,j,bi,bj) = 1.
220			angleSinC(i,j,bi,bj) = 0.
221	cnh	1.6	ENDDO
222			ENDDO
223	cnh	1.1
224	adcroft	1.16	C-- Cosine(lat) scaling
225	jmc	1.22	DO j=1-OLy,sNy+OLy
226			cosFacU(j,bi,bj)=1.
227			cosFacV(j,bi,bj)=1.
228			sqcosFacU(j,bi,bj)=1.
229			sqcosFacV(j,bi,bj)=1.
230	adcroft	1.16	ENDDO
231
232	jmc	1.21	C-- end bi,bj loops
233			ENDDO
234			ENDDO
235	adcroft	1.16
236	jmc	1.19	C-- Set default (=whole domain) for where relaxation to climatology applies
237	jmc	1.20	_BEGIN_MASTER(myThid)
238	jmc	1.19	IF ( latBandClimRelax.EQ.UNSET_RL ) THEN
239			latBandClimRelax = 0.
240			DO j=1,Ny
241			latBandClimRelax = latBandClimRelax + delY(j)
242			ENDDO
243			latBandClimRelax = latBandClimRelax*3. _d 0
244	jmc	1.21	ENDIF
245	jmc	1.20	_END_MASTER(myThid)
246	jmc	1.19
247	cnh	1.1	RETURN
248			END