model/src/ini_cartesian_grid.F

C $Header: /u/gcmpack/MITgcm/model/src/ini_cartesian_grid.F,v 1.18 2005/07/13 00:36:01 jmc Exp $
C $Name:  $

#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     !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)
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--   First find coordinate of tile corner (meaning outer corner of halo)
        xG0 = 0.
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 = 0.
        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*(
     &     xGloc(I,J)+xGloc(I+1,J)+xGloc(I,J+1)+xGloc(I+1,J+1) )
          yC(I,J,bi,bj) = 0.25*(
     &     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.5D0*(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*(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*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
          dYU(I,J,bi,bj) = 0.5*(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

       ENDDO ! bi
      ENDDO ! bj

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

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