pkg/layers/layers_fluxcalc.F

aa668e01f1 Gael*0001 #include "LAYERS_OPTIONS.h"
                0002 #ifdef ALLOW_GMREDI
                0003 #include "GMREDI_OPTIONS.h"
                0004 #endif
                0005
28ff9bcfe6 Jean*0006 C--  File layers_fluxcalc.F:
9479dbe556 Mart*0007 C--   Contents
                0008 C--   o LAYERS_FLUXCALC
cf336ab6c5 Ryan*0009 C--   o LAYERS_DIAPYCNAL
9479dbe556 Mart*0010 C--   o LAYERS_LOCATE
                0011
aa668e01f1 Gael*0012 CBOP 0
9479dbe556 Mart*0013 C     !ROUTINE: LAYERS_FLUXCALC
                0014 C     !INTERFACE:
aa668e01f1 Gael*0015       SUBROUTINE LAYERS_FLUXCALC(
406891c1a3 Gael*0016      I                  uVel,vVel,tracer,iLa,
50d8304171 Ryan*0017      O                  UH,VH,Hw,Hs,PIw,PIs,Uw,Vs,
aa668e01f1 Gael*0018      I                  myThid )
                0019
9479dbe556 Mart*0020 C     !DESCRIPTION: \bv
                0021 C     *==========================================================*
                0022 C     | SUBROUTINE LAYERS_FLUXCALC
                0023 C     | Calculate the transport in isotracer layers, for a chosen
                0024 C     | tracer. This is the meat of the LAYERS package.
                0025 C     *==========================================================*
                0026 C     \ev
aa668e01f1 Gael*0027
                0028 C !USES:
                0029       IMPLICIT NONE
9479dbe556 Mart*0030 C     == Global variables ===
aa668e01f1 Gael*0031 #include "SIZE.h"
                0032 #include "EEPARAMS.h"
                0033 #include "PARAMS.h"
                0034 #include "GRID.h"
                0035 #include "LAYERS_SIZE.h"
                0036 #include "LAYERS.h"
                0037 #ifdef ALLOW_GMREDI
                0038 # include "GMREDI.h"
                0039 #endif
                0040
                0041 C !INPUT PARAMETERS:
                0042 C     myThid    :: my Thread Id number
                0043 C     uVel  :: zonal velocity (m/s, i=1 held at western face)
                0044 C     vVel  :: meridional velocity (m/s, j=1 held at southern face)
d746ff23b0 Jean*0045 C     tracer :: potential temperature, salt or potential density prho
aa668e01f1 Gael*0046 C      UH   :: U integrated over layer (m^2/s)
                0047 C      VH   :: V integrated over layer (m^2/s)
                0048 C      Hw   :: Layer thickness at the U point (m)
                0049 C      Hs   :: Layer thickness at the V point (m)
14bb2ae926 Ryan*0050 C      PIw  :: 1 if layer exists, 0 otherwise (at U point)
                0051 C      PIs  :: 1 if layer exists, 0 otherwise (at V point)
50d8304171 Ryan*0052 C      Uw   :: average U over layer (m/s)
                0053 C      Vs   :: average V over layer (m/s)
406891c1a3 Gael*0054 C      iLa  :: layer coordinate index
aa668e01f1 Gael*0055       INTEGER myThid
28ff9bcfe6 Jean*0056       _RL uVel   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,     nSx,nSy)
                0057       _RL vVel   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,     nSx,nSy)
                0058       _RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,     nSx,nSy)
6f9adcbc06 Jean*0059 #ifdef LAYERS_UFLUX
28ff9bcfe6 Jean*0060       _RL UH     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
6f9adcbc06 Jean*0061 # ifdef LAYERS_THICKNESS
28ff9bcfe6 Jean*0062       _RL Hw     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
                0063       _RL PIw    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
50d8304171 Ryan*0064       _RL Uw     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
6f9adcbc06 Jean*0065 # else
50d8304171 Ryan*0066       _RL Hw(1), PIw(1), Uw(1)
6f9adcbc06 Jean*0067 # endif
                0068 #else
06e3ddabf4 Jean*0069       _RL UH(1), Hw(1), PIw(1), Uw(1)
6f9adcbc06 Jean*0070 #endif
                0071 #ifdef LAYERS_VFLUX
                0072       _RL VH     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
                0073 # ifdef LAYERS_THICKNESS
                0074       _RL Hs     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
                0075       _RL PIs    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
50d8304171 Ryan*0076       _RL Vs     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
6f9adcbc06 Jean*0077 # else
50d8304171 Ryan*0078       _RL Hs(1), PIs(1), Vs(1)
6f9adcbc06 Jean*0079 # endif
                0080 #else
50d8304171 Ryan*0081       _RL VH(1), Hs(1), PIs(1), Vs(1)
6f9adcbc06 Jean*0082 #endif
406891c1a3 Gael*0083       INTEGER iLa
aa668e01f1 Gael*0084 CEOP
                0085
                0086 #ifdef ALLOW_LAYERS
                0087
                0088 C !LOCAL VARIABLES:
                0089 C     bi, bj   :: tile indices
                0090 C     i,j      :: horizontal indices
                0091 C     k        :: vertical index for model grid
                0092 C     kci      :: index from CellIndex
406891c1a3 Gael*0093 C     kg       :: index for looping though layers_bounds
aa668e01f1 Gael*0094 C     kk       :: vertical index for ZZ (fine) grid
                0095 C     kgu,kgv  :: vertical index for isopycnal grid
9479dbe556 Mart*0096 C     kloc     :: local copy of kgu/v to reduce accesses to index arrays
                0097 C     mSteps   :: maximum number of steps for bisection method
6f9adcbc06 Jean*0098 C     prho     :: pot. density (less 1000) referenced to layers_krho pressure
5b31d86392 Mart*0099 C     TatU     :: temperature at U point
aa668e01f1 Gael*0100 C     TatV     :: temperature at V point
06e3ddabf4 Jean*0101 C     dzfac    :: temporary sublayer thickness
50d8304171 Ryan*0102 C     Tloc,Tp1 :: horizontally interpolated tracer values
aa668e01f1 Gael*0103
                0104       INTEGER bi, bj
