model/src/ini_spherical_polar_grid.F

1dbaea09ee Chri*0001 #include "CPP_OPTIONS.h"
924557e60a Chri*0002
fd5464a82d Jean*0003 #undef USE_BACKWARD_COMPATIBLE_GRID
fb481a83c2 Alis*0004
9366854e02 Chri*0005 CBOP
                0006 C     !ROUTINE: INI_SPHERICAL_POLAR_GRID
                0007 C     !INTERFACE:
924557e60a Chri*0008       SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )
758c85d317 Jean*0009
9366854e02 Chri*0010 C     !DESCRIPTION: \bv
758c85d317 Jean*0011 C     *==========================================================*
                0012 C     | SUBROUTINE INI_SPHERICAL_POLAR_GRID
                0013 C     | o Initialise model coordinate system arrays
                0014 C     *==========================================================*
                0015 C     | These arrays are used throughout the code in evaluating
                0016 C     | gradients, integrals and spatial avarages. This routine
                0017 C     | is called separately by each thread and initialise only
                0018 C     | the region of the domain it is "responsible" for.
                0019 C     | Under the spherical polar grid mode primitive distances
                0020 C     | in X and Y are in degrees. Distance in Z are in m or Pa
                0021 C     | depending on the vertical gridding mode.
                0022 C     *==========================================================*
9366854e02 Chri*0023 C     \ev
924557e60a Chri*0024
9366854e02 Chri*0025 C     !USES:
                0026       IMPLICIT NONE
924557e60a Chri*0027 C     === Global variables ===
                0028 #include "SIZE.h"
                0029 #include "EEPARAMS.h"
                0030 #include "PARAMS.h"
                0031 #include "GRID.h"
                0032
9366854e02 Chri*0033 C     !INPUT/OUTPUT PARAMETERS:
924557e60a Chri*0034 C     == Routine arguments ==
758c85d317 Jean*0035 C     myThid  :: my Thread Id Number
924557e60a Chri*0036       INTEGER myThid
                0037
9366854e02 Chri*0038 C     !LOCAL VARIABLES:
924557e60a Chri*0039 C     == Local variables ==
758c85d317 Jean*0040 C     bi,bj   :: tile indices
                0041 C     i, j    :: loop counters
06399abf63 Jean*0042 C     lat     :: Temporary variables used to hold latitude values.
                0043 C     dlat    :: Temporary variables used to hold latitudes increment.
                0044 C     dlon    :: Temporary variables used to hold longitude increment.
                0045 C     delXloc :: mesh spacing in X direction
                0046 C     delYloc :: mesh spacing in Y direction
                0047 C     xGloc   :: mesh corner-point location (local "Long" real array type)
                0048 C     yGloc   :: mesh corner-point location (local "Long" real array type)
1f3a025cae Jean*0049       LOGICAL skipCalcAngleC
924557e60a Chri*0050       INTEGER bi, bj
758c85d317 Jean*0051       INTEGER i,  j
06399abf63 Jean*0052       INTEGER gridNx, gridNy
                0053       _RL lat, dlat, dlon
                0054 C NOTICE the extended range of indices!!
                0055       _RL delXloc(0-OLx:sNx+OLx)
                0056       _RL delYloc(0-OLy:sNy+OLy)
                0057 C NOTICE the extended range of indices!!
                0058       _RL xGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
                0059       _RL yGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
9366854e02 Chri*0060 CEOP
fb481a83c2 Alis*0061
06399abf63 Jean*0062 C--   For each tile ...
924557e60a Chri*0063       DO bj = myByLo(myThid), myByHi(myThid)
                0064        DO bi = myBxLo(myThid), myBxHi(myThid)
fb481a83c2 Alis*0065
06399abf63 Jean*0066 C--     set tile local mesh (same units as delX,deY)
                0067 C       corresponding to coordinates of cell corners for N+1 grid-lines
                0068         CALL INI_LOCAL_GRID(
                0069      O                       xGloc, yGloc,
                0070      O                       delXloc, delYloc,
                0071      O                       gridNx, gridNy,
                0072      I                       bi, bj, myThid )
fb481a83c2 Alis*0073
                0074 C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
06399abf63 Jean*0075         DO j=1-OLy,sNy+OLy
                0076          DO i=1-OLx,sNx+OLx
758c85d317 Jean*0077           xG(i,j,bi,bj) = xGloc(i,j)
                0078           yG(i,j,bi,bj) = yGloc(i,j)
924557e60a Chri*0079          ENDDO
                0080         ENDDO
fb481a83c2 Alis*0081
                0082 C--     Calculate [xC,yC], coordinates of cell centers
