/[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.10 by cnh, Fri Nov 6 22:44:46 1998 UTC revision 1.23 by jmc, Mon Dec 12 19:01:01 2011 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  CStartOfInterface  CBOP
7    C     !ROUTINE: INI_CARTESIAN_GRID
8    C     !INTERFACE:
9        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     \==========================================================/  
10    
11    C     !DESCRIPTION: \bv
12    C     *==========================================================*
13    C     | SUBROUTINE INI_CARTESIAN_GRID
14    C     | o Initialise model coordinate system
15    C     *==========================================================*
16    C     | The grid arrays, initialised here, are used throughout
17    C     | the code in evaluating gradients, integrals and spatial
18    C     | avarages. This routine is called separately by each
19    C     | thread and initialises only the region of the domain
20    C     | it is "responsible" for.
21    C     | Under the cartesian grid mode primitive distances
22    C     | in X and Y are in metres. Distance in Z are in m or Pa
23    C     | depending on the vertical gridding mode.
24    C     *==========================================================*
25    C     \ev
26    
27    C     !USES:
28          IMPLICIT NONE
29  C     === Global variables ===  C     === Global variables ===
30  #include "SIZE.h"  #include "SIZE.h"
31  #include "EEPARAMS.h"  #include "EEPARAMS.h"
32  #include "PARAMS.h"  #include "PARAMS.h"
33  #include "GRID.h"  #include "GRID.h"
34    
35    C     !INPUT/OUTPUT PARAMETERS:
36  C     == Routine arguments ==  C     == Routine arguments ==
37  C     myThid -  Number of this instance of INI_CARTESIAN_GRID  C     myThid  :: my Thread Id Number
38        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
39    
40    C     !LOCAL VARIABLES:
41  C     == Local variables ==  C     == Local variables ==
42  C     xG, yG - Global coordinate location.  C     bi,bj   :: tile indices
43  C     zG  C     i, j    :: loop counters
44  C     xBase  - South-west corner location for process.  C     lat     :: Temporary variables used to hold latitude values.
45  C     yBase  C     dlat    :: Temporary variables used to hold latitudes increment.
46  C     zUpper - Work arrays for upper and lower  C     dlon    :: Temporary variables used to hold longitude increment.
47  C     zLower   cell-face heights.  C     delXloc :: mesh spacing in X direction
48  C     phi    - Temporary scalar  C     delYloc :: mesh spacing in Y direction
49  C     xBase  - Temporaries for lower corner coordinate  C     xGloc   :: mesh corner-point location (local "Long" real array type)
50  C     yBase  C     yGloc   :: mesh corner-point location (local "Long" real array type)
 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     I,J,K  
       _RL    xG, yG, zG  
       _RL    phi  
       _RL    zUpper(Nr), zLower(Nr)  
       _RL    xBase, yBase  
       INTEGER iG, jG  
51        INTEGER bi, bj        INTEGER bi, bj
52        INTEGER  I,  J, K        INTEGER i,  j
53          INTEGER gridNx, gridNy
54          _RL lat, dlat, dlon
55    C NOTICE the extended range of indices!!
56          _RL delXloc(0-OLx:sNx+OLx)
57          _RL delYloc(0-OLy:sNy+OLy)
58    C NOTICE the extended range of indices!!
59          _RL xGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
60          _RL yGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
61    CEOP
62    
63  C--   Simple example of inialisation on cartesian grid  C--   For each tile ...
 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  
       xC0 = 0. _d 0  
       yC0 = 0. _d 0  
64        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
        jG = myYGlobalLo + (bj-1)*sNy  
65         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
66          iG = myXGlobalLo + (bi-1)*sNx  
67          yBase = 0. _d 0  C--     set tile local mesh (same units as delX,deY)
68          xBase = 0. _d 0  C       corresponding to coordinates of cell corners for N+1 grid-lines
69          DO i=1,iG-1  
70           xBase = xBase + delX(i)          CALL INI_LOCAL_GRID(
71          ENDDO       O                       xGloc, yGloc,
72          DO j=1,jG-1       O                       delXloc, delYloc,
73           yBase = yBase + delY(j)       O                       gridNx, gridNy,
74          ENDDO       I                       bi, bj, myThid )
75          yG = yBase  
76          DO J=1,sNy  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
77           xG = xBase          DO j=1-OLy,sNy+OLy
78           DO I=1,sNx           DO i=1-OLx,sNx+OLx
79            xc(I,J,bi,bj)  = xG + delX(iG+i-1)*0.5 _d 0            xG(i,j,bi,bj) = xGloc(i,j)
80            yc(I,J,bi,bj)  = yG + delY(jG+j-1)*0.5 _d 0            yG(i,j,bi,bj) = yGloc(i,j)
           xG = xG + delX(iG+I-1)  
           dxF(I,J,bi,bj) = delX(iG+i-1)  
           dyF(I,J,bi,bj) = delY(jG+j-1)  
81           ENDDO           ENDDO
          yG = yG + delY(jG+J-1)  
82          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 )  
83    
84  C--   Calculate separation between other points  C--     Calculate [xC,yC], coordinates of cell centers
85  C     dxG, dyG are separations between cell corners along cell faces.          DO j=1-OLy,sNy+OLy
86        DO bj = myByLo(myThid), myByHi(myThid)           DO i=1-OLx,sNx+OLx
87         DO bi = myBxLo(myThid), myBxHi(myThid)  C         by averaging
88          DO J=1,sNy            xC(i,j,bi,bj) = 0.25 _d 0*(
89           DO I=1,sNx       &     xGloc(i,j)+xGloc(i+1,j)+xGloc(i,j+1)+xGloc(i+1,j+1) )
90            dxG(I,J,bi,bj) = (dxF(I,J,bi,bj)+dxF(I,J-1,bi,bj))*0.5 _d 0            yC(i,j,bi,bj) = 0.25 _d 0*(
91            dyG(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I-1,J,bi,bj))*0.5 _d 0       &     yGloc(i,j)+yGloc(i+1,j)+yGloc(i,j+1)+yGloc(i+1,j+1) )
92           ENDDO           ENDDO
93          ENDDO          ENDDO
94         ENDDO  
95        ENDDO  C--     Calculate [dxF,dyF], lengths between cell faces (through center)
96        _EXCH_XY_R4(dxG, myThid )          DO j=1-OLy,sNy+OLy
97        _EXCH_XY_R4(dyG, myThid )           DO i=1-OLx,sNx+OLx
98  C     dxV, dyU are separations between velocity points along cell faces.            dxF(i,j,bi,bj) = delXloc(i)
99        DO bj = myByLo(myThid), myByHi(myThid)            dyF(i,j,bi,bj) = delYloc(j)
        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  
           dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,bi,bj))*0.5 _d 0  
