pkg/layers/layers_calc.F

C $Header: /u/gcmpack/MITgcm/pkg/layers/layers_calc.F,v 1.8 2011/05/12 15:09:54 dfer Exp $
C $Name:  $

#include "LAYERS_OPTIONS.h"
#ifdef ALLOW_GMREDI
#include "GMREDI_OPTIONS.h"
#endif

CBOP 0
C !ROUTINE: LAYERS_CALC

C !INTERFACE:
      SUBROUTINE LAYERS_CALC(
     I                  myTime, myIter, myThid )

C !DESCRIPTION:
C ===================================================================
C     Calculate the transport in isopycnal layers.
C     This is the meat of the LAYERS package.
C ===================================================================

C !USES:
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "LAYERS_SIZE.h"
#include "LAYERS.h"
#ifdef ALLOW_GMREDI
# include "GMREDI.h"
#endif

C !INPUT PARAMETERS:
C     myTime :: Current time in simulation
C     myIter :: Current iteration number
C     myThid :: my Thread Id number
      _RL     myTime
      INTEGER myIter
      INTEGER myThid
CEOP

#ifdef ALLOW_LAYERS

C !LOCAL VARIABLES:
C     bi, bj   :: tile indices
C     i,j      :: horizontal indices
C     k        :: vertical index for model grid
C     kci      :: index from CellIndex
C     kg       :: index for looping though layers_G
C     kk       :: vertical index for ZZ (fine) grid
C     kgu,kgv  :: vertical index for isopycnal grid
C     TatV     :: temperature at U point
C     TatV     :: temperature at V point

      INTEGER bi, bj
      INTEGER i,j,k,kk,kg,kci,kp1
      INTEGER kgu(sNx+1,sNy+1), kgv(sNx+1,sNy+1)
      _RL TatU, TatV
      CHARACTER*(MAX_LEN_MBUF) msgBuf
#if (defined ALLOW_GMREDI) && (defined GM_BOLUS_ADVEC)
      INTEGER kcip1
      _RL delPsi, maskp1
#endif

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

C --- The tile loops
      DO bj=myByLo(myThid),myByHi(myThid)
      DO bi=myBxLo(myThid),myBxHi(myThid)

C     Initialize the search indices
      DO j = 1,sNy+1
        DO i = 1,sNx+1
C       The temperature index (layer_G) goes from cold to warm.
C       The water column goes from warm (k=1) to cold (k=Nr).
C       So initialize the search with the warmest value.
          kgu(i,j) = Nlayers
          kgv(i,j) = Nlayers
        ENDDO
      ENDDO

C     Reset the arrays
      DO kg=1,Nlayers
       DO j = 1,sNy+1
        DO i = 1,sNx+1
#ifdef LAYERS_UFLUX
         layers_UFlux(i,j,kg,bi,bj) = 0. _d 0
#ifdef LAYERS_THICKNESS
         layers_HU(i,j,kg,bi,bj) = 0. _d 0
#endif /* LAYERS_THICKNESS */
#endif /* LAYERS_UFLUX */
#ifdef LAYERS_VFLUX
         layers_VFlux(i,j,kg,bi,bj) = 0. _d 0
#ifdef LAYERS_THICKNESS
         layers_HV(i,j,kg,bi,bj) = 0. _d 0
#endif /* LAYERS_THICKNESS */
#endif /* LAYERS_VFLUX */
        ENDDO
       ENDDO
      ENDDO

C      _RL  theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
C     Sometimes it is done this way
C      DO j=1-Oly+1,sNy+Oly-1
C       DO i=1-Olx+1,sNx+Olx-1
      DO kk=1,NZZ
       k = MapIndex(kk)
       kci = CellIndex(kk)
       DO j = 1,sNy+1
        DO i = 1,sNx+1

#ifdef LAYERS_UFLUX
C ------ Find theta at the U point (west) on the fine Z grid
         kp1=k+1
         IF (hFacW(i,j,kp1,bi,bj) .EQ. 0.) kp1=k
         IF (LAYER_nb .EQ. 1) THEN
         TatU = MapFact(kk) *
     &    0.5 _d 0 * (theta(i-1,j,k,bi,bj)+theta(i,j,k,bi,bj)) +
     &    (1-MapFact(kk)) *
     &    0.5 _d 0 * (theta(i-1,j,kp1,bi,bj)+theta(i,j,kp1,bi,bj))
         ELSEIF (LAYER_nb .EQ. 2) THEN
         TatU = MapFact(kk) *
     &    0.5 _d 0 * (salt(i-1,j,k,bi,bj)+salt(i,j,k,bi,bj)) +
     &    (1-MapFact(kk)) *
     &    0.5 _d 0 * (salt(i-1,j,kp1,bi,bj)+salt(i,j,kp1,bi,bj))
         ENDIF

