/[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.17 by cnh, Wed Sep 26 18:09:15 2001 UTC revision 1.24 by jmc, Thu Dec 22 00:11:11 2011 UTC
# Line 7  CBOP Line 7  CBOP
7  C     !ROUTINE: INI_CARTESIAN_GRID  C     !ROUTINE: INI_CARTESIAN_GRID
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE INI_CARTESIAN_GRID( myThid )        SUBROUTINE INI_CARTESIAN_GRID( myThid )
10    
11  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
12  C     *==========================================================*  C     *==========================================================*
13  C     | SUBROUTINE INI_CARTESIAN_GRID                              C     | SUBROUTINE INI_CARTESIAN_GRID
14  C     | o Initialise model coordinate system                        C     | o Initialise model coordinate system
15  C     *==========================================================*  C     *==========================================================*
16  C     | The grid arrays, initialised here, are used throughout  C     | The grid arrays, initialised here, are used throughout
17  C     | the code in evaluating gradients, integrals and spatial  C     | the code in evaluating gradients, integrals and spatial
18  C     | avarages. This routine    C     | avarages. This routine is called separately by each
19  C     | is called separately by each thread and initialises only    C     | thread and initialises only the region of the domain
20  C     | the region of the domain it is "responsible" for.          C     | it is "responsible" for.
21  C     | Notes:                                                      C     | Under the cartesian grid mode primitive distances
22  C     | Two examples are included. One illustrates the              C     | in X and Y are in metres. Distance in Z are in m or Pa
23  C     | initialisation of a cartesian grid (this routine).  C     | depending on the vertical gridding mode.
 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.                    
24  C     *==========================================================*  C     *==========================================================*
25  C     \ev  C     \ev
26    
# Line 51  C     === Global variables === Line 34  C     === Global variables ===
34    
35  C     !INPUT/OUTPUT PARAMETERS:  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
39    
40  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
41  C     == Local variables ==  C     == Local variables ==
42        INTEGER iG, jG, bi, bj, I,  J  C     bi,bj   :: tile indices
43        _RL xG0, yG0  C     i, j    :: loop counters
44  C     "Long" real for temporary coordinate calculation  C     delXloc :: mesh spacing in X direction
45  C      NOTICE the extended range of indices!!  C     delYloc :: mesh spacing in Y direction
46        _RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)  C     xGloc   :: mesh corner-point location (local "Long" real array type)
47        _RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)  C     yGloc   :: mesh corner-point location (local "Long" real array type)
48  C     These functions return the "global" index with valid values beyond        INTEGER bi, bj
49  C     halo regions        INTEGER i,  j
50        INTEGER iGl,jGl        INTEGER gridNx, gridNy
51        iGl(I,bi) = 1+mod(myXGlobalLo-1+(bi-1)*sNx+I+Olx*Nx-1,Nx)  C NOTICE the extended range of indices!!
52        jGl(J,bj) = 1+mod(myYGlobalLo-1+(bj-1)*sNy+J+Oly*Ny-1,Ny)        _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  CEOP
58    
59  C     For each tile ...  C--   For each tile ...
60        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
61         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
62    
63  C--     "Global" index (place holder)  C--     set tile local mesh (same units as delX,deY)
64          jG = myYGlobalLo + (bj-1)*sNy  C       corresponding to coordinates of cell corners for N+1 grid-lines
65          iG = myXGlobalLo + (bi-1)*sNx  
66            CALL INI_LOCAL_GRID(
67  C--   First find coordinate of tile corner (meaning outer corner of halo)       O                       xGloc, yGloc,
68          xG0 = 0.       O                       delXloc, delYloc,
69  C       Find the X-coordinate of the outer grid-line of the "real" tile       O                       gridNx, gridNy,
70          DO i=1, iG-1       I                       bi, bj, myThid )
          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  
