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