pkg/kl10/kl10_calc.F

d8d1486ca1 Jean*0001 #include "KL10_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: KL10_CALC
                0005
                0006 C !INTERFACE: =======================================================
26b1a7a333 Jean*0007       SUBROUTINE KL10_CALC(
                0008      I                bi, bj, sigmaR, myTime, myIter, myThid )
d8d1486ca1 Jean*0009
                0010 C !DESCRIPTION: \bv
26b1a7a333 Jean*0011 C     *==========================================================*
d8d1486ca1 Jean*0012 C     | SUBROUTINE KL10_CALC                                     |
                0013 C     | o Compute all KL10 fields defined in KL10.h              |
26b1a7a333 Jean*0014 C     *==========================================================*
d8d1486ca1 Jean*0015 C     | This subroutine is based on SPEM code                    |
26b1a7a333 Jean*0016 C     *==========================================================*
                0017 C \ev
                0018
d8d1486ca1 Jean*0019 C--------------------------------------------------------------------
                0020
                0021 C JMK
                0022 C global parameters updated by kl_calc
26b1a7a333 Jean*0023 C     KLviscAz  :: KL eddy viscosity coefficient              (m^2/s)
                0024 C     KLdiffKzT :: KL diffusion coefficient for temperature   (m^2/s)
d8d1486ca1 Jean*0025
                0026 C !USES: ============================================================
26b1a7a333 Jean*0027       IMPLICIT NONE
d8d1486ca1 Jean*0028 #include "SIZE.h"
                0029 #include "EEPARAMS.h"
                0030 #include "PARAMS.h"
                0031 #include "EOS.h"
dc4a6ae782 Jean*0032 #include "GRID.h"
d8d1486ca1 Jean*0033 #include "DYNVARS.h"
                0034 #include "FFIELDS.h"
dc4a6ae782 Jean*0035 #include "KL10.h"
d5a47d279f Jean*0036 c#ifdef ALLOW_AUTODIFF_TAMC
                0037 c# include "tamc.h"
                0038 c#endif /* ALLOW_AUTODIFF_TAMC */
d8d1486ca1 Jean*0039
                0040 C !INPUT PARAMETERS: ===================================================
                0041 c Routine arguments
26b1a7a333 Jean*0042 C     bi, bj :: Current tile indices
                0043 C     sigmaR :: Vertical gradient of iso-neutral density
                0044 C     myTime :: Current time in simulation
                0045 C     myIter :: Current time-step number
                0046 C     myThid :: My Thread Id number
                0047       INTEGER bi, bj
                0048       _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0049       _RL     myTime
                0050       INTEGER myIter
                0051       INTEGER myThid
d8d1486ca1 Jean*0052
                0053 #ifdef ALLOW_KL10
                0054 C !LOCAL VARIABLES: ====================================================
                0055 c Local constants
26b1a7a333 Jean*0056 C     iMin, iMax, jMin, jMax :: array computation indices
d8d1486ca1 Jean*0057       INTEGER I, J, K, Km1, JJ
                0058       INTEGER   iMin ,iMax ,jMin ,jMax,di
28038c27d5 Lois*0059       _RL     KLviscTmp, tempu
d8d1486ca1 Jean*0060       _RL     b0, buoyFreqf, buoyFreqc, KLviscold,zsum,zsums
28038c27d5 Lois*0061 C note: rhoS(0) is never used but could be reached when evaluating
                0062 C       the "DO WHILE" expression around line 106
                0063       _RL     rhoS(0:Nr), RS(1:Nr)
d5a47d279f Jean*0064       _RL     dzp,ec,ep,es,epss(-1:0),epsw(-1:0),dz,KTemp
26b1a7a333 Jean*0065 c     _RL     bF(1:Nr)
                0066 c     _RL     theta_mcb(1:Nr),theta_mcb3(1:Nr)
