Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:39:02 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
c19aee9e8e Jean*0001 #include "DIAG_OPTIONS.h"
                0002 
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 
                0005 CBOP 0
                0006 C     !ROUTINE: DIAGNOSTICS_SET_CALC
                0007 
                0008 C     !INTERFACE:
                0009       SUBROUTINE DIAGNOSTICS_SET_CALC( myThid )
                0010 
                0011 C     !DESCRIPTION:
                0012 C     *==========================================================*
                0013 C     | S/R DIAGNOSTICS_SET_CALC
                0014 C     |  Set parameters and variables used in post-processing
                0015 C     |      diagnostics
                0016 C     *==========================================================*
                0017 
                0018 C     !USES:
                0019       IMPLICIT NONE
                0020 #include "EEPARAMS.h"
                0021 #include "SIZE.h"
20f8314a73 Jean*0022 #include "PARAMS.h"
c19aee9e8e Jean*0023 #include "GRID.h"
                0024 #include "DIAGNOSTICS_CALC.h"
20f8314a73 Jean*0025 #ifdef ALLOW_OBCS
                0026 # include "OBCS_GRID.h"
                0027 #endif /* ALLOW_OBCS */
c19aee9e8e Jean*0028 
                0029 C     !INPUT PARAMETERS:
                0030 C     myThid     ::  my thread Id number
                0031       INTEGER      myThid
                0032 CEOP
                0033 
                0034 C     !LOCAL VARIABLES:
                0035       INTEGER bi, bj
                0036       INTEGER i, j
                0037       INTEGER biG, bjG
                0038       _RL     dxLoc, dyLoc, d2Loc, d2Min
                0039       _RL     xLine, xy0, xyLoc, xyMin
                0040       CHARACTER*(MAX_LEN_MBUF) msgBuf
20f8314a73 Jean*0041 #ifdef ALLOW_OBCS
                0042       LOGICAL kPsi(1:sNx+1,1:sNy+1,nSx,nSy)
                0043 #endif /* ALLOW_OBCS */
c19aee9e8e Jean*0044 
                0045 C--   Set indices of grid-point location where Psi == 0
                0046       IF ( xPsi0.EQ.UNSET_RS .OR. yPsi0.EQ.UNSET_RS ) THEN
                0047 C-    Set indices to (-1,0) = disabled value
                0048         DO bj=myByLo(myThid),myByHi(myThid)
                0049          DO bi=myBxLo(myThid),myBxHi(myThid)
                0050            iPsi0(bi,bj) = -1
                0051            jPsi0(bi,bj) =  0
                0052          ENDDO
                0053         ENDDO
                0054       ELSE
20f8314a73 Jean*0055 #ifdef ALLOW_OBCS
                0056 C-     set flag where Psi is computed
                0057         DO bj=myByLo(myThid),myByHi(myThid)
                0058          DO bi=myBxLo(myThid),myBxHi(myThid)
                0059           DO j = 1,sNy+1
                0060            DO i = 1,sNx+1
                0061              kPsi(i,j,bi,bj) = .TRUE.
                0062            ENDDO
                0063           ENDDO
                0064           IF ( useOBCS ) THEN
                0065            DO j = 1,sNy+1
                0066             DO i = 1,sNx+1
                0067              kPsi(i,j,bi,bj) = OBCS_insideMask( i , j ,bi,bj).EQ.oneRS
                0068      &                    .OR. OBCS_insideMask(i-1, j ,bi,bj).EQ.oneRS
                0069      &                    .OR. OBCS_insideMask( i ,j-1,bi,bj).EQ.oneRS
                0070      &                    .OR. OBCS_insideMask(i-1,j-1,bi,bj).EQ.oneRS
                0071             ENDDO
                0072            ENDDO
                0073           ENDIF
                0074          ENDDO
                0075         ENDDO
                0076 #endif /* ALLOW_OBCS */
c19aee9e8e Jean*0077 C-      find minimum distance:
                0078         d2Min = -1. _d 0
                0079         DO bj=myByLo(myThid),myByHi(myThid)
                0080          DO bi=myBxLo(myThid),myBxHi(myThid)
                0081           DO j = 1,sNy+1
                0082            DO i = 1,sNx+1
                0083              dxLoc = xG(i,j,bi,bj)-xPsi0
                0084              dyLoc = yG(i,j,bi,bj)-yPsi0
                0085              d2Loc = dxLoc*dxLoc + dyLoc*dyLoc
20f8314a73 Jean*0086 #ifdef ALLOW_OBCS
                0087              IF ((d2Loc.LT.d2Min .OR. d2Min.EQ.-1. _d 0)
                0088      &                            .AND. kPsi(i,j,bi,bj) ) d2Min = d2Loc
                0089 #else
c19aee9e8e Jean*0090              IF ( d2Loc.LT.d2Min .OR. d2Min.EQ.-1. _d 0 ) d2Min = d2Loc
20f8314a73 Jean*0091 #endif
c19aee9e8e Jean*0092            ENDDO
                0093           ENDDO
                0094          ENDDO
                0095         ENDDO
                0096         d2Min = -d2Min
                0097         _GLOBAL_MAX_RL( d2Min, myThid )
                0098         d2Min = -d2Min
                0099 C-      find list of grid-points at minimum distance:
                0100         xyMin = 0.
                0101         xLine = (sNx+1)*nSx*nPx
                0102         DO bj=myByLo(myThid),myByHi(myThid)
                0103          DO bi=myBxLo(myThid),myBxHi(myThid)
                0104           iPsi0(bi,bj) = 0
                0105           jPsi0(bi,bj) = 0
                0106           biG = bi-1+(myXGlobalLo-1)/sNx
                0107           bjG = bj-1+(myYGlobalLo-1)/sNy
                0108           xy0 = biG*(sNx+1) + bjG*(sNy+1)*xLine
                0109           DO j = 1,sNy+1
                0110            DO i = 1,sNx+1
                0111              dxLoc = xG(i,j,bi,bj)-xPsi0
                0112              dyLoc = yG(i,j,bi,bj)-yPsi0
                0113              d2Loc = dxLoc*dxLoc + dyLoc*dyLoc
