/[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.4 by cnh, Mon May 25 18:01:32 1998 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$
3    
4  #include "CPP_EEOPTIONS.h"  c#include "PACKAGES_CONFIG.h"
5    #include "CPP_OPTIONS.h"
6    
7  CStartOfInterface  #undef  USE_BACKWARD_COMPATIBLE_GRID
8    
9    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     \==========================================================/  
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     zG  C     iG, jG  :: Global coordinate index. Usually used to hold
51  C     xBase  - South-west corner location for process.  C             :: the south-west global coordinate of a tile.
52  C     yBase  C     lat     :: Temporary variables used to hold latitude values.
53  C     zUpper - Work arrays for upper and lower  C     bi,bj   :: tile indices
54  C     zLower   cell-face heights.  C     i, j    :: loop counters
 C     phi    - Temporary scalar  
 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  
       _RL    xG, yG, zG  
       _RL    phi  
       _RL    zUpper(Nz), zLower(Nz)  
       _RL    xBase, yBase  
55        INTEGER iG, jG        INTEGER iG, jG
56        INTEGER bi, bj        INTEGER bi, bj
57        INTEGER  I,  J, K        INTEGER i,  j
58        _RL lat, latS, latN        _RL lat, dlat, dlon, xG0, yG0
59    
60    C     "Long" real for temporary coordinate calculation
61    C      NOTICE the extended range of indices!!
62          _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)
64    
65    C     The functions iGl, jGl return the "global" index with valid values beyond
66    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
85          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)
87    c#ifdef ALLOW_EXCH2
88    c      INTEGER tN
89    c#endif /* ALLOW_EXCH2 */
90    CEOP
91    
92  C--   Example of inialisation for spherical polar 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  
 C     Note: In the spherical polar case delX and delY are given in  
 C           degrees and are relative to some starting latitude and  
 C           longitude - phiMin and thetaMin.  
93        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
        jG = myYGlobalLo + (bj-1)*sNy  
94         DO bi = myBxLo(myThid), myBxHi(myThid)         DO bi = myBxLo(myThid), myBxHi(myThid)
95    
96    C--     "Global" index (place holder)
97            jG = myYGlobalLo + (bj-1)*sNy
98          iG = myXGlobalLo + (bi-1)*sNx          iG = myXGlobalLo + (bi-1)*sNx
99          yBase = phiMin  c#ifdef ALLOW_EXCH2
100          xBase = thetaMin  c        IF ( W2_useE2ioLayOut ) THEN
101          DO i=1,iG-1  cC- note: does not work for non-uniform delX or delY
102           xBase = xBase + delX(i)  c          tN = W2_myTileList(bi,bj)
103          ENDDO  c          iG = exch2_txGlobalo(tN)
104          DO j=1,jG-1  c          jG = exch2_tyGlobalo(tN)
105           yBase = yBase + delY(j)  c        ENDIF
106          ENDDO  c#endif /* ALLOW_EXCH2 */
         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)*deg2rad*rSphere*COS(yc(I,J,bi,bj)*deg2rad)  
           dyF(I,J,bi,bj) = delY(jG+j-1)*deg2rad*rSphere  
          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 )  