d8d1486ca1 Jean*0067 C     === Local variables ===
26b1a7a333 Jean*0068 C     msgBuf     :: Informational/error message buffer
d5a47d279f Jean*0069 c     CHARACTER*(MAX_LEN_MBUF) msgBuf
d8d1486ca1 Jean*0070 CEOP
26b1a7a333 Jean*0071
d8d1486ca1 Jean*0072       iMin = 2-OLx
                0073       iMax = sNx+OLx-1
                0074       jMin = 2-OLy
                0075       jMax = sNy+OLy-1
                0076
                0077       DO J=jMin,jMax
                0078          DO I=iMin,iMax
                0079             K=1
28038c27d5 Lois*0080             rhoS(1)=rhoInSitu(I,J,K,bi,bj)
d8d1486ca1 Jean*0081             RS(1)=rC(1)
                0082
                0083             KLeps(I-1,J-1,1,bi,bj)=0.0
26b1a7a333 Jean*0084 c           eps(k-1) = (dz(k-1)*eps0(k-1) +dz(k)*eps0(k))/(dz(k-1)+dz(k))
d8d1486ca1 Jean*0085             ep = 0.0
                0086             dzp = 0.0
                0087
                0088             KLviscAr(I,J,1,bi,bj) = viscArNr(1)
                0089             KLviscold = KLviscAr(I,J,1,bi,bj) ! at previous cell center
28038c27d5 Lois*0090 C Fill in density profile using vertical gradient of iso-neutral density (SigmaR)
                0091 C      (valid also for nonlinear EOS + remove pressure effect)
                0092             DO K=2,Nr
                0093                rhoS(K)= rhoS(K-1) + rkSign*drC(K)*SigmaR(I,J,K)
                0094                RS(K)=rC(K)
                0095             ENDDO
d8d1486ca1 Jean*0096 C get sorted densities rhoS, and the array with the depths from which
                0097 C the density came from, RS.
                0098             DO K=2,Nr
26b1a7a333 Jean*0099 c$$$               WRITE(msgBuf, '(A,I10.10,A,E10.4,A,E10.4)') 'Hellok ', K
d5a47d279f Jean*0100 c$$$     &              -1,' ',theta(I,J,K,bi,bj),' ',rhot
                0101 c$$$               CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0102 c$$$     &                             SQUEEZE_RIGHT, myThid )
d8d1486ca1 Jean*0103                IF ( (rhoS(K).LT.rhoS(K-1)).AND.(maskC(I,J,K,bi
                0104      &              ,bj).GT.0)) THEN
                0105                   JJ=K-1
                0106                   DO WHILE ( (JJ.GT.0).AND.(rhoS(K).LT.rhoS(JJ)) )
26b1a7a333 Jean*0107 c                    write(*,*) K,JJ,rhoS(K),rhoS(JJ)
d8d1486ca1 Jean*0108                      JJ=JJ-1
                0109                   ENDDO
                0110                   rhoS(JJ+1:K)=cshift(rhoS(JJ+1:K),-1)
                0111                   RS(JJ+1:K)=cshift(RS(JJ+1:K),-1)
                0112                ENDIF
                0113             ENDDO
                0114
                0115 C RS-R is dz....
                0116 C recip_drC=inverse distanance between centers,
                0117 C first is between surface and first center
78524d1402 Jean*0118 C diffKrNrS(K) = viscArNr(K) = background value
d8d1486ca1 Jean*0119
                0120             KLdiffKr(I,J,1,bi,bj) = MAX(KLviscAr(I,J,1,bi,bj),
dc4a6ae782 Jean*0121 #ifdef ALLOW_3D_DIFFKR
                0122      &                                  diffKr(I,J,1,bi,bj) )
                0123 #else
78524d1402 Jean*0124      &                                  diffKrNrS(1) )
dc4a6ae782 Jean*0125 #endif
d8d1486ca1 Jean*0126 C N at surface = zero or uses gradient
                0127             b0 = MAX(-gravity*mass2rUnit*
                0128      &              (rhoS(1) - rhoS(2))*recip_drC(2),0. _d 0)