20f8314a73 Jean*0114 #ifdef ALLOW_OBCS
                0115              IF ( d2Loc.EQ.d2Min .AND. kPsi(i,j,bi,bj) ) THEN
                0116 #else
c19aee9e8e Jean*0117              IF ( d2Loc.EQ.d2Min ) THEN
20f8314a73 Jean*0118 #endif
c19aee9e8e Jean*0119                xyLoc = xy0 + i + (j-1)*xLine
                0120                IF ( xyMin.EQ.0. _d 0 ) THEN
                0121                  xyMin = xyLoc
                0122                ELSE
                0123                  xyMin = MIN( xyMin, xyLoc )
                0124                ENDIF
                0125                iPsi0(bi,bj) = i
                0126                jPsi0(bi,bj) = j
                0127              ENDIF
                0128            ENDDO
                0129           ENDDO
                0130          ENDDO
                0131         ENDDO
d7550daaf0 Jean*0132         xyLoc = (sNx+1)*(sNy+1)*nSx*nSy*nPx*nPy + 2.
                0133         IF ( xyMin.EQ.0. _d 0 ) xyMin = xyLoc
c19aee9e8e Jean*0134         xyMin = -xyMin
                0135         _GLOBAL_MAX_RL( xyMin, myThid )
                0136         xyMin = -xyMin
                0137 C-      select only one (based on minimum global-map index)
                0138         _BARRIER
                0139         _BEGIN_MASTER( myThid )
                0140         WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_SET_CALC: ',
                0141      &    'setting indices iPsi0,jPsi0 where Psi == 0 :'
                0142         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0143      &                      SQUEEZE_RIGHT, myThid )
                0144         WRITE(msgBuf,'(A,1P1E19.6,A,0PF16.3)')
                0145      &    'DIAGNOSTICS_SET_CALC: d2Min=',d2Min, ', ijMin=',xyMin
                0146         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0147      &                      SQUEEZE_RIGHT, myThid )
d7550daaf0 Jean*0148         IF ( xyMin.EQ.xyLoc ) THEN
c19aee9e8e Jean*0149           WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_SET_CALC: ',
                0150      &      'Fail to find the minimum distance => use Default'
                0151           CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0152      &                        SQUEEZE_RIGHT, myThid )
                0153           DO bj=1,nSy
                0154            DO bi=1,nSx
                0155              iPsi0(bi,bj) = -1
                0156              jPsi0(bi,bj) =  0
                0157            ENDDO
                0158           ENDDO
                0159         ELSE
                0160           DO bj=1,nSy
                0161            DO bi=1,nSx
                0162              IF ( iPsi0(bi,bj).GT.0 ) THEN
                0163               biG = bi-1+(myXGlobalLo-1)/sNx
                0164               bjG = bj-1+(myYGlobalLo-1)/sNy
                0165               xy0 = biG*(sNx+1) + bjG*(sNy+1)*xLine
                0166               xyLoc = xy0 + iPsi0(bi,bj) + (jPsi0(bi,bj)-1)*xLine
                0167               d2Loc = ABS( xyLoc - xyMin )
                0168               IF ( d2Loc.GE.0.5 _d 0 ) THEN
                0169                WRITE(msgBuf,'(2(A,2I5),A,F16.3)')
                0170      &          ' discard: bi,bj=',bi,bj,
                0171      &          ' ; i,j=',iPsi0(bi,bj),jPsi0(bi,bj),' ; ijLoc=',xyLoc
                0172                CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0173      &                             SQUEEZE_RIGHT, myThid )
                0174                iPsi0(bi,bj) = 0
                0175                jPsi0(bi,bj) = 0
                0176               ELSE
                0177                WRITE(msgBuf,'(2(A,2I5),A,F16.3)')
                0178      &          ' SELECT : bi,bj=',bi,bj,
                0179      &          ' ; i,j=',iPsi0(bi,bj),jPsi0(bi,bj),' ; ijLoc=',xyLoc
                0180                CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0181      &                             SQUEEZE_RIGHT, myThid )
                0182               ENDIF
                0183              ENDIF
20f8314a73 Jean*0184 c            WRITE(standardMessageUnit,'(2(A,2I5))')
                0185 c    &        ' bi,bj=',bi,bj,' i,jPsi0=', iPsi0(bi,bj),jPsi0(bi,bj)
c19aee9e8e Jean*0186            ENDDO
                0187           ENDDO
                0188         ENDIF
                0189         WRITE(msgBuf,'(2A)')
                0190      &   '------------------------------------------------------------'
                0191         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0192      &                      SQUEEZE_RIGHT, myThid )
                0193         _END_MASTER( myThid )
                0194         _BARRIER
                0195       ENDIF
                0196 
                0197       RETURN
                0198       END