107    
108  C--   Calculate separation between other points  C--   First find coordinate of tile corner (meaning outer corner of halo)
109  C     dxG, dyG are separations between cell corners along cell faces.          xG0 = xgOrigin
110        DO bj = myByLo(myThid), myByHi(myThid)  C       Find the X-coordinate of the outer grid-line of the "real" tile
111         DO bi = myBxLo(myThid), myBxHi(myThid)          DO i=1, iG-1
112          DO J=1,sNy           xG0 = xG0 + delX(i)
113           DO I=1,sNx          ENDDO
114            jG = myYGlobalLo + (bj-1)*sNy + J-1  C       Back-step to the outer grid-line of the "halo" region
115            iG = myXGlobalLo + (bi-1)*sNx + I-1          DO i=1, Olx
116            lat = yc(I,J,bi,bj)-delY(jG) * 0.5 _d 0           xG0 = xG0 - delX( 1+MOD(Olx*Nx-1+iG-i,Nx) )
117            dxG(I,J,bi,bj) = rSphere*COS(lat*deg2rad)*delX(iG)*deg2rad          ENDDO
118            dyG(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I-1,J,bi,bj))*0.5 _d 0  C       Find the Y-coordinate of the outer grid-line of the "real" tile
119            yG0 = ygOrigin
120            DO j=1, jG-1
121             yG0 = yG0 + delY(j)
122            ENDDO
123    C       Back-step to the outer grid-line of the "halo" region
124            DO j=1, Oly
125             yG0 = yG0 - delY( 1+MOD(Oly*Ny-1+jG-j,Ny) )
126            ENDDO
127    
128    C--     Calculate coordinates of cell corners for N+1 grid-lines
129            DO j=1-Oly,sNy+Oly +1
130             xGloc(1-Olx,j) = xG0
131             DO i=1-Olx,sNx+Olx
132    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) )
134           ENDDO           ENDDO
135          ENDDO          ENDDO
136         ENDDO          DO i=1-Olx,sNx+Olx +1
137        ENDDO           yGloc(i,1-Oly) = yG0
138        _EXCH_XY_R4(dxG, myThid )           DO j=1-Oly,sNy+Oly
139        _EXCH_XY_R4(dyG, myThid )  c         yGloc(i,j+1) = yGloc(i,j) + delY(1+mod(Ny-1+jG-1+j,Ny))
140  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  
141           ENDDO           ENDDO
142          ENDDO          ENDDO
143         ENDDO  
144        ENDDO  C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
145        _EXCH_XY_R4(dxV, myThid )          DO j=1-Oly,sNy+Oly
146        _EXCH_XY_R4(dyU, myThid )           DO i=1-Olx,sNx+Olx
147  C     dxC, dyC is separation between cell centers            xG(i,j,bi,bj) = xGloc(i,j)
148        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  
149           ENDDO           ENDDO
150          ENDDO          ENDDO
151         ENDDO  
152        ENDDO  C--     Calculate [xC,yC], coordinates of cell centers
153        _EXCH_XY_R4(dxC, myThid )          DO j=1-Oly,sNy+Oly
154        _EXCH_XY_R4(dyC, myThid )           DO i=1-Olx,sNx+Olx
155  C     Calculate recipricols  C         by averaging
156        DO bj = myByLo(myThid), myByHi(myThid)            xC(i,j,bi,bj) = 0.25 _d 0*(
157         DO bi = myBxLo(myThid), myBxHi(myThid)       &     xGloc(i,j)+xGloc(i+1,j)+xGloc(i,j+1)+xGloc(i+1,j+1) )
158          DO J=1,sNy            yC(i,j,bi,bj) = 0.25 _d 0*(
159           DO I=1,sNx       &     yGloc(i,j)+yGloc(i+1,j)+yGloc(i,j+1)+yGloc(i+1,j+1) )
160            rDxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj)           ENDDO
161            rDyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj)          ENDDO
162            rDxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj)  
163            rDyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj)  C--     Calculate [dxF,dyF], lengths between cell faces (through center)
164            rDxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj)          DO j=1-Oly,sNy+Oly
165            rDyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj)           DO i=1-Olx,sNx+Olx
166            rDxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj)  C         by averaging
167            rDyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj)  c         dxF(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i,j+1,bi,bj))
168           ENDDO  c         dyF(i,j,bi,bj) = 0.5*(dyG(i,j,bi,bj)+dyG(i+1,j,bi,bj))
169          ENDDO  C         by formula
170         ENDDO            lat = yC(i,j,bi,bj)
171        ENDDO            dlon = delX( iGl(i,bi) )
172        _EXCH_XY_R4(rDxG, myThid )            dlat = delY( jGl(j,bj) )
173        _EXCH_XY_R4(rDyG, myThid )            dxF(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
174        _EXCH_XY_R4(rDxC, myThid )  #ifdef    USE_BACKWARD_COMPATIBLE_GRID
175        _EXCH_XY_R4(rDyC, myThid )            dxF(i,j,bi,bj) = delX(iGl(i,bi))*deg2rad*rSphere*
176        _EXCH_XY_R4(rDxF, myThid )       &                     COS(yC(i,j,bi,bj)*deg2rad)
177        _EXCH_XY_R4(rDyF, myThid )  #endif    /* USE_BACKWARD_COMPATIBLE_GRID */
178        _EXCH_XY_R4(rDxV, myThid )            dyF(i,j,bi,bj) = rSphere*dlat*deg2rad
179        _EXCH_XY_R4(rDyU, myThid )           ENDDO
180  C     Calculate vertical face area          ENDDO
181        DO bj = myByLo(myThid), myByHi(myThid)  
182         DO bi = myBxLo(myThid), myBxHi(myThid)  C--     Calculate [dxG,dyG], lengths along cell boundaries
183          DO J=1,sNy          DO j=1-Oly,sNy+Oly
184           DO I=1,sNx           DO i=1-Olx,sNx+Olx
185            jG = myYGlobalLo + (bj-1)*sNy + J-1  C         by averaging
186            latS = yc(i,j,bi,bj)-delY(jG)*0.5 _d 0  c         dxG(i,j,bi,bj) = 0.5*(dxF(i,j,bi,bj)+dxF(i,j-1,bi,bj))
187            latN = yc(i,j,bi,bj)+delY(jG)*0.5 _d 0  c         dyG(i,j,bi,bj) = 0.5*(dyF(i,j,bi,bj)+dyF(i-1,j,bi,bj))
188            zA(I,J,bi,bj) = dyF(I,J,bi,bj)  C         by formula
189       &    *rSphere*(SIN(latN*deg2rad)-SIN(latS*deg2rad))            lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
190              dlon = delX( iGl(i,bi) )
191              dlat = delY( jGl(j,bj) )
192              dxG(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
193              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
196            ENDDO
197    
198    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
200    C       referred to, especially in the denominator then it is a mistake!
201            DO j=1-Oly,sNy+Oly
202             DO i=1-Olx,sNx+Olx
203              dxC(i,j,bi,bj) = 0.
204              dyC(i,j,bi,bj) = 0.
205              dxV(i,j,bi,bj) = 0.
206              dyU(i,j,bi,bj) = 0.
207              rAw(i,j,bi,bj) = 0.
208              rAs(i,j,bi,bj) = 0.
209             ENDDO
210            ENDDO
211    
212    C--     Calculate [dxC], zonal length between cell centers
213            DO j=1-Oly,sNy+Oly
214             DO i=1-Olx+1,sNx+Olx ! NOTE range
215    C         by averaging
216              dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
217    C         by formula
218    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) ))
220    c         dxC(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
221    C         by difference
222    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))
224    c         dxC(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
225             ENDDO
226            ENDDO
227    
228    C--     Calculate [dyC], meridional length between cell centers
229            DO j=1-Oly+1,sNy+Oly ! NOTE range
230             DO i=1-Olx,sNx+Olx
231    C         by averaging
232              dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
233    C         by formula
234    c         dlat = 0.5*(delY( jGl(j,bj) ) + delY( jGl(j-1,bj) ))
235    c         dyC(i,j,bi,bj) = rSphere*dlat*deg2rad
236    C         by difference
237    c         dlat = (yC(i,j,bi,bj)+yC(i,j-1,bi,bj))
238    c         dyC(i,j,bi,bj) = rSphere*dlat*deg2rad
239             ENDDO
240            ENDDO
241    
242    C--     Calculate [dxV,dyU], length between velocity points (through corners)
243            DO j=1-Oly+1,sNy+Oly ! NOTE range
244             DO i=1-Olx+1,sNx+Olx ! NOTE range
245    C         by averaging (method I)
246              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 _d 0*(dyG(i,j,bi,bj)+dyG(i,j-1,bi,bj))
248    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))
250    c         dyU(i,j,bi,bj) = 0.5*(dyC(i,j,bi,bj)+dyC(i-1,j,bi,bj))
251             ENDDO
252            ENDDO
253    
254    C--     Calculate vertical face area (tracer cells)
255            DO j=1-Oly,sNy+Oly
256             DO i=1-Olx,sNx+Olx
257              lat=0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
258              dlon=delX( iGl(i,bi) )
259              dlat=delY( jGl(j,bj) )
260              rA(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
261         &        *ABS( SIN((lat+dlat)*deg2rad)-SIN(lat*deg2rad) )
262    #ifdef    USE_BACKWARD_COMPATIBLE_GRID
263              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
265              rA(i,j,bi,bj) = dyF(i,j,bi,bj)
266         &    *rSphere*(SIN(lon*deg2rad)-SIN(lat*deg2rad))
267    #endif    /* USE_BACKWARD_COMPATIBLE_GRID */
268             ENDDO
269            ENDDO
270    
271    C--     Calculate vertical face area (u cells)
272            DO j=1-Oly,sNy+Oly
273             DO i=1-Olx+1,sNx+Olx ! NOTE range
274    C         by averaging
275              rAw(i,j,bi,bj) = 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj))
276    C         by formula
277    c         lat=yGloc(i,j)
278    c         dlon=0.5*( delX( iGl(i,bi) ) + delX( iGl(i-1,bi) ) )
279    c         dlat=delY( jGl(j,bj) )
280    c         rAw(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
281    c    &        *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )
282             ENDDO
283            ENDDO
284    
285    C--     Calculate vertical face area (v cells)
286            DO j=1-Oly,sNy+Oly
287             DO i=1-Olx,sNx+Olx
288    C         by formula
289              lat=yC(i,j,bi,bj)
290              dlon=delX( iGl(i,bi) )
291              dlat=0.5 _d 0*( delY( jGl(j,bj) ) + delY( jGl(j-1,bj) ) )
292              rAs(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
293         &        *ABS( SIN(lat*deg2rad)-SIN((lat-dlat)*deg2rad) )
294    #ifdef    USE_BACKWARD_COMPATIBLE_GRID
295              lon=yC(i,j,bi,bj)-delY( jGl(j,bj) )
296              lat=yC(i,j,bi,bj)
297              rAs(i,j,bi,bj) = rSphere*delX(iGl(i,bi))*deg2rad
298         &    *rSphere*(SIN(lat*deg2rad)-SIN(lon*deg2rad))
299    #endif    /* USE_BACKWARD_COMPATIBLE_GRID */
300              IF (ABS(lat).GT.90..OR.ABS(lat-dlat).GT.90.) rAs(i,j,bi,bj)=0.
301           ENDDO           ENDDO
302          ENDDO          ENDDO
303         ENDDO  
304        ENDDO  C--     Calculate vertical face area (vorticity points)
305        _EXCH_XY_R4(zA, myThid )          DO j=1-Oly,sNy+Oly
306  C           DO i=1-Olx,sNx+Olx
307    C         by formula
308              lat =0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
309              dlon=0.5 _d 0*( delX( iGl(i,bi) ) + delX( iGl(i-1,bi) ) )
310              dlat=0.5 _d 0*( delY( jGl(j,bj) ) + delY( jGl(j-1,bj) ) )
311              rAz(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
312         &     *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.
314             ENDDO
315            ENDDO
316    
317    C--     Calculate trigonometric terms & grid orientation:
318            DO j=1-Oly,sNy+Oly
319             DO i=1-Olx,sNx+Olx
320              lat=0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
321              tanPhiAtU(i,j,bi,bj)=TAN(lat*deg2rad)
322              lat=0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
323              tanPhiAtV(i,j,bi,bj)=TAN(lat*deg2rad)
324              angleCosC(i,j,bi,bj) = 1.
325              angleSinC(i,j,bi,bj) = 0.
326             ENDDO
327            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
351          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.4  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22