C ------ Now that we know T everywhere, determine the binning.

         IF (TatU .GE. layers_G(Nlayers)) THEN
C        the point is in the hottest bin or hotter
          kgu(i,j) = Nlayers
         ELSE IF (TatU .LT. layers_G(2)) THEN
C        the point is in the coldest bin or colder
          kgu(i,j) = 1
         ELSE IF ( (TatU .GE. layers_G(kgu(i,j)))
     &    .AND. (TatU .LT. layers_G(kgu(i,j)+1)) ) THEN
C        already on the right bin -- do nothing
         ELSE IF (TatU .GE. layers_G(kgu(i,j))) THEN
C        have to hunt for the right bin by getting hotter
          DO WHILE (TatU .GE. layers_G(kgu(i,j)+1))
           kgu(i,j) = kgu(i,j) + 1
          ENDDO
C         now layers_G(kgu(i,j)+1) < TatU <= layers_G(kgu(i,j)+1)
         ELSE IF (TatU .LT. layers_G(kgu(i,j)+1)) THEN
C        have to hunt for the right bin by getting colder
          DO WHILE (TatU .LT. layers_G(kgu(i,j)))
           kgu(i,j) = kgu(i,j) - 1
          ENDDO
C         now layers_G(kgu(i,j)+1) <= TatU < layers_G(kgu(i,j)+1)
         ELSE
C        that should have covered all the options
          WRITE(msgBuf,'(A,1E14.6)')
     &     'S/R LAYERS_CALC: Couldnt find a bin in layers_G for TatU=',
     &     TatU
          CALL PRINT_ERROR( msgBuf, myThid )
          STOP 'ABNORMAL END: S/R LAYERS_INIT_FIXED'
         END IF

C ------ Augment the bin values
         layers_UFlux(i,j,kgu(i,j),bi,bj) =
     &    layers_UFlux(i,j,kgu(i,j),bi,bj) +
     &    dZZf(kk) * uVel(i,j,kci,bi,bj) * hFacW(i,j,kci,bi,bj)

#if (defined ALLOW_GMREDI) && (defined GM_BOLUS_ADVEC)
         IF ( GM_AdvForm .AND. useBOLUS ) THEN
           kcip1 = MIN(kci+1,Nr)
           maskp1 = 1.
           IF (kci.GE.Nr) maskp1 = 0.
           delPsi = GM_PsiX(i,j,kcip1,bi,bj)*maskp1
     &            - GM_PsiX(i,j, kci, bi,bj)
           layers_UFlux(i,j,kgu(i,j),bi,bj) =
     &      layers_UFlux(i,j,kgu(i,j),bi,bj)
     &      + delPsi*recip_drF(kci)*_recip_hFacW(i,j,kci,bi,bj)
     &      * dZZf(kk)*hFacW(i,j,kci,bi,bj)
         ENDIF
#endif

#ifdef LAYERS_THICKNESS
         layers_HU(i,j,kgu(i,j),bi,bj) = layers_HU(i,j,kgu(i,j),bi,bj)
     &    + dZZf(kk) * hFacW(i,j,kci,bi,bj)
#endif /* LAYERS_THICKNESS */

#endif /* LAYERS_UFLUX */

#ifdef LAYERS_VFLUX
C ------ Find theta at the V point (south) on the fine Z grid
         kp1=k+1
         IF (hFacS(i,j,kp1,bi,bj) .EQ. 0.) kp1=k
         IF (LAYER_nb .EQ. 1) THEN
         TatV = MapFact(kk) *
     &    0.5 _d 0 * (theta(i,j-1,k,bi,bj)+theta(i,j,k,bi,bj)) +
     &    (1-MapFact(kk)) *
     &    0.5 _d 0 * (theta(i,j-1,kp1,bi,bj)+theta(i,j,kp1,bi,bj))
         ELSEIF (LAYER_nb .EQ. 2) THEN
         TatV = MapFact(kk) *
     &    0.5 _d 0 * (salt(i,j-1,k,bi,bj)+salt(i,j,k,bi,bj)) +
     &    (1-MapFact(kk)) *
     &    0.5 _d 0 * (salt(i,j-1,kp1,bi,bj)+salt(i,j,kp1,bi,bj))
         ENDIF

