model/src/ini_local_grid.F

06399abf63 Jean*0001 #include "PACKAGES_CONFIG.h"
                0002 #include "CPP_OPTIONS.h"
                0003 #ifdef ALLOW_EXCH2
                0004 # include "W2_OPTIONS.h"
                0005 #endif /* ALLOW_EXCH2 */
                0006
                0007 CBOP
                0008 C     !ROUTINE: INI_LOCAL_GRID
                0009 C     !INTERFACE:
                0010       SUBROUTINE INI_LOCAL_GRID(
                0011      O                           xGloc, yGloc,
                0012      O                           delXloc, delYloc,
                0013      O                           gridNx, gridNy,
                0014      I                           bi, bj, myThid )
                0015
                0016 C     !DESCRIPTION: \bv
                0017 C     *==========================================================*
                0018 C     | SUBROUTINE INI_LOCAL_GRID
                0019 C     | o Initialise model tile-local horizontal grid
                0020 C     *==========================================================*
                0021 C     | Set local grid-point location (xGloc & yGloc) and
                0022 C     |  local grid-point spacing (delXloc,delYloc) keeping the
                0023 C     |  same units as grid-spacing input parameter (delX,delY
                0024 C     |  and xgOrigin,ygOrigin).
                0025 C     | This tile-local mesh setting will be used to build
                0026 C     |  the horizontal model grid, according to the selected
                0027 C     |  grid option (cartesian, spherical_polar or cylindrical)
                0028 C     *==========================================================*
                0029 C     \ev
                0030
                0031 C     !USES:
                0032       IMPLICIT NONE
                0033 C     === Global variables ===
                0034 #include "SIZE.h"
                0035 #include "EEPARAMS.h"
                0036 #include "PARAMS.h"
                0037 #ifdef ALLOW_EXCH2
                0038 # include "W2_EXCH2_SIZE.h"
                0039 # include "W2_EXCH2_TOPOLOGY.h"
                0040 #endif /* ALLOW_EXCH2 */
                0041 #include "SET_GRID.h"
                0042
                0043 C     !INPUT/OUTPUT PARAMETERS:
                0044 C     == Routine arguments ==
                0045 C     xGloc   :: mesh corner-point location (local "Long" real array type)
                0046 C     yGloc   :: mesh corner-point location (local "Long" real array type)
                0047 C     delXloc :: mesh spacing in X direction
                0048 C     delYloc :: mesh spacing in Y direction
                0049 C     gridNx  :: mesh total grid-point number in X direction
                0050 C     gridNy  :: mesh total grid-point number in Y direction
                0051 C     bi, bj  :: tile indices
                0052 C     myThid  :: my Thread Id Number
                0053
                0054 C NOTICE the extended range of indices!!
                0055       _RL xGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
                0056       _RL yGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
                0057 C NOTICE the extended range of indices!!
                0058       _RL delXloc(0-OLx:sNx+OLx)
                0059       _RL delYloc(0-OLy:sNy+OLy)
                0060       INTEGER gridNx, gridNy
                0061       INTEGER bi, bj
                0062       INTEGER myThid
                0063
                0064 C     !LOCAL VARIABLES:
                0065 C     == Local variables ==
                0066 C     xG0,yG0 :: coordinate of South-West tile-corner
                0067 C     iG0,jG0 :: Tile base X and Y indices within global index space
                0068 C     i, j    :: loop counters
                0069       INTEGER iG0, jG0
                0070       INTEGER i, j
                0071       _RL xG0, yG0
                0072 #ifdef ALLOW_EXCH2
                0073       INTEGER tN
                0074 #endif /* ALLOW_EXCH2 */
                0075
                0076 C     The functions iGl, jGl return the "global" index with valid values beyond
                0077 C     halo regions
                0078 C     cnh wrote:
                0079 C     >    I dont understand why we would ever have to multiply the
                0080 C     >    overlap by the total domain size e.g
                0081 C     >    OLx*Nx, OLy*Ny.
                0082 C     >    Can anybody explain? Lines are in ini_spherical_polar_grid.F.
                0083 C     >    Surprised the code works if its wrong, so I am puzzled.
                0084 C     jmc replied:
                0085 C     Yes, I can explain this since I put this modification to work
                0086 C     with small domain (where OLy > Ny, as for instance, zonal-average
                0087 C     case):
                0088 C     This has no effect on the acuracy of the evaluation of iGl(I,bi)
                0089 C     and jGl(j,bj) since we take mod(a+OLx*Nx,Nx) and mod(b+OLy*Ny,Ny).
                0090 C     But in case a or b is negative, then the FORTRAN function "mod"
                0091 C     does not return the matematical value of the "modulus" function,
                0092 C     and this is not good for your purpose.
                0093 C     This is why I add +OLx*Nx and +OLy*Ny to be sure that the 1rst
                0094 C     argument of the mod function is positive.
                0095 c     INTEGER iGl,jGl
                0096 c     iGl(i,bi) = 1+MOD(iG0+i-1+OLx*gridNx,gridNx)
                0097 c     jGl(j,bj) = 1+MOD(jG0+j-1+OLy*gridNy,gridNy)
                0098 CEOP
                0099
                0100 #ifdef ALLOW_EXCH2
                0101       gridNx = exch2_mydNx(1)
                0102       gridNy = exch2_mydNy(1)
                0103 #else /* ALLOW_EXCH2 */
                0104       gridNx = Nx
                0105       gridNy = Ny
                0106 #endif /* ALLOW_EXCH2 */
                0107
                0108 c     DO bj = myByLo(myThid), myByHi(myThid)
                0109 c      DO bi = myBxLo(myThid), myBxHi(myThid)
                0110 C     For this tile ...
                0111
                0112 C--   Set current tile base X and base Y indices within global index space
                0113 C     e.g., local indices of tile south-west corner grid point are (1,1)
                0114 C           and global indices are 1+iG0, 1+jG0
                0115 #ifdef ALLOW_EXCH2
                0116         tN = W2_myTileList(bi,bj)
                0117         iG0 = exch2_tBasex(tN)
                0118         jG0 = exch2_tBasey(tN)
                0119 #else  /* ALLOW_EXCH2 */
                0120         iG0 = myXGlobalLo - 1 + (bi-1)*sNx
                0121         jG0 = myYGlobalLo - 1 + (bj-1)*sNy
                0122 #endif /* ALLOW_EXCH2 */
                0123
                0124 C--   First find coordinate of tile corner (meaning outer corner of halo)
                0125         xG0 = xgOrigin
                0126 C       Find the X-coordinate of the outer grid-line of the "real" tile
                0127         DO i=1, iG0
                0128          xG0 = xG0 + delX(i)
                0129         ENDDO
                0130 C       Back-step to the outer grid-line of the "halo" region
                0131         DO i=1, OLx
                0132          xG0 = xG0 - delX( 1+MOD(iG0-i+OLx*gridNx,gridNx) )
                0133         ENDDO
                0134 C       Find the Y-coordinate of the outer grid-line of the "real" tile
                0135         yG0 = ygOrigin
                0136         DO j=1, jG0
                0137          yG0 = yG0 + delY(j)
                0138         ENDDO
                0139 C       Back-step to the outer grid-line of the "halo" region
                0140         DO j=1, OLy
                0141          yG0 = yG0 - delY( 1+MOD(jG0-j+OLy*gridNy,gridNy) )
                0142         ENDDO
                0143
                0144 C--     Make a local copy of current-tile grid-spacing
                0145         DO i=0-OLx,sNx+OLx
                0146           delXloc(i) = delX( 1+MOD(iG0+i-1+OLx*gridNx,gridNx) )
                0147         ENDDO
                0148         DO j=0-OLy,sNy+OLy
                0149           delYloc(j) = delY( 1+MOD(jG0+j-1+OLy*gridNy,gridNy) )
                0150         ENDDO
                0151
                0152 C--     Calculate coordinates of cell corners for N+1 grid-lines
                0153         DO j=1-OLy,sNy+OLy +1
                0154          xGloc(1-OLx,j) = xG0
                0155          DO i=1-OLx,sNx+OLx
                0156           xGloc(i+1,j) = xGloc(i,j) + delXloc(i)
                0157          ENDDO
                0158         ENDDO
                0159         DO i=1-OLx,sNx+OLx +1
                0160          yGloc(i,1-OLy) = yG0
                0161          DO j=1-OLy,sNy+OLy
                0162           yGloc(i,j+1) = yGloc(i,j) + delYloc(j)
                0163          ENDDO
                0164         ENDDO
                0165
                0166 C--   end bi,bj loops
                0167 c      ENDDO
                0168 c     ENDDO
                0169
                0170       RETURN
                0171       END