pkg/layers/layers_locate.F

f0f170d54b Mart*0001 #include "LAYERS_OPTIONS.h"
                0002
                0003 CBOP
                0004       SUBROUTINE LAYERS_LOCATE(
                0005      I                          xx, n, m, sNx, sNy, x,
                0006      O                          k,
                0007      I                          myThid )
                0008
                0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
                0011 C     | Find the index(-array) k such that x is between xx(k)
                0012 C     | and xx(k+1) by bisection, following Press et al.,
                0013 C     | Numerical Recipes in Fortran. xx must be monotonic.
                0014 C     *==========================================================*
                0015 C     \ev
                0016
                0017 C !USES:
                0018       IMPLICIT NONE
                0019 C !INPUT PARAMETERS:
                0020 C     xx        :: array of bin-boundaries (layers_boundaries)
                0021 C     n         :: length of xx
                0022 C     m         :: int(log2(n)) + 1 = length of bisection loop
                0023 C     sNx,sNy   :: size of index array and input x
                0024 C     x         :: input array of values
                0025 C     k         :: index array (output)
                0026 C     myThid    :: my Thread Id number
                0027       INTEGER n, m, sNx, sNy
                0028       _RL     xx(1:n+1)
                0029       _RL     x(sNx+1,sNy+1)
                0030       INTEGER k(sNx+1,sNy+1)
                0031       INTEGER myThid
                0032
                0033 C !LOCAL VARIABLES:
                0034 C     i,j      :: horizontal indices
                0035 C     l        :: bisection loop index
                0036 C     kl,ku,km :: work arrays and variables
                0037       INTEGER i, j
                0038 CEOP
                0039 #ifdef TARGET_NEC_SX
                0040       INTEGER l, km
                0041       INTEGER kl(sNx+1,sNy+1), ku(sNx+1,sNy+1)
                0042
                0043 C     bisection, following Press et al., Numerical Recipes in Fortran,
                0044 C     mostly, because it can be vectorized
                0045       DO j = 1,sNy+1
                0046        DO i = 1,sNx+1
                0047         kl(i,j)=1
                0048         ku(i,j)=n+1
                0049        ENDDO
                0050       ENDDO
                0051       DO l = 1,m
                0052        DO j = 1,sNy+1
                0053         DO i = 1,sNx+1
                0054          IF (ku(i,j)-kl(i,j).GT.1) THEN
                0055           km=(ku(i,j)+kl(i,j))/2
                0056 CML       IF ((xx(n).GE.xx(1)).EQV.(x(i,j).GE.xx(km))) THEN
                0057           IF ( ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))).OR.
                0058      &         ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))) ) THEN
                0059            kl(i,j)=km
                0060           ELSE
                0061            ku(i,j)=km
                0062           ENDIF
                0063          ENDIF
                0064         ENDDO
                0065        ENDDO
                0066       ENDDO
                0067       DO j = 1,sNy+1
                0068        DO i = 1,sNx+1
                0069         IF ( x(i,j).LT.xx(2) ) THEN
                0070          k(i,j)=1
                0071         ELSEIF ( x(i,j).GE.xx(n) ) THEN
                0072          k(i,j)=n
                0073         ELSE
                0074          k(i,j)=kl(i,j)
                0075         ENDIF
                0076        ENDDO
                0077       ENDDO
                0078 #else /* TARGET_NEC_SX */
                0079 C     the old way
                0080       DO j = 1,sNy+1
                0081        DO i = 1,sNx+1
                0082         IF (x(i,j) .GE. xx(n)) THEN
                0083 C     the point is in the hottest bin or hotter
                0084          k(i,j) = n
                0085         ELSEIF (x(i,j) .LT. xx(2)) THEN
                0086 C        the point is in the coldest bin or colder
                0087          k(i,j) = 1
                0088         ELSEIF ( (x(i,j) .GE. xx(k(i,j)))
                0089      &    .AND.   (x(i,j) .LT. xx(k(i,j)+1)) ) THEN
                0090 C     already on the right bin -- do nothing
                0091         ELSEIF (x(i,j) .GE. xx(k(i,j))) THEN
                0092 C     have to hunt for the right bin by getting hotter
                0093          DO WHILE (x(i,j) .GE. xx(k(i,j)+1))
                0094           k(i,j) = k(i,j) + 1
                0095          ENDDO
                0096 C     now xx(k) < x <= xx(k+1)
                0097         ELSEIF (x(i,j) .LT. xx(k(i,j)+1)) THEN
                0098 C     have to hunt for the right bin by getting colder
                0099          DO WHILE (x(i,j) .LT. xx(k(i,j)))
                0100           k(i,j) = k(i,j) - 1
                0101          ENDDO
                0102 C     now xx(k) <= x < xx(k+1)
                0103         ELSE
                0104 C     that should have covered all the options
                0105          k(i,j) = -1
                0106         ENDIF
                0107
                0108        ENDDO
                0109       ENDDO
                0110 #endif /* TARGET_NEC_SX */
                0111
                0112 #endif /* ALLOW_LAYERS */
                0113
                0114       RETURN
                0115       END