71    
72  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
73          DO J=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
74           DO I=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
75            xG(I,J,bi,bj) = xGloc(I,J)            xG(i,j,bi,bj) = xGloc(i,j)
76            yG(I,J,bi,bj) = yGloc(I,J)            yG(i,j,bi,bj) = yGloc(i,j)
77           ENDDO           ENDDO
78          ENDDO          ENDDO
79    
80  C--     Calculate [xC,yC], coordinates of cell centers  C--     Calculate [xC,yC], coordinates of cell centers
81          DO J=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
82           DO I=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
83  C         by averaging  C         by averaging
84            xC(I,J,bi,bj) = 0.25*(            xC(i,j,bi,bj) = 0.25 _d 0*(
85       &     xGloc(I,J)+xGloc(I+1,J)+xGloc(I,J+1)+xGloc(I+1,J+1) )       &     xGloc(i,j)+xGloc(i+1,j)+xGloc(i,j+1)+xGloc(i+1,j+1) )
86            yC(I,J,bi,bj) = 0.25*(            yC(i,j,bi,bj) = 0.25 _d 0*(
87       &     yGloc(I,J)+yGloc(I+1,J)+yGloc(I,J+1)+yGloc(I+1,J+1) )       &     yGloc(i,j)+yGloc(i+1,j)+yGloc(i,j+1)+yGloc(i+1,j+1) )
88           ENDDO           ENDDO
89          ENDDO          ENDDO
90    
91  C--     Calculate [dxF,dyF], lengths between cell faces (through center)  C--     Calculate [dxF,dyF], lengths between cell faces (through center)
92          DO J=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
93           DO I=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
94            dXF(I,J,bi,bj) = delX( iGl(I,bi) )            dxF(i,j,bi,bj) = delXloc(i)
95            dYF(I,J,bi,bj) = delY( jGl(J,bj) )            dyF(i,j,bi,bj) = delYloc(j)
96           ENDDO           ENDDO
97          ENDDO          ENDDO
98    
99  C--     Calculate [dxG,dyG], lengths along cell boundaries  C--     Calculate [dxG,dyG], lengths along cell boundaries
100          DO J=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
101           DO I=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
102            dXG(I,J,bi,bj) = delX( iGl(I,bi) )            dxG(i,j,bi,bj) = delXloc(i)
103            dYG(I,J,bi,bj) = delY( jGl(J,bj) )            dyG(i,j,bi,bj) = delYloc(j)
104           ENDDO           ENDDO
105          ENDDO          ENDDO
106    
107  C--     The following arrays are not defined in some parts of the halo  C--     The following arrays are not defined in some parts of the halo
108  C       region. We set them to zero here for safety. If they are ever  C       region. We set them to zero here for safety. If they are ever
109  C       referred to, especially in the denominator then it is a mistake!  C       referred to, especially in the denominator then it is a mistake!
110          DO J=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
111           DO I=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
112            dXC(I,J,bi,bj) = 0.            dxC(i,j,bi,bj) = 0.
113            dYC(I,J,bi,bj) = 0.            dyC(i,j,bi,bj) = 0.
114            dXV(I,J,bi,bj) = 0.            dxV(i,j,bi,bj) = 0.
115            dYU(I,J,bi,bj) = 0.            dyU(i,j,bi,bj) = 0.
116            rAw(I,J,bi,bj) = 0.            rAw(i,j,bi,bj) = 0.
117            rAs(I,J,bi,bj) = 0.            rAs(i,j,bi,bj) = 0.
118           ENDDO           ENDDO
119          ENDDO          ENDDO
120    
121  C--     Calculate [dxC], zonal length between cell centers  C--     Calculate [dxC], zonal length between cell centers
122          DO J=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
123           DO I=1-Olx+1,sNx+Olx ! NOTE range           DO i=1-OLx+1,sNx+OLx ! NOTE range
124            dXC(I,J,bi,bj) = 0.5D0*(dXF(I,J,bi,bj)+dXF(I-1,J,bi,bj))            dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
125           ENDDO           ENDDO
126          ENDDO          ENDDO
127    
128  C--     Calculate [dyC], meridional length between cell centers  C--     Calculate [dyC], meridional length between cell centers
129          DO J=1-Oly+1,sNy+Oly ! NOTE range          DO j=1-OLy+1,sNy+OLy ! NOTE range
130           DO I=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
131            dYC(I,J,bi,bj) = 0.5*(dYF(I,J,bi,bj)+dYF(I,J-1,bi,bj))            dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
132           ENDDO           ENDDO
133          ENDDO          ENDDO
134    
135  C--     Calculate [dxV,dyU], length between velocity points (through corners)  C--     Calculate [dxV,dyU], length between velocity points (through corners)
136          DO J=1-Oly+1,sNy+Oly ! NOTE range          DO j=1-OLy+1,sNy+OLy ! NOTE range
137           DO I=1-Olx+1,sNx+Olx ! NOTE range           DO i=1-OLx+1,sNx+OLx ! NOTE range
138  C         by averaging (method I)  C         by averaging (method I)
139            dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))            dxV(i,j,bi,bj) = 0.5 _d 0*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
140            dYU(I,J,bi,bj) = 0.5*(dYG(I,J,bi,bj)+dYG(I,J-1,bi,bj))            dyU(i,j,bi,bj) = 0.5 _d 0*(dyG(i,j,bi,bj)+dyG(i,j-1,bi,bj))
141  C         by averaging (method II)  C         by averaging (method II)
142  c         dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))  c         dxV(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
143  c         dYU(I,J,bi,bj) = 0.5*(dYC(I,J,bi,bj)+dYC(I-1,J,bi,bj))  c         dyU(i,j,bi,bj) = 0.5*(dyC(i,j,bi,bj)+dyC(i-1,j,bi,bj))
144           ENDDO           ENDDO
145          ENDDO          ENDDO
146    
147  C     Calculate vertical face area  C--     Calculate vertical face area
148          DO J=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
149           DO I=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
150            rA (I,J,bi,bj) = dxF(I,J,bi,bj)*dyF(I,J,bi,bj)            rA (i,j,bi,bj) = dxF(i,j,bi,bj)*dyF(i,j,bi,bj)
151            rAw(I,J,bi,bj) = dxC(I,J,bi,bj)*dyG(I,J,bi,bj)            rAw(i,j,bi,bj) = dxC(i,j,bi,bj)*dyG(i,j,bi,bj)
152            rAs(I,J,bi,bj) = dxG(I,J,bi,bj)*dyC(I,J,bi,bj)            rAs(i,j,bi,bj) = dxG(i,j,bi,bj)*dyC(i,j,bi,bj)
153            rAz(I,J,bi,bj) = dxV(I,J,bi,bj)*dyU(I,J,bi,bj)            rAz(i,j,bi,bj) = dxV(i,j,bi,bj)*dyU(i,j,bi,bj)
154            tanPhiAtU(I,J,bi,bj) = 0.  C--     Set trigonometric terms & grid orientation:
155            tanPhiAtV(I,J,bi,bj) = 0.            tanPhiAtU(i,j,bi,bj) = 0.
156              tanPhiAtV(i,j,bi,bj) = 0.
157              angleCosC(i,j,bi,bj) = 1.
158              angleSinC(i,j,bi,bj) = 0.
159           ENDDO           ENDDO
160          ENDDO          ENDDO
161    
162  C--     Cosine(lat) scaling  C--     Cosine(lat) scaling
163          DO J=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
164           cosFacU(J,bi,bj)=1.           cosFacU(j,bi,bj) = 1.
165           cosFacV(J,bi,bj)=1.           cosFacV(j,bi,bj) = 1.
166           sqcosFacU(J,bi,bj)=1.           sqcosFacU(j,bi,bj)=1.
167           sqcosFacV(J,bi,bj)=1.           sqcosFacV(j,bi,bj)=1.
168          ENDDO          ENDDO
169    
170         ENDDO ! bi  C--   end bi,bj loops
171        ENDDO ! bj         ENDDO
172          ENDDO
173    
174        RETURN        RETURN
175        END        END
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22