C $Header: /u/gcmpack/MITgcm/model/src/rotate_uv2en.F,v 1.1 2010/05/07 18:35:31 gforget Exp $
C $Name:  $

#include "CPP_OPTIONS.h"

C--  File rotate_uv2en.F: Routines to handle a vector coordinate system rotation.
C--   Contents
C--   o ROTATE_UV2EN_RL
C--   o ROTATE_UV2EN_RS

      subroutine ROTATE_UV2EN_RL(
     U          uFldX, vFldY,
     U          uFldE, vFldN,
     I          xy2en, switchGrid, applyMask, kSize, mythid
     &                     )

c     ==================================================================
c     SUBROUTINE rotate_uv2en_rl
c     ==================================================================
c
c     o uFldX/vFldY are in the model grid directions. 
c     o uFldE/vFldN are eastward/northward.
c     o This routine goes from uFldX/vFldY to uFldE/vFldN (for xy2en=.TRUE.) 
c         or vice versa (for xy2en=.FALSE.).
c     o uFldX/vFldY may be at the C grid velocity points, or at the A grid
c         velocity points (i.e. the C grid cell center). uFldE/vFldN are always 
c         at the cell center. If switchGrid=.TRUE. we go from C grid uFldX/vFldY
c         to A grid uFldE/vFldN, or vice versa. If switchGrid=.FALSE. we go
c         from A grid uFldX/vFldY to A grid uFldE/vFldN, or vice versa.
c     o If applyMask=.TRUE. then masks are applied to the output.
c         If kSize=1 (resp. nr) we then use the surface (resp. 3D) masks.
c     o In any case it is assumed that exchanges are done on the outside.
c
c     ==================================================================
c     SUBROUTINE rotate_uv2en_rl
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"

c     == routine arguments ==

      integer kSize 
      logical xy2en, switchGrid, applyMask
      _RL     uFldX(1-olx:snx+olx,1-oly:sny+oly,kSize,nsx,nsy)
      _RL     vFldY(1-olx:snx+olx,1-oly:sny+oly,kSize,nsx,nsy)
      _RL     uFldE(1-olx:snx+olx,1-oly:sny+oly,kSize,nsx,nsy)
      _RL     vFldN(1-olx:snx+olx,1-oly:sny+oly,kSize,nsx,nsy)

      integer mythid

c     == local variables ==

      integer bi,bj
      integer i,j,k,kk
      _RL     tmpU(1-olx:snx+olx,1-oly:sny+oly)
      _RL     tmpV(1-olx:snx+olx,1-oly:sny+oly)
      CHARACTER*(MAX_LEN_MBUF) msgBuf

c     == end of interface ==

      if ( (kSize.NE.1).AND.(kSize.NE.nr)
     &   .AND.(applyMask) ) then
        WRITE(msgBuf,'(2A,I4,A)') ' ROTATE_UV2EN: ',
     &       'no mask has ',kSize,' levels'
        CALL PRINT_ERROR(msgBuf, myThid)      
        STOP 'ABNROMAL END: S/R ROTATE_UV2EN' 
      endif

      do bj = mybylo(mythid),mybyhi(mythid)
       do bi = mybxlo(mythid),mybxhi(mythid)
        do k = 1,kSize

        if ( (kSize.EQ.1).AND.(usingPCoords) ) then
          kk=nr
        else
          kk=k
        endif

        if ( xy2en ) then
c go from uFldX/vFldY to uFldE/vFldN
        if ( switchGrid ) then
C 1a) go from C grid velocity points to A grid velocity points
        do i = 1-olx,snx+olx
          tmpU(i,sny+Oly)=0.
          tmpV(i,sny+Oly)=0.
        enddo
        do j = 1-oly,sny+oly-1
          tmpU(snx+Olx,j)=0.
          tmpV(snx+Olx,j)=0.
        do i = 1-olx,snx+olx-1
          tmpU(i,j) = 0.5 _d 0
     &          *( uFldX(i+1,j,k,bi,bj) + uFldX(i,j,k,bi,bj) )
          tmpV(i,j) = 0.5 _d 0
     &          *( vFldY(i,j+1,k,bi,bj) + vFldY(i,j,k,bi,bj) )
          if (applyMask) then
            tmpU(i,j) = tmpU(i,j) * maskC(i,j,kk,bi,bj)
            tmpV(i,j) = tmpV(i,j) * maskC(i,j,kk,bi,bj)
          endif
          enddo
          enddo
        else
C 1b) stay at A grid velocity points (i.e. at the C grid cell center)
          do j = 1-oly,sny+oly
          do i = 1-olx,snx+olx
            tmpU(i,j) = uFldX(i,j,k,bi,bj)
            tmpV(i,j) = vFldY(i,j,k,bi,bj)
            if (applyMask) then
              tmpU(i,j) = tmpU(i,j) * maskC(i,j,kk,bi,bj)
              tmpV(i,j) = tmpV(i,j) * maskC(i,j,kk,bi,bj)
            endif
          enddo
          enddo
        endif!if ( switchGrid ) then

