/[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.1 by cnh, Wed Apr 22 19:15:30 1998 UTC revision 1.16 by adcroft, Tue May 29 14:01:37 2001 UTC
# Line 1  Line 1 
1  C $Id$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  CStartOfInterface  CStartOfInterface
7        SUBROUTINE INI_CARTESIAN_GRID( myThid )        SUBROUTINE INI_CARTESIAN_GRID( myThid )
# Line 33  C     | Under the cartesian grid mode pr Line 34  C     | Under the cartesian grid mode pr
34  C     | and Y are in metres. Disktance in Z are in m or Pa       |  C     | and Y are in metres. Disktance in Z are in m or Pa       |
35  C     | depending on the vertical gridding mode.                 |  C     | depending on the vertical gridding mode.                 |
36  C     \==========================================================/  C     \==========================================================/
37          IMPLICIT NONE
38    
39  C     === Global variables ===  C     === Global variables ===
40  #include "SIZE.h"  #include "SIZE.h"
# Line 46  C     myThid -  Number of this instance Line 48  C     myThid -  Number of this instance
48  CEndOfInterface  CEndOfInterface
49    
50  C     == Local variables ==  C     == Local variables ==
51  C     xG, yG - Global coordinate location.        INTEGER iG, jG, bi, bj, I,  J
52  C     zG        _RL xG0, yG0
53  C     xBase  - South-west corner location for process.  
54  C     yBase  C "Long" real for temporary coordinate calculation
55  C     zUpper - Work arrays for upper and lower  C  NOTICE the extended range of indices!!
56  C     zLower   cell-face heights.        _RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
57  C     phi    - Temporary scalar        _RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
58  C     xBase  - Temporaries for lower corner coordinate  
59  C     yBase  C These functions return the "global" index with valid values beyond
60  C     iG, jG - Global coordinate index. Usually used to hold  C halo regions
61  C              the south-west global coordinate of a tile.        INTEGER iGl,jGl
62  C     bi,bj  - Loop counters        iGl(I,bi) = 1+mod(myXGlobalLo-1+(bi-1)*sNx+I+Olx*Nx-1,Nx)
63  C     zUpper - Temporary arrays holding z coordinates of        jGl(J,bj) = 1+mod(myYGlobalLo-1+(bj-1)*sNy+J+Oly*Ny-1,Ny)
64  C     zLower   upper and lower faces.  
65  C     I,J,K  C     For each tile ...
       _RL    xG, yG, zG  
       _RL    phi  
       _RL    zUpper(Nz), zLower(Nz)  
       _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  
66        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
        jG = myYGlobalLo + (bj-1)*sNy  
67         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
68    
69    C--     "Global" index (place holder)
70            jG = myYGlobalLo + (bj-1)*sNy
71          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 )  
72    
73  C--   Calculate separation between other points  C--   First find coordinate of tile corner (meaning outer corner of halo)
74  C     dxG, dyG are separations between cell corners along cell faces.          xG0 = 0.
75        DO bj = myByLo(myThid), myByHi(myThid)  C       Find the X-coordinate of the outer grid-line of the "real" tile
76         DO bi = myBxLo(myThid), myBxHi(myThid)          DO i=1, iG-1
77          DO J=1,sNy           xG0 = xG0 + delX(i)
78           DO I=1,sNx          ENDDO
79            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
80            dyG(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I-1,J,bi,bj))*0.5 _d 0          DO i=1, Olx
81             xG0 = xG0 - delX( 1+mod(Olx*Nx-1+iG-i,Nx) )
82            ENDDO
83    C       Find the Y-coordinate of the outer grid-line of the "real" tile
84            yG0 = 0.
85            DO j=1, jG-1
86             yG0 = yG0 + delY(j)
87            ENDDO
88    C       Back-step to the outer grid-line of the "halo" region
89            DO j=1, Oly
90             yG0 = yG0 - delY( 1+mod(Oly*Ny-1+jG-j,Ny) )
91            ENDDO
92    
93    C--     Calculate coordinates of cell corners for N+1 grid-lines
94            DO J=1-Oly,sNy+Oly +1
95             xGloc(1-Olx,J) = xG0
96             DO I=1-Olx,sNx+Olx
97    c         xGloc(I+1,J) = xGloc(I,J) + delX(1+mod(Nx-1+iG-1+i,Nx))
98              xGloc(I+1,J) = xGloc(I,J) + delX( iGl(I,bi) )
99           ENDDO           ENDDO
100          ENDDO          ENDDO
101         ENDDO          DO I=1-Olx,sNx+Olx +1
102        ENDDO           yGloc(I,1-Oly) = yG0
103        _EXCH_XY_R4(dxG, myThid )           DO J=1-Oly,sNy+Oly
104        _EXCH_XY_R4(dyG, myThid )  c         yGloc(I,J+1) = yGloc(I,J) + delY(1+mod(Ny-1+jG-1+j,Ny))
105  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  
106           ENDDO           ENDDO
107          ENDDO          ENDDO
108         ENDDO  
109        ENDDO  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
110        _EXCH_XY_R4(dxV, myThid )          DO J=1-Oly,sNy+Oly
111        _EXCH_XY_R4(dyU, myThid )           DO I=1-Olx,sNx+Olx
112  C     dxC, dyC is separation between cell centers            xG(I,J,bi,bj) = xGloc(I,J)
113        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 D0  
           dyC(I,J,bi,bj)    = (dyF(I,J,bi,bj)+dyF(I,J-1,bi,bj))*0.5 D0  
