/[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.13.2.1 by jmc, Fri Jan 26 16:57:19 2001 UTC revision 1.25 by jmc, Sun Feb 17 02:18:16 2013 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     \==========================================================/  
       IMPLICIT NONE  
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     xBase  - South-west corner location for process.  C     i, j    :: loop counters
44  C     yBase  C     delXloc :: mesh spacing in X direction
45  C     zUpper - Work arrays for upper and lower  C     delYloc :: mesh spacing in Y direction
46  C     zLower   cell-face heights.  C     xGloc   :: mesh corner-point location (local "Long" real array type)
47  C     phi    - Temporary scalar  C     yGloc   :: mesh corner-point location (local "Long" real array type)
 C     xBase  - Temporaries for lower corner coordinate  
 C     yBase  
 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    xGloc, yGloc  
       _RL    xBase, yBase  
       INTEGER iG, jG  
48        INTEGER bi, bj        INTEGER bi, bj
49        INTEGER  I,  J        INTEGER i,  j
50          INTEGER gridNx, gridNy
51    C NOTICE the extended range of indices!!
52          _RL delXloc(0-OLx:sNx+OLx)
53          _RL delYloc(0-OLy:sNy+OLy)
54    C NOTICE the extended range of indices!!
55          _RL xGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
56          _RL yGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
57    CEOP
58    
59  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  
60        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
        jG = myYGlobalLo + (bj-1)*sNy  
61         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
62          iG = myXGlobalLo + (bi-1)*sNx  
63          yBase = 0. _d 0  C--     set tile local mesh (same units as delX,deY)
64          xBase = 0. _d 0  C       corresponding to coordinates of cell corners for N+1 grid-lines
65          DO i=1,iG-1  
66           xBase = xBase + delX(i)          CALL INI_LOCAL_GRID(
67          ENDDO       O                       xGloc, yGloc,
68          DO j=1,jG-1       O                       delXloc, delYloc,
69           yBase = yBase + delY(j)       O                       gridNx, gridNy,
70          ENDDO       I                       bi, bj, myThid )
71          yGloc = yBase  
72          DO J=1,sNy  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
73           xGloc = xBase          DO j=1-OLy,sNy+OLy
74           DO I=1,sNx           DO i=1-OLx,sNx+OLx
75            xG(I,J,bi,bj)  = xGloc            xG(i,j,bi,bj) = xGloc(i,j)
76            yG(I,J,bi,bj)  = yGloc            yG(i,j,bi,bj) = yGloc(i,j)
           xc(I,J,bi,bj)  = xGloc + delX(iG+i-1)*0.5 _d 0  
           yc(I,J,bi,bj)  = yGloc + delY(jG+j-1)*0.5 _d 0  
           xGloc = xGloc + delX(iG+I-1)  
           dxF(I,J,bi,bj) = delX(iG+i-1)  
           dyF(I,J,bi,bj) = delY(jG+j-1)  
77           ENDDO           ENDDO
          yGloc = yGloc + delY(jG+J-1)  
78          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 )  
79    
80  C--   Calculate separation between other points  C--     Calculate [xC,yC], coordinates of cell centers
81  C     dxG, dyG are separations between cell corners along cell faces.          DO j=1-OLy,sNy+OLy
82        DO bj = myByLo(myThid), myByHi(myThid)           DO i=1-OLx,sNx+OLx
83         DO bi = myBxLo(myThid), myBxHi(myThid)  C         by averaging
84          DO J=1,sNy            xC(i,j,bi,bj) = 0.25 _d 0*(
85           DO I=1,sNx       &     xGloc(i,j)+xGloc(i+1,j)+xGloc(i,j+1)+xGloc(i+1,j+1) )
86            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*(
87            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) )
88           ENDDO           ENDDO
89          ENDDO          ENDDO
90         ENDDO  
91        ENDDO  C--     Calculate [dxF,dyF], lengths between cell faces (through center)
92        _EXCH_XY_R4(dxG, myThid )          DO j=1-OLy,sNy+OLy
93        _EXCH_XY_R4(dyG, myThid )           DO i=1-OLx,sNx+OLx
94  C     dxV, dyU are separations between velocity points along cell faces.            dxF(i,j,bi,bj) = delXloc(i)
95        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  
96           ENDDO           ENDDO
97          ENDDO          ENDDO
98         ENDDO  
99        ENDDO  C--     Calculate [dxG,dyG], lengths along cell boundaries
100        _EXCH_XY_R4(dxV, myThid )          DO j=1-OLy,sNy+OLy
101        _EXCH_XY_R4(dyU, myThid )           DO i=1-OLx,sNx+OLx
102  C     dxC, dyC is separation between cell centers            dxG(i,j,bi,bj) = delXloc(i)
103        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  
104           ENDDO           ENDDO
105          ENDDO          ENDDO
106         ENDDO  
107        ENDDO  C--     The following arrays are not defined in some parts of the halo
108        _EXCH_XY_R4(dxC, myThid )  C       region. We set them to zero here for safety.
109        _EXCH_XY_R4(dyC, myThid )  C       Note: this is now done earlier in main S/R INI_GRID
110  C     Calculate vertical face area  
111        DO bj = myByLo(myThid), myByHi(myThid)  C--     Calculate [dxC], zonal length between cell centers
112         DO bi = myBxLo(myThid), myBxHi(myThid)          DO j=1-OLy,sNy+OLy
113          DO J=1,sNy           DO i=1-OLx+1,sNx+OLx ! NOTE range
114           DO I=1,sNx            dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
           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)  
           tanPhiAtU(I,J,bi,bj) = 0. _d 0  
           tanPhiAtV(I,J,bi,bj) = 0. _d 0  