9479dbe556 Mart*0105       INTEGER i,j,k,kk,kg,kci,kp1,kloc
                0106       INTEGER mSteps
aa668e01f1 Gael*0107       INTEGER kgu(sNx+1,sNy+1), kgv(sNx+1,sNy+1)
5b31d86392 Mart*0108       _RL TatU(sNx+1,sNy+1), TatV(sNx+1,sNy+1)
f6012f0463 Ryan*0109       _RL dzfac
d746ff23b0 Jean*0110 #ifdef ALLOW_GMREDI
aa668e01f1 Gael*0111       INTEGER kcip1
                0112       _RL delPsi, maskp1
                0113 #endif
9479dbe556 Mart*0114       LOGICAL errorFlag
                0115       CHARACTER*(MAX_LEN_MBUF) msgBuf
aa668e01f1 Gael*0116
                0117 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0118
28ff9bcfe6 Jean*0119 C     compute maximum number of steps for bisection method (approx.
9479dbe556 Mart*0120 C     log2(Nlayers)) as log2(Nlayers) + 1 for safety
                0121       mSteps = int(log10(dble(Nlayers))/log10(2. _d 0))+1
                0122
aa668e01f1 Gael*0123 C --- The tile loops
                0124       DO bj=myByLo(myThid),myByHi(myThid)
                0125       DO bi=myBxLo(myThid),myBxHi(myThid)
                0126
                0127 C     Initialize the search indices
                0128       DO j = 1,sNy+1
                0129         DO i = 1,sNx+1
                0130 C       The temperature index (layer_G) goes from cold to warm.
                0131 C       The water column goes from warm (k=1) to cold (k=Nr).
                0132 C       So initialize the search with the warmest value.
                0133           kgu(i,j) = Nlayers
                0134           kgv(i,j) = Nlayers
                0135         ENDDO
                0136       ENDDO
                0137
                0138 C     Reset the arrays
                0139       DO kg=1,Nlayers
28ff9bcfe6 Jean*0140        DO j = 1-OLy,sNy+OLy
                0141         DO i = 1-OLx,sNx+OLx
aa668e01f1 Gael*0142 #ifdef LAYERS_UFLUX
1377839540 Jean*0143          UH (i,j,kg,bi,bj) = 0. _d 0
aa668e01f1 Gael*0144 #ifdef LAYERS_THICKNESS
50d8304171 Ryan*0145          Hw(i,j,kg,bi,bj) = 0. _d 0
1377839540 Jean*0146          PIw(i,j,kg,bi,bj) = 0. _d 0
50d8304171 Ryan*0147          Uw(i,j,kg,bi,bj) = 0. _d 0
aa668e01f1 Gael*0148 #endif /* LAYERS_THICKNESS */
                0149 #endif /* UH */
                0150 #ifdef LAYERS_VFLUX
1377839540 Jean*0151          VH (i,j,kg,bi,bj) = 0. _d 0
aa668e01f1 Gael*0152 #ifdef LAYERS_THICKNESS
50d8304171 Ryan*0153          Hs(i,j,kg,bi,bj) = 0. _d 0
1377839540 Jean*0154          PIs(i,j,kg,bi,bj) = 0. _d 0
50d8304171 Ryan*0155          Vs(i,j,kg,bi,bj) = 0. _d 0
aa668e01f1 Gael*0156 #endif /* LAYERS_THICKNESS */
                0157 #endif /* VH */
                0158         ENDDO
                0159        ENDDO
                0160       ENDDO
                0161
                0162       DO kk=1,NZZ
                0163        k = MapIndex(kk)
                0164        kci = CellIndex(kk)
d746ff23b0 Jean*0165 #ifdef ALLOW_GMREDI
                0166        kcip1 = MIN(kci+1,Nr)
                0167        maskp1 = 1.
                0168        IF (kci.GE.Nr) maskp1 = 0.
                0169 #endif /* ALLOW_GMREDI */
5b31d86392 Mart*0170 #ifdef LAYERS_UFLUX
aa668e01f1 Gael*0171        DO j = 1,sNy+1
                0172         DO i = 1,sNx+1
                0173
                0174 C ------ Find theta at the U point (west) on the fine Z grid
                0175          kp1=k+1
f6012f0463 Ryan*0176          IF (maskW(i,j,kp1,bi,bj).EQ.zeroRS) kp1=k
                0177          TatU(i,j) = MapFact(kk) *
                0178      &    0.5 _d 0 * (tracer(i-1,j,k,bi,bj)+tracer(i,j,k,bi,bj)) +
                0179      &    (1. _d 0 -MapFact(kk)) *
                0180      &    0.5 _d 0 * (tracer(i-1,j,kp1,bi,bj)+tracer(i,j,kp1,bi,bj))
aa668e01f1 Gael*0181
5b31d86392 Mart*0182         ENDDO
                0183        ENDDO
aa668e01f1 Gael*0184 C ------ Now that we know T everywhere, determine the binning.
9479dbe556 Mart*0185 C        find the layer indices kgu
                0186        CALL LAYERS_LOCATE(
28ff9bcfe6 Jean*0187      I      layers_bounds(1,iLa),Nlayers,mSteps,sNx,sNy,TatU,
9479dbe556 Mart*0188      O      kgu,
                0189      I      myThid )
28ff9bcfe6 Jean*0190 #ifndef TARGET_NEC_SX
9479dbe556 Mart*0191 C     check for failures
                0192        IF ( debugLevel .GE. debLevC ) THEN
                0193         errorFlag = .FALSE.
                0194         DO j = 1,sNy+1
                0195          DO i = 1,sNx+1
                0196           IF ( kgu(i,j) .LE. 0 ) THEN
                0197            WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)')
                0198      &          'S/R LAYERS_LOCATE: Could not find a bin in ',
                0199      &          'layers_bounds for TatU(',i,',',j,',)=',TatU(i,j)
                0200            CALL PRINT_ERROR( msgBuf, myThid )
                0201            errorFlag = .TRUE.
                0202           ENDIF
                0203          ENDDO
5b31d86392 Mart*0204         ENDDO
9479dbe556 Mart*0205         IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_FLUXCALC'
                0206        ENDIF
                0207 #endif /* ndef TARGET_NEC_SX */
