/[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.17 by cnh, Sun Feb 4 16:46:44 2001 UTC revision 1.27 by jmc, Sat Apr 17 18:25:12 2010 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    c#include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  #undef  USE_BACKWARD_COMPATIBLE_GRID  #undef  USE_BACKWARD_COMPATIBLE_GRID
8    
9  CStartOfInterface  CBOP
10    C     !ROUTINE: INI_SPHERICAL_POLAR_GRID
11    C     !INTERFACE:
12        SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )        SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )
 C     /==========================================================\  
 C     | SUBROUTINE INI_SPHERICAL_POLAR_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 spherical polar grid mode primitive distances  |  
 C     | in X and Y are in degrees. Distance in Z are in m or Pa  |  
 C     | depending on the vertical gridding mode.                 |  
 C     \==========================================================/  
       IMPLICIT NONE  
13    
14    C     !DESCRIPTION: \bv
15    C     *==========================================================*
16    C     | SUBROUTINE INI_SPHERICAL_POLAR_GRID
17    C     | o Initialise model coordinate system arrays
18    C     *==========================================================*
19    C     | These arrays are used throughout the code in evaluating
20    C     | gradients, integrals and spatial avarages. This routine
21    C     | is called separately by each thread and initialise only
22    C     | the region of the domain it is "responsible" for.
23    C     | Under the spherical polar grid mode primitive distances
24    C     | in X and Y are in degrees. Distance in Z are in m or Pa
25    C     | depending on the vertical gridding mode.
26    C     *==========================================================*
27    C     \ev
28    
29    C     !USES:
30          IMPLICIT NONE
31  C     === Global variables ===  C     === Global variables ===
32  #include "SIZE.h"  #include "SIZE.h"
33  #include "EEPARAMS.h"  #include "EEPARAMS.h"
34  #include "PARAMS.h"  #include "PARAMS.h"
35  #include "GRID.h"  #include "GRID.h"
36    c#ifdef ALLOW_EXCH2
37    c#include "W2_EXCH2_SIZE.h"
38    c#include "W2_EXCH2_TOPOLOGY.h"
39    c#include "W2_EXCH2_PARAMS.h"
40    c#endif /* ALLOW_EXCH2 */
41    
42    C     !INPUT/OUTPUT PARAMETERS:
43  C     == Routine arguments ==  C     == Routine arguments ==
44  C     myThid -  Number of this instance of INI_CARTESIAN_GRID  C     myThid  :: my Thread Id Number
45        INTEGER myThid        INTEGER myThid
 CEndOfInterface  
46    
47    C     !LOCAL VARIABLES:
48  C     == Local variables ==  C     == Local variables ==
49  C     xG, yG - Global coordinate location.  C     xG0,yG0 :: coordinate of South-West tile-corner
50  C     xBase  - South-west corner location for process.  C     iG, jG  :: Global coordinate index. Usually used to hold
51  C     yBase  C             :: the south-west global coordinate of a tile.
52  C     zUpper - Work arrays for upper and lower  C     lat     :: Temporary variables used to hold latitude values.
53  C     zLower   cell-face heights.  C     bi,bj   :: tile indices
54  C     phi    - Temporary scalar  C     i, j    :: loop counters
 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     xBase  - Lower coordinate for this threads cells  
 C     yBase  
 C     lat, latN, - Temporary variables used to hold latitude  
 C     latS         values.  
 C     I,J,K  
55        INTEGER iG, jG        INTEGER iG, jG
56        INTEGER bi, bj        INTEGER bi, bj
57        INTEGER  I,  J        INTEGER i,  j
58        _RL lat, dlat, dlon, xG0, yG0        _RL lat, dlat, dlon, xG0, yG0
59    
60    C     "Long" real for temporary coordinate calculation
61  C "Long" real for temporary coordinate calculation  C      NOTICE the extended range of indices!!
 C  NOTICE the extended range of indices!!  
