/[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.9 by cnh, Mon Nov 2 03:34:11 1998 UTC revision 1.20 by jmc, Tue Oct 17 18:52:34 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.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 )
10  C     /==========================================================\  C     !DESCRIPTION: \bv
11  C     | SUBROUTINE INI_CARTESIAN_GRID                            |  C     *==========================================================*
12  C     | o Initialise model coordinate system                     |  C     | SUBROUTINE INI_CARTESIAN_GRID                            
13  C     |==========================================================|  C     | o Initialise model coordinate system                      
14  C     | These arrays are used throughout the code in evaluating  |  C     *==========================================================*
15  C     | gradients, integrals and spatial avarages. This routine  |  C     | The grid arrays, initialised here, are used throughout
16  C     | is called separately by each thread and initialise only  |  C     | the code in evaluating gradients, integrals and spatial
17  C     | the region of the domain it is "responsible" for.        |  C     | avarages. This routine  
18  C     | Notes:                                                   |  C     | is called separately by each thread and initialises only  
19  C     | Two examples are included. One illustrates the           |  C     | the region of the domain it is "responsible" for.        
20  C     | initialisation of a cartesian grid. The other shows the  |  C     | Notes:                                                    
21  C     | inialisation of a spherical polar grid. Other orthonormal|  C     | Two examples are included. One illustrates the            
22  C     | grids can be fitted into this design. In this case       |  C     | initialisation of a cartesian grid (this routine).
23  C     | custom metric terms also need adding to account for the  |  C     | The other shows the  
24  C     | projections of velocity vectors onto these grids.        |  C     | inialisation of a spherical polar grid. Other orthonormal
25  C     | The structure used here also makes it possible to        |  C     | grids can be fitted into this design. In this case        
26  C     | implement less regular grid mappings. In particular      |  C     | custom metric terms also need adding to account for the  
27  C     | o Schemes which leave out blocks of the domain that are  |  C     | projections of velocity vectors onto these grids.        
28  C     |   all land could be supported.                           |  C     | The structure used here also makes it possible to        
29  C     | o Multi-level schemes such as icosohedral or cubic       |  C     | implement less regular grid mappings. In particular      
30  C     |   grid projections onto a sphere can also be fitted      |  C     | o Schemes which leave out blocks of the domain that are  
31  C     |   within the strategy we use.                            |  C     |   all land could be supported.                            
32  C     |   Both of the above also require modifying the support   |  C     | o Multi-level schemes such as icosohedral or cubic        
33  C     |   routines that map computational blocks to simulation   |  C     |   grid projections onto a sphere can also be fitted      
34  C     |   domain blocks.                                         |  C     |   within the strategy we use.                            
35  C     | Under the cartesian grid mode primitive distances in X   |  C     |   Both of the above also require modifying the support    
36  C     | and Y are in metres. Disktance in Z are in m or Pa       |  C     |   routines that map computational blocks to simulation    
37  C     | depending on the vertical gridding mode.                 |  C     |   domain blocks.                                          
38  C     \==========================================================/  C     | Under the cartesian grid mode primitive distances in X    
39    C     | and Y are in metres. Disktance in Z are in m or Pa        
40    C     | depending on the vertical gridding mode.                  
41    C     *==========================================================*
42    C     \ev
43    
44    C     !USES:
45          IMPLICIT NONE
46  C     === Global variables ===  C     === Global variables ===
47  #include "SIZE.h"  #include "SIZE.h"
48  #include "EEPARAMS.h"  #include "EEPARAMS.h"
49  #include "PARAMS.h"  #include "PARAMS.h"
50  #include "GRID.h"  #include "GRID.h"
51    
52    C     !INPUT/OUTPUT PARAMETERS:
53  C     == Routine arguments ==  C     == Routine arguments ==
54  C     myThid -  Number of this instance of INI_CARTESIAN_GRID  C     myThid -  Number of this instance of INI_CARTESIAN_GRID
55        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
56    
57    C     !LOCAL VARIABLES:
58  C     == Local variables ==  C     == Local variables ==
59  C     xG, yG - Global coordinate location.        INTEGER iG, jG, bi, bj, I,  J
60  C     zG        _RL xG0, yG0
61  C     xBase  - South-west corner location for process.  C     "Long" real for temporary coordinate calculation
62  C     yBase  C      NOTICE the extended range of indices!!
63  C     zUpper - Work arrays for upper and lower        _RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
64  C     zLower   cell-face heights.        _RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
65  C     phi    - Temporary scalar  C     These functions return the "global" index with valid values beyond
66  C     xBase  - Temporaries for lower corner coordinate  C     halo regions
67  C     yBase        INTEGER iGl,jGl
68  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)
69  C              the south-west global coordinate of a tile.        jGl(J,bj) = 1+mod(myYGlobalLo-1+(bj-1)*sNy+J+Oly*Ny-1,Ny)
70  C     bi,bj  - Loop counters  CEOP
71  C     zUpper - Temporary arrays holding z coordinates of  
72  C     zLower   upper and lower faces.  C     For each tile ...
 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  
   
 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  
       xC0 = 0. _d 0  
       yC0 = 0. _d 0  
