Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:38:54 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
bc85e7ab92 Jean*0001 #include "DIAG_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C     !ROUTINE: DIAG_CALC_PSIVEL
                0005 C     !INTERFACE:
                0006       SUBROUTINE DIAG_CALC_PSIVEL(
                0007      I                k, iPsi0, jPsi0, uTrans, vTrans,
                0008      O                psiVel, psiLoc,
fb7847904e Jean*0009      I                prtMsg, myTime, myIter, myThid )
bc85e7ab92 Jean*0010 C     !DESCRIPTION: \bv
                0011 C     *==========================================================*
                0012 C     | SUBROUTINE DIAG_CALC_PSIVEL
                0013 C     | o Calculate horizontal transport stream-function
                0014 C     |   from non-divergent horizontal transport.
                0015 C     *==========================================================*
                0016 C     *==========================================================*
                0017 C     \ev
                0018 
                0019 C     !USES:
                0020       IMPLICIT NONE
                0021 C     === Global data ===
                0022 #include "SIZE.h"
                0023 #include "EEPARAMS.h"
                0024 #include "PARAMS.h"
fb7847904e Jean*0025 #ifdef ALLOW_OBCS
                0026 # include "GRID.h"
                0027 # include "OBCS_GRID.h"
                0028 #endif /* ALLOW_OBCS */
bc85e7ab92 Jean*0029 
                0030 C     !INPUT/OUTPUT PARAMETERS:
                0031 C     === Routine arguments ===
                0032 C     k       :: current level
                0033 C     i,jPsi0 :: indices of grid-point location where Psi == 0
                0034 C     uTrans  :: horizontal transport, u-component
                0035 C     vTrans  :: horizontal transport, u-component
                0036 C     psiVel  :: horizontal stream-function
                0037 C     psiLoc  :: horizontal stream-function at special location
fb7847904e Jean*0038 C     prtMsg  :: do print message to standard-output
bc85e7ab92 Jean*0039 C     myTime  :: current time of simulation (s)
                0040 C     myIter  :: current iteration number
                0041 C     myThid  :: my Thread Id number
                0042       INTEGER k
                0043       INTEGER iPsi0(nSx,nSy)
                0044       INTEGER jPsi0(nSx,nSy)
                0045       _RL  uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0046       _RL  vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0047       _RL  psiVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0048       _RL  psiLoc(2)
fb7847904e Jean*0049       LOGICAL prtMsg
bc85e7ab92 Jean*0050       _RL  myTime
                0051       INTEGER myIter
                0052       INTEGER myThid
                0053 
                0054 C     !LOCAL VARIABLES:
                0055 C     === Local variables ====
                0056 C     bi, bj  :: tile indices
                0057 C     i, j    :: loop indices
                0058 C     dPsiX   :: tile stream-function increment along X-dir
                0059 C     dPsiY   :: tile stream-function increment along Y-dir
                0060 C     psiOri  :: stream-function value at tile origin
                0061       INTEGER bi, bj
                0062       INTEGER i, j
                0063       _RL    dPsiX (nSx,nSy)
                0064       _RL    dPsiY (nSx,nSy)
                0065       _RL    psiOri(nSx,nSy)
                0066       _RL    offSet
                0067       LOGICAL zeroPsi
                0068 c     CHARACTER*(MAX_LEN_MBUF) msgBuf
fb7847904e Jean*0069 #ifdef ALLOW_OBCS
                0070       INTEGER is, js, ix, jx, iy, jy, ijCnt
                0071       INTEGER npass, nPts, prev_nPts
                0072       LOGICAL kPsi(1:sNx+1,1:sNy+1)
                0073 #endif /* ALLOW_OBCS */
bc85e7ab92 Jean*0074 CEOP
                0075 
5ce5806717 Jean*0076 #ifdef ALLOW_DEBUG
                0077       IF (debugMode) CALL DEBUG_ENTER('DIAG_CALC_PSIVEL',myThid)
                0078 #endif
                0079 
bc85e7ab92 Jean*0080 C--   Initialise
                0081       zeroPsi = iPsi0(1,1).GE.0
                0082       psiLoc(1) = 0.
                0083       psiLoc(2) = 0.
                0084 
                0085 C--   step.1 : compute Psi over each tile separately
                0086       DO bj=myByLo(myThid),myByHi(myThid)
                0087        DO bi=myBxLo(myThid),myBxHi(myThid)