28ff9bcfe6 Jean*0208 C
5b31d86392 Mart*0209        DO j = 1,sNy+1
                0210         DO i = 1,sNx+1
                0211
9479dbe556 Mart*0212          kloc = kgu(i,j)
f6012f0463 Ryan*0213          dzfac = dZZf(kk) * hFacW(i,j,kci,bi,bj)
                0214
aa668e01f1 Gael*0215 C ------ Augment the bin values
9479dbe556 Mart*0216          UH(i,j,kloc,bi,bj) =
                0217      &    UH(i,j,kloc,bi,bj) +
f6012f0463 Ryan*0218      &    dzfac * uVel(i,j,kci,bi,bj)
aa668e01f1 Gael*0219
                0220 #ifdef ALLOW_GMREDI
d746ff23b0 Jean*0221          IF ( layers_bolus(iLa)  ) THEN
                0222            IF ( .NOT.GM_AdvForm ) THEN
9ecea62500 Jean*0223              delPsi = 0.25 _d 0 *(
                0224      &              ( rA(i-1,j,bi,bj)*Kwx(i-1,j,kcip1,bi,bj)
                0225      &               +rA( i ,j,bi,bj)*Kwx( i ,j,kcip1,bi,bj)
                0226      &              ) * maskW(i,j,kcip1,bi,bj) * maskp1
                0227      &            - ( rA(i-1,j,bi,bj)*Kwx(i-1,j, kci ,bi,bj)
                0228      &               +rA( i ,j,bi,bj)*Kwx( i ,j, kci ,bi,bj)
                0229      &              ) * maskW(i,j, kci ,bi,bj)
                0230      &                           ) * recip_rAw(i,j,bi,bj)
d746ff23b0 Jean*0231 #ifdef GM_BOLUS_ADVEC
                0232            ELSE
                0233              delPsi = GM_PsiX(i,j,kcip1,bi,bj)*maskp1
                0234      &              - GM_PsiX(i,j, kci, bi,bj)
                0235 #endif
                0236            ENDIF
9479dbe556 Mart*0237            UH(i,j,kloc,bi,bj) = UH(i,j,kloc,bi,bj)
aa668e01f1 Gael*0238      &      + delPsi*recip_drF(kci)*_recip_hFacW(i,j,kci,bi,bj)
f6012f0463 Ryan*0239      &      * dzfac
aa668e01f1 Gael*0240          ENDIF
d746ff23b0 Jean*0241 #endif /* ALLOW_GMREDI */
aa668e01f1 Gael*0242
                0243 #ifdef LAYERS_THICKNESS
f6012f0463 Ryan*0244          Hw(i,j,kloc,bi,bj) = Hw(i,j,kloc,bi,bj) + dzfac
aa668e01f1 Gael*0245 #endif /* LAYERS_THICKNESS */
                0246
5b31d86392 Mart*0247         ENDDO
                0248        ENDDO
aa668e01f1 Gael*0249 #endif /* LAYERS_UFLUX */
                0250
                0251 #ifdef LAYERS_VFLUX
5b31d86392 Mart*0252        DO j = 1,sNy+1
                0253         DO i = 1,sNx+1
aa668e01f1 Gael*0254 C ------ Find theta at the V point (south) on the fine Z grid
                0255          kp1=k+1
f6012f0463 Ryan*0256          IF (maskS(i,j,kp1,bi,bj).EQ.zeroRS) kp1=k
                0257          TatV(i,j) = MapFact(kk) *
                0258      &    0.5 _d 0 * (tracer(i,j-1,k,bi,bj)+tracer(i,j,k,bi,bj)) +
                0259      &    (1. _d 0 -MapFact(kk)) *
                0260      &    0.5 _d 0 * (tracer(i,j-1,kp1,bi,bj)+tracer(i,j,kp1,bi,bj))
aa668e01f1 Gael*0261
5b31d86392 Mart*0262         ENDDO
                0263        ENDDO
9479dbe556 Mart*0264 C ------ Now that we know T everywhere, determine the binning.
                0265 C        find the layer indices kgv
                0266        CALL LAYERS_LOCATE(
28ff9bcfe6 Jean*0267      I      layers_bounds(1,iLa),Nlayers,mSteps,sNx,sNy,TatV,
9479dbe556 Mart*0268      O      kgv,
                0269      I      myThid )
                0270 #ifndef TARGET_NEC_SX
                0271        IF ( debugLevel .GE. debLevC ) THEN
                0272 C     check for failures
                0273         errorFlag = .FALSE.
                0274         DO j = 1,sNy+1
                0275          DO i = 1,sNx+1
                0276           IF ( kgv(i,j) .LE. 0 ) THEN
                0277            WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)')
                0278      &          'S/R LAYERS_LOCATE: Could not find a bin in ',
                0279      &          'layers_bounds for TatV(',i,',',j,',)=',TatV(i,j)
                0280            CALL PRINT_ERROR( msgBuf, myThid )
                0281            errorFlag = .TRUE.
                0282           ENDIF
                0283          ENDDO
5b31d86392 Mart*0284         ENDDO
9479dbe556 Mart*0285         IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_FLUXCALC'
                0286        ENDIF
                0287 #endif /* ndef TARGET_NEC_SX */
                0288 C
5b31d86392 Mart*0289        DO j = 1,sNy+1
                0290         DO i = 1,sNx+1
                0291
9479dbe556 Mart*0292          kloc = kgv(i,j)
50d8304171 Ryan*0293          dzfac = dZZf(kk) * hFacS(i,j,kci,bi,bj)
                0294
                0295 C ------ debugging stuff
                0296 C         IF (i.EQ.10 .AND. j.EQ.10) THEN
                0297 C           WRITE(msgBuf,'(A,I3,A,F10.2,A,F6.2)')
                0298 C     &          '    kloc=', kloc,
                0299 C     &          ', TatV=',TatV(i,j),
                0300 C     &          ', dzfac=',dzfac
                0301 C           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0302 C     &                         SQUEEZE_RIGHT, myThid )
                0303 C         ENDIF
                0304
