/[MITgcm]/MITgcm/model/src/ini_spherical_polar_grid.F
ViewVC logotype

Diff of /MITgcm/model/src/ini_spherical_polar_grid.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.3 by cnh, Fri Apr 24 02:10:20 1998 UTC revision 1.14 by adcroft, Mon Mar 27 22:25:44 2000 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    
3  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
4    
5  CStartOfInterface  CStartOfInterface
6        SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )        SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )
# Line 33  C     | Under the spherical polar grid m Line 33  C     | Under the spherical polar grid m
33  C     | in X and Y are in degrees. Distance in Z are in m or Pa  |  C     | in X and Y are in degrees. Distance in Z are in m or Pa  |
34  C     | depending on the vertical gridding mode.                 |  C     | depending on the vertical gridding mode.                 |
35  C     \==========================================================/  C     \==========================================================/
36          IMPLICIT NONE
37    
38  C     === Global variables ===  C     === Global variables ===
39  #include "SIZE.h"  #include "SIZE.h"
# Line 47  CEndOfInterface Line 48  CEndOfInterface
48    
49  C     == Local variables ==  C     == Local variables ==
50  C     xG, yG - Global coordinate location.  C     xG, yG - Global coordinate location.
 C     zG  
51  C     xBase  - South-west corner location for process.  C     xBase  - South-west corner location for process.
52  C     yBase  C     yBase
 C     zUpper - Work arrays for upper and lower  
 C     zLower   cell-face heights.  
 C     phi    - Temporary scalar  
53  C     iG, jG - Global coordinate index. Usually used to hold  C     iG, jG - Global coordinate index. Usually used to hold
54  C              the south-west global coordinate of a tile.  C              the south-west global coordinate of a tile.
55  C     bi,bj  - Loop counters  C     bi,bj  - Loop counters
# Line 63  C     yBase Line 60  C     yBase
60  C     lat, latN, - Temporary variables used to hold latitude  C     lat, latN, - Temporary variables used to hold latitude
61  C     latS         values.  C     latS         values.
62  C     I,J,K  C     I,J,K
63        _RL    xG, yG, zG        _RL    xG, yG
       _RL    phi  
       _RL    zUpper(Nz), zLower(Nz)  