114           ENDDO           ENDDO
115          ENDDO          ENDDO
116         ENDDO  
117        ENDDO  C--     Calculate [xC,yC], coordinates of cell centers
118        _EXCH_XY_R4(dxC, myThid )          DO J=1-Oly,sNy+Oly
119        _EXCH_XY_R4(dyC, myThid )           DO I=1-Olx,sNx+Olx
120  C     Calculate recipricols  C         by averaging
121        DO bj = myByLo(myThid), myByHi(myThid)            xC(I,J,bi,bj) = 0.25*(
122         DO bi = myBxLo(myThid), myBxHi(myThid)       &     xGloc(I,J)+xGloc(I+1,J)+xGloc(I,J+1)+xGloc(I+1,J+1) )
123          DO J=1,sNy            yC(I,J,bi,bj) = 0.25*(
124           DO I=1,sNx       &     yGloc(I,J)+yGloc(I+1,J)+yGloc(I,J+1)+yGloc(I+1,J+1) )
125            rDxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)           ENDDO
126            rDyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)          ENDDO
127            rDxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)  
128            rDyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)  C--     Calculate [dxF,dyF], lengths between cell faces (through center)
129            rDxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)          DO J=1-Oly,sNy+Oly
130            rDyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)           DO I=1-Olx,sNx+Olx
131            rDxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)            dXF(I,J,bi,bj) = delX( iGl(I,bi) )
132            rDyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)            dYF(I,J,bi,bj) = delY( jGl(J,bj) )
133           ENDDO           ENDDO
134          ENDDO          ENDDO
135         ENDDO  
136        ENDDO  C--     Calculate [dxG,dyG], lengths along cell boundaries
137        _EXCH_XY_R4(rDxG, myThid )          DO J=1-Oly,sNy+Oly
138        _EXCH_XY_R4(rDyG, myThid )           DO I=1-Olx,sNx+Olx
139        _EXCH_XY_R4(rDxC, myThid )            dXG(I,J,bi,bj) = delX( iGl(I,bi) )
140        _EXCH_XY_R4(rDyC, myThid )            dYG(I,J,bi,bj) = delY( jGl(J,bj) )
141        _EXCH_XY_R4(rDxF, myThid )           ENDDO
142        _EXCH_XY_R4(rDyF, myThid )          ENDDO
143        _EXCH_XY_R4(rDxV, myThid )  
144        _EXCH_XY_R4(rDyU, myThid )  C--     The following arrays are not defined in some parts of the halo
145    C       region. We set them to zero here for safety. If they are ever
146    C       referred to, especially in the denominator then it is a mistake!
147            DO J=1-Oly,sNy+Oly
148             DO I=1-Olx,sNx+Olx
149              dXC(I,J,bi,bj) = 0.
150              dYC(I,J,bi,bj) = 0.
151              dXV(I,J,bi,bj) = 0.
152              dYU(I,J,bi,bj) = 0.
153              rAw(I,J,bi,bj) = 0.
154              rAs(I,J,bi,bj) = 0.
155             ENDDO
156            ENDDO
157    
158    C--     Calculate [dxC], zonal length between cell centers
159            DO J=1-Oly,sNy+Oly
160             DO I=1-Olx+1,sNx+Olx ! NOTE range
161              dXC(I,J,bi,bj) = 0.5D0*(dXF(I,J,bi,bj)+dXF(I-1,J,bi,bj))
162             ENDDO
163            ENDDO
164    
165    C--     Calculate [dyC], meridional length between cell centers
166            DO J=1-Oly+1,sNy+Oly ! NOTE range
167             DO I=1-Olx,sNx+Olx
168              dYC(I,J,bi,bj) = 0.5*(dYF(I,J,bi,bj)+dYF(I,J-1,bi,bj))
169             ENDDO
170            ENDDO
171    
172    C--     Calculate [dxV,dyU], length between velocity points (through corners)
173            DO J=1-Oly+1,sNy+Oly ! NOTE range
174             DO I=1-Olx+1,sNx+Olx ! NOTE range
175    C         by averaging (method I)
176              dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
177              dYU(I,J,bi,bj) = 0.5*(dYG(I,J,bi,bj)+dYG(I,J-1,bi,bj))
178    C         by averaging (method II)
179    c         dXV(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I-1,J,bi,bj))
180    c         dYU(I,J,bi,bj) = 0.5*(dYC(I,J,bi,bj)+dYC(I-1,J,bi,bj))
181             ENDDO
182            ENDDO
183    
184  C     Calculate vertical face area  C     Calculate vertical face area
185        DO bj = myByLo(myThid), myByHi(myThid)          DO J=1-Oly,sNy+Oly
186         DO bi = myBxLo(myThid), myBxHi(myThid)           DO I=1-Olx,sNx+Olx
187          DO J=1,sNy            rA (I,J,bi,bj) = dxF(I,J,bi,bj)*dyF(I,J,bi,bj)
188           DO I=1,sNx            rAw(I,J,bi,bj) = dxC(I,J,bi,bj)*dyG(I,J,bi,bj)
189            zA(I,J,bi,bj) = dxF(I,J,bi,bj)*dyF(I,J,bi,bj)            rAs(I,J,bi,bj) = dxG(I,J,bi,bj)*dyC(I,J,bi,bj)
190              rAz(I,J,bi,bj) = dxV(I,J,bi,bj)*dyU(I,J,bi,bj)
191              tanPhiAtU(I,J,bi,bj) = 0.
192              tanPhiAtV(I,J,bi,bj) = 0.
193           ENDDO           ENDDO
194          ENDDO          ENDDO
        ENDDO  
       ENDDO  