73        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
        jG = myYGlobalLo + (bj-1)*sNy  
74         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
75    
76    C--     "Global" index (place holder)
77            jG = myYGlobalLo + (bj-1)*sNy
78          iG = myXGlobalLo + (bi-1)*sNx          iG = myXGlobalLo + (bi-1)*sNx
         yBase = 0. _d 0  
         xBase = 0. _d 0  
         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)  
           dyF(I,J,bi,bj) = delY(jG+j-1)  
          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 )  
79    
80  C--   Calculate separation between other points  C--   First find coordinate of tile corner (meaning outer corner of halo)
81  C     dxG, dyG are separations between cell corners along cell faces.          xG0 = 0.
82        DO bj = myByLo(myThid), myByHi(myThid)  C       Find the X-coordinate of the outer grid-line of the "real" tile
83         DO bi = myBxLo(myThid), myBxHi(myThid)          DO i=1, iG-1
84          DO J=1,sNy           xG0 = xG0 + delX(i)
85           DO I=1,sNx          ENDDO
86            dxG(I,J,bi,bj) = (dxF(I,J,bi,bj)+dxF(I,J-1,bi,bj))*0.5 _d 0  C       Back-step to the outer grid-line of the "halo" region
87            dyG(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I-1,J,bi,bj))*0.5 _d 0          DO i=1, Olx
88             xG0 = xG0 - delX( 1+mod(Olx*Nx-1+iG-i,Nx) )
89            ENDDO
90    C       Find the Y-coordinate of the outer grid-line of the "real" tile
91            yG0 = 0.
92            DO j=1, jG-1
93             yG0 = yG0 + delY(j)
94            ENDDO
95    C       Back-step to the outer grid-line of the "halo" region
96            DO j=1, Oly
97             yG0 = yG0 - delY( 1+mod(Oly*Ny-1+jG-j,Ny) )
98            ENDDO
99    
100    C--     Calculate coordinates of cell corners for N+1 grid-lines
101            DO J=1-Oly,sNy+Oly +1
102             xGloc(1-Olx,J) = xG0
103             DO I=1-Olx,sNx+Olx
104    c         xGloc(I+1,J) = xGloc(I,J) + delX(1+mod(Nx-1+iG-1+i,Nx))
105              xGloc(I+1,J) = xGloc(I,J) + delX( iGl(I,bi) )
106           ENDDO           ENDDO
107          ENDDO          ENDDO
108         ENDDO          DO I=1-Olx,sNx+Olx +1
109        ENDDO           yGloc(I,1-Oly) = yG0
110        _EXCH_XY_R4(dxG, myThid )           DO J=1-Oly,sNy+Oly
111        _EXCH_XY_R4(dyG, myThid )  c         yGloc(I,J+1) = yGloc(I,J) + delY(1+mod(Ny-1+jG-1+j,Ny))
112  C     dxV, dyU are separations between velocity points along cell faces.            yGloc(I,J+1) = yGloc(I,J) + delY( jGl(J,bj) )
       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  
           dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,bi,bj))*0.5 _d 0  