fb7847904e Jean*0088         dPsiX (bi,bj) = 0.
                0089         dPsiY (bi,bj) = 0.
                0090         psiOri(bi,bj) = 0.
                0091         IF ( useOBCS ) THEN
                0092 C-    Case with OBC:
                0093 #ifdef ALLOW_OBCS
                0094 C- Note: OB may introduce discontinuity in domain & tile stream-function map;
                0095 C     within a tile: define a local "is-set" mask (=kPsi) and
                0096 C                    propagate stream-function value without assumption.
                0097 C     between tiles: present code is not "general", likely to work with
                0098 C                    simple OB setting and/or simple tile connection (no exch2).
                0099 C     A truly general algorithm requires to change CUMULSUM_Z_TILE (adding 1
                0100 C     more input dPsi/tile) and to account for disabled tile-connection due
                0101 C     to OB when setting cumsum tile-mapping (W2_SET_MAP_CUMSUM).
                0102          DO j=1,sNy+1
                0103           DO i=1,sNx+1
                0104             kPsi(i,j) = .FALSE.
                0105             psiVel(i,j,bi,bj) = 0.
                0106           ENDDO
                0107          ENDDO
                0108 C-    select starting point
                0109          ijCnt = sNx+sNy+1
                0110          is = 0
                0111          js = 0
                0112          DO j=1,sNy
                0113           DO i=1,sNx
                0114            IF ( OBCS_insideMask(i,j,bi,bj).EQ.1.
                0115      &                      .AND. (i+j).LE.ijCnt ) THEN
                0116              is = i
                0117              js = j
                0118              ijCnt = i+j
                0119            ENDIF
                0120           ENDDO
                0121          ENDDO
                0122          IF ( is.EQ.0 ) THEN
                0123            nPts = 0
                0124          ELSE
                0125            kPsi(is,js) = .TRUE.
                0126            nPts = 1
                0127          ENDIF
                0128          npass = 0
                0129          prev_nPts = 0
                0130          DO WHILE ( nPts.GT.prev_nPts )
                0131            prev_nPts = nPts
                0132            npass = npass + 1
                0133            DO j=1,sNy+1
                0134             DO i=1,sNx
                0135              IF ( OBCS_insideMask(i,j-1,bi,bj).EQ.1. .OR.
                0136      &            OBCS_insideMask(i, j ,bi,bj).EQ.1. ) THEN
                0137                IF ( kPsi(i,j) .AND. .NOT.kPsi(i+1,j) ) THEN
                0138                  nPts = nPts + 1
                0139                  kPsi(i+1,j) = .TRUE.
                0140                  psiVel(i+1,j,bi,bj) = psiVel(i,j,bi,bj)
                0141      &                               + vTrans(i,j,bi,bj)
                0142                ENDIF
                0143                IF ( .NOT.kPsi(i,j) .AND. kPsi(i+1,j) ) THEN
                0144                  nPts = nPts + 1
                0145                  kPsi(i,j) = .TRUE.
                0146                  psiVel(i,j,bi,bj) = psiVel(i+1,j,bi,bj)
                0147      &                               - vTrans(i,j,bi,bj)
                0148                ENDIF
                0149              ENDIF
                0150             ENDDO
                0151            ENDDO
                0152            DO j=1,sNy
                0153             DO i=1,sNx+1
                0154              IF ( OBCS_insideMask(i-1,j,bi,bj).EQ.1. .OR.
                0155      &            OBCS_insideMask( i ,j,bi,bj).EQ.1. ) THEN
                0156                IF ( kPsi(i,j) .AND. .NOT.kPsi(i,j+1) ) THEN
                0157                  nPts = nPts + 1
                0158                  kPsi(i,j+1) = .TRUE.
                0159                  psiVel(i,j+1,bi,bj) = psiVel(i,j,bi,bj)
                0160      &                               - uTrans(i,j,bi,bj)
                0161                ENDIF
                0162                IF ( .NOT.kPsi(i,j) .AND. kPsi(i,j+1) ) THEN
                0163                  nPts = nPts + 1
                0164                  kPsi(i,j) = .TRUE.
                0165                  psiVel(i,j,bi,bj) = psiVel(i,j+1,bi,bj)
                0166      &                               + uTrans(i,j,bi,bj)
                0167                ENDIF
                0168              ENDIF
                0169             ENDDO
                0170            ENDDO
                0171            IF ( prtMsg .AND. nPts.GT.prev_nPts ) THEN
                0172              _BEGIN_MASTER( myThid )
                0173              WRITE(standardMessageUnit,'(A,2I4,A,I6,I8)')
                0174      &          ' diag_calc_psivel: bi,bj=', bi, bj,
                0175      &          ' : npass,nPts=', npass, nPts
                0176              _END_MASTER( myThid )
                0177            ENDIF
                0178 C-    end do while (npass count)
                0179          ENDDO
                0180 C-    set tile increments
                0181          ix = 0
                0182          jx = 0
                0183          DO i=sNx+1,1,-1
                0184           DO j=1,sNy
                0185            IF ( kPsi(i,j) .AND. jx.EQ.0 ) THEN
                0186              ix = i
                0187              jx = j
                0188            ENDIF
                0189           ENDDO
                0190          ENDDO
                0191          IF ( jx.NE.0 ) dPsiX (bi,bj) = psiVel(ix,jx,bi,bj)
                0192          iy = 0
                0193          jy = 0
                0194          DO j=sNy+1,1,-1
                0195           DO i=1,sNx
                0196            IF ( kPsi(i,j) .AND. iy.EQ.0 ) THEN
                0197              iy = i
                0198              jy = j
                0199            ENDIF
                0200           ENDDO
                0201          ENDDO
                0202          IF ( iy.NE.0 ) dPsiY (bi,bj) = psiVel(iy,jy,bi,bj)
                0203          IF ( prtMsg ) THEN
                0204            _BEGIN_MASTER( myThid )
                0205            WRITE(standardMessageUnit,'(3(A,2I4))')
                0206      &          ' diag_calc_psivel:            is,js=', is,js,
                0207      &                 ' ; ix,jx=', ix,jx, ' ; iy,jy=', iy,jy
