C $Header: /u/gcmpack/MITgcm/pkg/kl10/kl10_calc.F,v 1.5 2015/02/23 21:20:15 jmc Exp $
C $Name: $
#include "KL10_OPTIONS.h"
CBOP
C !ROUTINE: KL10_CALC
C !INTERFACE: =======================================================
SUBROUTINE KL10_CALC(
I bi, bj, sigmaR, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE KL10_CALC |
C | o Compute all KL10 fields defined in KL10.h |
C *==========================================================*
C | This subroutine is based on SPEM code |
C *==========================================================*
C \ev
C--------------------------------------------------------------------
C JMK
C global parameters updated by kl_calc
C KLviscAz :: KL eddy viscosity coefficient (m^2/s)
C KLdiffKzT :: KL diffusion coefficient for temperature (m^2/s)
C !USES: ============================================================
IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "EOS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
#include "KL10.h"
c#ifdef ALLOW_AUTODIFF_TAMC
c# include "tamc.h"
c# include "tamc_keys.h"
c#endif /* ALLOW_AUTODIFF_TAMC */
C !INPUT PARAMETERS: ===================================================
c Routine arguments
C bi, bj :: Current tile indices
C sigmaR :: Vertical gradient of iso-neutral density
C myTime :: Current time in simulation
C myIter :: Current time-step number
C myThid :: My Thread Id number
INTEGER bi, bj
_RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL myTime
INTEGER myIter
INTEGER myThid
#ifdef ALLOW_KL10
C !LOCAL VARIABLES: ====================================================
c Local constants
C iMin, iMax, jMin, jMax :: array computation indices
INTEGER I, J, K, Km1, JJ
INTEGER iMin ,iMax ,jMin ,jMax,di
_RL KLviscTmp, rhot, tempu
_RL b0, buoyFreqf, buoyFreqc, KLviscold,zsum,zsums
_RL rhoS(1:Nr),RS(1:Nr)
_RL dzp,ec,ep,es,epss(-1:0),epsw(-1:0),dz,KTemp
c _RL bF(1:Nr)
c _RL theta_mcb(1:Nr),theta_mcb3(1:Nr)
C === Local variables ===
C msgBuf :: Informational/error message buffer
c CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP
iMin = 2-OLx
iMax = sNx+OLx-1
jMin = 2-OLy
jMax = sNy+OLy-1
DO J=jMin,jMax
DO I=iMin,iMax
K=1
CALL FIND_RHO_SCALAR(theta(I,J,K,bi,bj), salt(I,J,K,bi,bj),
& totPhiHyd(I,J,K,bi,bj),rhot,myThid )
rhoS(1)=rhot
RS(1)=rC(1)
KLeps(I-1,J-1,1,bi,bj)=0.0
c eps(k-1) = (dz(k-1)*eps0(k-1) +dz(k)*eps0(k))/(dz(k-1)+dz(k))
ep = 0.0
dzp = 0.0
KLviscAr(I,J,1,bi,bj) = viscArNr(1)
KLviscold = KLviscAr(I,J,1,bi,bj) ! at previous cell center
C get sorted densities rhoS, and the array with the depths from which
C the density came from, RS.
DO K=2,Nr
CALL FIND_RHO_SCALAR(theta(I,J,K,bi,bj), salt(I,J,K,bi
& ,bj),totPhiHyd(I,J,K,bi,bj),rhot,myThid )
rhoS(K)=rhot
RS(K)=rC(K)
c$$$ WRITE(msgBuf, '(A,I10.10,A,E10.4,A,E10.4)') 'Hellok ', K
c$$$ & -1,' ',theta(I,J,K,bi,bj),' ',rhot
c$$$ CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
c$$$ & SQUEEZE_RIGHT, myThid )
IF ( (rhoS(K).LT.rhoS(K-1)).AND.(maskC(I,J,K,bi
& ,bj).GT.0)) THEN
JJ=K-1
DO WHILE ( (JJ.GT.0).AND.(rhoS(K).LT.rhoS(JJ)) )
c write(*,*) K,JJ,rhoS(K),rhoS(JJ)
JJ=JJ-1
ENDDO
rhoS(JJ+1:K)=cshift(rhoS(JJ+1:K),-1)
RS(JJ+1:K)=cshift(RS(JJ+1:K),-1)
ENDIF
ENDDO
C RS-R is dz....
C recip_drC=inverse distanance between centers,
C first is between surface and first center
C diffKrNrS(K) = viscArNr(K) = background value
KLdiffKr(I,J,1,bi,bj) = MAX(KLviscAr(I,J,1,bi,bj),
#ifdef ALLOW_3D_DIFFKR
& diffKr(I,J,1,bi,bj) )
#else
& diffKrNrS(1) )
#endif
C N at surface = zero or uses gradient
b0 = MAX(-gravity*mass2rUnit*
& (rhoS(1) - rhoS(2))*recip_drC(2),0. _d 0)
c b0 = 0.
DO di=-1,0
epss(di)=0.0
epsw(di)=0.0
ENDDO
DO K=1,Nr
IF (K.LT.Nr) THEN
buoyFreqf = -gravity*mass2rUnit*
& (rhoS(K) - rhoS(K+1))*recip_drC(K+1)
ELSE
C N zero OR not zero near bottom (at the end of array)
buoyFreqf = -gravity*mass2rUnit*
& (rhoS(K-1) - rhoS(K))*recip_drC(K)
C buoyFreqf = 0.
ENDIF
buoyFreqf = MAX(buoyFreqf,0. _d 0) ! not < 0
buoyFreqc = (buoyFreqf + b0)*0.5 ! mean at cell center
C viscosity at cell center at K
C = 0.2*dz^2*N. 0.2 is mixing efficiency.
C to derive, note K = 0.2\epsilon/N^2. Then note that
C \epsilon = dz^2N^3 (Ozmidov scale)
KLviscTmp = MAX( viscArNr(K), 0.2*(RS(K)-rC(K))*
& (RS(K)-rC(K))*sqrt(buoyFreqc))
IF (K.GT.1) THEN
Km1=K-1
C viscosity at cell face above center at K
KTemp = 0.5*(KLviscTmp+KLviscold)
C Put an upper limit on viscosity to prevent instability when
C explicit viscosity is C used (e.g. for nonhydrostatic case) SAL
KTemp = MIN(KLviscMax,KTemp)
KLviscAr(I,J,K,bi,bj) = MAX(KTemp,viscArNr(K))
KLdiffKr(I,J,K,bi,bj) = MAX(KTemp,
#ifdef ALLOW_3D_DIFFKR
& diffKr(I,J,K,bi,bj) )
#else
& diffKrNrS(K) )
#endif
C Compute Epsilon for diagnostics:
C
C need to caclulate Im1 and Jm1 epsilon unfortunately... Here at
C i-1,j-1 we average the west nu(du/dz)^2 at i-1 and i, and the south
C nu(dv/dv)^2 at j-1 and j, and then add them
C
C dz is calculated from the face distances, with the cells assumed to be
C half way. Note the use of hfacW and hfacS to make these correct near
C bathy.
zsum=0.
ec=0.0
zsums=0.
es=0.
DO di=-1,0
IF (hfacW(I+di,J-1,K,bi,bj).GT.0.000001) THEN
dz = 0.5*(drF(K)*hfacW(I+di,J-1,K,bi,bj)
& +drF(Km1)*hfacW(I+di,J-1,Km1,bi,bj))
IF (dz.GT.0.00001) THEN
tempu = (uVel(I+di,J-1,Km1,bi,bj)-uVel(I+di,J
& -1,K,bi,bj))/dz
epsw(di)=tempu*tempu*KLviscAr(I+di,J-1,K,bi
& ,bj)
ec=ec+epsw(di)*dz
zsum = zsum+dz
ENDIF
ELSE
C This face is on the seafloor. set epsilon=the
C previous and dz = half the face.
dz=0.5*(drF(Km1)*hfacW(I+di,J-1,Km1,bi ,bj))
ec=ec+epsw(di)*dz
zsum = zsum+dz
ENDIF
C Now do the v-component
IF (hfacS(I-1,J+di,K,bi,bj).GT.0.000001) THEN
dz = 0.5*(drF(K)*hfacS(I-1,J+di,K,bi,bj)
& +drF(Km1)*hfacS(I-1,J+di,Km1,bi,bj))
IF (dz.GT.0.00001) THEN
tempu = (vVel(I-1,J+di,Km1,bi,bj)-vVel(I-1,J
& +di,K,bi,bj))/dz
epss(di)=tempu*tempu*KLviscAr(I-1,J+di,K,bi
& ,bj)
es = es+epss(di)*dz
zsums = zsums+dz
ENDIF
ELSE
C This face is on the seafloor. set epsilon=the
C previous and dz = half the face.
dz=+0.5*(drF(Km1)*hfacS(I-1,J+di,Km1 ,bi,bj))
es = es+epss(di)*dz
zsums = zsums+dz
ENDIF
ENDDO
C take the average of the du/dz terms
IF (zsum.GT.0.00001) THEN
ec=ec/zsum
ENDIF
C take the average of the dv/dz terms
IF (zsums.GT.0.00001) THEN
es=es/zsums
ENDIF
C add the u and v dissipations:
ec=es+ec
C Note this ec is defined on cell faces K=2..NR at the center of the
C cells (i.e. at XC), so its above the density variables.
C
C So to get at the center of the cells, just average this one and the previous one.
C And its a true average because the
KLeps(I-1,J-1,Km1,bi,bj) = 0.5*(ep+ec)
IF (Km1.EQ.1) THEN
KLeps(I-1,J-1,Km1,bi,bj) = ec
ENDIF
ep=ec
ENDIF
c$$$ WRITE(msgBuf, '(A,I10.10,A,E10.4,A,E10.4)') 'Hellok ', K
c$$$ & -1,' ',theta(I,J,K,bi,bj),' ',KLeps(I-1,J-1,Km1,bi
c$$$ & ,bj)
c$$$ CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
c$$$ & SQUEEZE_RIGHT, myThid )
b0 = buoyFreqf ! at previous cell face
KLviscold = KLviscTmp ! at previous cell center
ENDDO
C ENDDO K
C set on K=Nr
KLeps(I-1,J-1,Nr,bi,bj) =ep
ENDDO
C ENDDO J
ENDDO
C ENDDO I
#endif /* ALLOW_KL10 */
RETURN
END