C ------ Now that we know T everywhere, determine the binning
         IF (TatV .GE. layers_G(Nlayers)) THEN
C         the point is in the hottest bin or hotter
          kgv(i,j) = Nlayers
         ELSE IF (TatV .LT. layers_G(2)) THEN
C         the point is in the coldest bin or colder
          kgv(i,j) = 1
         ELSE IF ( (TatV .GE. layers_G(kgv(i,j)))
     &    .AND. (TatV .LT. layers_G(kgv(i,j)+1)) ) THEN
C         already on the right bin -- do nothing
         ELSE IF (TatV .GE. layers_G(kgv(i,j))) THEN
C         have to hunt for the right bin by getting hotter
          DO WHILE (TatV .GE. layers_G(kgv(i,j)+1))
           kgv(i,j) = kgv(i,j) + 1
          ENDDO
C         now layers_G(kgv(i,j)+1) < TatV <= layers_G(kgv(i,j)+1)
         ELSE IF (TatV .LT. layers_G(kgv(i,j)+1)) THEN
C         have to hunt for the right bin by getting colder
          DO WHILE (TatV .LT. layers_G(kgv(i,j)))
           kgv(i,j) = kgv(i,j) - 1
          ENDDO
C         now layers_G(kgv(i,j)+1) <= TatV < layers_G(kgv(i,j)+1)
         ELSE
C         that should have covered all the options
          WRITE(msgBuf,'(A,1E14.6)')
     &     'S/R LAYERS_CALC: Couldnt find a bin in layers_G for TatV=',
     &     TatV
          CALL PRINT_ERROR( msgBuf, myThid )
          STOP 'ABNORMAL END: S/R LAYERS_INIT_FIXED'
         END IF

C ------ Augment the bin values
         layers_VFlux(i,j,kgv(i,j),bi,bj) =
     &    layers_VFlux(i,j,kgv(i,j),bi,bj)
     &    + dZZf(kk) * vVel(i,j,kci,bi,bj) * hFacS(i,j,kci,bi,bj)

#if (defined ALLOW_GMREDI) && (defined GM_BOLUS_ADVEC)
         IF ( GM_AdvForm .AND. useBOLUS ) THEN
           kcip1 = MIN(kci+1,Nr)
           maskp1 = 1.
           IF (kci.GE.Nr) maskp1 = 0.
           delPsi = GM_PsiY(i,j,kcip1,bi,bj)*maskp1
     &            - GM_PsiY(i,j, kci, bi,bj)
           layers_VFlux(i,j,kgv(i,j),bi,bj) =
     &      layers_VFlux(i,j,kgv(i,j),bi,bj)
     &      + delPsi*recip_drF(kci)*_recip_hFacS(i,j,kci,bi,bj)
     &      * dZZf(kk)*hFacS(i,j,kci,bi,bj)
         ENDIF
#endif

#ifdef LAYERS_THICKNESS
         layers_HV(i,j,kgv(i,j),bi,bj) = layers_HV(i,j,kgv(i,j),bi,bj)
     &    + dZZf(kk) * hFacS(i,j,kci,bi,bj)
#endif /* LAYERS_THICKNESS */

#endif /* LAYERS_VFLUX */

        ENDDO
       ENDDO
      ENDDO

#ifdef ALLOW_TIMEAVE
C--   Time-average
      IF ( layers_taveFreq.GT.0. ) THEN

#ifdef LAYERS_UFLUX
         CALL TIMEAVE_CUMULATE( layers_UFlux_T, layers_UFlux, Nlayers,
     &                          deltaTclock, bi, bj, myThid )
#ifdef LAYERS_THICKNESS
         CALL TIMEAVE_CUMULATE( layers_HU_T, layers_HU, Nlayers,
     &                          deltaTclock, bi, bj, myThid )
