C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diag_calc_psivel.F,v 1.5 2014/11/26 18:31:34 jmc Exp $ C $Name: $ #include "DIAG_OPTIONS.h" CBOP C !ROUTINE: DIAG_CALC_PSIVEL C !INTERFACE: SUBROUTINE DIAG_CALC_PSIVEL( I k, iPsi0, jPsi0, uTrans, vTrans, O psiVel, psiLoc, I prtMsg, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE DIAG_CALC_PSIVEL C | o Calculate horizontal transport stream-function C | from non-divergent horizontal transport. C *==========================================================* C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global data === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #ifdef ALLOW_OBCS # include "GRID.h" # include "OBCS_GRID.h" #endif /* ALLOW_OBCS */ C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C k :: current level C i,jPsi0 :: indices of grid-point location where Psi == 0 C uTrans :: horizontal transport, u-component C vTrans :: horizontal transport, u-component C psiVel :: horizontal stream-function C psiLoc :: horizontal stream-function at special location C prtMsg :: do print message to standard-output C myTime :: current time of simulation (s) C myIter :: current iteration number C myThid :: my Thread Id number INTEGER k INTEGER iPsi0(nSx,nSy) INTEGER jPsi0(nSx,nSy) _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL psiVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) _RL psiLoc(2) LOGICAL prtMsg _RL myTime INTEGER myIter INTEGER myThid C !LOCAL VARIABLES: C === Local variables ==== C bi, bj :: tile indices C i, j :: loop indices C dPsiX :: tile stream-function increment along X-dir C dPsiY :: tile stream-function increment along Y-dir C psiOri :: stream-function value at tile origin INTEGER bi, bj INTEGER i, j _RL dPsiX (nSx,nSy) _RL dPsiY (nSx,nSy) _RL psiOri(nSx,nSy) _RL offSet LOGICAL zeroPsi c CHARACTER*(MAX_LEN_MBUF) msgBuf #ifdef ALLOW_OBCS INTEGER is, js, ix, jx, iy, jy, ijCnt INTEGER npass, nPts, prev_nPts LOGICAL kPsi(1:sNx+1,1:sNy+1) #endif /* ALLOW_OBCS */ CEOP #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_ENTER('DIAG_CALC_PSIVEL',myThid) #endif C-- Initialise zeroPsi = iPsi0(1,1).GE.0 psiLoc(1) = 0. psiLoc(2) = 0. C-- step.1 : compute Psi over each tile separately DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) dPsiX (bi,bj) = 0. dPsiY (bi,bj) = 0. psiOri(bi,bj) = 0. IF ( useOBCS ) THEN C- Case with OBC: #ifdef ALLOW_OBCS C- Note: OB may introduce discontinuity in domain & tile stream-function map; C within a tile: define a local "is-set" mask (=kPsi) and C propagate stream-function value without assumption. C between tiles: present code is not "general", likely to work with C simple OB setting and/or simple tile connection (no exch2). C A truly general algorithm requires to change CUMULSUM_Z_TILE (adding 1 C more input dPsi/tile) and to account for disabled tile-connection due C to OB when setting cumsum tile-mapping (W2_SET_MAP_CUMSUM). DO j=1,sNy+1 DO i=1,sNx+1 kPsi(i,j) = .FALSE. psiVel(i,j,bi,bj) = 0. ENDDO ENDDO C- select starting point ijCnt = sNx+sNy+1 is = 0 js = 0 DO j=1,sNy DO i=1,sNx IF ( OBCS_insideMask(i,j,bi,bj).EQ.1. & .AND. (i+j).LE.ijCnt ) THEN is = i js = j ijCnt = i+j ENDIF ENDDO ENDDO IF ( is.EQ.0 ) THEN nPts = 0 ELSE kPsi(is,js) = .TRUE. nPts = 1 ENDIF npass = 0 prev_nPts = 0 DO WHILE ( nPts.GT.prev_nPts ) prev_nPts = nPts npass = npass + 1 DO j=1,sNy+1 DO i=1,sNx IF ( OBCS_insideMask(i,j-1,bi,bj).EQ.1. .OR. & OBCS_insideMask(i, j ,bi,bj).EQ.1. ) THEN IF ( kPsi(i,j) .AND. .NOT.kPsi(i+1,j) ) THEN nPts = nPts + 1 kPsi(i+1,j) = .TRUE. psiVel(i+1,j,bi,bj) = psiVel(i,j,bi,bj) & + vTrans(i,j,bi,bj) ENDIF IF ( .NOT.kPsi(i,j) .AND. kPsi(i+1,j) ) THEN nPts = nPts + 1 kPsi(i,j) = .TRUE. psiVel(i,j,bi,bj) = psiVel(i+1,j,bi,bj) & - vTrans(i,j,bi,bj) ENDIF ENDIF ENDDO ENDDO DO j=1,sNy DO i=1,sNx+1 IF ( OBCS_insideMask(i-1,j,bi,bj).EQ.1. .OR. & OBCS_insideMask( i ,j,bi,bj).EQ.1. ) THEN IF ( kPsi(i,j) .AND. .NOT.kPsi(i,j+1) ) THEN nPts = nPts + 1 kPsi(i,j+1) = .TRUE. psiVel(i,j+1,bi,bj) = psiVel(i,j,bi,bj) & - uTrans(i,j,bi,bj) ENDIF IF ( .NOT.kPsi(i,j) .AND. kPsi(i,j+1) ) THEN nPts = nPts + 1 kPsi(i,j) = .TRUE. psiVel(i,j,bi,bj) = psiVel(i,j+1,bi,bj) & + uTrans(i,j,bi,bj) ENDIF ENDIF ENDDO ENDDO IF ( prtMsg .AND. nPts.GT.prev_nPts ) THEN _BEGIN_MASTER( myThid ) WRITE(standardMessageUnit,'(A,2I4,A,I6,I8)') & ' diag_calc_psivel: bi,bj=', bi, bj, & ' : npass,nPts=', npass, nPts _END_MASTER( myThid ) ENDIF C- end do while (npass count) ENDDO C- set tile increments ix = 0 jx = 0 DO i=sNx+1,1,-1 DO j=1,sNy IF ( kPsi(i,j) .AND. jx.EQ.0 ) THEN ix = i jx = j ENDIF ENDDO ENDDO IF ( jx.NE.0 ) dPsiX (bi,bj) = psiVel(ix,jx,bi,bj) iy = 0 jy = 0 DO j=sNy+1,1,-1 DO i=1,sNx IF ( kPsi(i,j) .AND. iy.EQ.0 ) THEN iy = i jy = j ENDIF ENDDO ENDDO IF ( iy.NE.0 ) dPsiY (bi,bj) = psiVel(iy,jy,bi,bj) IF ( prtMsg ) THEN _BEGIN_MASTER( myThid ) WRITE(standardMessageUnit,'(3(A,2I4))') & ' diag_calc_psivel: is,js=', is,js, & ' ; ix,jx=', ix,jx, ' ; iy,jy=', iy,jy c IF ( iPsi0(bi,bj)*jPsi0(bi,bj).GT.0 ) c & WRITE(standardMessageUnit,'(A,L5))') c & ' diag_calc_psivel: kPsi @ i,jPsi0 =', c & kPsi(iPsi0(bi,bj),jPsi0(bi,bj)) _END_MASTER( myThid ) ENDIF #endif /* ALLOW_OBCS */ ELSE C- Case without OBC: psiVel(1,1,bi,bj) = psiOri(bi,bj) j = 1 DO i=1,sNx psiVel(i+1,j,bi,bj) = psiVel(i,j,bi,bj) & + vTrans(i,j,bi,bj) ENDDO C- note: can vectorise inner loop DO j=1,sNy DO i=1,sNx+1 psiVel(i,j+1,bi,bj) = psiVel(i,j,bi,bj) & - uTrans(i,j,bi,bj) ENDDO ENDDO dPsiX (bi,bj) = psiVel(1+sNx,1,bi,bj) dPsiY (bi,bj) = psiVel(1,1+sNy,bi,bj) C- end with/without OBC cases ENDIF ENDDO ENDDO CALL CUMULSUM_Z_TILE_RL( O psiOri, psiLoc, I dPsiX, dPsiY, myThid ) C-- step.2 : account for Psi @ tile origin offSet = 0. DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy+1 DO i=1,sNx+1 psiVel(i,j,bi,bj) = psiVel(i,j,bi,bj) + psiOri(bi,bj) ENDDO ENDDO IF ( iPsi0(bi,bj)*jPsi0(bi,bj).GT.0 ) & offSet = -psiVel(iPsi0(bi,bj),jPsi0(bi,bj),bi,bj) ENDDO ENDDO IF ( zeroPsi ) THEN C-- step.3 : shift stream-function to satisfy Psi == 0 @ a particular location _GLOBAL_SUM_RL( offSet, myThid ) psiLoc(1) = psiLoc(1) + offSet psiLoc(2) = psiLoc(2) + offSet DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO j=1,sNy+1 DO i=1,sNx+1 psiVel(i,j,bi,bj) = psiVel(i,j,bi,bj) + offSet ENDDO ENDDO ENDDO ENDDO ENDIF #ifdef ALLOW_DEBUG IF (debugMode) CALL DEBUG_LEAVE('DIAG_CALC_PSIVEL',myThid) #endif RETURN END