C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_set_calc.F,v 1.3 2014/11/26 20:48:32 jmc Exp $
C $Name:  $

#include "DIAG_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

CBOP 0
C     !ROUTINE: DIAGNOSTICS_SET_CALC

C     !INTERFACE:
      SUBROUTINE DIAGNOSTICS_SET_CALC( myThid )

C     !DESCRIPTION:
C     *==========================================================*
C     | S/R DIAGNOSTICS_SET_CALC
C     |  Set parameters and variables used in post-processing
C     |      diagnostics
C     *==========================================================*

C     !USES:
      IMPLICIT NONE
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DIAGNOSTICS_CALC.h"
#ifdef ALLOW_OBCS
# include "OBCS_GRID.h"
#endif /* ALLOW_OBCS */

C     !INPUT PARAMETERS:
C     myThid     ::  my thread Id number
      INTEGER      myThid
CEOP

C     !LOCAL VARIABLES:
      INTEGER bi, bj
      INTEGER i, j
      INTEGER biG, bjG
      _RL     dxLoc, dyLoc, d2Loc, d2Min
      _RL     xLine, xy0, xyLoc, xyMin
      CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef ALLOW_OBCS
      LOGICAL kPsi(1:sNx+1,1:sNy+1,nSx,nSy)
#endif /* ALLOW_OBCS */

C--   Set indices of grid-point location where Psi == 0
      IF ( xPsi0.EQ.UNSET_RS .OR. yPsi0.EQ.UNSET_RS ) THEN
C-    Set indices to (-1,0) = disabled value
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
           iPsi0(bi,bj) = -1
           jPsi0(bi,bj) =  0
         ENDDO
        ENDDO
      ELSE
#ifdef ALLOW_OBCS
C-     set flag where Psi is computed
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO j = 1,sNy+1
           DO i = 1,sNx+1
             kPsi(i,j,bi,bj) = .TRUE.
           ENDDO
          ENDDO
          IF ( useOBCS ) THEN
           DO j = 1,sNy+1
            DO i = 1,sNx+1
             kPsi(i,j,bi,bj) = OBCS_insideMask( i , j ,bi,bj).EQ.oneRS
     &                    .OR. OBCS_insideMask(i-1, j ,bi,bj).EQ.oneRS
     &                    .OR. OBCS_insideMask( i ,j-1,bi,bj).EQ.oneRS
     &                    .OR. OBCS_insideMask(i-1,j-1,bi,bj).EQ.oneRS
            ENDDO
           ENDDO
          ENDIF
         ENDDO
        ENDDO
#endif /* ALLOW_OBCS */
C-      find minimum distance:
        d2Min = -1. _d 0
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO j = 1,sNy+1
           DO i = 1,sNx+1
             dxLoc = xG(i,j,bi,bj)-xPsi0
             dyLoc = yG(i,j,bi,bj)-yPsi0
             d2Loc = dxLoc*dxLoc + dyLoc*dyLoc
#ifdef ALLOW_OBCS
             IF ((d2Loc.LT.d2Min .OR. d2Min.EQ.-1. _d 0)
     &                            .AND. kPsi(i,j,bi,bj) ) d2Min = d2Loc
#else
             IF ( d2Loc.LT.d2Min .OR. d2Min.EQ.-1. _d 0 ) d2Min = d2Loc
#endif
           ENDDO
          ENDDO
         ENDDO
        ENDDO
        d2Min = -d2Min
        _GLOBAL_MAX_RL( d2Min, myThid )
        d2Min = -d2Min
C-      find list of grid-points at minimum distance:
        xyMin = 0.
        xLine = (sNx+1)*nSx*nPx
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          iPsi0(bi,bj) = 0
          jPsi0(bi,bj) = 0
          biG = bi-1+(myXGlobalLo-1)/sNx
          bjG = bj-1+(myYGlobalLo-1)/sNy
          xy0 = biG*(sNx+1) + bjG*(sNy+1)*xLine
          DO j = 1,sNy+1
           DO i = 1,sNx+1
             dxLoc = xG(i,j,bi,bj)-xPsi0
             dyLoc = yG(i,j,bi,bj)-yPsi0
             d2Loc = dxLoc*dxLoc + dyLoc*dyLoc
