model/src/ini_spherical_polar_grid.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_spherical_polar_grid.F,v 1.10 1998/11/06 22:44:47 cnh Exp $

#include "CPP_OPTIONS.h"

CStartOfInterface
      SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )
C     /==========================================================\
C     | SUBROUTINE INI_SPHERICAL_POLAR_GRID                      |
C     | o Initialise model coordinate system                     |
C     |==========================================================|
C     | These arrays are used throughout the code in evaluating  |
C     | gradients, integrals and spatial avarages. This routine  |
C     | is called separately by each thread and initialise 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. 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 spherical polar grid mode primitive distances  |
C     | in X and Y are in degrees. Distance in Z are in m or Pa  |
C     | depending on the vertical gridding mode.                 |
C     \==========================================================/

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"

C     == Routine arguments ==
C     myThid -  Number of this instance of INI_CARTESIAN_GRID
      INTEGER myThid
CEndOfInterface

C     == Local variables ==
C     xG, yG - Global coordinate location.
C     zG
C     xBase  - South-west corner location for process.
C     yBase
C     zUpper - Work arrays for upper and lower 
C     zLower   cell-face heights.
C     phi    - Temporary scalar
C     iG, jG - Global coordinate index. Usually used to hold
C              the south-west global coordinate of a tile.
C     bi,bj  - Loop counters
C     zUpper - Temporary arrays holding z coordinates of
C     zLower   upper and lower faces.
C     xBase  - Lower coordinate for this threads cells
C     yBase
C     lat, latN, - Temporary variables used to hold latitude
C     latS         values.
C     I,J,K
      _RL    xG, yG, zG
      _RL    phi
      _RL    zUpper(Nr), zLower(Nr)
      _RL    xBase, yBase
      INTEGER iG, jG
      INTEGER bi, bj
      INTEGER  I,  J, K
      _RL lat, latS, latN

C--   Example of inialisation for spherical polar grid
C--   First set coordinates of cell centers
C     This operation is only performed at start up so for more
C     complex configurations it is usually OK to pass iG, jG to a custom 
C     function and have it return xG and yG.
C     Set up my local grid first
C     Note: In the spherical polar case delX and delY are given in
C           degrees and are relative to some starting latitude and
C           longitude - phiMin and thetaMin.
      xC0 = thetaMin
      yC0 = phiMin
      DO bj = myByLo(myThid), myByHi(myThid)
       jG = myYGlobalLo + (bj-1)*sNy
       DO bi = myBxLo(myThid), myBxHi(myThid)
        iG = myXGlobalLo + (bi-1)*sNx
        yBase = phiMin
        xBase = thetaMin
        DO i=1,iG-1
         xBase = xBase + delX(i)
        ENDDO
        DO j=1,jG-1
         yBase = yBase + delY(j)
        ENDDO
        yG = yBase
        DO J=1,sNy
         xG = xBase
         DO I=1,sNx
          xc(I,J,bi,bj)  = xG + delX(iG+i-1)*0.5 _d 0
          yc(I,J,bi,bj)  = yG + delY(jG+j-1)*0.5 _d 0
          xG = xG + delX(iG+I-1)
          dxF(I,J,bi,bj) = delX(iG+i-1)*deg2rad
     &    *rSphere*COS(yc(I,J,bi,bj)*deg2rad)
          dyF(I,J,bi,bj) = delY(jG+j-1)*deg2rad*rSphere
         ENDDO
         yG = yG + delY(jG+J-1)
        ENDDO
       ENDDO
      ENDDO
C     Now sync. and get edge regions from other threads and/or processes.
C     Note: We could just set the overlap regions ourselves here but
C           exchanging edges is safer and is good practice!
      _EXCH_XY_R4( xc, myThid )
      _EXCH_XY_R4( yc, myThid )
      _EXCH_XY_R4(dxF, myThid )
      _EXCH_XY_R4(dyF, myThid )

C--   Calculate separation between other points
C     dxG, dyG are separations between cell corners along cell faces.
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO J=1,sNy
         DO I=1,sNx
          jG = myYGlobalLo + (bj-1)*sNy + J-1
          iG = myXGlobalLo + (bi-1)*sNx + I-1
          lat = yc(I,J,bi,bj)-delY(jG) * 0.5 _d 0
          dxG(I,J,bi,bj) = rSphere*COS(lat*deg2rad)*delX(iG)*deg2rad
          dyG(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I-1,J,bi,bj))*0.5 _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XY_R4(dxG, myThid )
      _EXCH_XY_R4(dyG, myThid )
