/[MITgcm]/MITgcm/model/src/ini_cartesian_grid.F
ViewVC logotype

Diff of /MITgcm/model/src/ini_cartesian_grid.F

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

revision 1.4 by cnh, Mon May 25 18:01:32 1998 UTC revision 1.22 by jmc, Sat Apr 17 18:25:12 2010 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  c#include "PACKAGES_CONFIG.h"
5    #include "CPP_OPTIONS.h"
6    
7  CStartOfInterface  CBOP
8    C     !ROUTINE: INI_CARTESIAN_GRID
9    C     !INTERFACE:
10        SUBROUTINE INI_CARTESIAN_GRID( myThid )        SUBROUTINE INI_CARTESIAN_GRID( myThid )
 C     /==========================================================\  
 C     | SUBROUTINE INI_CARTESIAN_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 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     \==========================================================/  
11    
12    C     !DESCRIPTION: \bv
13    C     *==========================================================*
14    C     | SUBROUTINE INI_CARTESIAN_GRID
15    C     | o Initialise model coordinate system
16    C     *==========================================================*
17    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    C     *==========================================================*
44    C     \ev
45    
46    C     !USES:
47          IMPLICIT NONE
48  C     === Global variables ===  C     === Global variables ===
49  #include "SIZE.h"  #include "SIZE.h"
50  #include "EEPARAMS.h"  #include "EEPARAMS.h"
51  #include "PARAMS.h"  #include "PARAMS.h"
52  #include "GRID.h"  #include "GRID.h"
53    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    
59    C     !INPUT/OUTPUT PARAMETERS:
60  C     == Routine arguments ==  C     == Routine arguments ==
61  C     myThid -  Number of this instance of INI_CARTESIAN_GRID  C     myThid  ::  Number of this instance of INI_CARTESIAN_GRID
62        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
63    
64    C     !LOCAL VARIABLES:
65  C     == Local variables ==  C     == Local variables ==
66  C     xG, yG - Global coordinate location.        INTEGER iG, jG, bi, bj, i, j
67  C     zG        _RL xG0, yG0
68  C     xBase  - South-west corner location for process.  C     "Long" real for temporary coordinate calculation
69  C     yBase  C      NOTICE the extended range of indices!!
70  C     zUpper - Work arrays for upper and lower        _RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
71  C     zLower   cell-face heights.        _RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
72  C     phi    - Temporary scalar  C     These functions return the "global" index with valid values beyond
73  C     xBase  - Temporaries for lower corner coordinate  C     halo regions
74  C     yBase        INTEGER iGl,jGl
75  C     iG, jG - Global coordinate index. Usually used to hold        iGl(i,bi) = 1+MOD(myXGlobalLo-1+(bi-1)*sNx+i+Olx*Nx-1,Nx)
76  C              the south-west global coordinate of a tile.        jGl(j,bj) = 1+MOD(myYGlobalLo-1+(bj-1)*sNy+j+Oly*Ny-1,Ny)
77  C     bi,bj  - Loop counters  c#ifdef ALLOW_EXCH2
78  C     zUpper - Temporary arrays holding z coordinates of  c      INTEGER tN
79  C     zLower   upper and lower faces.  c#endif /* ALLOW_EXCH2 */
80  C     I,J,K  CEOP
81        _RL    xG, yG, zG  
82        _RL    phi  C     For each tile ...
       _RL    zUpper(Nz), zLower(Nz)  
       _RL    xBase, yBase  
       INTEGER iG, jG  
       INTEGER bi, bj  
       INTEGER  I,  J, K  
   
 C--   Simple example of inialisation on cartesian 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  
83        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
        jG = myYGlobalLo + (bj-1)*sNy  
84         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
85    
86    C--     "Global" index (place holder)
87            jG = myYGlobalLo + (bj-1)*sNy
88          iG = myXGlobalLo + (bi-1)*sNx          iG = myXGlobalLo + (bi-1)*sNx
89          yBase = 0. _d 0  c#ifdef ALLOW_EXCH2
90          xBase = 0. _d 0  c        IF ( W2_useE2ioLayOut ) THEN
91          DO i=1,iG-1  cC- note: does not work for non-uniform delX or delY
92           xBase = xBase + delX(i)  c          tN = W2_myTileList(bi,bj)
93          ENDDO  c          iG = exch2_txGlobalo(tN)
94          DO j=1,jG-1  c          jG = exch2_tyGlobalo(tN)
95           yBase = yBase + delY(j)  c        ENDIF
96          ENDDO  c#endif /* ALLOW_EXCH2 */
97          yG = yBase  
98          DO J=1,sNy  C--   First find coordinate of tile corner (meaning outer corner of halo)
99           xG = xBase          xG0 = xgOrigin
100           DO I=1,sNx  C       Find the X-coordinate of the outer grid-line of the "real" tile
101            xc(I,J,bi,bj)  = xG + delX(iG+i-1)*0.5 _d 0          DO i=1, iG-1
102            yc(I,J,bi,bj)  = yG + delY(jG+j-1)*0.5 _d 0           xG0 = xG0 + delX(i)
103            xG = xG + delX(iG+I-1)          ENDDO
104            dxF(I,J,bi,bj) = delX(iG+i-1)  C       Back-step to the outer grid-line of the "halo" region
105            dyF(I,J,bi,bj) = delY(jG+j-1)          DO i=1, Olx
106             xG0 = xG0 - delX( 1+MOD(Olx*Nx-1+iG-i,Nx) )
107            ENDDO
108    C       Find the Y-coordinate of the outer grid-line of the "real" tile
109            yG0 = ygOrigin
110            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             yG0 = yG0 - delY( 1+MOD(Oly*Ny-1+jG-j,Ny) )
116            ENDDO
117    
118    C--     Calculate coordinates of cell corners for N+1 grid-lines
119            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             ENDDO
125            ENDDO
126            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           ENDDO           ENDDO
          yG = yG + delY(jG+J-1)  