C 2) rotation
          do j = 1-oly,sny+oly
          do i = 1-olx,snx+olx
            uFldE(i,j,k,bi,bj) = 
     &         angleCosC(i,j,bi,bj)*tmpU(i,j)
     &        -angleSinC(i,j,bi,bj)*tmpV(i,j)
            vFldN(i,j,k,bi,bj) = 
     &         angleSinC(i,j,bi,bj)*tmpU(i,j)
     &        +angleCosC(i,j,bi,bj)*tmpV(i,j)
          enddo
          enddo

      else!if (xy2en) then
c go from uFldE/vFldN to uFldX/vFldY

C 1) rotation
          do j = 1-oly,sny+oly
          do i = 1-olx,snx+olx
            tmpU(i,j) = 
     &         angleCosC(i,j,bi,bj)*uFldE(i,j,k,bi,bj)
     &        +angleSinC(i,j,bi,bj)*vFldN(i,j,k,bi,bj)
            tmpV(i,j) =
     &        -angleSinC(i,j,bi,bj)*uFldE(i,j,k,bi,bj)
     &        +angleCosC(i,j,bi,bj)*vFldN(i,j,k,bi,bj)
          enddo
          enddo

        if ( switchGrid ) then
C 2a) go from A grid velocity points to C grid velocity points
          do i = 1-olx,snx+olx
            uFldX(i,1,k,bi,bj)=0.
            vFldY(i,1,k,bi,bj)=0.
          enddo
          do j = 1-oly+1,sny+oly
             uFldX(1,j,k,bi,bj)=0.
             vFldY(1,j,k,bi,bj)=0.
          do i = 1-olx+1,snx+olx
            uFldX(i,j,k,bi,bj) = 0.5 _d 0
     &         *( tmpU(i-1,j) + tmpU(i,j) )
            vFldY(i,j,k,bi,bj) = 0.5 _d 0
     &         *( tmpV(i,j-1) + tmpV(i,j) )
            if (applyMask) then
              uFldX(i,j,k,bi,bj)=uFldX(i,j,k,bi,bj)*maskW(i,j,kk,bi,bj)
              vFldY(i,j,k,bi,bj)=vFldY(i,j,k,bi,bj)*maskS(i,j,kk,bi,bj)
            endif
          enddo
          enddo
        else
C 2b) stay at A grid velocity points (i.e. at the C grid cell center)
          do j = 1-oly,sny+oly
          do i = 1-olx,snx+olx
            uFldX(i,j,k,bi,bj) = tmpU(i,j)
            vFldY(i,j,k,bi,bj) = tmpV(i,j)
            if (applyMask) then
              uFldX(i,j,k,bi,bj)=uFldX(i,j,k,bi,bj)*maskC(i,j,kk,bi,bj)
              vFldY(i,j,k,bi,bj)=vFldY(i,j,k,bi,bj)*maskC(i,j,kk,bi,bj)
            endif
          enddo
          enddo
        endif!if ( switchGrid ) then

        endif!if (xy2en) then

        enddo
       enddo
      enddo

      return
      end