115           ENDDO           ENDDO
116          ENDDO          ENDDO
117    
118    C--     Calculate [dyC], meridional length between cell centers
119            DO j=1-OLy+1,sNy+OLy ! NOTE range
120             DO i=1-OLx,sNx+OLx
121              dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
122             ENDDO
123            ENDDO
124    
125    C--     Calculate [dxV,dyU], length between velocity points (through corners)
126            DO j=1-OLy+1,sNy+OLy ! NOTE range
127             DO i=1-OLx+1,sNx+OLx ! NOTE range
128    C         by averaging (method I)
129              dxV(i,j,bi,bj) = 0.5 _d 0*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
130              dyU(i,j,bi,bj) = 0.5 _d 0*(dyG(i,j,bi,bj)+dyG(i,j-1,bi,bj))
131    C         by averaging (method II)
132    c         dxV(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
133    c         dyU(i,j,bi,bj) = 0.5*(dyC(i,j,bi,bj)+dyC(i-1,j,bi,bj))
134             ENDDO
135            ENDDO
136    
137    C--     Calculate vertical face area
138            DO j=1-OLy,sNy+OLy
139             DO i=1-OLx,sNx+OLx
140              rA (i,j,bi,bj) = dxF(i,j,bi,bj)*dyF(i,j,bi,bj)
141              rAw(i,j,bi,bj) = dxC(i,j,bi,bj)*dyG(i,j,bi,bj)
142              rAs(i,j,bi,bj) = dxG(i,j,bi,bj)*dyC(i,j,bi,bj)
143              rAz(i,j,bi,bj) = dxV(i,j,bi,bj)*dyU(i,j,bi,bj)
144    C--     Set trigonometric terms & grid orientation:
145    C       Note: this is now done earlier in main S/R INI_GRID
146    c         tanPhiAtU(i,j,bi,bj) = 0.
147    c         tanPhiAtV(i,j,bi,bj) = 0.
148    c         angleCosC(i,j,bi,bj) = 1.
149    c         angleSinC(i,j,bi,bj) = 0.
150             ENDDO
151            ENDDO
152    
153    C--   end bi,bj loops
154         ENDDO         ENDDO
155        ENDDO        ENDDO
       _EXCH_XY_R4 (rA       , myThid )  
       _EXCH_XY_R4 (rAw      , myThid )  
       _EXCH_XY_R4 (rAs      , myThid )  
       _EXCH_XY_R4 (tanPhiAtU , myThid )  
       _EXCH_XY_R4 (tanPhiAtV , myThid )  
156    
 C  
157        RETURN        RETURN
158        END        END
Legend:
Removed from v.1.13.2.1  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22