/[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.15 by adcroft, Fri Feb 2 21:04:48 2001 UTC revision 1.23 by jmc, Wed Jul 13 00:36:01 2005 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6  #undef  USE_BACKWARD_COMPATIBLE_GRID  #undef  USE_BACKWARD_COMPATIBLE_GRID
7    
8  CStartOfInterface  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        IMPLICIT NONE  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.
44  C     xBase  - South-west corner location for process.  C     xBase  - South-west corner location for process.
# Line 71  C     I,J,K Line 62  C     I,J,K
62        _RL lat, dlat, dlon, xG0, yG0        _RL lat, dlat, dlon, xG0, yG0
63    
64    
65  C "Long" real for temporary coordinate calculation  C     "Long" real for temporary coordinate calculation
66  C  NOTICE the extended range of indices!!  C      NOTICE the extended range of indices!!
67        _RL xGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)        _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)        _RL yGloc(1-Olx:sNx+Olx+1,1-Oly:sNy+Oly+1)
69    
70  C These functions return the "global" index with valid values beyond  C     The functions iGl, jGl return the "global" index with valid values beyond
71  C halo regions  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        INTEGER iGl,jGl
90        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)
91        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)
92    CEOP
93    
94    
95  C     For each tile ...  C     For each tile ...
96        DO bj = myByLo(myThid), myByHi(myThid)        DO bj = myByLo(myThid), myByHi(myThid)
# Line 175  C         by formula Line 185  C         by formula
185            dlon = delX( iGl(I,bi) )            dlon = delX( iGl(I,bi) )
186            dlat = delY( jGl(J,bj) )            dlat = delY( jGl(J,bj) )
187            dXG(I,J,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad            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            dYG(I,J,bi,bj) = rSphere*dlat*deg2rad
190           ENDDO           ENDDO
191          ENDDO          ENDDO
# Line 289  C--     Calculate vertical face area (vo Line 300  C--     Calculate vertical face area (vo
300          DO J=1-Oly,sNy+Oly          DO J=1-Oly,sNy+Oly
301           DO I=1-Olx,sNx+Olx           DO I=1-Olx,sNx+Olx
302  C         by formula  C         by formula
303            lat=yC(I,J,bi,bj)            lat =0.5 _d 0*(yGloc(I,J)+yGloc(I,J+1))
304            dlon=delX( iGl(I,bi) )            dlon=0.5 _d 0*( delX( iGl(I,bi) ) + delX( iGl(I-1,bi) ) )
305            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) ) )
306            rAz(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad            rAz(I,J,bi,bj) = rSphere*rSphere*dlon*deg2rad
307       &     *abs( sin(lat*deg2rad)-sin((lat-dlat)*deg2rad) )       &     *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.            IF (abs(lat).GT.90..OR.abs(lat-dlat).GT.90.) rAz(I,J,bi,bj)=0.
309           ENDDO           ENDDO
310          ENDDO          ENDDO
311    
312  C--     Calculate trigonometric terms  C--     Calculate trigonometric terms & grid orientation:
313          DO J=1-Oly,sNy+Oly          DO J=1-Oly,sNy+Oly
314           DO I=1-Olx,sNx+Olx           DO I=1-Olx,sNx+Olx
315            lat=0.5*(yGloc(I,J)+yGloc(I,J+1))            lat=0.5*(yGloc(I,J)+yGloc(I,J+1))
316            tanPhiAtU(i,j,bi,bj)=tan(lat*deg2rad)            tanPhiAtU(I,J,bi,bj)=tan(lat*deg2rad)
317            lat=0.5*(yGloc(I,J)+yGloc(I+1,J))            lat=0.5*(yGloc(I,J)+yGloc(I+1,J))
318            tanPhiAtV(i,j,bi,bj)=tan(lat*deg2rad)            tanPhiAtV(I,J,bi,bj)=tan(lat*deg2rad)
319              angleCosC(I,J,bi,bj) = 1.
320              angleSinC(I,J,bi,bj) = 0.
321           ENDDO           ENDDO
322          ENDDO          ENDDO
323    
324    C--     Cosine(lat) scaling
325            DO J=1-OLy,sNy+OLy
326             jG = myYGlobalLo + (bj-1)*sNy + J-1
327             jG = min(max(1,jG),Ny)
328             IF (cosPower.NE.0.) THEN
329              cosFacU(J,bi,bj)=COS(yC(1,J,bi,bj)*deg2rad)
330         &                    **cosPower
331              cosFacV(J,bi,bj)=COS((yC(1,J,bi,bj)-0.5*delY(jG))*deg2rad)
332         &                    **cosPower
333              cosFacU(J,bi,bj)=ABS(cosFacU(J,bi,bj))
334              cosFacV(J,bi,bj)=ABS(cosFacV(J,bi,bj))
335              sqcosFacU(J,bi,bj)=sqrt(cosFacU(J,bi,bj))
336              sqcosFacV(J,bi,bj)=sqrt(cosFacV(J,bi,bj))
337             ELSE
338              cosFacU(J,bi,bj)=1.
339              cosFacV(J,bi,bj)=1.
340              sqcosFacU(J,bi,bj)=1.
341              sqcosFacV(J,bi,bj)=1.
342             ENDIF
343            ENDDO
344    
345         ENDDO ! bi         ENDDO ! bi
346        ENDDO ! bj        ENDDO ! bj
347    
       write(0,*) ' yC=', (yC(1,j,1,1),j=1,sNy)  
       write(0,*) 'dxF=', (dXF(1,j,1,1),j=1,sNy)  
       write(0,*) 'dyF=', (dYF(1,j,1,1),j=1,sNy)  
       write(0,*) 'dxG=', (dXG(1,j,1,1),j=1,sNy)  
       write(0,*) 'dyG=', (dYG(1,j,1,1),j=1,sNy)  
       write(0,*) 'dxC=', (dXC(1,j,1,1),j=1,sNy)  
       write(0,*) 'dyC=', (dYC(1,j,1,1),j=1,sNy)  
       write(0,*) 'dxV=', (dXV(1,j,1,1),j=1,sNy)  
       write(0,*) 'dyU=', (dYU(1,j,1,1),j=1,sNy)  
       write(0,*) ' rA=', (rA(1,j,1,1),j=1,sNy)  
       write(0,*) 'rAw=', (rAw(1,j,1,1),j=1,sNy)  
       write(0,*) 'rAs=', (rAs(1,j,1,1),j=1,sNy)  
   
348        RETURN        RETURN
349        END        END
Legend:
Removed from v.1.15  
changed lines
  Added in v.1.23
  ViewVC Help
Powered by ViewVC 1.1.22