C     dxC, dyC is separation between cell centers
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO J=1,sNy
         DO I=1,sNx
          dxC(I,J,bi,bj)    = (dxF(I,J,bi,bj)+dxF(I-1,J,bi,bj))*0.5 _d 0
          dyC(I,J,bi,bj)    = (dyF(I,J,bi,bj)+dyF(I,J-1,bi,bj))*0.5 _d 0
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XY_R4(dxC, myThid )
      _EXCH_XY_R4(dyC, myThid )
C     dxV, dyU are separations between velocity points along cell faces.
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO J=1,sNy
         DO I=1,sNx
          dxV(I,J,bi,bj) = (dxG(I,J,bi,bj)+dxG(I-1,J,bi,bj))*0.5 _d 0
#ifdef OLD_UV_GEOMETRY
          dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,bi,bj))*0.5 _d 0
#else
          dyU(I,J,bi,bj) = (dyC(I,J,bi,bj)+dyC(I-1,J,bi,bj))*0.5 _d 0
#endif
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XY_R4(dxV, myThid )
      _EXCH_XY_R4(dyU, myThid )
C     Calculate vertical face area and trigonometric terms
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO J=1,sNy
         DO I=1,sNx
          jG = myYGlobalLo + (bj-1)*sNy + J-1
          latS = yc(i,j,bi,bj)-delY(jG)*0.5 _d 0
          latN = yc(i,j,bi,bj)+delY(jG)*0.5 _d 0
#ifdef OLD_UV_GEOMETRY
          rA(I,J,bi,bj) = dyF(I,J,bi,bj)
     &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))
#else
          rA(I,J,bi,bj) = rSphere*delX(iG)*deg2rad
     &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))
#endif
C         Area cannot be zero but sin can be if lat if < -90.
          IF ( rA(I,J,bi,bj) .LT. 0. ) rA(I,J,bi,bj) = -rA(I,J,bi,bj)
          tanPhiAtU(i,j,bi,bj)=tan(_yC(i,j,bi,bj)*deg2rad)
          tanPhiAtV(i,j,bi,bj)=tan(latS*deg2rad)
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XY_R4 (rA       , myThid )
      _EXCH_XY_R4 (tanPhiAtU , myThid )
      _EXCH_XY_R4 (tanPhiAtV , myThid )
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO J=1,sNy
         DO I=1,sNx
          iG = myXGlobalLo + (bi-1)*sNx + I-1
          jG = myYGlobalLo + (bj-1)*sNy + J-1
          latS = yc(i,j-1,bi,bj)
          latN = yc(i,j,bi,bj)
#ifdef OLD_UV_GEOMETRY
          rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))
          rAs(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I,J-1,bi,bj))
#else
          rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))
          rAs(I,J,bi,bj) = rSphere*delX(iG)*deg2rad
     &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))
#endif
         ENDDO
        ENDDO
       ENDDO
      ENDDO
      _EXCH_XY_R4 (rAw      , myThid )
      _EXCH_XY_R4 (rAs      , myThid )
C
      RETURN
      END