195    
196        DO bj = myByLo(myThid), myByHi(myThid)  C--     Cosine(lat) scaling
197         DO bi = myBxLo(myThid), myBxHi(myThid)          DO J=1-OLy,sNy+OLy
198          DO K=1,Nz           cosFacU(J,bi,bj)=1.
199           DO J=1,sNy           cosFacV(J,bi,bj)=1.
200            DO I=1,sNx           sqcosFacU(J,bi,bj)=1.
201             IF (HFacC(I,J,K,bi,bj) .NE. 0. D0 ) THEN           sqcosFacV(J,bi,bj)=1.
202              rHFacC(I,J,K,bi,bj) = 1. D0 / HFacC(I,J,K,bi,bj)          ENDDO
203             ELSE  
204              rHFacC(I,J,K,bi,bj) = 0. D0         ENDDO ! bi
205             ENDIF        ENDDO ! bj
            IF (HFacW(I,J,K,bi,bj) .NE. 0. D0 ) THEN  
             rHFacW(I,J,K,bi,bj) = 1. D0 / HFacW(I,J,K,bi,bj)  
             maskW(I,J,K,bi,bj) = 1. D0  
            ELSE  
             rHFacW(I,J,K,bi,bj) = 0. D0  
             maskW(I,J,K,bi,bj) = 0.0 D0  
            ENDIF  
            IF (HFacS(I,J,K,bi,bj) .NE. 0. D0 ) THEN  
             rHFacS(I,J,K,bi,bj) = 1. D0 / HFacS(I,J,K,bi,bj)  
             maskS(I,J,K,bi,bj) = 1. D0  
            ELSE  
             rHFacS(I,J,K,bi,bj) = 0. D0  
             maskS(I,J,K,bi,bj) = 0. D0  
            ENDIF  
           ENDDO  
          ENDDO  
         ENDDO  
        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 )  
206    
 C  
207        RETURN        RETURN
208        END        END
   
 C     $Id$  
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22