100           ENDDO           ENDDO
101          ENDDO          ENDDO
102         ENDDO  
103        ENDDO  C--     Calculate [dxG,dyG], lengths along cell boundaries
104        _EXCH_XY_R4(dxV, myThid )          DO j=1-OLy,sNy+OLy
105        _EXCH_XY_R4(dyU, myThid )           DO i=1-OLx,sNx+OLx
106  C     dxC, dyC is separation between cell centers            dxG(i,j,bi,bj) = delXloc(i)
107        DO bj = myByLo(myThid), myByHi(myThid)            dyG(i,j,bi,bj) = delYloc(j)
        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  
108           ENDDO           ENDDO
109          ENDDO          ENDDO
110         ENDDO  
111        ENDDO  C--     The following arrays are not defined in some parts of the halo
112        _EXCH_XY_R4(dxC, myThid )  C       region. We set them to zero here for safety. If they are ever
113        _EXCH_XY_R4(dyC, myThid )  C       referred to, especially in the denominator then it is a mistake!
114  C     Calculate vertical face area          DO j=1-OLy,sNy+OLy
115        DO bj = myByLo(myThid), myByHi(myThid)           DO i=1-OLx,sNx+OLx
116         DO bi = myBxLo(myThid), myBxHi(myThid)            dxC(i,j,bi,bj) = 0.
117          DO J=1,sNy            dyC(i,j,bi,bj) = 0.
118           DO I=1,sNx            dxV(i,j,bi,bj) = 0.
119            rA(I,J,bi,bj) = dxF(I,J,bi,bj)*dyF(I,J,bi,bj)            dyU(i,j,bi,bj) = 0.
120            tanPhiAtU(I,J,bi,bj) = 0. _d 0            rAw(i,j,bi,bj) = 0.
121            tanPhiAtV(I,J,bi,bj) = 0. _d 0            rAs(i,j,bi,bj) = 0.
122             ENDDO
123            ENDDO
124    
125    C--     Calculate [dxC], zonal length between cell centers
126            DO j=1-OLy,sNy+OLy
127             DO i=1-OLx+1,sNx+OLx ! NOTE range
128              dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
129             ENDDO
130            ENDDO
131    
132    C--     Calculate [dyC], meridional length between cell centers
133            DO j=1-OLy+1,sNy+OLy ! NOTE range
134             DO i=1-OLx,sNx+OLx
135              dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
136           ENDDO           ENDDO
137          ENDDO          ENDDO
138    
139    C--     Calculate [dxV,dyU], length between velocity points (through corners)
140            DO j=1-OLy+1,sNy+OLy ! NOTE range
141             DO i=1-OLx+1,sNx+OLx ! NOTE range
142    C         by averaging (method I)
143              dxV(i,j,bi,bj) = 0.5 _d 0*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
144              dyU(i,j,bi,bj) = 0.5 _d 0*(dyG(i,j,bi,bj)+dyG(i,j-1,bi,bj))
145    C         by averaging (method II)
146    c         dxV(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
147    c         dyU(i,j,bi,bj) = 0.5*(dyC(i,j,bi,bj)+dyC(i-1,j,bi,bj))
148             ENDDO
149            ENDDO
150    
151    C--     Calculate vertical face area
152            DO j=1-OLy,sNy+OLy
153             DO i=1-OLx,sNx+OLx
154              rA (i,j,bi,bj) = dxF(i,j,bi,bj)*dyF(i,j,bi,bj)
155              rAw(i,j,bi,bj) = dxC(i,j,bi,bj)*dyG(i,j,bi,bj)
156              rAs(i,j,bi,bj) = dxG(i,j,bi,bj)*dyC(i,j,bi,bj)
157              rAz(i,j,bi,bj) = dxV(i,j,bi,bj)*dyU(i,j,bi,bj)
158    C--     Set trigonometric terms & grid orientation:
159              tanPhiAtU(i,j,bi,bj) = 0.
160              tanPhiAtV(i,j,bi,bj) = 0.
161              angleCosC(i,j,bi,bj) = 1.
162              angleSinC(i,j,bi,bj) = 0.
163             ENDDO
164            ENDDO
165    
166    C--     Cosine(lat) scaling
167            DO j=1-OLy,sNy+OLy
168             cosFacU(j,bi,bj) = 1.
169             cosFacV(j,bi,bj) = 1.
170             sqcosFacU(j,bi,bj)=1.
171             sqcosFacV(j,bi,bj)=1.
172            ENDDO
173    
174    C--   end bi,bj loops
175         ENDDO         ENDDO
176        ENDDO        ENDDO
       _EXCH_XY_R4 (rA       , myThid )  
       _EXCH_XY_R4 (tanPhiAtU , myThid )  
       _EXCH_XY_R4 (tanPhiAtV , myThid )  
177    
 C  
178        RETURN        RETURN
179        END        END
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22