132          ENDDO          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 )  
133    
134  C--   Calculate separation between other points  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
135  C     dxG, dyG are separations between cell corners along cell faces.          DO j=1-Oly,sNy+Oly
136        DO bj = myByLo(myThid), myByHi(myThid)           DO i=1-Olx,sNx+Olx
137         DO bi = myBxLo(myThid), myBxHi(myThid)            xG(i,j,bi,bj) = xGloc(i,j)
138          DO J=1,sNy            yG(i,j,bi,bj) = yGloc(i,j)
          DO I=1,sNx  
           dxG(I,J,bi,bj) = (dxF(I,J,bi,bj)+dxF(I,J-1,bi,bj))*0.5 _d 0  
           dyG(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I-1,J,bi,bj))*0.5 _d 0  
139           ENDDO           ENDDO
140          ENDDO          ENDDO
141         ENDDO  
142        ENDDO  C--     Calculate [xC,yC], coordinates of cell centers
143        _EXCH_XY_R4(dxG, myThid )          DO j=1-Oly,sNy+Oly
144        _EXCH_XY_R4(dyG, myThid )           DO i=1-Olx,sNx+Olx
145  C     dxV, dyU are separations between velocity points along cell faces.  C         by averaging
146        DO bj = myByLo(myThid), myByHi(myThid)            xC(i,j,bi,bj) = 0.25 _d 0*(
147         DO bi = myBxLo(myThid), myBxHi(myThid)       &     xGloc(i,j)+xGloc(i+1,j)+xGloc(i,j+1)+xGloc(i+1,j+1) )
148          DO J=1,sNy            yC(i,j,bi,bj) = 0.25 _d 0*(
149           DO I=1,sNx       &     yGloc(i,j)+yGloc(i+1,j)+yGloc(i,j+1)+yGloc(i+1,j+1) )
           dxV(I,J,bi,bj) = (dxG(I,J,bi,bj)+dxG(I-1,J,bi,bj))*0.5 _d 0  
           dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,bi,bj))*0.5 _d 0  
150           ENDDO           ENDDO
151          ENDDO          ENDDO
152         ENDDO  
153        ENDDO  C--     Calculate [dxF,dyF], lengths between cell faces (through center)
154        _EXCH_XY_R4(dxV, myThid )          DO j=1-Oly,sNy+Oly
155        _EXCH_XY_R4(dyU, myThid )           DO i=1-Olx,sNx+Olx
156  C     dxC, dyC is separation between cell centers            dxF(i,j,bi,bj) = delX( iGl(i,bi) )
157        DO bj = myByLo(myThid), myByHi(myThid)            dyF(i,j,bi,bj) = delY( jGl(j,bj) )
        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 D0  
           dyC(I,J,bi,bj)    = (dyF(I,J,bi,bj)+dyF(I,J-1,bi,bj))*0.5 D0  
158           ENDDO           ENDDO
159          ENDDO          ENDDO
160         ENDDO  
161        ENDDO  C--     Calculate [dxG,dyG], lengths along cell boundaries
162        _EXCH_XY_R4(dxC, myThid )          DO j=1-Oly,sNy+Oly
163        _EXCH_XY_R4(dyC, myThid )           DO i=1-Olx,sNx+Olx
164  C     Calculate recipricols            dxG(i,j,bi,bj) = delX( iGl(i,bi) )
165        DO bj = myByLo(myThid), myByHi(myThid)            dyG(i,j,bi,bj) = delY( jGl(j,bj) )
        DO bi = myBxLo(myThid), myBxHi(myThid)  
         DO J=1,sNy  
          DO I=1,sNx  
           rDxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)  
           rDyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)  
           rDxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)  
           rDyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)  
           rDxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)  
           rDyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)  
           rDxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)  
           rDyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)  
