/[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.5 by cnh, Mon Jun 8 21:43:01 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_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 )
 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(Nz), zLower(Nz)  
       _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 D0  
           dyC(I,J,bi,bj)    = (dyF(I,J,bi,bj)+dyF(I,J-1,bi,bj))*0.5 D0  
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 recipricols          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            rDxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)            dyU(i,j,bi,bj) = 0.
120            rDyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)            rAw(i,j,bi,bj) = 0.
121            rDxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)            rAs(i,j,bi,bj) = 0.
           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)  
122           ENDDO           ENDDO
123          ENDDO          ENDDO
124         ENDDO  
125        ENDDO  C--     Calculate [dxC], zonal length between cell centers
126        _EXCH_XY_R4(rDxG, myThid )          DO j=1-OLy,sNy+OLy
127        _EXCH_XY_R4(rDyG, myThid )           DO i=1-OLx+1,sNx+OLx ! NOTE range
128        _EXCH_XY_R4(rDxC, myThid )            dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
       _EXCH_XY_R4(rDyC, myThid )  
       _EXCH_XY_R4(rDxF, myThid )  
       _EXCH_XY_R4(rDyF, myThid )  
       _EXCH_XY_R4(rDxV, myThid )  
       _EXCH_XY_R4(rDyU, myThid )  
 C     Calculate vertical face area  
       DO bj = myByLo(myThid), myByHi(myThid)  
        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)  
129           ENDDO           ENDDO
130          ENDDO          ENDDO
        ENDDO  
       ENDDO  
131    
132        DO bj = myByLo(myThid), myByHi(myThid)  C--     Calculate [dyC], meridional length between cell centers
133         DO bi = myBxLo(myThid), myBxHi(myThid)          DO j=1-OLy+1,sNy+OLy ! NOTE range
134          DO K=1,Nz           DO i=1-OLx,sNx+OLx
135           DO J=1,sNy            dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
136            DO I=1,sNx           ENDDO
137             IF (HFacC(I,J,K,bi,bj) .NE. 0. D0 ) THEN          ENDDO
138              rHFacC(I,J,K,bi,bj) = 1. D0 / HFacC(I,J,K,bi,bj)  
139             ELSE  C--     Calculate [dxV,dyU], length between velocity points (through corners)
140              rHFacC(I,J,K,bi,bj) = 0. D0          DO j=1-OLy+1,sNy+OLy ! NOTE range
141             ENDIF           DO i=1-OLx+1,sNx+OLx ! NOTE range
142             IF (HFacW(I,J,K,bi,bj) .NE. 0. D0 ) THEN  C         by averaging (method I)
143              rHFacW(I,J,K,bi,bj) = 1. D0 / HFacW(I,J,K,bi,bj)            dxV(i,j,bi,bj) = 0.5 _d 0*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
144              maskW(I,J,K,bi,bj) = 1. D0            dyU(i,j,bi,bj) = 0.5 _d 0*(dyG(i,j,bi,bj)+dyG(i,j-1,bi,bj))
145             ELSE  C         by averaging (method II)
146              rHFacW(I,J,K,bi,bj) = 0. D0  c         dxV(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
147              maskW(I,J,K,bi,bj) = 0.0 D0  c         dyU(i,j,bi,bj) = 0.5*(dyC(i,j,bi,bj)+dyC(i-1,j,bi,bj))
148             ENDIF           ENDDO
149             IF (HFacS(I,J,K,bi,bj) .NE. 0. D0 ) THEN          ENDDO
150              rHFacS(I,J,K,bi,bj) = 1. D0 / HFacS(I,J,K,bi,bj)  
151              maskS(I,J,K,bi,bj) = 1. D0  C--     Calculate vertical face area
152             ELSE          DO j=1-OLy,sNy+OLy
153              rHFacS(I,J,K,bi,bj) = 0. D0           DO i=1-OLx,sNx+OLx
154              maskS(I,J,K,bi,bj) = 0. D0            rA (i,j,bi,bj) = dxF(i,j,bi,bj)*dyF(i,j,bi,bj)
155             ENDIF            rAw(i,j,bi,bj) = dxC(i,j,bi,bj)*dyG(i,j,bi,bj)
156            ENDDO            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           ENDDO
164          ENDDO          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
 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 )  
   
 CcnhDebugStarts  
       tanPhiAtU = 0. _d 0  
       tanPhiAtV = 0. _d 0  
       _EXCH_XY_R4 (tanPhiAtU , myThid )  
       _EXCH_XY_R4 (tanPhiAtV , myThid )  
 CcnhDebugEnds  
177    
 C  
178        RETURN        RETURN
179        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22