62        _RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)        _RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
63        _RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)        _RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
64    
65  C These functions return the "global" index with valid values beyond  C     The functions iGl, jGl return the "global" index with valid values beyond
66  C halo regions  C     halo regions
67    C     cnh wrote:
68    C     >    I dont understand why we would ever have to multiply the
69    C     >    overlap by the total domain size e.g
70    C     >    OLx*Nx, OLy*Ny.
71    C     >    Can anybody explain? Lines are in ini_spherical_polar_grid.F.
72    C     >    Surprised the code works if its wrong, so I am puzzled.
73    C     jmc replied:
74    C     Yes, I can explain this since I put this modification to work
75    C     with small domain (where Oly > Ny, as for instance, zonal-average
76    C     case):
77    C     This has no effect on the acuracy of the evaluation of iGl(I,bi)
78    C     and jGl(j,bj) since we take mod(a+OLx*Nx,Nx) and mod(b+OLy*Ny,Ny).
79    C     But in case a or b is negative, then the FORTRAN function "mod"
80    C     does not return the matematical value of the "modulus" function,
81    C     and this is not good for your purpose.
82    C     This is why I add +OLx*Nx and +OLy*Ny to be sure that the 1rst
83    C     argument of the mod function is positive.
84        INTEGER iGl,jGl        INTEGER iGl,jGl
85        iGl(I,bi) = 1+mod(myXGlobalLo-1+(bi-1)*sNx+I+Olx*Nx-1,Nx)        iGl(i,bi) = 1+MOD(myXGlobalLo-1+(bi-1)*sNx+i+Olx*Nx-1,Nx)
86        jGl(J,bj) = 1+mod(myYGlobalLo-1+(bj-1)*sNy+J+Oly*Ny-1,Ny)        jGl(j,bj) = 1+MOD(myYGlobalLo-1+(bj-1)*sNy+j+Oly*Ny-1,Ny)
87    c#ifdef ALLOW_EXCH2
88    c      INTEGER tN
89    c#endif /* ALLOW_EXCH2 */
90    CEOP
91    
92  C     For each tile ...  C     For each tile ...
93        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
# Line 90  C     For each tile ... Line 96  C     For each tile ...
96  C--     "Global" index (place holder)  C--     "Global" index (place holder)
97          jG = myYGlobalLo + (bj-1)*sNy          jG = myYGlobalLo + (bj-1)*sNy
98          iG = myXGlobalLo + (bi-1)*sNx          iG = myXGlobalLo + (bi-1)*sNx
99    c#ifdef ALLOW_EXCH2
100    c        IF ( W2_useE2ioLayOut ) THEN
101    cC- note: does not work for non-uniform delX or delY
102    c          tN = W2_myTileList(bi,bj)
103    c          iG = exch2_txGlobalo(tN)
104    c          jG = exch2_tyGlobalo(tN)
105    c        ENDIF
106    c#endif /* ALLOW_EXCH2 */
107    
108  C--   First find coordinate of tile corner (meaning outer corner of halo)  C--   First find coordinate of tile corner (meaning outer corner of halo)
109          xG0 = thetaMin          xG0 = xgOrigin
110  C       Find the X-coordinate of the outer grid-line of the "real" tile  C       Find the X-coordinate of the outer grid-line of the "real" tile
111          DO i=1, iG-1          DO i=1, iG-1
112           xG0 = xG0 + delX(i)           xG0 = xG0 + delX(i)
113          ENDDO          ENDDO
114  C       Back-step to the outer grid-line of the "halo" region  C       Back-step to the outer grid-line of the "halo" region
115          DO i=1, Olx          DO i=1, Olx
116           xG0 = xG0 - delX( 1+mod(Olx*Nx-1+iG-i,Nx) )           xG0 = xG0 - delX( 1+MOD(Olx*Nx-1+iG-i,Nx) )
117          ENDDO          ENDDO
118  C       Find the Y-coordinate of the outer grid-line of the "real" tile  C       Find the Y-coordinate of the outer grid-line of the "real" tile
119          yG0 = phiMin          yG0 = ygOrigin
120          DO j=1, jG-1          DO j=1, jG-1
121           yG0 = yG0 + delY(j)           yG0 = yG0 + delY(j)
122          ENDDO          ENDDO
123  C       Back-step to the outer grid-line of the "halo" region  C       Back-step to the outer grid-line of the "halo" region
124          DO j=1, Oly          DO j=1, Oly
125           yG0 = yG0 - delY( 1+mod(Oly*Ny-1+jG-j,Ny) )           yG0 = yG0 - delY( 1+MOD(Oly*Ny-1+jG-j,Ny) )
126          ENDDO          ENDDO
127    
128  C--     Calculate coordinates of cell corners for N+1 grid-lines  C--     Calculate coordinates of cell corners for N+1 grid-lines
129          DO J=1-Oly,sNy+Oly +1          DO j=1-Oly,sNy+Oly +1
130           xGloc(1-Olx,J) = xG0           xGloc(1-Olx,j) = xG0
131           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
132  c         xGloc(I+1,J) = xGloc(I,J) + delX(1+mod(Nx-1+iG-1+i,Nx))  c         xGloc(i+1,j) = xGloc(i,j) + delX(1+mod(Nx-1+iG-1+i,Nx))
133            xGloc(I+1,J) = xGloc(I,J) + delX( iGl(I,bi) )            xGloc(i+1,j) = xGloc(i,j) + delX( iGl(i,bi) )
134           ENDDO           ENDDO
135          ENDDO          ENDDO
136          DO I=1-Olx,sNx+Olx +1          DO i=1-Olx,sNx+Olx +1
137           yGloc(I,1-Oly) = yG0           yGloc(i,1-Oly) = yG0
138           DO J=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
139  c         yGloc(I,J+1) = yGloc(I,J) + delY(1+mod(Ny-1+jG-1+j,Ny))  c         yGloc(i,j+1) = yGloc(i,j) + delY(1+mod(Ny-1+jG-1+j,Ny))
140            yGloc(I,J+1) = yGloc(I,J) + delY( jGl(J,bj) )            yGloc(i,j+1) = yGloc(i,j) + delY( jGl(j,bj) )
141           ENDDO           ENDDO
142          ENDDO          ENDDO
143    
144  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
145          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
146           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
147            xG(I,J,bi,bj) = xGloc(I,J)            xG(i,j,bi,bj) = xGloc(i,j)
148            yG(I,J,bi,bj) = yGloc(I,J)            yG(i,j,bi,bj) = yGloc(i,j)
149           ENDDO           ENDDO
150          ENDDO          ENDDO
151    
152  C--     Calculate [xC,yC], coordinates of cell centers  C--     Calculate [xC,yC], coordinates of cell centers
153          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
154           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
155  C         by averaging  C         by averaging
156            xC(I,J,bi,bj) = 0.25*(            xC(i,j,bi,bj) = 0.25 _d 0*(
157       &     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) )
158            yC(I,J,bi,bj) = 0.25*(            yC(i,j,bi,bj) = 0.25 _d 0*(
159       &     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) )
160           ENDDO           ENDDO
161          ENDDO          ENDDO
162    
163  C--     Calculate [dxF,dyF], lengths between cell faces (through center)  C--     Calculate [dxF,dyF], lengths between cell faces (through center)
164          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
165           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
166  C         by averaging  C         by averaging
167  c         dXF(I,J,bi,bj) = 0.5*(dXG(I,J,bi,bj)+dXG(I,J+1,bi,bj))  c         dxF(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i,j+1,bi,bj))
168  c         dYF(I,J,bi,bj) = 0.5*(dYG(I,J,bi,bj)+dYG(I+1,J,bi,bj))  c         dyF(i,j,bi,bj) = 0.5*(dyG(i,j,bi,bj)+dyG(i+1,j,bi,bj))
169  C         by formula  C         by formula
170            lat = yC(I,J,bi,bj)            lat = yC(i,j,bi,bj)
171            dlon = delX( iGl(I,bi) )            dlon = delX( iGl(i,bi) )
172            dlat = delY( jGl(J,bj) )            dlat = delY( jGl(j,bj) )
173            dXF(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad            dxF(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
174  #ifdef    USE_BACKWARD_COMPATIBLE_GRID  #ifdef    USE_BACKWARD_COMPATIBLE_GRID
175            dXF(I,J,bi,bj) = delX(iGl(I,bi))*deg2rad*rSphere*            dxF(i,j,bi,bj) = delX(iGl(i,bi))*deg2rad*rSphere*
176       &                     COS(yc(I,J,bi,bj)*deg2rad)       &                     COS(yC(i,j,bi,bj)*deg2rad)
177  #endif    /* USE_BACKWARD_COMPATIBLE_GRID */  #endif    /* USE_BACKWARD_COMPATIBLE_GRID */
178            dYF(I,J,bi,bj) = rSphere*dlat*deg2rad            dyF(i,j,bi,bj) = rSphere*dlat*deg2rad
179           ENDDO           ENDDO
180          ENDDO          ENDDO
181    
182  C--     Calculate [dxG,dyG], lengths along cell boundaries  C--     Calculate [dxG,dyG], lengths along cell boundaries
183          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
184           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
185  C         by averaging  C         by averaging
186  c         dXG(I,J,bi,bj) = 0.5*(dXF(I,J,bi,bj)+dXF(I,J-1,bi,bj))  c         dxG(i,j,bi,bj) = 0.5*(dxF(i,j,bi,bj)+dxF(i,j-1,bi,bj))
187  c         dYG(I,J,bi,bj) = 0.5*(dYF(I,J,bi,bj)+dYF(I-1,J,bi,bj))  c         dyG(i,j,bi,bj) = 0.5*(dyF(i,j,bi,bj)+dyF(i-1,j,bi,bj))
188  C         by formula  C         by formula
189            lat = 0.5*(yGloc(I,J)+yGloc(I+1,J))            lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
190            dlon = delX( iGl(I,bi) )            dlon = delX( iGl(i,bi) )
191            dlat = delY( jGl(J,bj) )            dlat = delY( jGl(j,bj) )
192            dXG(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad            dxG(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
193            dYG(I,J,bi,bj) = rSphere*dlat*deg2rad            if (dxG(i,j,bi,bj).LT.1.) dxG(i,j,bi,bj)=0.
194              dyG(i,j,bi,bj) = rSphere*dlat*deg2rad
195           ENDDO           ENDDO
196          ENDDO          ENDDO
197    
198  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
199  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
200  C       referred to, especially in the denominator then it is a mistake!  C       referred to, especially in the denominator then it is a mistake!
201          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
202           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
203            dXC(I,J,bi,bj) = 0.            dxC(i,j,bi,bj) = 0.
204            dYC(I,J,bi,bj) = 0.            dyC(i,j,bi,bj) = 0.
205            dXV(I,J,bi,bj) = 0.            dxV(i,j,bi,bj) = 0.
206            dYU(I,J,bi,bj) = 0.            dyU(i,j,bi,bj) = 0.
207            rAw(I,J,bi,bj) = 0.            rAw(i,j,bi,bj) = 0.
208            rAs(I,J,bi,bj) = 0.            rAs(i,j,bi,bj) = 0.
209           ENDDO           ENDDO
210          ENDDO          ENDDO
211    
212  C--     Calculate [dxC], zonal length between cell centers  C--     Calculate [dxC], zonal length between cell centers
213          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
214           DO I=1-Olx+1,sNx+Olx ! NOTE range           DO i=1-Olx+1,sNx+Olx ! NOTE range
215  C         by averaging  C         by averaging
216            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))
217  C         by formula  C         by formula
218  c         lat = 0.5*(yC(I,J,bi,bj)+yC(I-1,J,bi,bj))  c         lat = 0.5*(yC(i,j,bi,bj)+yC(i-1,j,bi,bj))
219  c         dlon = 0.5*(delX( iGl(I,bi) ) + delX( iGl(I-1,bi) ))  c         dlon = 0.5*(delX( iGl(i,bi) ) + delX( iGl(i-1,bi) ))
220  c         dXC(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad  c         dxC(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
221  C         by difference  C         by difference
222  c         lat = 0.5*(yC(I,J,bi,bj)+yC(I-1,J,bi,bj))  c         lat = 0.5*(yC(i,j,bi,bj)+yC(i-1,j,bi,bj))
223  c         dlon = (xC(I,J,bi,bj)+xC(I-1,J,bi,bj))  c         dlon = (xC(i,j,bi,bj)+xC(i-1,j,bi,bj))
224  c         dXC(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad  c         dxC(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
225           ENDDO           ENDDO
226          ENDDO          ENDDO
227    
228  C--     Calculate [dyC], meridional length between cell centers  C--     Calculate [dyC], meridional length between cell centers
229          DO J=1-Oly+1,sNy+Oly ! NOTE range          DO j=1-Oly+1,sNy+Oly ! NOTE range
230           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
231  C         by averaging  C         by averaging
232            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))
233  C         by formula  C         by formula
234  c         dlat = 0.5*(delY( jGl(J,bj) ) + delY( jGl(J-1,bj) ))  c         dlat = 0.5*(delY( jGl(j,bj) ) + delY( jGl(j-1,bj) ))
235  c         dYC(I,J,bi,bj) = rSphere*dlat*deg2rad  c         dyC(i,j,bi,bj) = rSphere*dlat*deg2rad
236  C         by difference  C         by difference
237  c         dlat = (yC(I,J,bi,bj)+yC(I,J-1,bi,bj))  c         dlat = (yC(i,j,bi,bj)+yC(i,j-1,bi,bj))
238  c         dYC(I,J,bi,bj) = rSphere*dlat*deg2rad  c         dyC(i,j,bi,bj) = rSphere*dlat*deg2rad
239           ENDDO           ENDDO
240          ENDDO          ENDDO
241    
242  C--     Calculate [dxV,dyU], length between velocity points (through corners)  C--     Calculate [dxV,dyU], length between velocity points (through corners)
243          DO J=1-Oly+1,sNy+Oly ! NOTE range          DO j=1-Oly+1,sNy+Oly ! NOTE range
244           DO I=1-Olx+1,sNx+Olx ! NOTE range           DO i=1-Olx+1,sNx+Olx ! NOTE range
245  C         by averaging (method I)  C         by averaging (method I)
246            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))
247            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))
248  C         by averaging (method II)  C         by averaging (method II)
249  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))
250  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))
251           ENDDO           ENDDO
252          ENDDO          ENDDO
253    
254  C--     Calculate vertical face area (tracer cells)  C--     Calculate vertical face area (tracer cells)
255          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
256           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
257            lat=0.5*(yGloc(I,J)+yGloc(I+1,J))            lat=0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
258            dlon=delX( iGl(I,bi) )            dlon=delX( iGl(i,bi) )
259            dlat=delY( jGl(J,bj) )            dlat=delY( jGl(j,bj) )
260            rA(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad            rA(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
261       &        *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )       &        *ABS( SIN((lat+dlat)*deg2rad)-SIN(lat*deg2rad) )
262  #ifdef    USE_BACKWARD_COMPATIBLE_GRID  #ifdef    USE_BACKWARD_COMPATIBLE_GRID
263            lat=yC(I,J,bi,bj)-delY( jGl(J,bj) )*0.5 _d 0            lat=yC(i,j,bi,bj)-delY( jGl(j,bj) )*0.5 _d 0
264            lon=yC(I,J,bi,bj)+delY( jGl(J,bj) )*0.5 _d 0            lon=yC(i,j,bi,bj)+delY( jGl(j,bj) )*0.5 _d 0
265            rA(I,J,bi,bj) = dyF(I,J,bi,bj)            rA(i,j,bi,bj) = dyF(i,j,bi,bj)
266       &    *rSphere*(SIN(lon*deg2rad)-SIN(lat*deg2rad))       &    *rSphere*(SIN(lon*deg2rad)-SIN(lat*deg2rad))
267  #endif    /* USE_BACKWARD_COMPATIBLE_GRID */  #endif    /* USE_BACKWARD_COMPATIBLE_GRID */
268           ENDDO           ENDDO
269          ENDDO          ENDDO
270    
271  C--     Calculate vertical face area (u cells)  C--     Calculate vertical face area (u cells)
272          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
273           DO I=1-Olx+1,sNx+Olx ! NOTE range           DO i=1-Olx+1,sNx+Olx ! NOTE range
274  C         by averaging  C         by averaging
275            rAw(I,J,bi,bj) = 0.5*(rA(I,J,bi,bj)+rA(I-1,J,bi,bj))            rAw(i,j,bi,bj) = 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj))
276  C         by formula  C         by formula
277  c         lat=yGloc(I,J)  c         lat=yGloc(i,j)
278  c         dlon=0.5*( delX( iGl(I,bi) ) + delX( iGl(I-1,bi) ) )  c         dlon=0.5*( delX( iGl(i,bi) ) + delX( iGl(i-1,bi) ) )
279  c         dlat=delY( jGl(J,bj) )  c         dlat=delY( jGl(j,bj) )
280  c         rAw(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad  c         rAw(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
281  c    &        *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )  c    &        *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )
282           ENDDO           ENDDO
283          ENDDO          ENDDO
284    
285  C--     Calculate vertical face area (v cells)  C--     Calculate vertical face area (v cells)
286          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
287           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
288  C         by formula  C         by formula
289            lat=yC(I,J,bi,bj)            lat=yC(i,j,bi,bj)
290            dlon=delX( iGl(I,bi) )            dlon=delX( iGl(i,bi) )
291            dlat=0.5*( delY( jGl(J,bj) ) + delY( jGl(J-1,bj) ) )            dlat=0.5 _d 0*( delY( jGl(j,bj) ) + delY( jGl(j-1,bj) ) )
292            rAs(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad            rAs(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
293       &        *abs( sin(lat*deg2rad)-sin((lat-dlat)*deg2rad) )       &        *ABS( SIN(lat*deg2rad)-SIN((lat-dlat)*deg2rad) )
294  #ifdef    USE_BACKWARD_COMPATIBLE_GRID  #ifdef    USE_BACKWARD_COMPATIBLE_GRID
295            lon=yC(I,J,bi,bj)-delY( jGl(J,bj) )            lon=yC(i,j,bi,bj)-delY( jGl(j,bj) )
296            lat=yC(I,J,bi,bj)            lat=yC(i,j,bi,bj)
297            rAs(I,J,bi,bj) = rSphere*delX(iGl(I,bi))*deg2rad            rAs(i,j,bi,bj) = rSphere*delX(iGl(i,bi))*deg2rad
298       &    *rSphere*(SIN(lat*deg2rad)-SIN(lon*deg2rad))       &    *rSphere*(SIN(lat*deg2rad)-SIN(lon*deg2rad))
299  #endif    /* USE_BACKWARD_COMPATIBLE_GRID */  #endif    /* USE_BACKWARD_COMPATIBLE_GRID */
300            IF (abs(lat).GT.90..OR.abs(lat-dlat).GT.90.) rAs(I,J,bi,bj)=0.            IF (ABS(lat).GT.90..OR.ABS(lat-dlat).GT.90.) rAs(i,j,bi,bj)=0.
301           ENDDO           ENDDO
302          ENDDO          ENDDO
303    
304  C--     Calculate vertical face area (vorticity points)  C--     Calculate vertical face area (vorticity points)
305          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
306           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
307  C         by formula  C         by formula
308            lat=yC(I,J,bi,bj)            lat =0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
309            dlon=delX( iGl(I,bi) )            dlon=0.5 _d 0*( delX( iGl(i,bi) ) + delX( iGl(i-1,bi) ) )
310            dlat=0.5*( delY( jGl(J,bj) ) + delY( jGl(J-1,bj) ) )            dlat=0.5 _d 0*( delY( jGl(j,bj) ) + delY( jGl(j-1,bj) ) )
311            rAz(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad            rAz(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
312       &     *abs( sin(lat*deg2rad)-sin((lat-dlat)*deg2rad) )       &     *ABS( SIN(lat*deg2rad)-SIN((lat-dlat)*deg2rad) )
313            IF (abs(lat).GT.90..OR.abs(lat-dlat).GT.90.) rAz(I,J,bi,bj)=0.            IF (ABS(lat).GT.90..OR.ABS(lat-dlat).GT.90.) rAz(i,j,bi,bj)=0.
314           ENDDO           ENDDO
315          ENDDO          ENDDO
316    
317  C--     Calculate trigonometric terms  C--     Calculate trigonometric terms & grid orientation:
318          DO J=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
319           DO I=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
320            lat=0.5*(yGloc(I,J)+yGloc(I,J+1))            lat=0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
321            tanPhiAtU(i,j,bi,bj)=tan(lat*deg2rad)            tanPhiAtU(i,j,bi,bj)=TAN(lat*deg2rad)
322            lat=0.5*(yGloc(I,J)+yGloc(I+1,J))            lat=0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
323            tanPhiAtV(i,j,bi,bj)=tan(lat*deg2rad)            tanPhiAtV(i,j,bi,bj)=TAN(lat*deg2rad)
324              angleCosC(i,j,bi,bj) = 1.
325              angleSinC(i,j,bi,bj) = 0.
326           ENDDO           ENDDO
327          ENDDO          ENDDO
328    
329    C--     Cosine(lat) scaling
330            DO j=1-OLy,sNy+OLy
331             jG = myYGlobalLo + (bj-1)*sNy + j-1
332             jG = MIN(MAX(1,jG),Ny)
333             IF (cosPower.NE.0.) THEN
334              cosFacU(j,bi,bj)=COS(yC(1,j,bi,bj)*deg2rad)
335         &                    **cosPower
336              cosFacV(j,bi,bj)=COS((yC(1,j,bi,bj)-0.5*delY(jG))*deg2rad)
337         &                    **cosPower
338              cosFacU(j,bi,bj)=ABS(cosFacU(j,bi,bj))
339              cosFacV(j,bi,bj)=ABS(cosFacV(j,bi,bj))
340              sqcosFacU(j,bi,bj)=SQRT(cosFacU(j,bi,bj))
341              sqcosFacV(j,bi,bj)=SQRT(cosFacV(j,bi,bj))
342             ELSE
343              cosFacU(j,bi,bj)=1.
344              cosFacV(j,bi,bj)=1.
345              sqcosFacU(j,bi,bj)=1.
346              sqcosFacV(j,bi,bj)=1.
347             ENDIF
348            ENDDO
349    
350         ENDDO ! bi         ENDDO ! bi
351        ENDDO ! bj        ENDDO ! bj
352    
353          IF ( rotateGrid ) THEN
354           CALL ROTATE_SPHERICAL_POLAR_GRID( xC, yC, myThid )
355           CALL ROTATE_SPHERICAL_POLAR_GRID( xG, yG, myThid )
356           CALL CALC_ANGLES( myThid )
357          ENDIF
358    
359        RETURN        RETURN
360        END        END
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22