166           ENDDO           ENDDO
167          ENDDO          ENDDO
168         ENDDO  
169        ENDDO  C--     The following arrays are not defined in some parts of the halo
170        _EXCH_XY_R4(rDxG, myThid )  C       region. We set them to zero here for safety. If they are ever
171        _EXCH_XY_R4(rDyG, myThid )  C       referred to, especially in the denominator then it is a mistake!
172        _EXCH_XY_R4(rDxC, myThid )          DO j=1-Oly,sNy+Oly
173        _EXCH_XY_R4(rDyC, myThid )           DO i=1-Olx,sNx+Olx
174        _EXCH_XY_R4(rDxF, myThid )            dxC(i,j,bi,bj) = 0.
175        _EXCH_XY_R4(rDyF, myThid )            dyC(i,j,bi,bj) = 0.
176        _EXCH_XY_R4(rDxV, myThid )            dxV(i,j,bi,bj) = 0.
177        _EXCH_XY_R4(rDyU, myThid )            dyU(i,j,bi,bj) = 0.
178  C     Calculate vertical face area            rAw(i,j,bi,bj) = 0.
179        DO bj = myByLo(myThid), myByHi(myThid)            rAs(i,j,bi,bj) = 0.
        DO bi = myBxLo(myThid), myBxHi(myThid)  
         DO J=1,sNy  
          DO I=1,sNx  
           zA(I,J,bi,bj) = dxF(I,J,bi,bj)*dyF(I,J,bi,bj)  
180           ENDDO           ENDDO
181          ENDDO          ENDDO
        ENDDO  
       ENDDO  
182    
183        DO bj = myByLo(myThid), myByHi(myThid)  C--     Calculate [dxC], zonal length between cell centers
184         DO bi = myBxLo(myThid), myBxHi(myThid)          DO j=1-Oly,sNy+Oly
185          DO K=1,Nz           DO i=1-Olx+1,sNx+Olx ! NOTE range
186           DO J=1,sNy            dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
187            DO I=1,sNx           ENDDO
188             IF (HFacC(I,J,K,bi,bj) .NE. 0. D0 ) THEN          ENDDO
189              rHFacC(I,J,K,bi,bj) = 1. D0 / HFacC(I,J,K,bi,bj)  
190             ELSE  C--     Calculate [dyC], meridional length between cell centers
191              rHFacC(I,J,K,bi,bj) = 0. D0          DO j=1-Oly+1,sNy+Oly ! NOTE range
192             ENDIF           DO i=1-Olx,sNx+Olx
193             IF (HFacW(I,J,K,bi,bj) .NE. 0. D0 ) THEN            dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
             rHFacW(I,J,K,bi,bj) = 1. D0 / HFacW(I,J,K,bi,bj)  
             maskW(I,J,K,bi,bj) = 1. D0  
            ELSE  
             rHFacW(I,J,K,bi,bj) = 0. D0  
             maskW(I,J,K,bi,bj) = 0.0 D0  
            ENDIF  
            IF (HFacS(I,J,K,bi,bj) .NE. 0. D0 ) THEN  
             rHFacS(I,J,K,bi,bj) = 1. D0 / HFacS(I,J,K,bi,bj)  
             maskS(I,J,K,bi,bj) = 1. D0  
            ELSE  
             rHFacS(I,J,K,bi,bj) = 0. D0  
             maskS(I,J,K,bi,bj) = 0. D0  
            ENDIF  
           ENDDO  
194           ENDDO           ENDDO
195          ENDDO          ENDDO
196    
197    C--     Calculate [dxV,dyU], length between velocity points (through corners)
198            DO j=1-Oly+1,sNy+Oly ! NOTE range
199             DO i=1-Olx+1,sNx+Olx ! NOTE range
200    C         by averaging (method I)
201              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    C         by averaging (method II)
204    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             ENDDO
207            ENDDO
208    
209    C--     Calculate vertical face area
210            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    C--     Set trigonometric terms & grid orientation:
217              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             ENDDO
222            ENDDO
223    
224    C--     Cosine(lat) scaling
225            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            ENDDO
231    
232    C--   end bi,bj loops
233         ENDDO         ENDDO
234        ENDDO        ENDDO
 C     Now sync. and get/send edge regions that are shared with  
 C     other threads.  
       _EXCH_XYZ_R4(rHFacC    , myThid )  
       _EXCH_XYZ_R4(rHFacW    , myThid )  
       _EXCH_XYZ_R4(rHFacS    , myThid )  
       _EXCH_XYZ_R4(maskW    , myThid )  
       _EXCH_XYZ_R4(maskS    , myThid )  
       _EXCH_XY_R4 (zA       , myThid )  
235    
236  C  C--   Set default (=whole domain) for where relaxation to climatology applies
237          _BEGIN_MASTER(myThid)
238          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          ENDIF
245          _END_MASTER(myThid)
246    
247        RETURN        RETURN
248        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.22
  ViewVC Help
Powered by ViewVC 1.1.22