#include "CPP_OPTIONS.h"

#undef USE_BACKWARD_COMPATIBLE_GRID

CBOP
C     !ROUTINE: INI_SPHERICAL_POLAR_GRID
C     !INTERFACE:
      SUBROUTINE INI_SPHERICAL_POLAR_GRID( myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE INI_SPHERICAL_POLAR_GRID
C     | o Initialise model coordinate system arrays
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     | 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     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid  :: my Thread Id Number
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
C     bi,bj   :: tile indices
C     i, j    :: loop counters
C     lat     :: Temporary variables used to hold latitude values.
C     dlat    :: Temporary variables used to hold latitudes increment.
C     dlon    :: Temporary variables used to hold longitude increment.
C     delXloc :: mesh spacing in X direction
C     delYloc :: mesh spacing in Y direction
C     xGloc   :: mesh corner-point location (local "Long" real array type)
C     yGloc   :: mesh corner-point location (local "Long" real array type)
      LOGICAL skipCalcAngleC
      INTEGER bi, bj
      INTEGER i,  j
      INTEGER gridNx, gridNy
      _RL lat, dlat, dlon
C NOTICE the extended range of indices!!
      _RL delXloc(0-OLx:sNx+OLx)
      _RL delYloc(0-OLy:sNy+OLy)
C NOTICE the extended range of indices!!
      _RL xGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
      _RL yGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
CEOP

C--   For each tile ...
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)

C--     set tile local mesh (same units as delX,deY)
C       corresponding to coordinates of cell corners for N+1 grid-lines
        CALL INI_LOCAL_GRID(
     O                       xGloc, yGloc,
     O                       delXloc, delYloc,
     O                       gridNx, gridNy,
     I                       bi, bj, myThid )

C--     Make a permanent copy of [xGloc,yGloc] in [xG,yG]
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          xG(i,j,bi,bj) = xGloc(i,j)
          yG(i,j,bi,bj) = yGloc(i,j)
         ENDDO
        ENDDO

C--     Calculate [xC,yC], coordinates of cell centers
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
C         by averaging
          xC(i,j,bi,bj) = 0.25 _d 0*(
     &     xGloc(i,j)+xGloc(i+1,j)+xGloc(i,j+1)+xGloc(i+1,j+1) )
          yC(i,j,bi,bj) = 0.25 _d 0*(
     &     yGloc(i,j)+yGloc(i+1,j)+yGloc(i,j+1)+yGloc(i+1,j+1) )
         ENDDO
        ENDDO

C--     Calculate [dxF,dyF], lengths between cell faces (through center)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
C         by averaging
c         dxF(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i,j+1,bi,bj))
c         dyF(i,j,bi,bj) = 0.5*(dyG(i,j,bi,bj)+dyG(i+1,j,bi,bj))
C         by formula
          lat = yC(i,j,bi,bj)
          dlon = delXloc(i)
          dlat = delYloc(j)
          dxF(i,j,bi,bj) = rSphere*COS(lat*deg2rad)*dlon*deg2rad
          dyF(i,j,bi,bj) = rSphere*dlat*deg2rad
         ENDDO
        ENDDO