#ifdef ALLOW_OBCS
             IF ( d2Loc.EQ.d2Min .AND. kPsi(i,j,bi,bj) ) THEN
#else
             IF ( d2Loc.EQ.d2Min ) THEN
#endif
               xyLoc = xy0 + i + (j-1)*xLine
               IF ( xyMin.EQ.0. _d 0 ) THEN
                 xyMin = xyLoc
               ELSE
                 xyMin = MIN( xyMin, xyLoc )
               ENDIF
               iPsi0(bi,bj) = i
               jPsi0(bi,bj) = j
             ENDIF
           ENDDO
          ENDDO
         ENDDO
        ENDDO
        xyLoc = (sNx+1)*(sNy+1)*nSx*nSy*nPx*nPy + 2.
        IF ( xyMin.EQ.0. _d 0 ) xyMin = xyLoc
        xyMin = -xyMin
        _GLOBAL_MAX_RL( xyMin, myThid )
        xyMin = -xyMin
C-      select only one (based on minimum global-map index)
        _BARRIER
        _BEGIN_MASTER( myThid )
        WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_SET_CALC: ',
     &    'setting indices iPsi0,jPsi0 where Psi == 0 :'
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
        WRITE(msgBuf,'(A,1P1E19.6,A,0PF16.3)')
     &    'DIAGNOSTICS_SET_CALC: d2Min=',d2Min, ', ijMin=',xyMin
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
        IF ( xyMin.EQ.xyLoc ) THEN
          WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_SET_CALC: ',
     &      'Fail to find the minimum distance => use Default'
          CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                        SQUEEZE_RIGHT, myThid )
          DO bj=1,nSy
           DO bi=1,nSx
             iPsi0(bi,bj) = -1
             jPsi0(bi,bj) =  0
           ENDDO
          ENDDO
        ELSE
          DO bj=1,nSy
           DO bi=1,nSx
             IF ( iPsi0(bi,bj).GT.0 ) THEN
              biG = bi-1+(myXGlobalLo-1)/sNx
              bjG = bj-1+(myYGlobalLo-1)/sNy
              xy0 = biG*(sNx+1) + bjG*(sNy+1)*xLine
              xyLoc = xy0 + iPsi0(bi,bj) + (jPsi0(bi,bj)-1)*xLine
              d2Loc = ABS( xyLoc - xyMin )
              IF ( d2Loc.GE.0.5 _d 0 ) THEN
               WRITE(msgBuf,'(2(A,2I5),A,F16.3)')
     &          ' discard: bi,bj=',bi,bj,
     &          ' ; i,j=',iPsi0(bi,bj),jPsi0(bi,bj),' ; ijLoc=',xyLoc
               CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                             SQUEEZE_RIGHT, myThid )
               iPsi0(bi,bj) = 0
               jPsi0(bi,bj) = 0
              ELSE
               WRITE(msgBuf,'(2(A,2I5),A,F16.3)')
     &          ' SELECT : bi,bj=',bi,bj,
     &          ' ; i,j=',iPsi0(bi,bj),jPsi0(bi,bj),' ; ijLoc=',xyLoc
               CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                             SQUEEZE_RIGHT, myThid )
              ENDIF
             ENDIF
c            WRITE(standardMessageUnit,'(2(A,2I5))')
c    &        ' bi,bj=',bi,bj,' i,jPsi0=', iPsi0(bi,bj),jPsi0(bi,bj)
           ENDDO
          ENDDO
        ENDIF
        WRITE(msgBuf,'(2A)')
     &   '------------------------------------------------------------'
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
        _END_MASTER( myThid )
        _BARRIER
      ENDIF

      RETURN
      END