aa668e01f1 Gael*0305 C ------ Augment the bin values
9479dbe556 Mart*0306          VH(i,j,kloc,bi,bj) =
50d8304171 Ryan*0307      &    VH(i,j,kloc,bi,bj) + dzfac * vVel(i,j,kci,bi,bj)
aa668e01f1 Gael*0308
                0309 #ifdef ALLOW_GMREDI
d746ff23b0 Jean*0310          IF ( layers_bolus(iLa) ) THEN
                0311            IF ( .NOT.GM_AdvForm ) THEN
9ecea62500 Jean*0312              delPsi = 0.25 _d 0 *(
                0313      &              ( rA(i,j-1,bi,bj)*Kwy(i,j-1,kcip1,bi,bj)
                0314      &               +rA(i, j ,bi,bj)*Kwy(i, j ,kcip1,bi,bj)
                0315      &              ) * maskS(i,j,kcip1,bi,bj) * maskp1
                0316      &            - ( rA(i,j-1,bi,bj)*Kwy(i,j-1, kci ,bi,bj)
                0317      &               +rA(i, j ,bi,bj)*Kwy(i, j , kci ,bi,bj)
                0318      &              ) * maskS(i,j, kci ,bi,bj)
                0319      &                           ) * recip_rAs(i,j,bi,bj)
d746ff23b0 Jean*0320 #ifdef GM_BOLUS_ADVEC
                0321            ELSE
                0322              delPsi = GM_PsiY(i,j,kcip1,bi,bj)*maskp1
                0323      &              - GM_PsiY(i,j, kci, bi,bj)
                0324 #endif
                0325            ENDIF
9479dbe556 Mart*0326            VH(i,j,kloc,bi,bj) = VH(i,j,kloc,bi,bj)
aa668e01f1 Gael*0327      &      + delPsi*recip_drF(kci)*_recip_hFacS(i,j,kci,bi,bj)
50d8304171 Ryan*0328      &      * dzfac
aa668e01f1 Gael*0329          ENDIF
d746ff23b0 Jean*0330 #endif /* ALLOW_GMREDI */
aa668e01f1 Gael*0331
                0332 #ifdef LAYERS_THICKNESS
50d8304171 Ryan*0333          Hs(i,j,kloc,bi,bj) = Hs(i,j,kloc,bi,bj) + dzfac
aa668e01f1 Gael*0334 #endif /* LAYERS_THICKNESS */
                0335
                0336         ENDDO
                0337        ENDDO
5b31d86392 Mart*0338 #endif /* LAYERS_VFLUX */
aa668e01f1 Gael*0339       ENDDO
d746ff23b0 Jean*0340
14bb2ae926 Ryan*0341 C--   Now that we know the thicknesses, compute the heaviside function
                0342 C--   (Needs another loop through Ng)
                0343 #ifdef LAYERS_THICKNESS
                0344       DO kg=1,Nlayers
                0345        DO j = 1,sNy+1
                0346         DO i = 1,sNx+1
                0347 #ifdef LAYERS_UFLUX
                0348          IF (Hw(i,j,kg,bi,bj) .GT. 0.) THEN
                0349           PIw(i,j,kg,bi,bj) = 1. _d 0
50d8304171 Ryan*0350           Uw(i,j,kg,bi,bj) =
14bb2ae926 Ryan*0351      &        UH(i,j,kg,bi,bj) / Hw(i,j,kg,bi,bj)
                0352          ENDIF
                0353 #endif /* LAYERS_UFLUX */
                0354 #ifdef LAYERS_VFLUX
                0355          IF (Hs(i,j,kg,bi,bj) .GT. 0.) THEN
                0356           PIs(i,j,kg,bi,bj) = 1. _d 0
50d8304171 Ryan*0357           Vs(i,j,kg,bi,bj) =
14bb2ae926 Ryan*0358      &        VH(i,j,kg,bi,bj) / Hs(i,j,kg,bi,bj)
                0359          ENDIF
                0360 #endif /* LAYERS_VFLUX */
                0361         ENDDO
                0362        ENDDO
                0363       ENDDO
                0364 #endif /* LAYERS_THICKNESS */
aa668e01f1 Gael*0365
                0366 C --- End bi,bj loop
                0367       ENDDO
                0368       ENDDO
                0369
9479dbe556 Mart*0370       RETURN
                0371       END
28ff9bcfe6 Jean*0372
                0373 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
cf336ab6c5 Ryan*0374 CBOP 0
                0375 C     !ROUTINE: LAYERS_DIAPYCNAL
                0376 C     !INTERFACE:
                0377       SUBROUTINE LAYERS_DIAPYCNAL(
                0378      I                  tracer,iLa,
                0379      O                  TtendSurf, TtendDiffh, TtendDiffr,
50d8304171 Ryan*0380      O                  TtendAdvh, TtendAdvr, Ttendtot,
cf336ab6c5 Ryan*0381      O                  StendSurf, StendDiffh, StendDiffr,
50d8304171 Ryan*0382      O                  StendAdvh, StendAdvr, Stendtot,
cf336ab6c5 Ryan*0383      O                  Hc, PIc,
                0384      I                  myThid )
                0385
                0386 C     !DESCRIPTION: \bv
                0387 C     *==========================================================*
                0388 C     | SUBROUTINE LAYERS_DIAPYCNAL
                0389 C     | Calculate the diapycnal velocity in isotracer layers, for a chosen
                0390 C     | tracer.
                0391 C     *==========================================================*
                0392 C     \ev
                0393       IMPLICIT NONE
                0394 #include "SIZE.h"
                0395 #include "EEPARAMS.h"
                0396 #include "PARAMS.h"
                0397 #include "GRID.h"
                0398 #include "LAYERS_SIZE.h"
                0399 #include "LAYERS.h"
                0400
                0401 C !INPUT PARAMETERS:
                0402 C     myThid    :: my Thread Id number
                0403 C     tracer    :: potential temperature, salt or potential density prho
                0404 C     iLa       :: layer coordinate index
1377839540 Jean*0405 C     TtendSurf :: temperature tendency due to surface forcing times thickness
cf336ab6c5 Ryan*0406 C     TtendDiffh:: temperature tendency due to horiz. diffusion times thickness
                0407 C     TtendDiffr:: temperature tendency due to vert. diffusion times thickness