1	adcroft	1.11	C $Header: /u/gcmpack/models/MITgcmUV/model/src/ini_spherical_polar_grid.F,v 1.10 1998/11/06 22:44:47 cnh Exp $
2	cnh	1.1
3	cnh	1.10	#include "CPP_OPTIONS.h"
4	cnh	1.1
5			CStartOfInterface
6			SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )
7			C /==========================================================\
8			C \| SUBROUTINE INI_SPHERICAL_POLAR_GRID \|
9			C \| o Initialise model coordinate system \|
10			C \|==========================================================\|
11			C \| These arrays are used throughout the code in evaluating \|
12			C \| gradients, integrals and spatial avarages. This routine \|
13			C \| is called separately by each thread and initialise only \|
14			C \| the region of the domain it is "responsible" for. \|
15			C \| Notes: \|
16			C \| Two examples are included. One illustrates the \|
17			C \| initialisation of a cartesian grid. The other shows the \|
18			C \| inialisation of a spherical polar grid. Other orthonormal\|
19			C \| grids can be fitted into this design. In this case \|
20			C \| custom metric terms also need adding to account for the \|
21			C \| projections of velocity vectors onto these grids. \|
22			C \| The structure used here also makes it possible to \|
23			C \| implement less regular grid mappings. In particular \|
24			C \| o Schemes which leave out blocks of the domain that are \|
25			C \| all land could be supported. \|
26			C \| o Multi-level schemes such as icosohedral or cubic \|
27			C \| grid projections onto a sphere can also be fitted \|
28			C \| within the strategy we use. \|
29			C \| Both of the above also require modifying the support \|
30			C \| routines that map computational blocks to simulation \|
31			C \| domain blocks. \|
32			C \| Under the spherical polar grid mode primitive distances \|
33			C \| in X and Y are in degrees. Distance in Z are in m or Pa \|
34			C \| depending on the vertical gridding mode. \|
35			C \==========================================================/
36
37			C === Global variables ===
38			#include "SIZE.h"
39			#include "EEPARAMS.h"
40			#include "PARAMS.h"
41			#include "GRID.h"
42
43			C == Routine arguments ==
44			C myThid - Number of this instance of INI_CARTESIAN_GRID
45			INTEGER myThid
46			CEndOfInterface
47
48			C == Local variables ==
49			C xG, yG - Global coordinate location.
50			C zG
51			C xBase - South-west corner location for process.
52			C yBase
53			C zUpper - Work arrays for upper and lower
54			C zLower cell-face heights.
55			C phi - Temporary scalar
56			C iG, jG - Global coordinate index. Usually used to hold
57			C the south-west global coordinate of a tile.
58			C bi,bj - Loop counters
59			C zUpper - Temporary arrays holding z coordinates of
60			C zLower upper and lower faces.
61			C xBase - Lower coordinate for this threads cells
62			C yBase
63			C lat, latN, - Temporary variables used to hold latitude
64			C latS values.
65			C I,J,K
66			_RL xG, yG, zG
67			_RL phi
68	cnh	1.8	_RL zUpper(Nr), zLower(Nr)
69	cnh	1.1	_RL xBase, yBase
70			INTEGER iG, jG
71			INTEGER bi, bj
72			INTEGER I, J, K
73			_RL lat, latS, latN
74
75			C-- Example of inialisation for spherical polar grid
76			C-- First set coordinates of cell centers
77			C This operation is only performed at start up so for more
78			C complex configurations it is usually OK to pass iG, jG to a custom
79			C function and have it return xG and yG.
80			C Set up my local grid first
81			C Note: In the spherical polar case delX and delY are given in
82			C degrees and are relative to some starting latitude and
83			C longitude - phiMin and thetaMin.
84	cnh	1.5	xC0 = thetaMin
85			yC0 = phiMin
86	cnh	1.1	DO bj = myByLo(myThid), myByHi(myThid)
87			jG = myYGlobalLo + (bj-1)*sNy
88			DO bi = myBxLo(myThid), myBxHi(myThid)
89			iG = myXGlobalLo + (bi-1)*sNx
90			yBase = phiMin
91			xBase = thetaMin
92			DO i=1,iG-1
93			xBase = xBase + delX(i)
94			ENDDO
95			DO j=1,jG-1
96			yBase = yBase + delY(j)
97			ENDDO
98			yG = yBase
99			DO J=1,sNy
100			xG = xBase
101			DO I=1,sNx
102			xc(I,J,bi,bj) = xG + delX(iG+i-1)*0.5 _d 0
103			yc(I,J,bi,bj) = yG + delY(jG+j-1)*0.5 _d 0
104			xG = xG + delX(iG+I-1)
105	cnh	1.10	dxF(I,J,bi,bj) = delX(iG+i-1)*deg2rad
106			& rSphereCOS(yc(I,J,bi,bj)*deg2rad)
107	cnh	1.1	dyF(I,J,bi,bj) = delY(jG+j-1)deg2radrSphere
108			ENDDO
109			yG = yG + delY(jG+J-1)
110			ENDDO
111			ENDDO
112			ENDDO
113			C Now sync. and get edge regions from other threads and/or processes.
114			C Note: We could just set the overlap regions ourselves here but
115			C exchanging edges is safer and is good practice!
116			_EXCH_XY_R4( xc, myThid )
117			_EXCH_XY_R4( yc, myThid )
118			_EXCH_XY_R4(dxF, myThid )
119			_EXCH_XY_R4(dyF, myThid )
120
121			C-- Calculate separation between other points
122			C dxG, dyG are separations between cell corners along cell faces.
123			DO bj = myByLo(myThid), myByHi(myThid)
124			DO bi = myBxLo(myThid), myBxHi(myThid)
125			DO J=1,sNy
126			DO I=1,sNx
127			jG = myYGlobalLo + (bj-1)*sNy + J-1
128			iG = myXGlobalLo + (bi-1)*sNx + I-1
129			lat = yc(I,J,bi,bj)-delY(jG) * 0.5 _d 0
130			dxG(I,J,bi,bj) = rSphereCOS(latdeg2rad)delX(iG)deg2rad
131			dyG(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I-1,J,bi,bj))*0.5 _d 0
132			ENDDO
133			ENDDO
134			ENDDO
135			ENDDO
136			_EXCH_XY_R4(dxG, myThid )
137			_EXCH_XY_R4(dyG, myThid )
138	adcroft	1.11	C dxC, dyC is separation between cell centers
139	cnh	1.1	DO bj = myByLo(myThid), myByHi(myThid)
140			DO bi = myBxLo(myThid), myBxHi(myThid)
141			DO J=1,sNy
142			DO I=1,sNx
143	adcroft	1.11	dxC(I,J,bi,bj) = (dxF(I,J,bi,bj)+dxF(I-1,J,bi,bj))*0.5 _d 0
144			dyC(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I,J-1,bi,bj))*0.5 _d 0
145	cnh	1.1	ENDDO
146			ENDDO
147			ENDDO
148			ENDDO
149	adcroft	1.11	_EXCH_XY_R4(dxC, myThid )
150			_EXCH_XY_R4(dyC, myThid )
151			C dxV, dyU are separations between velocity points along cell faces.
152	cnh	1.1	DO bj = myByLo(myThid), myByHi(myThid)
153			DO bi = myBxLo(myThid), myBxHi(myThid)
154			DO J=1,sNy
155			DO I=1,sNx
156	adcroft	1.11	dxV(I,J,bi,bj) = (dxG(I,J,bi,bj)+dxG(I-1,J,bi,bj))*0.5 _d 0
157			#ifdef OLD_UV_GEOMETRY
158			dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,bi,bj))*0.5 _d 0
159			#else
160			dyU(I,J,bi,bj) = (dyC(I,J,bi,bj)+dyC(I-1,J,bi,bj))*0.5 _d 0
161			#endif
162	cnh	1.1	ENDDO
163			ENDDO
164			ENDDO
165			ENDDO
166	adcroft	1.11	_EXCH_XY_R4(dxV, myThid )
167			_EXCH_XY_R4(dyU, myThid )
168	adcroft	1.6	C Calculate vertical face area and trigonometric terms
169	cnh	1.1	DO bj = myByLo(myThid), myByHi(myThid)
170			DO bi = myBxLo(myThid), myBxHi(myThid)
171			DO J=1,sNy
172			DO I=1,sNx
173			jG = myYGlobalLo + (bj-1)*sNy + J-1
174			latS = yc(i,j,bi,bj)-delY(jG)*0.5 _d 0
175			latN = yc(i,j,bi,bj)+delY(jG)*0.5 _d 0
176	adcroft	1.11	#ifdef OLD_UV_GEOMETRY
177	cnh	1.8	rA(I,J,bi,bj) = dyF(I,J,bi,bj)
178	cnh	1.1	& rSphere(SIN(latNdeg2rad)-SIN(latSdeg2rad))
179	adcroft	1.11	#else
180			rA(I,J,bi,bj) = rSpheredelX(iG)deg2rad
181			& rSphere(SIN(latNdeg2rad)-SIN(latSdeg2rad))
182			#endif
183	cnh	1.10	C Area cannot be zero but sin can be if lat if < -90.
184			IF ( rA(I,J,bi,bj) .LT. 0. ) rA(I,J,bi,bj) = -rA(I,J,bi,bj)
185	adcroft	1.6	tanPhiAtU(i,j,bi,bj)=tan(_yC(i,j,bi,bj)*deg2rad)
186			tanPhiAtV(i,j,bi,bj)=tan(latS*deg2rad)
187	cnh	1.1	ENDDO
188			ENDDO
189			ENDDO
190			ENDDO
191	cnh	1.8	_EXCH_XY_R4 (rA , myThid )
192	cnh	1.5	_EXCH_XY_R4 (tanPhiAtU , myThid )
193			_EXCH_XY_R4 (tanPhiAtV , myThid )
194	adcroft	1.11	DO bj = myByLo(myThid), myByHi(myThid)
195			DO bi = myBxLo(myThid), myBxHi(myThid)
196			DO J=1,sNy
197			DO I=1,sNx
198			iG = myXGlobalLo + (bi-1)*sNx + I-1
199			jG = myYGlobalLo + (bj-1)*sNy + J-1
200			latS = yc(i,j-1,bi,bj)
201			latN = yc(i,j,bi,bj)
202			#ifdef OLD_UV_GEOMETRY
203			rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))
204			rAs(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I,J-1,bi,bj))
205			#else
206			rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))
207			rAs(I,J,bi,bj) = rSpheredelX(iG)deg2rad
208			& rSphere(SIN(latNdeg2rad)-SIN(latSdeg2rad))
209			#endif
210			ENDDO
211			ENDDO
212			ENDDO
213			ENDDO
214			_EXCH_XY_R4 (rAw , myThid )
215			_EXCH_XY_R4 (rAs , myThid )
216	cnh	1.1	C
217			RETURN
218			END