pkg/mom_fluxform/mom_uv_boundary.F

229b6d70e7 Jean*0001 #include "MOM_FLUXFORM_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: MOM_UV_BOUNDARY
                0005
                0006 C !INTERFACE: ==========================================================
                0007       SUBROUTINE MOM_UV_BOUNDARY (
                0008      I               bi,bj,k,
                0009      I               uFld, vFld,
                0010      O               uBnd, vBnd,
                0011      I               myTime, myIter, myThid )
                0012
                0013 C !DESCRIPTION:
                0014 C Set velocity at a boundary for a momentum conserving advection
                0015 C  Note: really conserve momentum when "steps" (vertical plane)
                0016 C  or coastline (horizontal plane) are only 1 grid-point wide.
                0017
                0018 C !USES: ===============================================================
                0019 C     == Global variables ==
                0020       IMPLICIT NONE
                0021 #include "SIZE.h"
                0022 #include "EEPARAMS.h"
                0023 #include "PARAMS.h"
                0024 #include "GRID.h"
                0025
                0026 C !INPUT PARAMETERS: ===================================================
                0027 C  bi,bj          :: tile indices
                0028 C  k              :: vertical level
                0029 C  uFld           :: zonal      velocity
                0030 C  vFld           :: meridional velocity
                0031 C  myTime         :: current time
                0032 C  myIter         :: current iteration number
                0033 C  myThid         :: My Thread Id. number
                0034       INTEGER bi,bj
                0035       INTEGER k
                0036       _RL uFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0037       _RL vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
                0038       _RL     myTime
                0039       INTEGER myIter
                0040       INTEGER myThid
                0041
                0042 C !OUTPUT PARAMETERS: ==================================================
                0043 C  uBnd           :: zonal      velocity extended to boundaries
                0044 C  vBnd           :: meridional velocity extended to boundaries
                0045       _RL uBnd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0046       _RL vBnd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0047
                0048 #ifdef MOM_BOUNDARY_CONSERVE
                0049 C !LOCAL VARIABLES: ====================================================
                0050 C  i,j            :: loop indices
                0051       INTEGER i,j
                0052       INTEGER km1,kp1
                0053       _RL maskM1, maskP1
                0054       _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0055       _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0056       _RL  aBndU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0057       _RL  aBndV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0058       _RL  aBndW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0059       _RL tmpVar
                0060       LOGICAL useMomBndConserve
                0061       PARAMETER ( useMomBndConserve = .TRUE. )
                0062 CEOP
                0063
                0064 C     Initialise output array
                0065       DO j=1-OLy,sNy+OLy
                0066        DO i=1-OLx,sNx+OLx
                0067          uBnd(i,j) = uFld(i,j,k,bi,bj)
                0068          vBnd(i,j) = vFld(i,j,k,bi,bj)
                0069        ENDDO
                0070       ENDDO
                0071
                0072       IF ( useMomBndConserve ) THEN
                0073
                0074 C-    Initialise intermediate arrays:
                0075         km1 = MAX( k-1, 1 )
                0076         kp1 = MIN( k+1, Nr )
                0077         maskM1 = 1.
                0078         maskP1 = 1.
                0079         IF ( k.EQ.1  ) maskM1 = 0.
                0080         IF ( k.EQ.Nr ) maskP1 = 0.
                0081         DO j=1-OLy,sNy+OLy
                0082          DO i=1-OLx,sNx+OLx
                0083           aBndU(i,j) = 0.
                0084           aBndV(i,j) = 0.
                0085           aBndW(i,j) = 0.
                0086          ENDDO
                0087         ENDDO
                0088
                0089 C-      Calculate Divergence in 3 directions:
                0090         DO j=1-OLy,sNy+OLy
                0091          DO i=1-OLx,sNx+OLx
                0092           uTrans(i,j) = uFld(i,j,k,bi,bj)
                0093      &                * dyG(i,j,bi,bj)*deepFacC(k)
                0094      &                * drF(k)*hFacW(i,j,k,bi,bj)*rhoFacC(k)
                0095           vTrans(i,j) = vFld(i,j,k,bi,bj)
                0096      &                * dxG(i,j,bi,bj)*deepFacC(k)
                0097      &                * drF(k)*hFacS(i,j,k,bi,bj)*rhoFacC(k)
                0098          ENDDO
                0099         ENDDO
                0100         DO j=1-OLy,sNy+OLy-1
                0101          DO i=1-OLx,sNx+OLx-1
                0102           aBndU(i,j) = uTrans(i+1,j)-uTrans(i,j)
                0103           aBndV(i,j) = vTrans(i,j+1)-vTrans(i,j)
                0104           aBndW(i,j) = ABS(aBndU(i,j)+aBndV(i,j))
                0105           aBndU(i,j) = ABS(aBndU(i,j))
                0106           aBndV(i,j) = ABS(aBndV(i,j))
                0107          ENDDO
                0108         ENDDO
                0109 C-      Normalise by the sum:
                0110         DO j=1-OLy,sNy+OLy-1
                0111          DO i=1-OLx,sNx+OLx-1
                0112           tmpVar = aBndU(i,j)+aBndV(i,j)+aBndW(i,j)
                0113           IF ( tmpVar.GT.0. ) THEN
                0114             tmpVar = 1. _d 0 / tmpVar
                0115             aBndU(i,j) = aBndU(i,j)*tmpVar
                0116             aBndV(i,j) = aBndV(i,j)*tmpVar
                0117             aBndW(i,j) = aBndW(i,j)*tmpVar
                0118           ENDIF
                0119          ENDDO
                0120         ENDDO
                0121
                0122 C-      At a boundary, replace uFld,vFld by a weighted average
                0123 C       Note: multiply by 2 to cancel the 1/2 factor in advections S/R
                0124         DO j=1-OLy+1,sNy+OLy-1
                0125          DO i=1-OLx+1,sNx+OLx-1
                0126           IF (maskW(i,j,k,bi,bj).EQ.0.) THEN
                0127 C       Note: only 1 set of aBnd_U,V,W is non-zero (either i-1 or i)
                0128 C        and  only 1 uFld is non-zero (either i-1 or i+1)
                0129 C        and  only 1 uFld is non-zero (either k-1 or k+1)
                0130             uBnd(i,j) = (
                0131      &        (aBndU(i-1,j)+aBndU(i,j))
                0132      &                     *(uFld(i-1,j,k,bi,bj)+uFld(i+1,j,k,bi,bj))
                0133      &       +(aBndV(i-1,j)+aBndV(i,j))
                0134      &                     *(uFld(i,j-1,k,bi,bj)+uFld(i,j+1,k,bi,bj))
                0135      &       +(aBndW(i-1,j)+aBndW(i,j))
                0136      &                     *(uFld(i,j,km1,bi,bj)*maskM1
                0137      &                      +uFld(i,j,kp1,bi,bj)*maskP1)
                0138      &                  )*2. _d 0
                0139           ENDIF
                0140           IF (maskS(i,j,k,bi,bj).EQ.0.) THEN
                0141 C       Note: only 1 set of aBnd_U,V,W is non-zero (either j-1 or j)
                0142 C        and  only 1 vFld is non-zero (either j-1 or j+1)
                0143 C        and  only 1 vFld is non-zero (either k-1 or k+1)
                0144             vBnd(i,j) = (
                0145      &        (aBndU(i,j-1)+aBndU(i,j))
                0146      &                     *(vFld(i-1,j,k,bi,bj)+vFld(i+1,j,k,bi,bj))
                0147      &       +(aBndV(i,j-1)+aBndV(i,j))
                0148      &                     *(vFld(i,j-1,k,bi,bj)+vFld(i,j+1,k,bi,bj))
                0149      &       +(aBndW(i,j-1)+aBndW(i,j))
                0150      &                     *(vFld(i,j,km1,bi,bj)*maskM1
                0151      &                      +vFld(i,j,kp1,bi,bj)*maskP1)
                0152      &                  )*2. _d 0
                0153           ENDIF
                0154          ENDDO
                0155         ENDDO
                0156       ENDIF
                0157 #endif /* MOM_BOUNDARY_CONSERVE */
                0158
                0159       RETURN
                0160       END