06399abf63 Jean*0083         DO j=1-OLy,sNy+OLy
                0084          DO i=1-OLx,sNx+OLx
fb481a83c2 Alis*0085 C         by averaging
758c85d317 Jean*0086           xC(i,j,bi,bj) = 0.25 _d 0*(
                0087      &     xGloc(i,j)+xGloc(i+1,j)+xGloc(i,j+1)+xGloc(i+1,j+1) )
                0088           yC(i,j,bi,bj) = 0.25 _d 0*(
                0089      &     yGloc(i,j)+yGloc(i+1,j)+yGloc(i,j+1)+yGloc(i+1,j+1) )
924557e60a Chri*0090          ENDDO
                0091         ENDDO
fb481a83c2 Alis*0092
                0093 C--     Calculate [dxF,dyF], lengths between cell faces (through center)
06399abf63 Jean*0094         DO j=1-OLy,sNy+OLy
                0095          DO i=1-OLx,sNx+OLx
fb481a83c2 Alis*0096 C         by averaging
758c85d317 Jean*0097 c         dxF(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i,j+1,bi,bj))
                0098 c         dyF(i,j,bi,bj) = 0.5*(dyG(i,j,bi,bj)+dyG(i+1,j,bi,bj))
fb481a83c2 Alis*0099 C         by formula
758c85d317 Jean*0100           lat = yC(i,j,bi,bj)
06399abf63 Jean*0101           dlon = delXloc(i)
                0102           dlat = delYloc(j)
fd5464a82d Jean*0103           dxF(i,j,bi,bj) = rSphere*COS(lat*deg2rad)*dlon*deg2rad
758c85d317 Jean*0104           dyF(i,j,bi,bj) = rSphere*dlat*deg2rad
fb481a83c2 Alis*0105          ENDDO
                0106         ENDDO
                0107
                0108 C--     Calculate [dxG,dyG], lengths along cell boundaries
06399abf63 Jean*0109         DO j=1-OLy,sNy+OLy
                0110          DO i=1-OLx,sNx+OLx
fb481a83c2 Alis*0111 C         by averaging
758c85d317 Jean*0112 c         dxG(i,j,bi,bj) = 0.5*(dxF(i,j,bi,bj)+dxF(i,j-1,bi,bj))
                0113 c         dyG(i,j,bi,bj) = 0.5*(dyF(i,j,bi,bj)+dyF(i-1,j,bi,bj))
fb481a83c2 Alis*0114 C         by formula
758c85d317 Jean*0115           lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
06399abf63 Jean*0116           dlon = delXloc(i)
                0117           dlat = delYloc(j)