#endif /* LAYERS_THICKNESS */
#endif /* LAYERS_UFLUX */
#ifdef LAYERS_VFLUX
         CALL TIMEAVE_CUMULATE( layers_VFlux_T, layers_VFlux, Nlayers,
     &                          deltaTclock, bi, bj, myThid )
#ifdef LAYERS_THICKNESS
         CALL TIMEAVE_CUMULATE( layers_HV_T, layers_HV, Nlayers,
     &                          deltaTclock, bi, bj, myThid )
#endif /* LAYERS_THICKNESS */
#endif /* LAYERS_VFLUX */

         layers_TimeAve(bi,bj)=layers_TimeAve(bi,bj)+deltaTclock

      ENDIF
#endif /* ALLOW_TIMEAVE */

C --- End bi,bj loop
      ENDDO
      ENDDO

#endif /* ALLOW_LAYERS */

      RETURN
      END
1	dfer	1.9	C $Header: /u/gcmpack/MITgcm/pkg/layers/layers_calc.F,v 1.8 2011/05/12 15:09:54 dfer Exp $
2	rpa	1.1	C $Name: $
3
4			#include "LAYERS_OPTIONS.h"
5	jmc	1.6	#ifdef ALLOW_GMREDI
6	dfer	1.5	#include "GMREDI_OPTIONS.h"
7	jmc	1.6	#endif
8	rpa	1.1
9	jmc	1.4	CBOP 0
10			C !ROUTINE: LAYERS_CALC
11	rpa	1.1
12	jmc	1.4	C !INTERFACE:
13	rpa	1.1	SUBROUTINE LAYERS_CALC(
14	jmc	1.4	I myTime, myIter, myThid )
15	rpa	1.1
16	jmc	1.4	C !DESCRIPTION:
17	rpa	1.1	C ===================================================================
18			C Calculate the transport in isopycnal layers.
19			C This is the meat of the LAYERS package.
20			C ===================================================================
21
22	jmc	1.4	C !USES:
23	rpa	1.1	IMPLICIT NONE
24			#include "SIZE.h"
25	jmc	1.4	#include "EEPARAMS.h"
26			#include "PARAMS.h"
27	rpa	1.1	#include "GRID.h"
28			#include "DYNVARS.h"
29			#include "LAYERS_SIZE.h"
30			#include "LAYERS.h"
31	jmc	1.6	#ifdef ALLOW_GMREDI
32			# include "GMREDI.h"
33	dfer	1.5	#endif
34	rpa	1.1
35	jmc	1.4	C !INPUT PARAMETERS:
36			C myTime :: Current time in simulation
37			C myIter :: Current iteration number
38			C myThid :: my Thread Id number
39			_RL myTime
40			INTEGER myIter
41			INTEGER myThid
42			CEOP
43	rpa	1.1
44			#ifdef ALLOW_LAYERS
45
46	jmc	1.4	C !LOCAL VARIABLES:
47			C bi, bj :: tile indices
48	rpa	1.1	C i,j :: horizontal indices
49			C k :: vertical index for model grid
50			C kci :: index from CellIndex
51			C kg :: index for looping though layers_G
52			C kk :: vertical index for ZZ (fine) grid
53			C kgu,kgv :: vertical index for isopycnal grid
54			C TatV :: temperature at U point
55			C TatV :: temperature at V point
56
57	jmc	1.4	INTEGER bi, bj
58	dfer	1.9	INTEGER i,j,k,kk,kg,kci,kp1
59	rpa	1.1	INTEGER kgu(sNx+1,sNy+1), kgv(sNx+1,sNy+1)
60	jmc	1.6	_RL TatU, TatV
61	rpa	1.1	CHARACTER*(MAX_LEN_MBUF) msgBuf
62	jmc	1.6	#if (defined ALLOW_GMREDI) && (defined GM_BOLUS_ADVEC)
63			INTEGER kcip1
64			_RL delPsi, maskp1
65			#endif
66	rpa	1.1
67	jmc	1.4	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
68
69			C --- The tile loops
70	rpa	1.1	DO bj=myByLo(myThid),myByHi(myThid)
71	jmc	1.4	DO bi=myBxLo(myThid),myBxHi(myThid)
72	rpa	1.1
73	dfer	1.5	C Initialize the search indices
74	rpa	1.1	DO j = 1,sNy+1
75			DO i = 1,sNx+1
76			C The temperature index (layer_G) goes from cold to warm.
77			C The water column goes from warm (k=1) to cold (k=Nr).
78			C So initialize the search with the warmest value.
79			kgu(i,j) = Nlayers
80			kgv(i,j) = Nlayers
81			ENDDO
82			ENDDO
83
84			C Reset the arrays
85			DO kg=1,Nlayers
86			DO j = 1,sNy+1
87			DO i = 1,sNx+1
88			#ifdef LAYERS_UFLUX
89			layers_UFlux(i,j,kg,bi,bj) = 0. _d 0
90			#ifdef LAYERS_THICKNESS
91			layers_HU(i,j,kg,bi,bj) = 0. _d 0
92			#endif /* LAYERS_THICKNESS */
93			#endif /* LAYERS_UFLUX */
94			#ifdef LAYERS_VFLUX
95			layers_VFlux(i,j,kg,bi,bj) = 0. _d 0
96			#ifdef LAYERS_THICKNESS
97			layers_HV(i,j,kg,bi,bj) = 0. _d 0
98			#endif /* LAYERS_THICKNESS */
99			#endif /* LAYERS_VFLUX */
100			ENDDO
101			ENDDO
102			ENDDO
103
104			C _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
105			C Sometimes it is done this way
106			C DO j=1-Oly+1,sNy+Oly-1
107			C DO i=1-Olx+1,sNx+Olx-1
108			DO kk=1,NZZ
109			k = MapIndex(kk)
110			kci = CellIndex(kk)
111			DO j = 1,sNy+1
112			DO i = 1,sNx+1
113
114			#ifdef LAYERS_UFLUX
115			C ------ Find theta at the U point (west) on the fine Z grid
116	dfer	1.9	kp1=k+1
117			IF (hFacW(i,j,kp1,bi,bj) .EQ. 0.) kp1=k
118	dfer	1.5	IF (LAYER_nb .EQ. 1) THEN
119	rpa	1.1	TatU = MapFact(kk) *
120			& 0.5 _d 0 * (theta(i-1,j,k,bi,bj)+theta(i,j,k,bi,bj)) +
121			& (1-MapFact(kk)) *
122	dfer	1.9	& 0.5 _d 0 * (theta(i-1,j,kp1,bi,bj)+theta(i,j,kp1,bi,bj))
123	dfer	1.5	ELSEIF (LAYER_nb .EQ. 2) THEN
124			TatU = MapFact(kk) *
125			& 0.5 _d 0 * (salt(i-1,j,k,bi,bj)+salt(i,j,k,bi,bj)) +
126			& (1-MapFact(kk)) *
127	dfer	1.9	& 0.5 _d 0 * (salt(i-1,j,kp1,bi,bj)+salt(i,j,kp1,bi,bj))
128	dfer	1.5	ENDIF
129	rpa	1.1
130			C ------ Now that we know T everywhere, determine the binning.
131
132			IF (TatU .GE. layers_G(Nlayers)) THEN
133			C the point is in the hottest bin or hotter
134			kgu(i,j) = Nlayers
135			ELSE IF (TatU .LT. layers_G(2)) THEN
136			C the point is in the coldest bin or colder
137			kgu(i,j) = 1
138			ELSE IF ( (TatU .GE. layers_G(kgu(i,j)))
139			& .AND. (TatU .LT. layers_G(kgu(i,j)+1)) ) THEN
140			C already on the right bin -- do nothing
141			ELSE IF (TatU .GE. layers_G(kgu(i,j))) THEN
142			C have to hunt for the right bin by getting hotter
143			DO WHILE (TatU .GE. layers_G(kgu(i,j)+1))
144			kgu(i,j) = kgu(i,j) + 1
145			ENDDO
146			C now layers_G(kgu(i,j)+1) < TatU <= layers_G(kgu(i,j)+1)
147			ELSE IF (TatU .LT. layers_G(kgu(i,j)+1)) THEN
148			C have to hunt for the right bin by getting colder
149			DO WHILE (TatU .LT. layers_G(kgu(i,j)))
150			kgu(i,j) = kgu(i,j) - 1
151			ENDDO
152			C now layers_G(kgu(i,j)+1) <= TatU < layers_G(kgu(i,j)+1)
153			ELSE
154			C that should have covered all the options
155			WRITE(msgBuf,'(A,1E14.6)')
156			& 'S/R LAYERS_CALC: Couldnt find a bin in layers_G for TatU=',
157			& TatU
158			CALL PRINT_ERROR( msgBuf, myThid )
159			STOP 'ABNORMAL END: S/R LAYERS_INIT_FIXED'
160			END IF
161
162			C ------ Augment the bin values
163			layers_UFlux(i,j,kgu(i,j),bi,bj) =
164			& layers_UFlux(i,j,kgu(i,j),bi,bj) +
165	rpa	1.2	& dZZf(kk) * uVel(i,j,kci,bi,bj) * hFacW(i,j,kci,bi,bj)
166	rpa	1.1
167	dfer	1.5	#if (defined ALLOW_GMREDI) && (defined GM_BOLUS_ADVEC)
168	dfer	1.7	IF ( GM_AdvForm .AND. useBOLUS ) THEN
169	dfer	1.5	kcip1 = MIN(kci+1,Nr)
170			maskp1 = 1.
171			IF (kci.GE.Nr) maskp1 = 0.
172			delPsi = GM_PsiX(i,j,kcip1,bi,bj)*maskp1
173			& - GM_PsiX(i,j, kci, bi,bj)
174			layers_UFlux(i,j,kgu(i,j),bi,bj) =
175	jmc	1.6	& layers_UFlux(i,j,kgu(i,j),bi,bj)
176	dfer	1.5	& + delPsirecip_drF(kci)_recip_hFacW(i,j,kci,bi,bj)
177			& * dZZf(kk)*hFacW(i,j,kci,bi,bj)
178			ENDIF
179			#endif
180
181	rpa	1.1	#ifdef LAYERS_THICKNESS
182			layers_HU(i,j,kgu(i,j),bi,bj) = layers_HU(i,j,kgu(i,j),bi,bj)
183	rpa	1.2	& + dZZf(kk) * hFacW(i,j,kci,bi,bj)
184	rpa	1.1	#endif /* LAYERS_THICKNESS */
185
186			#endif /* LAYERS_UFLUX */
187
188			#ifdef LAYERS_VFLUX
189			C ------ Find theta at the V point (south) on the fine Z grid
190	dfer	1.9	kp1=k+1
191			IF (hFacS(i,j,kp1,bi,bj) .EQ. 0.) kp1=k
192	dfer	1.5	IF (LAYER_nb .EQ. 1) THEN
193	rpa	1.1	TatV = MapFact(kk) *
194			& 0.5 _d 0 * (theta(i,j-1,k,bi,bj)+theta(i,j,k,bi,bj)) +
195			& (1-MapFact(kk)) *
196	dfer	1.9	& 0.5 _d 0 * (theta(i,j-1,kp1,bi,bj)+theta(i,j,kp1,bi,bj))
197	dfer	1.5	ELSEIF (LAYER_nb .EQ. 2) THEN
198			TatV = MapFact(kk) *
199			& 0.5 _d 0 * (salt(i,j-1,k,bi,bj)+salt(i,j,k,bi,bj)) +
200			& (1-MapFact(kk)) *
201	dfer	1.9	& 0.5 _d 0 * (salt(i,j-1,kp1,bi,bj)+salt(i,j,kp1,bi,bj))
202	dfer	1.5	ENDIF
203	rpa	1.1
204			C ------ Now that we know T everywhere, determine the binning
205			IF (TatV .GE. layers_G(Nlayers)) THEN
206			C the point is in the hottest bin or hotter
207			kgv(i,j) = Nlayers
208			ELSE IF (TatV .LT. layers_G(2)) THEN
209			C the point is in the coldest bin or colder
210			kgv(i,j) = 1
211			ELSE IF ( (TatV .GE. layers_G(kgv(i,j)))
212			& .AND. (TatV .LT. layers_G(kgv(i,j)+1)) ) THEN
213			C already on the right bin -- do nothing
214			ELSE IF (TatV .GE. layers_G(kgv(i,j))) THEN
215			C have to hunt for the right bin by getting hotter
216			DO WHILE (TatV .GE. layers_G(kgv(i,j)+1))
217			kgv(i,j) = kgv(i,j) + 1
218			ENDDO
219			C now layers_G(kgv(i,j)+1) < TatV <= layers_G(kgv(i,j)+1)
220			ELSE IF (TatV .LT. layers_G(kgv(i,j)+1)) THEN
221			C have to hunt for the right bin by getting colder
222			DO WHILE (TatV .LT. layers_G(kgv(i,j)))
223			kgv(i,j) = kgv(i,j) - 1
224			ENDDO
225			C now layers_G(kgv(i,j)+1) <= TatV < layers_G(kgv(i,j)+1)
226			ELSE
227			C that should have covered all the options
228			WRITE(msgBuf,'(A,1E14.6)')
229			& 'S/R LAYERS_CALC: Couldnt find a bin in layers_G for TatV=',
230			& TatV
231			CALL PRINT_ERROR( msgBuf, myThid )
232			STOP 'ABNORMAL END: S/R LAYERS_INIT_FIXED'
233			END IF
234
235			C ------ Augment the bin values
236			layers_VFlux(i,j,kgv(i,j),bi,bj) =
237			& layers_VFlux(i,j,kgv(i,j),bi,bj)
238	rpa	1.2	& + dZZf(kk) * vVel(i,j,kci,bi,bj) * hFacS(i,j,kci,bi,bj)
239	rpa	1.1
240	dfer	1.5	#if (defined ALLOW_GMREDI) && (defined GM_BOLUS_ADVEC)
241	dfer	1.7	IF ( GM_AdvForm .AND. useBOLUS ) THEN
242	dfer	1.5	kcip1 = MIN(kci+1,Nr)
243			maskp1 = 1.
244			IF (kci.GE.Nr) maskp1 = 0.
245	dfer	1.8	delPsi = GM_PsiY(i,j,kcip1,bi,bj)*maskp1
246			& - GM_PsiY(i,j, kci, bi,bj)
247			layers_VFlux(i,j,kgv(i,j),bi,bj) =
248			& layers_VFlux(i,j,kgv(i,j),bi,bj)
249			& + delPsirecip_drF(kci)_recip_hFacS(i,j,kci,bi,bj)
250			& * dZZf(kk)*hFacS(i,j,kci,bi,bj)
251	dfer	1.5	ENDIF
252			#endif
253
254	rpa	1.1	#ifdef LAYERS_THICKNESS
255			layers_HV(i,j,kgv(i,j),bi,bj) = layers_HV(i,j,kgv(i,j),bi,bj)
256	rpa	1.2	& + dZZf(kk) * hFacS(i,j,kci,bi,bj)
257	rpa	1.1	#endif /* LAYERS_THICKNESS */
258
259			#endif /* LAYERS_VFLUX */
260	jmc	1.4
261	rpa	1.1	ENDDO
262			ENDDO
263			ENDDO
264
265			#ifdef ALLOW_TIMEAVE
266			C-- Time-average
267	dfer	1.3	IF ( layers_taveFreq.GT.0. ) THEN
268	rpa	1.1
269			#ifdef LAYERS_UFLUX
270			CALL TIMEAVE_CUMULATE( layers_UFlux_T, layers_UFlux, Nlayers,
271			& deltaTclock, bi, bj, myThid )
272			#ifdef LAYERS_THICKNESS
273			CALL TIMEAVE_CUMULATE( layers_HU_T, layers_HU, Nlayers,
274			& deltaTclock, bi, bj, myThid )
275			#endif /* LAYERS_THICKNESS */
276			#endif /* LAYERS_UFLUX */
277			#ifdef LAYERS_VFLUX
278			CALL TIMEAVE_CUMULATE( layers_VFlux_T, layers_VFlux, Nlayers,
279			& deltaTclock, bi, bj, myThid )
280			#ifdef LAYERS_THICKNESS
281			CALL TIMEAVE_CUMULATE( layers_HV_T, layers_HV, Nlayers,
282			& deltaTclock, bi, bj, myThid )
283			#endif /* LAYERS_THICKNESS */
284			#endif /* LAYERS_VFLUX */
285
286	jmc	1.4	layers_TimeAve(bi,bj)=layers_TimeAve(bi,bj)+deltaTclock
287	rpa	1.1
288			ENDIF
289			#endif /* ALLOW_TIMEAVE */
290
291			C --- End bi,bj loop
292			ENDDO
293			ENDDO
294
295			#endif /* ALLOW_LAYERS */
296
297			RETURN
298			END