50d8304171 Ryan*0408 C     TtendAdvh:: salinity tendency due to horiz. advection times thickness
                0409 C     TtendAdvr:: salinity tendency due to vert. advection times thickness
06e3ddabf4 Jean*0410 C     StendSurf :: salinity tendency due to surface forcing times thickness
50d8304171 Ryan*0411 C     StendDiffh:: salinity tendency due to horiz. diffusion times thickness
                0412 C     StendDiffr:: salinity tendency due to vert. diffusion times thickness
                0413 C     StendAdvh :: salinity tendency due to horiz. advection times thickness
                0414 C     StendAdvr :: salinity tendency due to vert. advection times thickness
cf336ab6c5 Ryan*0415 C     Hc        :: Layer thickness at the tracer point (m)
                0416 C     PIw       :: 1 if layer exists, 0 otherwise (at tracer point)
                0417       INTEGER iLa, myThid
                0418       _RL tracer    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,     nSx,nSy)
                0419       _RL Hc        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
                0420       _RL PIc       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
50d8304171 Ryan*0421       _RL TtendSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0422       _RL TtendDiffh(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
06e3ddabf4 Jean*0423       _RL TtendDiffr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
50d8304171 Ryan*0424       _RL TtendAdvh (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
06e3ddabf4 Jean*0425       _RL TtendAdvr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0426       _RL Ttendtot  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
50d8304171 Ryan*0427       _RL StendSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0428       _RL StendDiffh(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
06e3ddabf4 Jean*0429       _RL StendDiffr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
50d8304171 Ryan*0430       _RL StendAdvh (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
06e3ddabf4 Jean*0431       _RL StendAdvr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
                0432       _RL Stendtot  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
cf336ab6c5 Ryan*0433 CEOP
                0434
                0435 #ifdef LAYERS_THERMODYNAMICS
                0436
                0437 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0438 C !LOCAL VARIABLES:
                0439 C     bi, bj   :: tile indices
                0440 C     i,j      :: horizontal indices
                0441 C     k        :: vertical index for model grid
f6012f0463 Ryan*0442 C     kp1      :: vertical index for model grid next cell
cf336ab6c5 Ryan*0443 C     kci      :: index from CellIndex
                0444 C     kg       :: index for looping though layers_bounds
                0445 C     kk       :: vertical index for ZZ (fine) grid
                0446 C     kloc     :: local copy of kgu/v to reduce accesses to index arrays
                0447 C     mSteps   :: maximum number of steps for bisection method
                0448 C     TatC     :: temperature at C point
da9f56e003 Jean*0449       _RL Hcw       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
cf336ab6c5 Ryan*0450       INTEGER bi, bj
da9f56e003 Jean*0451       INTEGER i,j,k,kk,kg,kci,kloc
cf336ab6c5 Ryan*0452       INTEGER mSteps
                0453       INTEGER kgc(sNx+1,sNy+1)
06e3ddabf4 Jean*0454       INTEGER kgcw(sNx+1,sNy+1)
f6012f0463 Ryan*0455       _RL TatC(sNx+1,sNy+1), dzfac, Tfac, Sfac
cf336ab6c5 Ryan*0456       LOGICAL errorFlag
                0457       CHARACTER*(MAX_LEN_MBUF) msgBuf
da9f56e003 Jean*0458 #ifdef LAYERS_FINEGRID_DIAPYCNAL
                0459       INTEGER kp1
                0460 #endif
cf336ab6c5 Ryan*0461
50d8304171 Ryan*0462 C --  constants for T and S forcing, gets reset later for rho
                0463       Tfac = 1. _d 0
                0464       Sfac = 1. _d 0
06e3ddabf4 Jean*0465
cf336ab6c5 Ryan*0466 C     compute maximum number of steps for bisection method (approx.
                0467 C     log2(Nlayers)) as log2(Nlayers) + 1 for safety
                0468       mSteps = int(log10(dble(Nlayers))/log10(2. _d 0))+1
                0469
                0470 C      STOP 'DEBUG END: S/R LAYERS_DIAPYCNAL'
                0471
                0472 C --- The tile loops
                0473       DO bj=myByLo(myThid),myByHi(myThid)
                0474       DO bi=myBxLo(myThid),myBxHi(myThid)
                0475
                0476 C     Initialize the search indices
                0477       DO j = 1,sNy+1
                0478         DO i = 1,sNx+1
                0479 C       The temperature index (layer_G) goes from cold to warm.
                0480 C       The water column goes from warm (k=1) to cold (k=Nr).
                0481 C       So initialize the search with the warmest value.
                0482           kgc(i,j) = Nlayers
50d8304171 Ryan*0483           kgcw(i,j) = Nlayers-1
cf336ab6c5 Ryan*0484         ENDDO
                0485       ENDDO
                0486
                0487 C     Reset the arrays
50d8304171 Ryan*0488 C --- These are at the w point
                0489       DO kg=1,Nlayers-1
cf336ab6c5 Ryan*0490        DO j = 1-OLy,sNy+OLy
                0491         DO i = 1-OLx,sNx+OLx
                0492          TtendSurf (i,j,kg,bi,bj) = 0. _d 0
                0493          TtendDiffh(i,j,kg,bi,bj) = 0. _d 0
                0494          TtendDiffr(i,j,kg,bi,bj) = 0. _d 0
50d8304171 Ryan*0495          TtendAdvh(i,j,kg,bi,bj)  = 0. _d 0
                0496          TtendAdvr(i,j,kg,bi,bj)  = 0. _d 0
                0497          Ttendtot(i,j,kg,bi,bj)   = 0. _d 0
cf336ab6c5 Ryan*0498          StendSurf (i,j,kg,bi,bj) = 0. _d 0
                0499          StendDiffh(i,j,kg,bi,bj) = 0. _d 0
06e3ddabf4 Jean*0500          StendDiffr(i,j,kg,bi,bj) = 0. _d 0
50d8304171 Ryan*0501          StendAdvh(i,j,kg,bi,bj)  = 0. _d 0
06e3ddabf4 Jean*0502          StendAdvr(i,j,kg,bi,bj)  = 0. _d 0
50d8304171 Ryan*0503          Stendtot(i,j,kg,bi,bj)   = 0. _d 0
06e3ddabf4 Jean*0504          Hcw(i,j,kg,bi,bj) = 0. _d 0
50d8304171 Ryan*0505         ENDDO
                0506        ENDDO
                0507       ENDDO
                0508 C --- These are at the c point
                0509       DO kg=1,Nlayers
                0510        DO j = 1-OLy,sNy+OLy
06e3ddabf4 Jean*0511         DO i = 1-OLx,sNx+OLx
cf336ab6c5 Ryan*0512          Hc(i,j,kg,bi,bj) = 0. _d 0
1377839540 Jean*0513          PIc(i,j,kg,bi,bj) = 0. _d 0
cf336ab6c5 Ryan*0514         ENDDO
                0515        ENDDO
                0516       ENDDO
                0517
cdcbace143 Ryan*0518 #ifdef LAYERS_FINEGRID_DIAPYCNAL
cf336ab6c5 Ryan*0519       DO kk=1,NZZ
                0520        k = MapIndex(kk)
                0521        kci = CellIndex(kk)
                0522        DO j = 1,sNy+1
                0523         DO i = 1,sNx+1
f6012f0463 Ryan*0524 C ------ Find theta at the V point (south) on the fine Z grid
                0525          kp1=k+1
                0526          IF (maskC(i,j,kp1,bi,bj).EQ.zeroRS) kp1=k
                0527          TatC(i,j) = MapFact(kk) * tracer(i,j,k,bi,bj) +
                0528      &    (1. _d 0 -MapFact(kk)) * tracer(i,j,kp1,bi,bj)
cdcbace143 Ryan*0529         ENDDO
                0530        ENDDO
                0531 #else
                0532       DO kk=1,Nr
                0533        k = kk
                0534        kci = kk
                0535        DO j = 1,sNy+1
                0536         DO i = 1,sNx+1
                0537          TatC(i,j) = tracer(i,j,k,bi,bj)
                0538         ENDDO
                0539        ENDDO
                0540 #endif /* LAYERS_FINEGRID_DIAPYCNAL */
da9f56e003 Jean*0541
f6012f0463 Ryan*0542 C ------ debugging stuff
                0543 c         IF (i.EQ.38 .AND. j.EQ.4 .AND. bi.EQ.1 .AND. bj.EQ.1) THEN
cdcbace143 Ryan*0544 c           i=38
                0545 c           j=4
f6012f0463 Ryan*0546 c           WRITE(msgBuf,
                0547 c     &       '(A,I3,A,I3,A,I3,A,F7.2,A,F7.2,A,F7.2,A,F7.2,A,F3.1)')
                0548 c     &          'LAYERS_DEBUG: iLa=', iLa,
                0549 c     &          ', kk=', kk,
                0550 c     &          ', k=', k,
                0551 c     &          ', tracer=', tracer(i,j,k,bi,bj),
                0552 c     &          ', TatC=',TatC(i,j),
                0553 c     &          ', hFacC=',hFacC(i,j,k,bi,bj)
                0554 c           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0555 c     &                         SQUEEZE_RIGHT, myThid )
                0556 c         ENDIF
cf336ab6c5 Ryan*0557 C ------ Now that we know T everywhere, determine the binning.
50d8304171 Ryan*0558 C        find the layer indices kgc for the center point
cf336ab6c5 Ryan*0559        CALL LAYERS_LOCATE(
                0560      I      layers_bounds(1,iLa),Nlayers,mSteps,sNx,sNy,TatC,
                0561      O      kgc,
                0562      I      myThid )
                0563 #ifndef TARGET_NEC_SX
                0564 C     check for failures
                0565        IF ( debugLevel .GE. debLevC ) THEN
                0566         errorFlag = .FALSE.
                0567         DO j = 1,sNy+1
                0568          DO i = 1,sNx+1
                0569           IF ( kgc(i,j) .LE. 0 ) THEN
                0570            WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)')
                0571      &          'S/R LAYERS_LOCATE: Could not find a bin in ',
                0572      &          'layers_bounds for TatC(',i,',',j,',)=',TatC(i,j)
                0573            CALL PRINT_ERROR( msgBuf, myThid )
                0574            errorFlag = .TRUE.
                0575           ENDIF
                0576          ENDDO
                0577         ENDDO
                0578         IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_DIAPYCNAL'
                0579        ENDIF
                0580 #endif /* ndef TARGET_NEC_SX */
50d8304171 Ryan*0581
                0582 C        find the layer indices kgcw for the w point
                0583        CALL LAYERS_LOCATE(
                0584      I      layers_bounds_w(1,iLa),Nlayers-1,mSteps,sNx,sNy,TatC,
                0585      O      kgcw,
                0586      I      myThid )
                0587 #ifndef TARGET_NEC_SX
                0588 C     check for failures
                0589        IF ( debugLevel .GE. debLevC ) THEN
                0590         errorFlag = .FALSE.
                0591         DO j = 1,sNy+1
                0592          DO i = 1,sNx+1
                0593           IF ( kgcw(i,j) .LE. 0 ) THEN
                0594            WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)')
                0595      &          'S/R LAYERS_LOCATE: Could not find a bin in ',
                0596      &          'layers_bounds for TatC(',i,',',j,',)=',TatC(i,j)
                0597            CALL PRINT_ERROR( msgBuf, myThid )
                0598            errorFlag = .TRUE.
                0599           ENDIF
                0600          ENDDO
                0601         ENDDO
                0602         IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_DIAPYCNAL'
                0603        ENDIF
                0604 #endif /* ndef TARGET_NEC_SX */
cf336ab6c5 Ryan*0605
1377839540 Jean*0606 C ------ Augment the bin values
cf336ab6c5 Ryan*0607        DO j = 1,sNy+1
                0608         DO i = 1,sNx+1
cdcbace143 Ryan*0609 #ifdef LAYERS_FINEGRID_DIAPYCNAL
06e3ddabf4 Jean*0610          dzfac = dZZf(kk) * hFacC(i,j,kci,bi,bj)
cdcbace143 Ryan*0611 #else
                0612          dzfac = dRf(kci) * hFacC(i,j,kci,bi,bj)
                0613 #endif /* LAYERS_FINEGRID_DIAPYCNAL */
50d8304171 Ryan*0614          kloc = kgcw(i,j)
                0615
                0616 C ------- Thickness at w point
                0617          Hcw(i,j,kloc,bi,bj) = Hcw(i,j,kloc,bi,bj)
                0618      &     + dzfac
                0619 C ------- Thickness at c point
                0620          Hc(i,j,kgc(i,j),bi,bj) = Hc(i,j,kgc(i,j),bi,bj)
                0621      &     + dzfac
06e3ddabf4 Jean*0622
50d8304171 Ryan*0623 C ------- Now rescale dzfac to include the layer coordinate spacing
                0624          dzfac = dzfac * layers_recip_delta(kloc,iLa)
                0625
84f265d4b3 dfer*0626 #ifdef LAYERS_PRHO_REF
50d8304171 Ryan*0627          IF ( layers_num(iLa) .EQ. 3 ) THEN
                0628            Tfac = layers_alpha(i,j,kci,bi,bj)
                0629            Sfac = layers_beta(i,j,kci,bi,bj)
                0630          ENDIF
84f265d4b3 dfer*0631 #endif
cf336ab6c5 Ryan*0632          IF (kci.EQ.1) THEN
                0633 C ------- We are in the surface layer
                0634           TtendSurf(i,j,kloc,bi,bj) =
06e3ddabf4 Jean*0635      &     TtendSurf(i,j,kloc,bi,bj) +
50d8304171 Ryan*0636      &     Tfac * dzfac * layers_surfflux(i,j,1,1,bi,bj)
cf336ab6c5 Ryan*0637           StendSurf(i,j,kloc,bi,bj) =
06e3ddabf4 Jean*0638      &     StendSurf(i,j,kloc,bi,bj) +
50d8304171 Ryan*0639      &     Sfac * dzfac * layers_surfflux(i,j,1,2,bi,bj)
cf336ab6c5 Ryan*0640          ENDIF
50d8304171 Ryan*0641
06e3ddabf4 Jean*0642 #ifdef SHORTWAVE_HEATING
50d8304171 Ryan*0643           TtendSurf(i,j,kloc,bi,bj) =
                0644      &     TtendSurf(i,j,kloc,bi,bj) +
                0645      &     Tfac * dzfac * layers_sw(i,j,kci,1,bi,bj)
                0646 #endif /* SHORTWAVE_HEATING */
                0647
                0648 C ------- Diffusion
cf336ab6c5 Ryan*0649          TtendDiffh(i,j,kloc,bi,bj) =
50d8304171 Ryan*0650      &     TtendDiffh(i,j,kloc,bi,bj) + dzfac * Tfac *
                0651      &    (layers_dfx(i,j,kci,1,bi,bj)+
                0652      &     layers_dfy(i,j,kci,1,bi,bj))
cf336ab6c5 Ryan*0653          TtendDiffr(i,j,kloc,bi,bj) =
06e3ddabf4 Jean*0654      &     TtendDiffr(i,j,kloc,bi,bj) +
50d8304171 Ryan*0655      &     dzfac * Tfac * layers_dfr(i,j,kci,1,bi,bj)
cf336ab6c5 Ryan*0656          StendDiffh(i,j,kloc,bi,bj) =
50d8304171 Ryan*0657      &     StendDiffh(i,j,kloc,bi,bj) + dzfac * Sfac *
                0658      &    (layers_dfx(i,j,kci,2,bi,bj)+
                0659      &     layers_dfy(i,j,kci,2,bi,bj))
cf336ab6c5 Ryan*0660          StendDiffr(i,j,kloc,bi,bj) =
06e3ddabf4 Jean*0661      &     StendDiffr(i,j,kloc,bi,bj) +
50d8304171 Ryan*0662      &     dzfac * Sfac * layers_dfr(i,j,kci,2,bi,bj)
                0663 C ------- Advection
                0664          TtendAdvh(i,j,kloc,bi,bj) =
                0665      &     TtendAdvh(i,j,kloc,bi,bj) + dzfac * Tfac *
                0666      &    (layers_afx(i,j,kci,1,bi,bj)+
                0667      &     layers_afy(i,j,kci,1,bi,bj))
                0668          TtendAdvr(i,j,kloc,bi,bj) =
06e3ddabf4 Jean*0669      &     TtendAdvr(i,j,kloc,bi,bj) +
50d8304171 Ryan*0670      &     dzfac * Tfac * layers_afr(i,j,kci,1,bi,bj)
                0671          StendAdvh(i,j,kloc,bi,bj) =
                0672      &     StendAdvh(i,j,kloc,bi,bj) + dzfac * Sfac *
                0673      &    (layers_afx(i,j,kci,2,bi,bj)+
                0674      &     layers_afy(i,j,kci,2,bi,bj))
                0675          StendAdvr(i,j,kloc,bi,bj) =
06e3ddabf4 Jean*0676      &     StendAdvr(i,j,kloc,bi,bj) +
50d8304171 Ryan*0677      &     dzfac * Sfac * layers_afr(i,j,kci,2,bi,bj)
                0678 C -------- Total Tendency
06e3ddabf4 Jean*0679          Ttendtot(i,j,kloc,bi,bj) =
                0680      &     Ttendtot(i,j,kloc,bi,bj) +
50d8304171 Ryan*0681      &     dzfac * Tfac * layers_tottend(i,j,kci,1,bi,bj)
06e3ddabf4 Jean*0682          Stendtot(i,j,kloc,bi,bj) =
                0683      &     Stendtot(i,j,kloc,bi,bj) +
50d8304171 Ryan*0684      &     dzfac * Sfac * layers_tottend(i,j,kci,2,bi,bj)
cf336ab6c5 Ryan*0685         ENDDO
                0686        ENDDO
                0687       ENDDO
1377839540 Jean*0688
cf336ab6c5 Ryan*0689 C--   Now that we know the thicknesses, compute the heaviside function
                0690 C--   (Needs another loop through Ng)
                0691       DO kg=1,Nlayers
                0692        DO j = 1,sNy+1
                0693         DO i = 1,sNx+1
                0694          IF (Hc(i,j,kg,bi,bj) .GT. 0.) THEN
                0695           PIc(i,j,kg,bi,bj) = 1. _d 0
                0696          ENDIF
                0697         ENDDO
                0698        ENDDO
                0699       ENDDO
                0700
                0701 C --- End bi,bj loop
                0702       ENDDO
                0703       ENDDO
                0704
                0705 #endif /* LAYERS_THERMODYNAMICS */
                0706
                0707       RETURN
                0708       END
                0709
                0710 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
9479dbe556 Mart*0711
                0712 CBOP
28ff9bcfe6 Jean*0713       SUBROUTINE LAYERS_LOCATE(
                0714      I                          xx,n,m,sNx,sNy,x,
                0715      O                          k,
                0716      I                          myThid )
9479dbe556 Mart*0717
                0718 C     !DESCRIPTION: \bv
                0719 C     *==========================================================*
28ff9bcfe6 Jean*0720 C     | Find the index(-array) k such that x is between xx(k)
                0721 C     | and xx(k+1) by bisection, following Press et al.,
9479dbe556 Mart*0722 C     | Numerical Recipes in Fortran. xx must be monotonic.
                0723 C     *==========================================================*
                0724 C     \ev
                0725
                0726 C !USES:
                0727       IMPLICIT NONE
                0728 C !INPUT PARAMETERS:
28ff9bcfe6 Jean*0729 C     xx        :: array of bin-boundaries (layers_boundaries)
9479dbe556 Mart*0730 C     n         :: length of xx
                0731 C     m         :: int(log2(n)) + 1 = length of bisection loop
28ff9bcfe6 Jean*0732 C     sNx,sNy   :: size of index array and input x
9479dbe556 Mart*0733 C     x         :: input array of values
                0734 C     k         :: index array (output)
28ff9bcfe6 Jean*0735 C     myThid    :: my Thread Id number
                0736       INTEGER n,m,sNx,sNy
9479dbe556 Mart*0737       _RL     xx(1:n+1)
28ff9bcfe6 Jean*0738       _RL     x(snx+1,sny+1)
                0739       INTEGER k(snx+1,sny+1)
                0740       INTEGER myThid
9479dbe556 Mart*0741
                0742 C !LOCAL VARIABLES:
                0743 C     i,j      :: horizontal indices
                0744 C     l        :: bisection loop index
28ff9bcfe6 Jean*0745 C     kl,ku,km :: work arrays and variables
                0746       INTEGER i,j
9479dbe556 Mart*0747 CEOP
                0748 #ifdef TARGET_NEC_SX
28ff9bcfe6 Jean*0749       INTEGER l, km
                0750       INTEGER kl(sNx+1,sNy+1), ku(sNx+1,sNy+1)
                0751
9479dbe556 Mart*0752 C     bisection, following Press et al., Numerical Recipes in Fortran,
                0753 C     mostly, because it can be vectorized
                0754       DO j = 1,sNy+1
                0755        DO i = 1,sNx+1
                0756         kl(i,j)=1
                0757         ku(i,j)=n+1
                0758        END DO
                0759       END DO
                0760       DO l = 1,m
                0761        DO j = 1,sNy+1
                0762         DO i = 1,sNx+1
                0763          IF (ku(i,j)-kl(i,j).GT.1) THEN
                0764           km=(ku(i,j)+kl(i,j))/2
                0765 CML       IF ((xx(n).GE.xx(1)).EQV.(x(i,j).GE.xx(km))) THEN
                0766           IF ( ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))).OR.
                0767      &         ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))) ) THEN
                0768            kl(i,j)=km
                0769           ELSE
                0770            ku(i,j)=km
                0771           END IF
                0772          END IF
                0773         END DO
                0774        END DO
                0775       END DO
                0776       DO j = 1,sNy+1
                0777        DO i = 1,sNx+1
                0778         IF ( x(i,j).LT.xx(2) ) THEN
                0779          k(i,j)=1
                0780         ELSE IF ( x(i,j).GE.xx(n) ) THEN
                0781          k(i,j)=n
                0782         ELSE
                0783          k(i,j)=kl(i,j)
                0784         END IF
                0785        END DO
                0786       END DO
                0787 #else
                0788 C     the old way
                0789       DO j = 1,sNy+1
                0790        DO i = 1,sNx+1
                0791         IF (x(i,j) .GE. xx(n)) THEN
                0792 C     the point is in the hottest bin or hotter
                0793          k(i,j) = n
                0794         ELSE IF (x(i,j) .LT. xx(2)) THEN
                0795 C        the point is in the coldest bin or colder
                0796          k(i,j) = 1
                0797         ELSE IF ( (x(i,j) .GE. xx(k(i,j)))
                0798      &    .AND.   (x(i,j) .LT. xx(k(i,j)+1)) ) THEN
                0799 C     already on the right bin -- do nothing
                0800         ELSE IF (x(i,j) .GE. xx(k(i,j))) THEN
                0801 C     have to hunt for the right bin by getting hotter
                0802          DO WHILE (x(i,j) .GE. xx(k(i,j)+1))
                0803           k(i,j) = k(i,j) + 1
                0804          ENDDO
                0805 C     now xx(k) < x <= xx(k+1)
                0806         ELSE IF (x(i,j) .LT. xx(k(i,j)+1)) THEN
                0807 C     have to hunt for the right bin by getting colder
                0808          DO WHILE (x(i,j) .LT. xx(k(i,j)))
                0809           k(i,j) = k(i,j) - 1
                0810          ENDDO
                0811 C     now xx(k) <= x < xx(k+1)
                0812         ELSE
                0813 C     that should have covered all the options
                0814          k(i,j) = -1
                0815         ENDIF
28ff9bcfe6 Jean*0816
9479dbe556 Mart*0817        ENDDO
                0818       ENDDO
                0819 #endif /* TARGET_NEC_SX */
                0820
aa668e01f1 Gael*0821 #endif /* ALLOW_LAYERS */
                0822
                0823       RETURN
                0824       END
cf336ab6c5 Ryan*0825