model/src/calc_grid_angles.F

c679ff8f37 Jean*0001 #include "CPP_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP
                0005 C     !ROUTINE: CALC_GRID_ANGLES
                0006 C     !INTERFACE:
294d039451 Jean*0007       SUBROUTINE CALC_GRID_ANGLES( skipCalcAngleC, myThid )
c679ff8f37 Jean*0008
                0009 C     !DESCRIPTION: \bv
                0010 C     *===================================================================*
                0011 C     | SUBROUTINE CALC_GRID_ANGLES
                0012 C     | o calculate the angle between geographical north and model grid
                0013 C     |   north, assuming that yG holds the geographical coordinates
                0014 C     *===================================================================*
                0015 C     \ev
                0016
                0017 C     !USES:
                0018       IMPLICIT NONE
                0019 C     === Global variables ===
                0020 #include "SIZE.h"
                0021 #include "EEPARAMS.h"
                0022 #include "PARAMS.h"
                0023 #include "GRID.h"
                0024
                0025 C     !INPUT/OUTPUT PARAMETERS:
                0026 C     == Routine arguments ==
294d039451 Jean*0027 C     skipCalcAngleC :: skip setting of grid-angle at cell-center location
c679ff8f37 Jean*0028 C     myThid  :: my Thread Id Number
294d039451 Jean*0029       LOGICAL skipCalcAngleC
c679ff8f37 Jean*0030       INTEGER myThid
                0031 CEOP
                0032
                0033 C     !LOCAL VARIABLES:
                0034 C     == Local variables ==
                0035 C     bi,bj  :: Tile indices
                0036 C     i, j   :: Loop counters
                0037       INTEGER bi, bj
                0038       INTEGER  i,  j
                0039 C     pseudo velocities
                0040       _RL uPseudo(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0041       _RL vPseudo(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
294d039451 Jean*0042       _RL uC, vC, uNorm, tmpVal
c679ff8f37 Jean*0043 CEOP
                0044
294d039451 Jean*0045 C-    For each tile ...
c679ff8f37 Jean*0046       DO bj = myByLo(myThid), myByHi(myThid)
                0047        DO bi = myBxLo(myThid), myBxHi(myThid)
                0048
294d039451 Jean*0049 C-    compute pseudo velocities from stream function psi = -yG*deg2rad,
c679ff8f37 Jean*0050 C     that is, zonal flow
                0051         DO j = 1-OLy,sNy+OLy-1
                0052          DO i = 1-OLx,sNx+OLx
                0053           IF ( _dyG(i,j,bi,bj).GT.0. ) THEN
                0054             uPseudo(i,j) =
                0055      &         - ( yG(i,j,bi,bj) - yG(i,j+1,bi,bj) )*deg2rad
                0056      &         / _dyG(i,j,bi,bj)
                0057           ELSE
                0058             uPseudo(i,j) = 0.
                0059           ENDIF
294d039451 Jean*0060           u2zonDir(i,j,bi,bj) = rSphere*uPseudo(i,j)
c679ff8f37 Jean*0061          ENDDO
                0062         ENDDO
                0063         DO j = 1-OLy,sNy+OLy
                0064          DO i = 1-OLx,sNx+OLx-1
                0065           IF ( _dxG(i,j,bi,bj).GT.0. ) THEN
                0066             vPseudo(i,j) =
                0067      &         + ( yG(i,j,bi,bj) - yG(i+1,j,bi,bj) )*deg2rad
                0068      &         / _dxG(i,j,bi,bj)
                0069           ELSE
                0070             vPseudo(i,j) = 0.
                0071           ENDIF
294d039451 Jean*0072           v2zonDir(i,j,bi,bj) = rSphere*vPseudo(i,j)
c679ff8f37 Jean*0073          ENDDO
                0074         ENDDO
294d039451 Jean*0075         IF ( .NOT.skipCalcAngleC ) THEN
                0076          DO j = 1-OLy,sNy+OLy-1
                0077           DO i = 1-OLx,sNx+OLx-1
                0078            uC = 0.5*(uPseudo(i,j) + uPseudo(i+1,j))
                0079            vC = 0.5*(vPseudo(i,j) + vPseudo(i,j+1))
                0080            uNorm = SQRT(uC*uC+vC*vC)
                0081            IF ( uNorm .NE. 0. _d 0 ) uNorm = 1. _d 0/uNorm
                0082            angleCosC(i,j,bi,bj) =  uC*uNorm
                0083            angleSinC(i,j,bi,bj) = -vC*uNorm
                0084           ENDDO
                0085          ENDDO
                0086         ENDIF
                0087
                0088 C-   To check angular momentum conservation, use an alternative definition
                0089 C    of grid-angles cosine (@ U pt) & sine (@ V pt) (consistent with
                0090 C    stream-function of solid-body velocity field).
c679ff8f37 Jean*0091         DO j = 1-OLy,sNy+OLy-1
294d039451 Jean*0092          DO i = 1-OLx,sNx+OLx
                0093 C- Note: most natural way would be to divide by dyG (as below); but scaling by
                0094 C        dxC/rAw ensures that u2zonDir is exactly =1 with current Lat-Lon grid
                0095 c         tmpVal = _dyG(i,j,bi,bj) * COS( deg2rad*
                0096           tmpVal = _rAw(i,j,bi,bj) * COS( deg2rad*
                0097      &             ( yG(i,j,bi,bj) + yG(i,j+1,bi,bj) )*halfRL )
                0098           IF ( tmpVal.GT.0. ) THEN
                0099             u2zonDir(i,j,bi,bj) =  rSphere
                0100      &          *( SIN( yG(i,j+1,bi,bj)*deg2rad )
                0101      &           - SIN( yG(i, j, bi,bj)*deg2rad )
                0102      &           )* _dxC(i,j,bi,bj)/tmpVal
                0103 c    &           )/tmpVal
                0104           ELSE
                0105             u2zonDir(i,j,bi,bj) = 1.
                0106           ENDIF
                0107          ENDDO
                0108         ENDDO
                0109         DO j = 1-OLy,sNy+OLy
c679ff8f37 Jean*0110          DO i = 1-OLx,sNx+OLx-1
294d039451 Jean*0111 C- Note: most natural way would be to divide by dxG (as below); for symetry
                0112 C        reason with u2zonDir expression, we use instead dyC/rAs scaling factor
                0113 c         tmpVal = _dxG(i,j,bi,bj) * COS( deg2rad*
                0114           tmpVal = _rAs(i,j,bi,bj) * COS( deg2rad*
                0115      &             ( yG(i,j,bi,bj) + yG(i+1,j,bi,bj) )*halfRL )
                0116           IF ( tmpVal.GT.0. ) THEN
                0117             v2zonDir(i,j,bi,bj) = -rSphere
                0118      &          *( SIN( yG(i+1,j,bi,bj)*deg2rad )
                0119      &           - SIN( yG(i,j,bi,bj)*deg2rad )
                0120      &           )* _dyC(i,j,bi,bj)/tmpVal
                0121 c    &           )/tmpVal
                0122           ELSE
                0123             v2zonDir(i,j,bi,bj) = 0.
                0124           ENDIF
c679ff8f37 Jean*0125          ENDDO
                0126         ENDDO
294d039451 Jean*0127
                0128 C-    bi,bj-loops
c679ff8f37 Jean*0129        ENDDO
                0130       ENDDO
                0131
                0132       RETURN
                0133       END