26b1a7a333 Jean*0129 c           b0 = 0.
d8d1486ca1 Jean*0130             DO di=-1,0
                0131                epss(di)=0.0
                0132                epsw(di)=0.0
                0133             ENDDO
                0134
                0135             DO K=1,Nr
                0136                IF (K.LT.Nr) THEN
                0137                   buoyFreqf = -gravity*mass2rUnit*
                0138      &              (rhoS(K) - rhoS(K+1))*recip_drC(K+1)
                0139                ELSE
                0140 C N zero OR not zero near bottom (at the end of array)
                0141                   buoyFreqf = -gravity*mass2rUnit*
                0142      &              (rhoS(K-1) - rhoS(K))*recip_drC(K)
                0143 C                  buoyFreqf = 0.
                0144                ENDIF
                0145                buoyFreqf = MAX(buoyFreqf,0. _d 0) ! not < 0
                0146                buoyFreqc = (buoyFreqf + b0)*0.5   ! mean at cell center
                0147
                0148 C viscosity at cell center at K
                0149 C = 0.2*dz^2*N.  0.2 is mixing efficiency.
                0150 C to derive, note K = 0.2\epsilon/N^2.  Then note that
                0151 C \epsilon = dz^2N^3 (Ozmidov scale)
                0152                KLviscTmp = MAX( viscArNr(K), 0.2*(RS(K)-rC(K))*
                0153      &                        (RS(K)-rC(K))*sqrt(buoyFreqc))
                0154
                0155                IF (K.GT.1) THEN
                0156                   Km1=K-1
                0157
                0158 C viscosity at cell face above center at K
                0159                   KTemp = 0.5*(KLviscTmp+KLviscold)
                0160 C Put an upper limit on viscosity to prevent instability when
                0161 C explicit viscosity is C used (e.g. for nonhydrostatic case) SAL
                0162                   KTemp = MIN(KLviscMax,KTemp)
                0163                   KLviscAr(I,J,K,bi,bj) = MAX(KTemp,viscArNr(K))
dc4a6ae782 Jean*0164                   KLdiffKr(I,J,K,bi,bj) = MAX(KTemp,
                0165 #ifdef ALLOW_3D_DIFFKR
                0166      &                                        diffKr(I,J,K,bi,bj) )
                0167 #else