C--     Calculate [dxG,dyG], lengths along cell boundaries
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
C         by averaging
c         dxG(i,j,bi,bj) = 0.5*(dxF(i,j,bi,bj)+dxF(i,j-1,bi,bj))
c         dyG(i,j,bi,bj) = 0.5*(dyF(i,j,bi,bj)+dyF(i-1,j,bi,bj))
C         by formula
          lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
          dlon = delXloc(i)
          dlat = delYloc(j)
          dxG(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
          IF (dxG(i,j,bi,bj).LT.1.) dxG(i,j,bi,bj)=0.
          dyG(i,j,bi,bj) = rSphere*dlat*deg2rad
         ENDDO
        ENDDO

C--     The following arrays are not defined in some parts of the halo
C       region. We set them to zero here for safety. If they are ever
C       referred to, especially in the denominator then it is a mistake!
C       Note: this is now done earlier in main S/R INI_GRID
c       DO j=1-OLy,sNy+OLy
c        DO i=1-OLx,sNx+OLx
c         dxC(i,j,bi,bj) = 0.
c         dyC(i,j,bi,bj) = 0.
c         dxV(i,j,bi,bj) = 0.
c         dyU(i,j,bi,bj) = 0.
c         rAw(i,j,bi,bj) = 0.
c         rAs(i,j,bi,bj) = 0.
c        ENDDO
c       ENDDO

C--     Calculate [dxC], zonal length between cell centers
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx+1,sNx+OLx ! NOTE range
C         by averaging
          dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
C         by formula
c         lat = 0.5*(yC(i,j,bi,bj)+yC(i-1,j,bi,bj))
c         dlon = 0.5*( delXloc(i) + delXloc(i-1) )
c         dxC(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
C         by difference
c         lat = 0.5*(yC(i,j,bi,bj)+yC(i-1,j,bi,bj))
c         dlon = (xC(i,j,bi,bj)-xC(i-1,j,bi,bj))
c         dxC(i,j,bi,bj) = rSphere*COS(deg2rad*lat)*dlon*deg2rad
         ENDDO
        ENDDO

C--     Calculate [dyC], meridional length between cell centers
        DO j=1-OLy+1,sNy+OLy ! NOTE range
         DO i=1-OLx,sNx+OLx
C         by averaging
          dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
C         by formula
c         dlat = 0.5*( delYloc(j) + delYloc(j-1) )
c         dyC(i,j,bi,bj) = rSphere*dlat*deg2rad
C         by difference
c         dlat = (yC(i,j,bi,bj)-yC(i,j-1,bi,bj))
c         dyC(i,j,bi,bj) = rSphere*dlat*deg2rad
         ENDDO
        ENDDO

C--     Calculate [dxV,dyU], length between velocity points (through corners)
        DO j=1-OLy+1,sNy+OLy ! NOTE range
         DO i=1-OLx+1,sNx+OLx ! NOTE range
C         by averaging (method I)
          dxV(i,j,bi,bj) = 0.5 _d 0*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
          dyU(i,j,bi,bj) = 0.5 _d 0*(dyG(i,j,bi,bj)+dyG(i,j-1,bi,bj))
C         by averaging (method II)
c         dxV(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
c         dyU(i,j,bi,bj) = 0.5*(dyC(i,j,bi,bj)+dyC(i-1,j,bi,bj))
         ENDDO
        ENDDO

C--     Calculate vertical face area (tracer cells)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          lat=0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
          dlon = delXloc(i)
          dlat = delYloc(j)
          rA(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
     &        *ABS( SIN((lat+dlat)*deg2rad)-SIN(lat*deg2rad) )
#ifdef USE_BACKWARD_COMPATIBLE_GRID
          lat = yC(i,j,bi,bj) - delYloc(j)*0.5 _d 0
          dlat= yC(i,j,bi,bj) + delYloc(j)*0.5 _d 0
          rA(i,j,bi,bj) = dyF(i,j,bi,bj)
     &    *rSphere*(SIN(dlat*deg2rad)-SIN(lat*deg2rad))
#endif /* USE_BACKWARD_COMPATIBLE_GRID */
         ENDDO
        ENDDO

C--     Calculate vertical face area (u cells)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx+1,sNx+OLx ! NOTE range
C         by averaging
          rAw(i,j,bi,bj) = 0.5 _d 0*(rA(i,j,bi,bj)+rA(i-1,j,bi,bj))
C         by formula
c         lat=yGloc(i,j)
c         dlon = 0.5*( delXloc(i) + delXloc(i-1) )
c         dlat = delYloc(j)
c         rAw(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
c    &        *abs( sin((lat+dlat)*deg2rad)-sin(lat*deg2rad) )
         ENDDO
        ENDDO

C--     Calculate vertical face area (v cells)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
C         by formula
          lat=yC(i,j,bi,bj)
          dlon = delXloc(i)
          dlat = 0.5 _d 0*( delYloc(j) + delYloc(j-1) )
#ifdef USE_BACKWARD_COMPATIBLE_GRID
          dlat= delYloc(j)
#endif /* USE_BACKWARD_COMPATIBLE_GRID */
          rAs(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
     &        *ABS( SIN(lat*deg2rad)-SIN((lat-dlat)*deg2rad) )
          IF (ABS(lat).GT.90..OR.ABS(lat-dlat).GT.90.) rAs(i,j,bi,bj)=0.
         ENDDO
        ENDDO

C--     Calculate vertical face area (vorticity points)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
C         by formula
          lat  = 0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
          dlon = 0.5 _d 0*( delXloc(i) + delXloc(i-1) )
          dlat = 0.5 _d 0*( delYloc(j) + delYloc(j-1) )
          rAz(i,j,bi,bj) = rSphere*rSphere*dlon*deg2rad
     &     *ABS( SIN(lat*deg2rad)-SIN((lat-dlat)*deg2rad) )
          IF (ABS(lat).GT.90..OR.ABS(lat-dlat).GT.90.) rAz(i,j,bi,bj)=0.
         ENDDO
        ENDDO

C--     Calculate trigonometric terms & grid orientation:
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          lat=0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
          tanPhiAtU(i,j,bi,bj)=TAN(lat*deg2rad)
          lat=0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
          tanPhiAtV(i,j,bi,bj)=TAN(lat*deg2rad)
C       Note: this is now done earlier in main S/R INI_GRID
c         angleCosC(i,j,bi,bj) = 1.
c         angleSinC(i,j,bi,bj) = 0.
         ENDDO
        ENDDO

C--     Cosine(lat) scaling
        DO j=1-OLy,sNy+OLy
         i = 1
         IF (cosPower.NE.0.) THEN
          lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i,j+1))
          cosFacU(j,bi,bj) = ABS( COS(lat*deg2rad) )**cosPower
          lat = 0.5 _d 0*(yGloc(i,j)+yGloc(i+1,j))
          cosFacV(j,bi,bj) = ABS( COS(lat*deg2rad) )**cosPower
          sqcosFacU(j,bi,bj) = SQRT(cosFacU(j,bi,bj))
          sqcosFacV(j,bi,bj) = SQRT(cosFacV(j,bi,bj))
         ELSE
          cosFacU(j,bi,bj) = 1.
          cosFacV(j,bi,bj) = 1.
          sqcosFacU(j,bi,bj)=1.
          sqcosFacV(j,bi,bj)=1.
         ENDIF
        ENDDO

C--   end bi,bj loops
       ENDDO
      ENDDO

      IF ( rotateGrid ) THEN
       CALL ROTATE_SPHERICAL_POLAR_GRID( xC, yC, myThid )
       CALL ROTATE_SPHERICAL_POLAR_GRID( xG, yG, myThid )
       skipCalcAngleC = .FALSE.
       CALL CALC_GRID_ANGLES( skipCalcAngleC, myThid )
      ENDIF

      RETURN
      END