758c85d317 Jean*0118           dxG(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
fd5464a82d Jean*0119           IF (dxG(i,j,bi,bj).LT.1.) dxG(i,j,bi,bj)=0.
758c85d317 Jean*0120           dyG(i,j,bi,bj) = rSphere*dlat*deg2rad
fb481a83c2 Alis*0121          ENDDO
                0122         ENDDO
                0123
                0124 C--     The following arrays are not defined in some parts of the halo
                0125 C       region. We set them to zero here for safety. If they are ever
                0126 C       referred to, especially in the denominator then it is a mistake!
1f3a025cae Jean*0127 C       Note: this is now done earlier in main S/R INI_GRID
                0128 c       DO j=1-OLy,sNy+OLy
                0129 c        DO i=1-OLx,sNx+OLx
                0130 c         dxC(i,j,bi,bj) = 0.
                0131 c         dyC(i,j,bi,bj) = 0.
                0132 c         dxV(i,j,bi,bj) = 0.
                0133 c         dyU(i,j,bi,bj) = 0.
                0134 c         rAw(i,j,bi,bj) = 0.
                0135 c         rAs(i,j,bi,bj) = 0.
                0136 c        ENDDO
                0137 c       ENDDO
fb481a83c2 Alis*0138
                0139 C--     Calculate [dxC], zonal length between cell centers
06399abf63 Jean*0140         DO j=1-OLy,sNy+OLy
                0141          DO i=1-OLx+1,sNx+OLx ! NOTE range
fb481a83c2 Alis*0142 C         by averaging
758c85d317 Jean*0143           dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
fb481a83c2 Alis*0144 C         by formula
758c85d317 Jean*0145 c         lat = 0.5*(yC(i,j,bi,bj)+yC(i-1,j,bi,bj))
06399abf63 Jean*0146 c         dlon = 0.5*( delXloc(i) + delXloc(i-1) )
758c85d317 Jean*0147 c         dxC(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
fb481a83c2 Alis*0148 C         by difference
758c85d317 Jean*0149 c         lat = 0.5*(yC(i,j,bi,bj)+yC(i-1,j,bi,bj))
06399abf63 Jean*0150 c         dlon = (xC(i,j,bi,bj)-xC(i-1,j,bi,bj))
758c85d317 Jean*0151 c         dxC(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
fb481a83c2 Alis*0152          ENDDO
                0153         ENDDO
                0154
                0155 C--     Calculate [dyC], meridional length between cell centers
06399abf63 Jean*0156         DO j=1-OLy+1,sNy+OLy ! NOTE range
                0157          DO i=1-OLx,sNx+OLx
fb481a83c2 Alis*0158 C         by averaging
758c85d317 Jean*0159           dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
fb481a83c2 Alis*0160 C         by formula
06399abf63 Jean*0161 c         dlat = 0.5*( delYloc(j) + delYloc(j-1) )
758c85d317 Jean*0162 c         dyC(i,j,bi,bj) = rSphere*dlat*deg2rad
fb481a83c2 Alis*0163 C         by difference
06399abf63 Jean*0164 c         dlat = (yC(i,j,bi,bj)-yC(i,j-1,bi,bj))
758c85d317 Jean*0165 c         dyC(i,j,bi,bj) = rSphere*dlat*deg2rad
fb481a83c2 Alis*0166          ENDDO
                0167         ENDDO
                0168
                0169 C--     Calculate [dxV,dyU], length between velocity points (through corners)
06399abf63 Jean*0170         DO j=1-OLy+1,sNy+OLy ! NOTE range
                0171          DO i=1-OLx+1,sNx+OLx ! NOTE range
fb481a83c2 Alis*0172 C         by averaging (method I)
758c85d317 Jean*0173           dxV(i,j,bi,bj) = 0.5 _d 0*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
                0174           dyU(i,j,bi,bj) = 0.5 _d 0*(dyG(i,j,bi,bj)+dyG(i,j-1,bi,bj))
fb481a83c2 Alis*0175 C         by averaging (method II)
758c85d317 Jean*0176 c         dxV(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
                0177 c         dyU(i,j,bi,bj) = 0.5*(dyC(i,j,bi,bj)+dyC(i-1,j,bi,bj))
fb481a83c2 Alis*0178          ENDDO
                0179         ENDDO
                0180
                0181 C--     Calculate vertical face area (tracer cells)
06399abf63 Jean*0182         DO j=1-OLy,sNy+OLy
                0183          DO i=1-OLx,sNx+OLx
758c85d317 Jean*0184           lat=0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
06399abf63 Jean*0185           dlon = delXloc(i)
                0186           dlat = delYloc(j)
758c85d317 Jean*0187           rA(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
                0188      &        *ABS( SIN((lat+dlat)*deg2rad)-SIN(lat*deg2rad) )
fd5464a82d Jean*0189 #ifdef USE_BACKWARD_COMPATIBLE_GRID
06399abf63 Jean*0190           lat = yC(i,j,bi,bj) - delYloc(j)*0.5 _d 0
                0191           dlat= yC(i,j,bi,bj) + delYloc(j)*0.5 _d 0
758c85d317 Jean*0192           rA(i,j,bi,bj) = dyF(i,j,bi,bj)
fd5464a82d Jean*0193      &    *rSphere*(SIN(dlat*deg2rad)-SIN(lat*deg2rad))
                0194 #endif /* USE_BACKWARD_COMPATIBLE_GRID */
924557e60a Chri*0195          ENDDO
                0196         ENDDO
fb481a83c2 Alis*0197
                0198 C--     Calculate vertical face area (u cells)
06399abf63 Jean*0199         DO j=1-OLy,sNy+OLy
                0200          DO i=1-OLx+1,sNx+OLx ! NOTE range
fb481a83c2 Alis*0201 C         by averaging
758c85d317 Jean*0202           rAw(i,j,bi,bj) = 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj))
fb481a83c2 Alis*0203 C         by formula
758c85d317 Jean*0204 c         lat=yGloc(i,j)
06399abf63 Jean*0205 c         dlon = 0.5*( delXloc(i) + delXloc(i-1) )
                0206 c         dlat = delYloc(j)
758c85d317 Jean*0207 c         rAw(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
fb481a83c2 Alis*0208 c    &        *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )
                0209          ENDDO
                0210         ENDDO
                0211
                0212 C--     Calculate vertical face area (v cells)
06399abf63 Jean*0213         DO j=1-OLy,sNy+OLy
                0214          DO i=1-OLx,sNx+OLx
fb481a83c2 Alis*0215 C         by formula
758c85d317 Jean*0216           lat=yC(i,j,bi,bj)
06399abf63 Jean*0217           dlon = delXloc(i)
                0218           dlat = 0.5 _d 0*( delYloc(j) + delYloc(j-1) )
fd5464a82d Jean*0219 #ifdef USE_BACKWARD_COMPATIBLE_GRID
06399abf63 Jean*0220           dlat= delYloc(j)
fd5464a82d Jean*0221 #endif /* USE_BACKWARD_COMPATIBLE_GRID */
758c85d317 Jean*0222           rAs(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
                0223      &        *ABS( SIN(lat*deg2rad)-SIN((lat-dlat)*deg2rad) )
                0224           IF (ABS(lat).GT.90..OR.ABS(lat-dlat).GT.90.) rAs(i,j,bi,bj)=0.
fb481a83c2 Alis*0225          ENDDO
                0226         ENDDO
                0227
                0228 C--     Calculate vertical face area (vorticity points)
06399abf63 Jean*0229         DO j=1-OLy,sNy+OLy
                0230          DO i=1-OLx,sNx+OLx
fb481a83c2 Alis*0231 C         by formula
06399abf63 Jean*0232           lat  = 0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
                0233           dlon = 0.5 _d 0*( delXloc(i) + delXloc(i-1) )
                0234           dlat = 0.5 _d 0*( delYloc(j) + delYloc(j-1) )
758c85d317 Jean*0235           rAz(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
                0236      &     *ABS( SIN(lat*deg2rad)-SIN((lat-dlat)*deg2rad) )
                0237           IF (ABS(lat).GT.90..OR.ABS(lat-dlat).GT.90.) rAz(i,j,bi,bj)=0.
fb481a83c2 Alis*0238          ENDDO
                0239         ENDDO
                0240
51910c84f4 Jean*0241 C--     Calculate trigonometric terms & grid orientation:
06399abf63 Jean*0242         DO j=1-OLy,sNy+OLy
                0243          DO i=1-OLx,sNx+OLx
758c85d317 Jean*0244           lat=0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
                0245           tanPhiAtU(i,j,bi,bj)=TAN(lat*deg2rad)
                0246           lat=0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
                0247           tanPhiAtV(i,j,bi,bj)=TAN(lat*deg2rad)
1f3a025cae Jean*0248 C       Note: this is now done earlier in main S/R INI_GRID
                0249 c         angleCosC(i,j,bi,bj) = 1.
                0250 c         angleSinC(i,j,bi,bj) = 0.
c1701ff971 Alis*0251          ENDDO
                0252         ENDDO
fb481a83c2 Alis*0253
aea29c8517 Alis*0254 C--     Cosine(lat) scaling
758c85d317 Jean*0255         DO j=1-OLy,sNy+OLy
fd5464a82d Jean*0256          i = 1
aea29c8517 Alis*0257          IF (cosPower.NE.0.) THEN
fd5464a82d Jean*0258           lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
                0259           cosFacU(j,bi,bj) = ABS( COS(lat*deg2rad) )**cosPower
                0260           lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
                0261           cosFacV(j,bi,bj) = ABS( COS(lat*deg2rad) )**cosPower
                0262           sqcosFacU(j,bi,bj) = SQRT(cosFacU(j,bi,bj))
                0263           sqcosFacV(j,bi,bj) = SQRT(cosFacV(j,bi,bj))
aea29c8517 Alis*0264          ELSE
fd5464a82d Jean*0265           cosFacU(j,bi,bj) = 1.
                0266           cosFacV(j,bi,bj) = 1.
758c85d317 Jean*0267           sqcosFacU(j,bi,bj)=1.
                0268           sqcosFacV(j,bi,bj)=1.
aea29c8517 Alis*0269          ENDIF
                0270         ENDDO
                0271
fd5464a82d Jean*0272 C--   end bi,bj loops
                0273        ENDDO
                0274       ENDDO
fb481a83c2 Alis*0275
7514c1bd55 Mart*0276       IF ( rotateGrid ) THEN
                0277        CALL ROTATE_SPHERICAL_POLAR_GRID( xC, yC, myThid )
                0278        CALL ROTATE_SPHERICAL_POLAR_GRID( xG, yG, myThid )
1f3a025cae Jean*0279        skipCalcAngleC = .FALSE.
                0280        CALL CALC_GRID_ANGLES( skipCalcAngleC, myThid )
7514c1bd55 Mart*0281       ENDIF
                0282
924557e60a Chri*0283       RETURN
                0284       END