113           ENDDO           ENDDO
114          ENDDO          ENDDO
115         ENDDO  
116        ENDDO  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
117        _EXCH_XY_R4(dxV, myThid )          DO J=1-Oly,sNy+Oly
118        _EXCH_XY_R4(dyU, myThid )           DO I=1-Olx,sNx+Olx
119  C     dxC, dyC is separation between cell centers            xG(I,J,bi,bj) = xGloc(I,J)
120        DO bj = myByLo(myThid), myByHi(myThid)            yG(I,J,bi,bj) = yGloc(I,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  
121           ENDDO           ENDDO
122          ENDDO          ENDDO
123         ENDDO  
124        ENDDO  C--     Calculate [xC,yC], coordinates of cell centers
125        _EXCH_XY_R4(dxC, myThid )          DO J=1-Oly,sNy+Oly
126        _EXCH_XY_R4(dyC, myThid )           DO I=1-Olx,sNx+Olx
127  C     Calculate vertical face area  C         by averaging
128        DO bj = myByLo(myThid), myByHi(myThid)            xC(I,J,bi,bj) = 0.25*(
129         DO bi = myBxLo(myThid), myBxHi(myThid)       &     xGloc(I,J)+xGloc(I+1,J)+xGloc(I,J+1)+xGloc(I+1,J+1) )
130          DO J=1,sNy            yC(I,J,bi,bj) = 0.25*(
131           DO I=1,sNx       &     yGloc(I,J)+yGloc(I+1,J)+yGloc(I,J+1)+yGloc(I+1,J+1) )
132            rA(I,J,bi,bj) = dxF(I,J,bi,bj)*dyF(I,J,bi,bj)           ENDDO
133            tanPhiAtU(I,J,bi,bj) = 0. _d 0          ENDDO
134            tanPhiAtV(I,J,bi,bj) = 0. _d 0  
135    C--     Calculate [dxF,dyF], lengths between cell faces (through center)
136            DO J=1-Oly,sNy+Oly
137             DO I=1-Olx,sNx+Olx
138              dXF(I,J,bi,bj) = delX( iGl(I,bi) )
139              dYF(I,J,bi,bj) = delY( jGl(J,bj) )
140             ENDDO
141            ENDDO
142    
143    C--     Calculate [dxG,dyG], lengths along cell boundaries
144            DO J=1-Oly,sNy+Oly
145             DO I=1-Olx,sNx+Olx
146              dXG(I,J,bi,bj) = delX( iGl(I,bi) )
147              dYG(I,J,bi,bj) = delY( jGl(J,bj) )
148           ENDDO           ENDDO
149          ENDDO          ENDDO
        ENDDO  
       ENDDO  
       _EXCH_XY_R4 (rA       , myThid )  
       _EXCH_XY_R4 (tanPhiAtU , myThid )  
       _EXCH_XY_R4 (tanPhiAtV , myThid )  
150    
151  C  C--     The following arrays are not defined in some parts of the halo
152    C       region. We set them to zero here for safety. If they are ever
153    C       referred to, especially in the denominator then it is a mistake!
154            DO J=1-Oly,sNy+Oly
155             DO I=1-Olx,sNx+Olx
156              dXC(I,J,bi,bj) = 0.
157              dYC(I,J,bi,bj) = 0.
158              dXV(I,J,bi,bj) = 0.
159              dYU(I,J,bi,bj) = 0.
160              rAw(I,J,bi,bj) = 0.
161              rAs(I,J,bi,bj) = 0.
162             ENDDO
163            ENDDO
164    
165    C--     Calculate [dxC], zonal length between cell centers
166            DO J=1-Oly,sNy+Oly
167             DO I=1-Olx+1,sNx+Olx ! NOTE range
168              dXC(I,J,bi,bj) = 0.5D0*(dXF(I,J,bi,bj)+dXF(I-1,J,bi,bj))
169             ENDDO
170            ENDDO
171    
172    C--     Calculate [dyC], meridional length between cell centers
173            DO J=1-Oly+1,sNy+Oly ! NOTE range
174             DO I=1-Olx,sNx+Olx
175              dYC(I,J,bi,bj) = 0.5*(dYF(I,J,bi,bj)+dYF(I,J-1,bi,bj))
176             ENDDO
177            ENDDO
178    
179    C--     Calculate [dxV,dyU], length between velocity points (through corners)
180            DO J=1-Oly+1,sNy+Oly ! NOTE range
181             DO I=1-Olx+1,sNx+Olx ! NOTE range
182    C         by averaging (method I)
183              dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
184              dYU(I,J,bi,bj) = 0.5*(dYG(I,J,bi,bj)+dYG(I,J-1,bi,bj))
185    C         by averaging (method II)
186    c         dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
187    c         dYU(I,J,bi,bj) = 0.5*(dYC(I,J,bi,bj)+dYC(I-1,J,bi,bj))
188             ENDDO
189            ENDDO
190    
191    C--     Calculate vertical face area
192            DO J=1-Oly,sNy+Oly
193             DO I=1-Olx,sNx+Olx
194              rA (I,J,bi,bj) = dxF(I,J,bi,bj)*dyF(I,J,bi,bj)
195              rAw(I,J,bi,bj) = dxC(I,J,bi,bj)*dyG(I,J,bi,bj)
196              rAs(I,J,bi,bj) = dxG(I,J,bi,bj)*dyC(I,J,bi,bj)
197              rAz(I,J,bi,bj) = dxV(I,J,bi,bj)*dyU(I,J,bi,bj)
198    C--     Set trigonometric terms & grid orientation:
199              tanPhiAtU(I,J,bi,bj) = 0.
200              tanPhiAtV(I,J,bi,bj) = 0.
201              angleCosC(I,J,bi,bj) = 1.
202              angleSinC(I,J,bi,bj) = 0.
203             ENDDO
204            ENDDO
205    
206    C--     Cosine(lat) scaling
207            DO J=1-OLy,sNy+OLy
208             cosFacU(J,bi,bj)=1.
209             cosFacV(J,bi,bj)=1.
210             sqcosFacU(J,bi,bj)=1.
211             sqcosFacV(J,bi,bj)=1.
212            ENDDO
213    
214           ENDDO ! bi
215          ENDDO ! bj
216    
217    C--   Set default (=whole domain) for where relaxation to climatology applies
218          _BEGIN_MASTER(myThid)
219          IF ( latBandClimRelax.EQ.UNSET_RL ) THEN
220            latBandClimRelax = 0.
221            DO j=1,Ny
222              latBandClimRelax = latBandClimRelax + delY(j)
223            ENDDO
224            latBandClimRelax = latBandClimRelax*3. _d 0
225          ENDIF
226          _END_MASTER(myThid)
227    
228        RETURN        RETURN
229        END        END
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22