C $Header: /u/gcmpack/MITgcm/pkg/layers/layers_fluxcalc.F,v 1.17 2015/06/15 21:34:07 jmc Exp $
C $Name: $
#include "LAYERS_OPTIONS.h"
#ifdef ALLOW_GMREDI
#include "GMREDI_OPTIONS.h"
#endif
C-- File layers_fluxcalc.F:
C-- Contents
C-- o LAYERS_FLUXCALC
C-- o LAYERS_DIAPYCNAL
C-- o LAYERS_LOCATE
CBOP 0
C !ROUTINE: LAYERS_FLUXCALC
C !INTERFACE:
SUBROUTINE LAYERS_FLUXCALC(
I uVel,vVel,tracer,iLa,
O UH,VH,Hw,Hs,PIw,PIs,Uw,Vs,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE LAYERS_FLUXCALC
C | Calculate the transport in isotracer layers, for a chosen
C | tracer. This is the meat of the LAYERS package.
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "LAYERS_SIZE.h"
#include "LAYERS.h"
#ifdef ALLOW_GMREDI
# include "GMREDI.h"
#endif
C !INPUT PARAMETERS:
C myThid :: my Thread Id number
C uVel :: zonal velocity (m/s, i=1 held at western face)
C vVel :: meridional velocity (m/s, j=1 held at southern face)
C tracer :: potential temperature, salt or potential density prho
C UH :: U integrated over layer (m^2/s)
C VH :: V integrated over layer (m^2/s)
C Hw :: Layer thickness at the U point (m)
C Hs :: Layer thickness at the V point (m)
C PIw :: 1 if layer exists, 0 otherwise (at U point)
C PIs :: 1 if layer exists, 0 otherwise (at V point)
C Uw :: average U over layer (m/s)
C Vs :: average V over layer (m/s)
C iLa :: layer coordinate index
INTEGER myThid
_RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr, nSx,nSy)
_RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr, nSx,nSy)
_RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr, nSx,nSy)
#ifdef LAYERS_UFLUX
_RL UH (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
# ifdef LAYERS_THICKNESS
_RL Hw (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
_RL PIw (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
_RL Uw (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
else
_RL Hw(1), PIw(1), Uw(1)
# endif
#else
_RL UH(1), Hw(1), PIw(1), Uw(1)
#endif
#ifdef LAYERS_VFLUX
_RL VH (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
# ifdef LAYERS_THICKNESS
_RL Hs (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
_RL PIs (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
_RL Vs (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
else
_RL Hs(1), PIs(1), Vs(1)
# endif
#else
_RL VH(1), Hs(1), PIs(1), Vs(1)
#endif
INTEGER iLa
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_bounds
C kk :: vertical index for ZZ (fine) grid
C kgu,kgv :: vertical index for isopycnal grid
C kloc :: local copy of kgu/v to reduce accesses to index arrays
C mSteps :: maximum number of steps for bisection method
C prho :: pot. density (less 1000) referenced to layers_krho pressure
C TatU :: temperature at U point
C TatV :: temperature at V point
C dzfac :: temporary sublayer thickness
C Tloc,Tp1 :: horizontally interpolated tracer values
INTEGER bi, bj
INTEGER i,j,k,kk,kg,kci,kp1,kloc
INTEGER mSteps
INTEGER kgu(sNx+1,sNy+1), kgv(sNx+1,sNy+1)
_RL TatU(sNx+1,sNy+1), TatV(sNx+1,sNy+1)
_RL dzfac
#ifdef ALLOW_GMREDI
INTEGER kcip1
_RL delPsi, maskp1
#endif
LOGICAL errorFlag
CHARACTER*(MAX_LEN_MBUF) msgBuf
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C compute maximum number of steps for bisection method (approx.
C log2(Nlayers)) as log2(Nlayers) + 1 for safety
mSteps = int(log10(dble(Nlayers))/log10(2. _d 0))+1
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-OLy,sNy+OLy
DO i = 1-OLx,sNx+OLx
#ifdef LAYERS_UFLUX
UH (i,j,kg,bi,bj) = 0. _d 0
#ifdef LAYERS_THICKNESS
Hw(i,j,kg,bi,bj) = 0. _d 0
PIw(i,j,kg,bi,bj) = 0. _d 0
Uw(i,j,kg,bi,bj) = 0. _d 0
#endif /* LAYERS_THICKNESS */
#endif /* UH */
#ifdef LAYERS_VFLUX
VH (i,j,kg,bi,bj) = 0. _d 0
#ifdef LAYERS_THICKNESS
Hs(i,j,kg,bi,bj) = 0. _d 0
PIs(i,j,kg,bi,bj) = 0. _d 0
Vs(i,j,kg,bi,bj) = 0. _d 0
#endif /* LAYERS_THICKNESS */
#endif /* VH */
ENDDO
ENDDO
ENDDO
DO kk=1,NZZ
k = MapIndex(kk)
kci = CellIndex(kk)
#ifdef ALLOW_GMREDI
kcip1 = MIN(kci+1,Nr)
maskp1 = 1.
IF (kci.GE.Nr) maskp1 = 0.
#endif /* ALLOW_GMREDI */
#ifdef LAYERS_UFLUX
DO j = 1,sNy+1
DO i = 1,sNx+1
C ------ Find theta at the U point (west) on the fine Z grid
kp1=k+1
IF (maskW(i,j,kp1,bi,bj).EQ.zeroRS) kp1=k
TatU(i,j) = MapFact(kk) *
& 0.5 _d 0 * (tracer(i-1,j,k,bi,bj)+tracer(i,j,k,bi,bj)) +
& (1. _d 0 -MapFact(kk)) *
& 0.5 _d 0 * (tracer(i-1,j,kp1,bi,bj)+tracer(i,j,kp1,bi,bj))
ENDDO
ENDDO
C ------ Now that we know T everywhere, determine the binning.
C find the layer indices kgu
CALL LAYERS_LOCATE(
I layers_bounds(1,iLa),Nlayers,mSteps,sNx,sNy,TatU,
O kgu,
I myThid )
#ifndef TARGET_NEC_SX
C check for failures
IF ( debugLevel .GE. debLevC ) THEN
errorFlag = .FALSE.
DO j = 1,sNy+1
DO i = 1,sNx+1
IF ( kgu(i,j) .LE. 0 ) THEN
WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)')
& 'S/R LAYERS_LOCATE: Could not find a bin in ',
& 'layers_bounds for TatU(',i,',',j,',)=',TatU(i,j)
CALL PRINT_ERROR( msgBuf, myThid )
errorFlag = .TRUE.
ENDIF
ENDDO
ENDDO
IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_FLUXCALC'
ENDIF
#endif /* ndef TARGET_NEC_SX */
C
DO j = 1,sNy+1
DO i = 1,sNx+1
kloc = kgu(i,j)
dzfac = dZZf(kk) * hFacW(i,j,kci,bi,bj)
C ------ Augment the bin values
UH(i,j,kloc,bi,bj) =
& UH(i,j,kloc,bi,bj) +
& dzfac * uVel(i,j,kci,bi,bj)
#ifdef ALLOW_GMREDI
IF ( layers_bolus(iLa) ) THEN
IF ( .NOT.GM_AdvForm ) THEN
delPsi = 0.25 _d 0 *(
& ( rA(i-1,j,bi,bj)*Kwx(i-1,j,kcip1,bi,bj)
& +rA( i ,j,bi,bj)*Kwx( i ,j,kcip1,bi,bj)
& ) * maskW(i,j,kcip1,bi,bj) * maskp1
& - ( rA(i-1,j,bi,bj)*Kwx(i-1,j, kci ,bi,bj)
& +rA( i ,j,bi,bj)*Kwx( i ,j, kci ,bi,bj)
& ) * maskW(i,j, kci ,bi,bj)
& ) * recip_rAw(i,j,bi,bj)
#ifdef GM_BOLUS_ADVEC
ELSE
delPsi = GM_PsiX(i,j,kcip1,bi,bj)*maskp1
& - GM_PsiX(i,j, kci, bi,bj)
#endif
ENDIF
UH(i,j,kloc,bi,bj) = UH(i,j,kloc,bi,bj)
& + delPsi*recip_drF(kci)*_recip_hFacW(i,j,kci,bi,bj)
& * dzfac
ENDIF
#endif /* ALLOW_GMREDI */
#ifdef LAYERS_THICKNESS
Hw(i,j,kloc,bi,bj) = Hw(i,j,kloc,bi,bj) + dzfac
#endif /* LAYERS_THICKNESS */
ENDDO
ENDDO
#endif /* LAYERS_UFLUX */
#ifdef LAYERS_VFLUX
DO j = 1,sNy+1
DO i = 1,sNx+1
C ------ Find theta at the V point (south) on the fine Z grid
kp1=k+1
IF (maskS(i,j,kp1,bi,bj).EQ.zeroRS) kp1=k
TatV(i,j) = MapFact(kk) *
& 0.5 _d 0 * (tracer(i,j-1,k,bi,bj)+tracer(i,j,k,bi,bj)) +
& (1. _d 0 -MapFact(kk)) *
& 0.5 _d 0 * (tracer(i,j-1,kp1,bi,bj)+tracer(i,j,kp1,bi,bj))
ENDDO
ENDDO
C ------ Now that we know T everywhere, determine the binning.
C find the layer indices kgv
CALL LAYERS_LOCATE(
I layers_bounds(1,iLa),Nlayers,mSteps,sNx,sNy,TatV,
O kgv,
I myThid )
#ifndef TARGET_NEC_SX
IF ( debugLevel .GE. debLevC ) THEN
C check for failures
errorFlag = .FALSE.
DO j = 1,sNy+1
DO i = 1,sNx+1
IF ( kgv(i,j) .LE. 0 ) THEN
WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)')
& 'S/R LAYERS_LOCATE: Could not find a bin in ',
& 'layers_bounds for TatV(',i,',',j,',)=',TatV(i,j)
CALL PRINT_ERROR( msgBuf, myThid )
errorFlag = .TRUE.
ENDIF
ENDDO
ENDDO
IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_FLUXCALC'
ENDIF
#endif /* ndef TARGET_NEC_SX */
C
DO j = 1,sNy+1
DO i = 1,sNx+1
kloc = kgv(i,j)
dzfac = dZZf(kk) * hFacS(i,j,kci,bi,bj)
C ------ debugging stuff
C IF (i.EQ.10 .AND. j.EQ.10) THEN
C WRITE(msgBuf,'(A,I3,A,F10.2,A,F6.2)')
C & ' kloc=', kloc,
C & ', TatV=',TatV(i,j),
C & ', dzfac=',dzfac
C CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
C & SQUEEZE_RIGHT, myThid )
C ENDIF
C ------ Augment the bin values
VH(i,j,kloc,bi,bj) =
& VH(i,j,kloc,bi,bj) + dzfac * vVel(i,j,kci,bi,bj)
#ifdef ALLOW_GMREDI
IF ( layers_bolus(iLa) ) THEN
IF ( .NOT.GM_AdvForm ) THEN
delPsi = 0.25 _d 0 *(
& ( rA(i,j-1,bi,bj)*Kwy(i,j-1,kcip1,bi,bj)
& +rA(i, j ,bi,bj)*Kwy(i, j ,kcip1,bi,bj)
& ) * maskS(i,j,kcip1,bi,bj) * maskp1
& - ( rA(i,j-1,bi,bj)*Kwy(i,j-1, kci ,bi,bj)
& +rA(i, j ,bi,bj)*Kwy(i, j , kci ,bi,bj)
& ) * maskS(i,j, kci ,bi,bj)
& ) * recip_rAs(i,j,bi,bj)
#ifdef GM_BOLUS_ADVEC
ELSE
delPsi = GM_PsiY(i,j,kcip1,bi,bj)*maskp1
& - GM_PsiY(i,j, kci, bi,bj)
#endif
ENDIF
VH(i,j,kloc,bi,bj) = VH(i,j,kloc,bi,bj)
& + delPsi*recip_drF(kci)*_recip_hFacS(i,j,kci,bi,bj)
& * dzfac
ENDIF
#endif /* ALLOW_GMREDI */
#ifdef LAYERS_THICKNESS
Hs(i,j,kloc,bi,bj) = Hs(i,j,kloc,bi,bj) + dzfac
#endif /* LAYERS_THICKNESS */
ENDDO
ENDDO
#endif /* LAYERS_VFLUX */
ENDDO
C-- Now that we know the thicknesses, compute the heaviside function
C-- (Needs another loop through Ng)
#ifdef LAYERS_THICKNESS
DO kg=1,Nlayers
DO j = 1,sNy+1
DO i = 1,sNx+1
#ifdef LAYERS_UFLUX
IF (Hw(i,j,kg,bi,bj) .GT. 0.) THEN
PIw(i,j,kg,bi,bj) = 1. _d 0
Uw(i,j,kg,bi,bj) =
& UH(i,j,kg,bi,bj) / Hw(i,j,kg,bi,bj)
ENDIF
#endif /* LAYERS_UFLUX */
#ifdef LAYERS_VFLUX
IF (Hs(i,j,kg,bi,bj) .GT. 0.) THEN
PIs(i,j,kg,bi,bj) = 1. _d 0
Vs(i,j,kg,bi,bj) =
& VH(i,j,kg,bi,bj) / Hs(i,j,kg,bi,bj)
ENDIF
#endif /* LAYERS_VFLUX */
ENDDO
ENDDO
ENDDO
#endif /* LAYERS_THICKNESS */
C --- End bi,bj loop
ENDDO
ENDDO
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP 0
C !ROUTINE: LAYERS_DIAPYCNAL
C !INTERFACE:
SUBROUTINE LAYERS_DIAPYCNAL(
I tracer,iLa,
O TtendSurf, TtendDiffh, TtendDiffr,
O TtendAdvh, TtendAdvr, Ttendtot,
O StendSurf, StendDiffh, StendDiffr,
O StendAdvh, StendAdvr, Stendtot,
O Hc, PIc,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE LAYERS_DIAPYCNAL
C | Calculate the diapycnal velocity in isotracer layers, for a chosen
C | tracer.
C *==========================================================*
C \ev
IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "LAYERS_SIZE.h"
#include "LAYERS.h"
C !INPUT PARAMETERS:
C myThid :: my Thread Id number
C tracer :: potential temperature, salt or potential density prho
C iLa :: layer coordinate index
C TtendSurf :: temperature tendency due to surface forcing times thickness
C TtendDiffh:: temperature tendency due to horiz. diffusion times thickness
C TtendDiffr:: temperature tendency due to vert. diffusion times thickness
C TtendAdvh:: salinity tendency due to horiz. advection times thickness
C TtendAdvr:: salinity tendency due to vert. advection times thickness
C StendSurf :: salinity tendency due to surface forcing times thickness
C StendDiffh:: salinity tendency due to horiz. diffusion times thickness
C StendDiffr:: salinity tendency due to vert. diffusion times thickness
C StendAdvh :: salinity tendency due to horiz. advection times thickness
C StendAdvr :: salinity tendency due to vert. advection times thickness
C Hc :: Layer thickness at the tracer point (m)
C PIw :: 1 if layer exists, 0 otherwise (at tracer point)
INTEGER iLa, myThid
_RL tracer (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr, nSx,nSy)
_RL Hc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
_RL PIc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers,nSx,nSy)
_RL TtendSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
_RL TtendDiffh(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
_RL TtendDiffr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
_RL TtendAdvh (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
_RL TtendAdvr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
_RL Ttendtot (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
_RL StendSurf (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
_RL StendDiffh(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
_RL StendDiffr(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
_RL StendAdvh (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
_RL StendAdvr (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
_RL Stendtot (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
CEOP
#ifdef LAYERS_THERMODYNAMICS
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C !LOCAL VARIABLES:
C bi, bj :: tile indices
C i,j :: horizontal indices
C k :: vertical index for model grid
C kp1 :: vertical index for model grid next cell
C kci :: index from CellIndex
C kg :: index for looping though layers_bounds
C kk :: vertical index for ZZ (fine) grid
C kloc :: local copy of kgu/v to reduce accesses to index arrays
C mSteps :: maximum number of steps for bisection method
C TatC :: temperature at C point
_RL Hcw (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nlayers-1,nSx,nSy)
INTEGER bi, bj
INTEGER i,j,k,kk,kg,kci,kloc
INTEGER mSteps
INTEGER kgc(sNx+1,sNy+1)
INTEGER kgcw(sNx+1,sNy+1)
_RL TatC(sNx+1,sNy+1), dzfac, Tfac, Sfac
LOGICAL errorFlag
CHARACTER*(MAX_LEN_MBUF) msgBuf
#ifdef LAYERS_FINEGRID_DIAPYCNAL
INTEGER kp1
#endif
C -- constants for T and S forcing, gets reset later for rho
Tfac = 1. _d 0
Sfac = 1. _d 0
C compute maximum number of steps for bisection method (approx.
C log2(Nlayers)) as log2(Nlayers) + 1 for safety
mSteps = int(log10(dble(Nlayers))/log10(2. _d 0))+1
C STOP 'DEBUG END: S/R LAYERS_DIAPYCNAL'
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.
kgc(i,j) = Nlayers
kgcw(i,j) = Nlayers-1
ENDDO
ENDDO
C Reset the arrays
C --- These are at the w point
DO kg=1,Nlayers-1
DO j = 1-OLy,sNy+OLy
DO i = 1-OLx,sNx+OLx
TtendSurf (i,j,kg,bi,bj) = 0. _d 0
TtendDiffh(i,j,kg,bi,bj) = 0. _d 0
TtendDiffr(i,j,kg,bi,bj) = 0. _d 0
TtendAdvh(i,j,kg,bi,bj) = 0. _d 0
TtendAdvr(i,j,kg,bi,bj) = 0. _d 0
Ttendtot(i,j,kg,bi,bj) = 0. _d 0
StendSurf (i,j,kg,bi,bj) = 0. _d 0
StendDiffh(i,j,kg,bi,bj) = 0. _d 0
StendDiffr(i,j,kg,bi,bj) = 0. _d 0
StendAdvh(i,j,kg,bi,bj) = 0. _d 0
StendAdvr(i,j,kg,bi,bj) = 0. _d 0
Stendtot(i,j,kg,bi,bj) = 0. _d 0
Hcw(i,j,kg,bi,bj) = 0. _d 0
ENDDO
ENDDO
ENDDO
C --- These are at the c point
DO kg=1,Nlayers
DO j = 1-OLy,sNy+OLy
DO i = 1-OLx,sNx+OLx
Hc(i,j,kg,bi,bj) = 0. _d 0
PIc(i,j,kg,bi,bj) = 0. _d 0
ENDDO
ENDDO
ENDDO
#ifdef LAYERS_FINEGRID_DIAPYCNAL
DO kk=1,NZZ
k = MapIndex(kk)
kci = CellIndex(kk)
DO j = 1,sNy+1
DO i = 1,sNx+1
C ------ Find theta at the V point (south) on the fine Z grid
kp1=k+1
IF (maskC(i,j,kp1,bi,bj).EQ.zeroRS) kp1=k
TatC(i,j) = MapFact(kk) * tracer(i,j,k,bi,bj) +
& (1. _d 0 -MapFact(kk)) * tracer(i,j,kp1,bi,bj)
ENDDO
ENDDO
#else
DO kk=1,Nr
k = kk
kci = kk
DO j = 1,sNy+1
DO i = 1,sNx+1
TatC(i,j) = tracer(i,j,k,bi,bj)
ENDDO
ENDDO
#endif /* LAYERS_FINEGRID_DIAPYCNAL */
C ------ debugging stuff
c IF (i.EQ.38 .AND. j.EQ.4 .AND. bi.EQ.1 .AND. bj.EQ.1) THEN
c i=38
c j=4
c WRITE(msgBuf,
c & '(A,I3,A,I3,A,I3,A,F7.2,A,F7.2,A,F7.2,A,F7.2,A,F3.1)')
c & 'LAYERS_DEBUG: iLa=', iLa,
c & ', kk=', kk,
c & ', k=', k,
c & ', tracer=', tracer(i,j,k,bi,bj),
c & ', TatC=',TatC(i,j),
c & ', hFacC=',hFacC(i,j,k,bi,bj)
c CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
c & SQUEEZE_RIGHT, myThid )
c ENDIF
C ------ Now that we know T everywhere, determine the binning.
C find the layer indices kgc for the center point
CALL LAYERS_LOCATE(
I layers_bounds(1,iLa),Nlayers,mSteps,sNx,sNy,TatC,
O kgc,
I myThid )
#ifndef TARGET_NEC_SX
C check for failures
IF ( debugLevel .GE. debLevC ) THEN
errorFlag = .FALSE.
DO j = 1,sNy+1
DO i = 1,sNx+1
IF ( kgc(i,j) .LE. 0 ) THEN
WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)')
& 'S/R LAYERS_LOCATE: Could not find a bin in ',
& 'layers_bounds for TatC(',i,',',j,',)=',TatC(i,j)
CALL PRINT_ERROR( msgBuf, myThid )
errorFlag = .TRUE.
ENDIF
ENDDO
ENDDO
IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_DIAPYCNAL'
ENDIF
#endif /* ndef TARGET_NEC_SX */
C find the layer indices kgcw for the w point
CALL LAYERS_LOCATE(
I layers_bounds_w(1,iLa),Nlayers-1,mSteps,sNx,sNy,TatC,
O kgcw,
I myThid )
#ifndef TARGET_NEC_SX
C check for failures
IF ( debugLevel .GE. debLevC ) THEN
errorFlag = .FALSE.
DO j = 1,sNy+1
DO i = 1,sNx+1
IF ( kgcw(i,j) .LE. 0 ) THEN
WRITE(msgBuf,'(2A,I3,A,I3,A,1E14.6)')
& 'S/R LAYERS_LOCATE: Could not find a bin in ',
& 'layers_bounds for TatC(',i,',',j,',)=',TatC(i,j)
CALL PRINT_ERROR( msgBuf, myThid )
errorFlag = .TRUE.
ENDIF
ENDDO
ENDDO
IF ( errorFlag ) STOP 'ABNORMAL END: S/R LAYERS_DIAPYCNAL'
ENDIF
#endif /* ndef TARGET_NEC_SX */
C ------ Augment the bin values
DO j = 1,sNy+1
DO i = 1,sNx+1
#ifdef LAYERS_FINEGRID_DIAPYCNAL
dzfac = dZZf(kk) * hFacC(i,j,kci,bi,bj)
#else
dzfac = dRf(kci) * hFacC(i,j,kci,bi,bj)
#endif /* LAYERS_FINEGRID_DIAPYCNAL */
kloc = kgcw(i,j)
C ------- Thickness at w point
Hcw(i,j,kloc,bi,bj) = Hcw(i,j,kloc,bi,bj)
& + dzfac
C ------- Thickness at c point
Hc(i,j,kgc(i,j),bi,bj) = Hc(i,j,kgc(i,j),bi,bj)
& + dzfac
C ------- Now rescale dzfac to include the layer coordinate spacing
dzfac = dzfac * layers_recip_delta(kloc,iLa)
IF ( layers_num(iLa) .EQ. 3 ) THEN
Tfac = layers_alpha(i,j,kci,bi,bj)
Sfac = layers_beta(i,j,kci,bi,bj)
ENDIF
IF (kci.EQ.1) THEN
C ------- We are in the surface layer
TtendSurf(i,j,kloc,bi,bj) =
& TtendSurf(i,j,kloc,bi,bj) +
& Tfac * dzfac * layers_surfflux(i,j,1,1,bi,bj)
StendSurf(i,j,kloc,bi,bj) =
& StendSurf(i,j,kloc,bi,bj) +
& Sfac * dzfac * layers_surfflux(i,j,1,2,bi,bj)
ENDIF
#ifdef SHORTWAVE_HEATING
TtendSurf(i,j,kloc,bi,bj) =
& TtendSurf(i,j,kloc,bi,bj) +
& Tfac * dzfac * layers_sw(i,j,kci,1,bi,bj)
#endif /* SHORTWAVE_HEATING */
C ------- Diffusion
TtendDiffh(i,j,kloc,bi,bj) =
& TtendDiffh(i,j,kloc,bi,bj) + dzfac * Tfac *
& (layers_dfx(i,j,kci,1,bi,bj)+
& layers_dfy(i,j,kci,1,bi,bj))
TtendDiffr(i,j,kloc,bi,bj) =
& TtendDiffr(i,j,kloc,bi,bj) +
& dzfac * Tfac * layers_dfr(i,j,kci,1,bi,bj)
StendDiffh(i,j,kloc,bi,bj) =
& StendDiffh(i,j,kloc,bi,bj) + dzfac * Sfac *
& (layers_dfx(i,j,kci,2,bi,bj)+
& layers_dfy(i,j,kci,2,bi,bj))
StendDiffr(i,j,kloc,bi,bj) =
& StendDiffr(i,j,kloc,bi,bj) +
& dzfac * Sfac * layers_dfr(i,j,kci,2,bi,bj)
C ------- Advection
TtendAdvh(i,j,kloc,bi,bj) =
& TtendAdvh(i,j,kloc,bi,bj) + dzfac * Tfac *
& (layers_afx(i,j,kci,1,bi,bj)+
& layers_afy(i,j,kci,1,bi,bj))
TtendAdvr(i,j,kloc,bi,bj) =
& TtendAdvr(i,j,kloc,bi,bj) +
& dzfac * Tfac * layers_afr(i,j,kci,1,bi,bj)
StendAdvh(i,j,kloc,bi,bj) =
& StendAdvh(i,j,kloc,bi,bj) + dzfac * Sfac *
& (layers_afx(i,j,kci,2,bi,bj)+
& layers_afy(i,j,kci,2,bi,bj))
StendAdvr(i,j,kloc,bi,bj) =
& StendAdvr(i,j,kloc,bi,bj) +
& dzfac * Sfac * layers_afr(i,j,kci,2,bi,bj)
C -------- Total Tendency
Ttendtot(i,j,kloc,bi,bj) =
& Ttendtot(i,j,kloc,bi,bj) +
& dzfac * Tfac * layers_tottend(i,j,kci,1,bi,bj)
Stendtot(i,j,kloc,bi,bj) =
& Stendtot(i,j,kloc,bi,bj) +
& dzfac * Sfac * layers_tottend(i,j,kci,2,bi,bj)
ENDDO
ENDDO
ENDDO
C-- Now that we know the thicknesses, compute the heaviside function
C-- (Needs another loop through Ng)
DO kg=1,Nlayers
DO j = 1,sNy+1
DO i = 1,sNx+1
IF (Hc(i,j,kg,bi,bj) .GT. 0.) THEN
PIc(i,j,kg,bi,bj) = 1. _d 0
ENDIF
ENDDO
ENDDO
ENDDO
C --- End bi,bj loop
ENDDO
ENDDO
#endif /* LAYERS_THERMODYNAMICS */
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
SUBROUTINE LAYERS_LOCATE(
I xx,n,m,sNx,sNy,x,
O k,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | Find the index(-array) k such that x is between xx(k)
C | and xx(k+1) by bisection, following Press et al.,
C | Numerical Recipes in Fortran. xx must be monotonic.
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C !INPUT PARAMETERS:
C xx :: array of bin-boundaries (layers_boundaries)
C n :: length of xx
C m :: int(log2(n)) + 1 = length of bisection loop
C sNx,sNy :: size of index array and input x
C x :: input array of values
C k :: index array (output)
C myThid :: my Thread Id number
INTEGER n,m,sNx,sNy
_RL xx(1:n+1)
_RL x(snx+1,sny+1)
INTEGER k(snx+1,sny+1)
INTEGER myThid
C !LOCAL VARIABLES:
C i,j :: horizontal indices
C l :: bisection loop index
C kl,ku,km :: work arrays and variables
INTEGER i,j
CEOP
#ifdef TARGET_NEC_SX
INTEGER l, km
INTEGER kl(sNx+1,sNy+1), ku(sNx+1,sNy+1)
C bisection, following Press et al., Numerical Recipes in Fortran,
C mostly, because it can be vectorized
DO j = 1,sNy+1
DO i = 1,sNx+1
kl(i,j)=1
ku(i,j)=n+1
END
DO
END
DO
DO l = 1,m
DO j = 1,sNy+1
DO i = 1,sNx+1
IF (ku(i,j)-kl(i,j).GT.1) THEN
km=(ku(i,j)+kl(i,j))/2
CML IF ((xx(n).GE.xx(1)).EQV.(x(i,j).GE.xx(km))) THEN
IF ( ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))).OR.
& ((xx(n).GE.xx(1)).AND.(x(i,j).GE.xx(km))) ) THEN
kl(i,j)=km
ELSE
ku(i,j)=km
END
IF
END
IF
END
DO
END
DO
END
DO
DO j = 1,sNy+1
DO i = 1,sNx+1
IF ( x(i,j).LT.xx(2) ) THEN
k(i,j)=1
ELSE IF ( x(i,j).GE.xx(n) ) THEN
k(i,j)=n
ELSE
k(i,j)=kl(i,j)
END
IF
END
DO
END
DO
#else
C the old way
DO j = 1,sNy+1
DO i = 1,sNx+1
IF (x(i,j) .GE. xx(n)) THEN
C the point is in the hottest bin or hotter
k(i,j) = n
ELSE IF (x(i,j) .LT. xx(2)) THEN
C the point is in the coldest bin or colder
k(i,j) = 1
ELSE IF ( (x(i,j) .GE. xx(k(i,j)))
& .AND. (x(i,j) .LT. xx(k(i,j)+1)) ) THEN
C already on the right bin -- do nothing
ELSE IF (x(i,j) .GE. xx(k(i,j))) THEN
C have to hunt for the right bin by getting hotter
DO WHILE (x(i,j) .GE. xx(k(i,j)+1))
k(i,j) = k(i,j) + 1
ENDDO
C now xx(k) < x <= xx(k+1)
ELSE IF (x(i,j) .LT. xx(k(i,j)+1)) THEN
C have to hunt for the right bin by getting colder
DO WHILE (x(i,j) .LT. xx(k(i,j)))
k(i,j) = k(i,j) - 1
ENDDO
C now xx(k) <= x < xx(k+1)
ELSE
C that should have covered all the options
k(i,j) = -1
ENDIF
ENDDO
ENDDO
#endif /* TARGET_NEC_SX */
#endif /* ALLOW_LAYERS */
RETURN
END