subroutine ROTATE_UV2EN_RS( U uFldX, vFldY, U uFldE, vFldN, I xy2en, switchGrid, applyMask, kSize, mythid & ) c ================================================================== c SUBROUTINE rotate_uv2en_rs c ================================================================== c c o uFldX/vFldY are in the model grid directions. c o uFldE/vFldN are eastward/northward. c o This routine goes from uFldX/vFldY to uFldE/vFldN (for xy2en=.TRUE.) c or vice versa (for xy2en=.FALSE.). c o uFldX/vFldY may be at the C grid velocity points, or at the A grid c velocity points (i.e. the C grid cell center). uFldE/vFldN are always c at the cell center. If switchGrid=.TRUE. we go from C grid uFldX/vFldY c to A grid uFldE/vFldN, or vice versa. If switchGrid=.FALSE. we go c from A grid uFldX/vFldY to A grid uFldE/vFldN, or vice versa. c o If applyMask=.TRUE. then masks are applied to the output. c If kSize=1 (resp. nr) we then use the surface (resp. 3D) masks. c o In any case it is assumed that exchanges are done on the outside. c c ================================================================== c SUBROUTINE rotate_uv2en_rs c ================================================================== implicit none c == global variables == #include "EEPARAMS.h" #include "SIZE.h" #include "PARAMS.h" #include "GRID.h" c == routine arguments == integer kSize logical xy2en, switchGrid, applyMask _RS uFldX(1-olx:snx+olx,1-oly:sny+oly,kSize,nsx,nsy) _RS vFldY(1-olx:snx+olx,1-oly:sny+oly,kSize,nsx,nsy) _RS uFldE(1-olx:snx+olx,1-oly:sny+oly,kSize,nsx,nsy) _RS vFldN(1-olx:snx+olx,1-oly:sny+oly,kSize,nsx,nsy) integer mythid c == local variables == integer bi,bj integer i,j,k,kk _RS tmpU(1-olx:snx+olx,1-oly:sny+oly) _RS tmpV(1-olx:snx+olx,1-oly:sny+oly) CHARACTER*(MAX_LEN_MBUF) msgBuf c == end of interface == if ( (kSize.NE.1).AND.(kSize.NE.nr) & .AND.(applyMask) ) then WRITE(msgBuf,'(2A,I4,A)') ' ROTATE_UV2EN: ', & 'no mask has ',kSize,' levels' CALL PRINT_ERROR(msgBuf, myThid) STOP 'ABNROMAL END: S/R ROTATE_UV2EN' endif do bj = mybylo(mythid),mybyhi(mythid) do bi = mybxlo(mythid),mybxhi(mythid) do k = 1,kSize if ( (kSize.EQ.1).AND.(usingPCoords) ) then kk=nr else kk=k endif if ( xy2en ) then c go from uFldX/vFldY to uFldE/vFldN if ( switchGrid ) then C 1a) go from C grid velocity points to A grid velocity points do i = 1-olx,snx+olx tmpU(i,sny+Oly)=0. tmpV(i,sny+Oly)=0. enddo do j = 1-oly,sny+oly-1 tmpU(snx+Olx,j)=0. tmpV(snx+Olx,j)=0. do i = 1-olx,snx+olx-1 tmpU(i,j) = 0.5 _d 0 & *( uFldX(i+1,j,k,bi,bj) + uFldX(i,j,k,bi,bj) ) tmpV(i,j) = 0.5 _d 0 & *( vFldY(i,j+1,k,bi,bj) + vFldY(i,j,k,bi,bj) ) if (applyMask) then tmpU(i,j) = tmpU(i,j) * maskC(i,j,kk,bi,bj) tmpV(i,j) = tmpV(i,j) * maskC(i,j,kk,bi,bj) endif enddo enddo else C 1b) stay at A grid velocity points (i.e. at the C grid cell center) do j = 1-oly,sny+oly do i = 1-olx,snx+olx tmpU(i,j) = uFldX(i,j,k,bi,bj) tmpV(i,j) = vFldY(i,j,k,bi,bj) if (applyMask) then tmpU(i,j) = tmpU(i,j) * maskC(i,j,kk,bi,bj) tmpV(i,j) = tmpV(i,j) * maskC(i,j,kk,bi,bj) endif enddo enddo endif!if ( switchGrid ) then C 2) rotation do j = 1-oly,sny+oly do i = 1-olx,snx+olx uFldE(i,j,k,bi,bj) = & angleCosC(i,j,bi,bj)*tmpU(i,j) & -angleSinC(i,j,bi,bj)*tmpV(i,j) vFldN(i,j,k,bi,bj) = & angleSinC(i,j,bi,bj)*tmpU(i,j) & +angleCosC(i,j,bi,bj)*tmpV(i,j) enddo enddo else!if (xy2en) then c go from uFldE/vFldN to uFldX/vFldY C 1) rotation do j = 1-oly,sny+oly do i = 1-olx,snx+olx tmpU(i,j) = & angleCosC(i,j,bi,bj)*uFldE(i,j,k,bi,bj) & +angleSinC(i,j,bi,bj)*vFldN(i,j,k,bi,bj) tmpV(i,j) = & -angleSinC(i,j,bi,bj)*uFldE(i,j,k,bi,bj) & +angleCosC(i,j,bi,bj)*vFldN(i,j,k,bi,bj) enddo enddo if ( switchGrid ) then C 2a) go from A grid velocity points to C grid velocity points do i = 1-olx,snx+olx uFldX(i,1,k,bi,bj)=0. vFldY(i,1,k,bi,bj)=0. enddo do j = 1-oly+1,sny+oly uFldX(1,j,k,bi,bj)=0. vFldY(1,j,k,bi,bj)=0. do i = 1-olx+1,snx+olx uFldX(i,j,k,bi,bj) = 0.5 _d 0 & *( tmpU(i-1,j) + tmpU(i,j) ) vFldY(i,j,k,bi,bj) = 0.5 _d 0 & *( tmpV(i,j-1) + tmpV(i,j) ) if (applyMask) then uFldX(i,j,k,bi,bj)=uFldX(i,j,k,bi,bj)*maskW(i,j,kk,bi,bj) vFldY(i,j,k,bi,bj)=vFldY(i,j,k,bi,bj)*maskS(i,j,kk,bi,bj) endif enddo enddo else C 2b) stay at A grid velocity points (i.e. at the C grid cell center) do j = 1-oly,sny+oly do i = 1-olx,snx+olx uFldX(i,j,k,bi,bj) = tmpU(i,j) vFldY(i,j,k,bi,bj) = tmpV(i,j) if (applyMask) then uFldX(i,j,k,bi,bj)=uFldX(i,j,k,bi,bj)*maskC(i,j,kk,bi,bj) vFldY(i,j,k,bi,bj)=vFldY(i,j,k,bi,bj)*maskC(i,j,kk,bi,bj) endif enddo enddo endif!if ( switchGrid ) then endif!if (xy2en) then enddo enddo enddo return end