78524d1402 Jean*0168      &                                        diffKrNrS(K) )
dc4a6ae782 Jean*0169 #endif
d8d1486ca1 Jean*0170
                0171 C Compute Epsilon for diagnostics:
                0172 C
                0173 C need to caclulate Im1 and Jm1 epsilon unfortunately...  Here at
                0174 C i-1,j-1 we average the west nu(du/dz)^2 at i-1 and i, and the south
                0175 C nu(dv/dv)^2 at j-1 and j, and then add them
                0176 C
                0177 C dz is calculated from the face distances, with the cells assumed to be
                0178 C half way.  Note the use of hfacW and hfacS to make these correct near
                0179 C bathy.
                0180                   zsum=0.
                0181                   ec=0.0
                0182                   zsums=0.
                0183                   es=0.
                0184                   DO di=-1,0
                0185                      IF (hfacW(I+di,J-1,K,bi,bj).GT.0.000001) THEN
                0186                         dz = 0.5*(drF(K)*hfacW(I+di,J-1,K,bi,bj)
d5a47d279f Jean*0187      &                       +drF(Km1)*hfacW(I+di,J-1,Km1,bi,bj))
d8d1486ca1 Jean*0188                         IF (dz.GT.0.00001) THEN
                0189                            tempu = (uVel(I+di,J-1,Km1,bi,bj)-uVel(I+di,J
d5a47d279f Jean*0190      &                          -1,K,bi,bj))/dz
d8d1486ca1 Jean*0191                            epsw(di)=tempu*tempu*KLviscAr(I+di,J-1,K,bi
d5a47d279f Jean*0192      &                          ,bj)
d8d1486ca1 Jean*0193                            ec=ec+epsw(di)*dz
                0194                            zsum = zsum+dz
                0195                         ENDIF
                0196                      ELSE
                0197 C                       This face is on the seafloor.  set epsilon=the
                0198 C                       previous and dz = half the face.
                0199                         dz=0.5*(drF(Km1)*hfacW(I+di,J-1,Km1,bi ,bj))
                0200                         ec=ec+epsw(di)*dz
                0201                         zsum = zsum+dz
                0202                      ENDIF
                0203 C Now do the v-component
                0204                      IF (hfacS(I-1,J+di,K,bi,bj).GT.0.000001) THEN
                0205                         dz = 0.5*(drF(K)*hfacS(I-1,J+di,K,bi,bj)
d5a47d279f Jean*0206      &                       +drF(Km1)*hfacS(I-1,J+di,Km1,bi,bj))
d8d1486ca1 Jean*0207                         IF (dz.GT.0.00001) THEN
                0208                            tempu = (vVel(I-1,J+di,Km1,bi,bj)-vVel(I-1,J
d5a47d279f Jean*0209      &                          +di,K,bi,bj))/dz
d8d1486ca1 Jean*0210                            epss(di)=tempu*tempu*KLviscAr(I-1,J+di,K,bi
d5a47d279f Jean*0211      &                          ,bj)
d8d1486ca1 Jean*0212                            es = es+epss(di)*dz
                0213                            zsums = zsums+dz
                0214                         ENDIF
                0215                      ELSE
                0216 C                       This face is on the seafloor.  set epsilon=the
                0217 C                       previous and dz = half the face.
                0218                         dz=+0.5*(drF(Km1)*hfacS(I-1,J+di,Km1 ,bi,bj))
                0219                         es = es+epss(di)*dz
                0220                         zsums = zsums+dz
                0221                      ENDIF
                0222                   ENDDO
                0223 C                 take the average of the du/dz terms
                0224                   IF (zsum.GT.0.00001) THEN
                0225                      ec=ec/zsum
                0226                   ENDIF
                0227 C                 take the average of the dv/dz terms
                0228                   IF (zsums.GT.0.00001) THEN
                0229                      es=es/zsums
                0230                   ENDIF
                0231 C add the u and v dissipations:
                0232                   ec=es+ec
                0233
                0234 C Note this ec is defined on cell faces K=2..NR at the center of the
                0235 C cells (i.e. at XC), so its above the density variables.
                0236 C
26b1a7a333 Jean*0237 C So to get at the center of the cells, just average this one and the previous one.
                0238 C  And its a true average because the
d8d1486ca1 Jean*0239
                0240                   KLeps(I-1,J-1,Km1,bi,bj) = 0.5*(ep+ec)
                0241                   IF (Km1.EQ.1) THEN
                0242                      KLeps(I-1,J-1,Km1,bi,bj) = ec
                0243                   ENDIF
                0244                   ep=ec
                0245                ENDIF
                0246 c$$$               WRITE(msgBuf, '(A,I10.10,A,E10.4,A,E10.4)') 'Hellok ', K
d5a47d279f Jean*0247 c$$$     &              -1,' ',theta(I,J,K,bi,bj),' ',KLeps(I-1,J-1,Km1,bi
                0248 c$$$     &              ,bj)
                0249 c$$$               CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0250 c$$$     &                             SQUEEZE_RIGHT, myThid )
d8d1486ca1 Jean*0251
                0252                b0        = buoyFreqf ! at previous cell face
                0253                KLviscold = KLviscTmp ! at previous cell center
                0254             ENDDO
                0255 C           ENDDO K
                0256 C     set on K=Nr
                0257             KLeps(I-1,J-1,Nr,bi,bj) =ep
                0258
                0259          ENDDO
                0260 C           ENDDO J
                0261       ENDDO
                0262 C           ENDDO I
                0263
                0264 #endif /* ALLOW_KL10 */
                0265
                0266       RETURN
                0267       END