64        _RL    xBase, yBase        _RL    xBase, yBase
65        INTEGER iG, jG        INTEGER iG, jG
66        INTEGER bi, bj        INTEGER bi, bj
67        INTEGER  I,  J, K        INTEGER  I,  J
68        _RL lat, latS, latN        _RL lat, latS, latN
69    
70  C--   Example of inialisation for spherical polar grid  C--   Example of inialisation for spherical polar grid
# Line 81  C     Set up my local grid first Line 76  C     Set up my local grid first
76  C     Note: In the spherical polar case delX and delY are given in  C     Note: In the spherical polar case delX and delY are given in
77  C           degrees and are relative to some starting latitude and  C           degrees and are relative to some starting latitude and
78  C           longitude - phiMin and thetaMin.  C           longitude - phiMin and thetaMin.
79          xC0 = thetaMin
80          yC0 = phiMin
81        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
82         jG = myYGlobalLo + (bj-1)*sNy         jG = myYGlobalLo + (bj-1)*sNy
83         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
# Line 100  C           longitude - phiMin and theta Line 97  C           longitude - phiMin and theta
97            xc(I,J,bi,bj)  = xG + delX(iG+i-1)*0.5 _d 0            xc(I,J,bi,bj)  = xG + delX(iG+i-1)*0.5 _d 0
98            yc(I,J,bi,bj)  = yG + delY(jG+j-1)*0.5 _d 0            yc(I,J,bi,bj)  = yG + delY(jG+j-1)*0.5 _d 0
99            xG = xG + delX(iG+I-1)            xG = xG + delX(iG+I-1)
100            dxF(I,J,bi,bj) = delX(iG+i-1)*deg2rad*rSphere*COS(yc(I,J,bi,bj)*deg2rad)            dxF(I,J,bi,bj) = delX(iG+i-1)*deg2rad
101         &    *rSphere*COS(yc(I,J,bi,bj)*deg2rad)
102            dyF(I,J,bi,bj) = delY(jG+j-1)*deg2rad*rSphere            dyF(I,J,bi,bj) = delY(jG+j-1)*deg2rad*rSphere
103           ENDDO           ENDDO
104           yG = yG + delY(jG+J-1)           yG = yG + delY(jG+J-1)
# Line 132  C     dxG, dyG are separations between c Line 130  C     dxG, dyG are separations between c
130        ENDDO        ENDDO
131        _EXCH_XY_R4(dxG, myThid )        _EXCH_XY_R4(dxG, myThid )
132        _EXCH_XY_R4(dyG, myThid )        _EXCH_XY_R4(dyG, myThid )
133  C     dxV, dyU are separations between velocity points along cell faces.  C     dxC, dyC is separation between cell centers
134        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
135         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
136          DO J=1,sNy          DO J=1,sNy
137           DO I=1,sNx           DO I=1,sNx
138            dxV(I,J,bi,bj) = (dxG(I,J,bi,bj)+dxG(I-1,J,bi,bj))*0.5 _d 0            dxC(I,J,bi,bj)    = (dxF(I,J,bi,bj)+dxF(I-1,J,bi,bj))*0.5 _d 0
139            dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,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
140           ENDDO           ENDDO
141          ENDDO          ENDDO
142         ENDDO         ENDDO
143        ENDDO        ENDDO
144        _EXCH_XY_R4(dxV, myThid )        _EXCH_XY_R4(dxC, myThid )
145        _EXCH_XY_R4(dyU, myThid )        _EXCH_XY_R4(dyC, myThid )
146  C     dxC, dyC is separation between cell centers  C     dxV, dyU are separations between velocity points along cell faces.
147        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
148         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
149          DO J=1,sNy          DO J=1,sNy
150           DO I=1,sNx           DO I=1,sNx
151            dxC(I,J,bi,bj)    = (dxF(I,J,bi,bj)+dxF(I-1,J,bi,bj))*0.5 _d 0            dxV(I,J,bi,bj) = (dxG(I,J,bi,bj)+dxG(I-1,J,bi,bj))*0.5 _d 0
152            dyC(I,J,bi,bj)    = (dyF(I,J,bi,bj)+dyF(I,J-1,bi,bj))*0.5 _d 0  #ifdef OLD_UV_GEOMETRY
153              dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,bi,bj))*0.5 _d 0
154    #else
155              dyU(I,J,bi,bj) = (dyC(I,J,bi,bj)+dyC(I-1,J,bi,bj))*0.5 _d 0
156    #endif
157           ENDDO           ENDDO
158          ENDDO          ENDDO
159         ENDDO         ENDDO
160        ENDDO        ENDDO
161        _EXCH_XY_R4(dxC, myThid )        _EXCH_XY_R4(dxV, myThid )
162        _EXCH_XY_R4(dyC, myThid )        _EXCH_XY_R4(dyU, myThid )
163  C     Calculate recipricols  C     Calculate vertical face area and trigonometric terms
164        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
165         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
166          DO J=1,sNy          DO J=1,sNy
167           DO I=1,sNx           DO I=1,sNx
168            rDxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)            jG = myYGlobalLo + (bj-1)*sNy + J-1
169            rDyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)            iG = myXGlobalLo + (bi-1)*sNx + I-1
170            rDxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)            latS = yc(i,j,bi,bj)-delY(jG)*0.5 _d 0
171            rDyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)            latN = yc(i,j,bi,bj)+delY(jG)*0.5 _d 0
172            rDxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)  #ifdef OLD_UV_GEOMETRY
173            rDyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)            rA(I,J,bi,bj) = dyF(I,J,bi,bj)
174            rDxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)       &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))
175            rDyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)  #else
176              rA(I,J,bi,bj) = rSphere*delX(iG)*deg2rad
177         &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))
178    #endif
179    C         Area cannot be zero but sin can be if lat if < -90.
180              IF ( rA(I,J,bi,bj) .LT. 0. ) rA(I,J,bi,bj) = -rA(I,J,bi,bj)
181              tanPhiAtU(i,j,bi,bj)=tan(_yC(i,j,bi,bj)*deg2rad)
182              tanPhiAtV(i,j,bi,bj)=tan(latS*deg2rad)
183           ENDDO           ENDDO
184          ENDDO          ENDDO
185         ENDDO         ENDDO
186        ENDDO        ENDDO
187        _EXCH_XY_R4(rDxG, myThid )        _EXCH_XY_R4 (rA       , myThid )
188        _EXCH_XY_R4(rDyG, myThid )        _EXCH_XY_R4 (tanPhiAtU , myThid )
189        _EXCH_XY_R4(rDxC, myThid )        _EXCH_XY_R4 (tanPhiAtV , myThid )
       _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  
190        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
191         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
192          DO J=1,sNy          DO J=1,sNy
193           DO I=1,sNx           DO I=1,sNx
194              iG = myXGlobalLo + (bi-1)*sNx + I-1
195            jG = myYGlobalLo + (bj-1)*sNy + J-1            jG = myYGlobalLo + (bj-1)*sNy + J-1
196            latS = yc(i,j,bi,bj)-delY(jG)*0.5 _d 0            latS = yc(i,j-1,bi,bj)
197            latN = yc(i,j,bi,bj)+delY(jG)*0.5 _d 0            latN = yc(i,j,bi,bj)
198            zA(I,J,bi,bj) = dyF(I,J,bi,bj)  #ifdef OLD_UV_GEOMETRY
199              rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))
200              rAs(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I,J-1,bi,bj))
201    #else
202              rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))
203              rAs(I,J,bi,bj) = rSphere*delX(iG)*deg2rad
204       &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))       &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))
205    #endif
206           ENDDO           ENDDO
207          ENDDO          ENDDO
208         ENDDO         ENDDO
209        ENDDO        ENDDO
210          _EXCH_XY_R4 (rAw      , myThid )
211          _EXCH_XY_R4 (rAs      , myThid )
212  C  C
213        RETURN        RETURN
214        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22