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