7c5612a70e Jean*0208 c          IF ( iPsi0(bi,bj)*jPsi0(bi,bj).GT.0 )
                0209 c    &     WRITE(standardMessageUnit,'(A,L5))')
                0210 c    &          ' diag_calc_psivel:          kPsi @ i,jPsi0 =',
                0211 c    &                          kPsi(iPsi0(bi,bj),jPsi0(bi,bj))
fb7847904e Jean*0212            _END_MASTER( myThid )
                0213          ENDIF
                0214 #endif /* ALLOW_OBCS */
                0215         ELSE
                0216 C-    Case without OBC:
bc85e7ab92 Jean*0217          psiVel(1,1,bi,bj) = psiOri(bi,bj)
                0218          j = 1
                0219          DO i=1,sNx
                0220            psiVel(i+1,j,bi,bj) = psiVel(i,j,bi,bj)
                0221      &                         + vTrans(i,j,bi,bj)
                0222          ENDDO
                0223 C-     note: can vectorise inner loop
                0224          DO j=1,sNy
                0225           DO i=1,sNx+1
                0226            psiVel(i,j+1,bi,bj) = psiVel(i,j,bi,bj)
                0227      &                         - uTrans(i,j,bi,bj)
                0228           ENDDO
                0229          ENDDO
                0230          dPsiX (bi,bj) = psiVel(1+sNx,1,bi,bj)
                0231          dPsiY (bi,bj) = psiVel(1,1+sNy,bi,bj)
fb7847904e Jean*0232 C-    end with/without OBC cases
                0233         ENDIF
bc85e7ab92 Jean*0234        ENDDO
                0235       ENDDO
                0236 
                0237       CALL CUMULSUM_Z_TILE_RL(
                0238      O                  psiOri, psiLoc,
                0239      I                  dPsiX, dPsiY, myThid )
                0240 
                0241 C--   step.2 : account for Psi @ tile origin
                0242       offSet = 0.
                0243       DO bj=myByLo(myThid),myByHi(myThid)
                0244        DO bi=myBxLo(myThid),myBxHi(myThid)
                0245         DO j=1,sNy+1
                0246          DO i=1,sNx+1
                0247            psiVel(i,j,bi,bj) = psiVel(i,j,bi,bj) + psiOri(bi,bj)
                0248          ENDDO
                0249         ENDDO
                0250         IF ( iPsi0(bi,bj)*jPsi0(bi,bj).GT.0 )
                0251      &   offSet = -psiVel(iPsi0(bi,bj),jPsi0(bi,bj),bi,bj)
                0252        ENDDO
                0253       ENDDO
                0254 
                0255       IF ( zeroPsi ) THEN
                0256 C--   step.3 : shift stream-function to satisfy Psi == 0 @ a particular location
                0257         _GLOBAL_SUM_RL( offSet, myThid )
236dd4d19a Jean*0258         psiLoc(1) = psiLoc(1) + offSet
                0259         psiLoc(2) = psiLoc(2) + offSet
bc85e7ab92 Jean*0260         DO bj=myByLo(myThid),myByHi(myThid)
                0261          DO bi=myBxLo(myThid),myBxHi(myThid)
                0262           DO j=1,sNy+1
                0263            DO i=1,sNx+1
                0264              psiVel(i,j,bi,bj) = psiVel(i,j,bi,bj) + offSet
                0265            ENDDO
                0266           ENDDO
                0267          ENDDO
                0268         ENDDO
                0269       ENDIF
                0270 
5ce5806717 Jean*0271 #ifdef ALLOW_DEBUG
                0272       IF (debugMode) CALL DEBUG_LEAVE('DIAG_CALC_PSIVEL',myThid)
                0273 #endif
                0274 
bc85e7ab92 Jean*0275       RETURN
                0276       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!