Back to home page

MITgcm
 
 

    

File indexing completed on 2023-06-14 05:10:21 UTC

view on githubraw file Latest commit 5b0716a6 on 2023-06-13 17:16:45 UTC
31a3206180 Mart*0001 #include "GGL90_OPTIONS.h"
0320e25227 Mart*0002 #ifdef ALLOW_GMREDI
                0003 # include "GMREDI_OPTIONS.h"
                0004 #endif
                0005 #undef GM_EG_PROGNOSTIC
                0006 
                0007 C--  File ggl90_idemix.F:
                0008 C--   Contents
                0009 C--   o GGL90_IDEMIX
                0010 C--   o IDEMIX_gofx2
                0011 C--   o IDEMIX_hofx1
                0012 
                0013 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
31a3206180 Mart*0014 
                0015 CBOP
                0016 C     !ROUTINE: GGL90_IDEMIX
                0017 C     !INTERFACE: ======================================================
                0018       SUBROUTINE GGL90_IDEMIX(
31f96e9372 Jean*0019      I     bi, bj, hFacI, recip_hFacI, sigmaR,
                0020      O     gTKE,
                0021      I     myTime, myIter, myThid )
31a3206180 Mart*0022 
                0023 C     !DESCRIPTION: \bv
                0024 C     *==========================================================*
                0025 C     | S/R GGL90_IDEMIX
                0026 C     |
57a618d609 Mart*0027 C     | IDEMIX1 model as described in
31a3206180 Mart*0028 C     | - Olbers, D. and Eden, C. (2013), JPO, doi:10.1175/JPO-D-12-0207.1
cdaa8119a1 Jean*0029 C     | in a nutshell:
31a3206180 Mart*0030 C     | computes contribution of internal wave field to vertical mixing
                0031 C     *==========================================================*
                0032 C     \ev
                0033 
                0034 C     !USES:
                0035       IMPLICIT NONE
                0036 C     === Global variables ===
                0037 #include "SIZE.h"
                0038 #include "EEPARAMS.h"
                0039 #include "PARAMS.h"
                0040 #include "DYNVARS.h"
                0041 #include "GGL90.h"
                0042 #include "FFIELDS.h"
                0043 #include "GRID.h"
                0044 
                0045 #ifdef ALLOW_GMREDI
0320e25227 Mart*0046 # include "GMREDI.h"
31a3206180 Mart*0047 #endif
                0048 
                0049 C     !INPUT/OUTPUT PARAMETERS:
                0050 C     === Routine arguments ===
                0051 C     bi, bj :: Current tile indices
cdafb98dea Mart*0052 C     hFacI  :: thickness factors for w-cells (interface)
                0053 C               with reciprocal of hFacI = recip_hFacI
31a3206180 Mart*0054 C     sigmaR :: Vertical gradient of iso-neutral density
5b0716a6b3 Mart*0055 C     gTKE   :: dissipation of IW energy (output of S/R GGL90_IDEMIX)
31a3206180 Mart*0056 C     myTime :: Current time in simulation
                0057 C     myIter :: Current time-step number
                0058 C     myThid :: My Thread Id number
                0059       INTEGER bi, bj
cdafb98dea Mart*0060       _RL       hFacI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0061       _RL recip_hFacI(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
31f96e9372 Jean*0062       _RL        gTKE(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0063       _RL     sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
31a3206180 Mart*0064       _RL     myTime
                0065       INTEGER myIter
                0066       INTEGER myThid
                0067 
                0068 #ifdef ALLOW_GGL90_IDEMIX
0320e25227 Mart*0069 C     !FUNCTIONS:
                0070       _RL  IDEMIX_gofx2, IDEMIX_hofx1
31f96e9372 Jean*0071 #ifdef ALLOW_DIAGNOSTICS
                0072       LOGICAL  DIAGNOSTICS_IS_ON
                0073       EXTERNAL DIAGNOSTICS_IS_ON
                0074 #endif /* ALLOW_DIAGNOSTICS */
0320e25227 Mart*0075 
                0076 C     !LOCAL VARIABLES:
31a3206180 Mart*0077 C     === Local variables ===
                0078       INTEGER iMin ,iMax ,jMin ,jMax
0320e25227 Mart*0079       INTEGER i, j, k, kl, kp1, km1
                0080       INTEGER kSrf, kTop, kBot
31a3206180 Mart*0081       INTEGER errCode
0320e25227 Mart*0082       _RL  deltaTloc
31f96e9372 Jean*0083 C     cstar :: vertical integral over N, eq (13) in Olbers+Eden (2013)
5b0716a6b3 Mart*0084       _RL  fxa, fxb, fxc, cstar, twoOverPi, pijstar, recip_pijstar
0320e25227 Mart*0085       _RL  coordFac, recip_coordFac
                0086       _RL  dfx        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0087       _RL  dfy        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31f96e9372 Jean*0088 c     bN0  :: vertically integrated N
0320e25227 Mart*0089       _RL  bN0        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31f96e9372 Jean*0090       _RL  delta      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0091       _RL  Nsquare    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0320e25227 Mart*0092       _RL  a3d        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0093       _RL  b3d        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0094       _RL  c3d        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
31f96e9372 Jean*0095 C     v0 :: mean lateral group velocity
                0096       _RL  v0         (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0097 C     c0 ::  mean vertical group velocity, defined at interfaces (wVel-like)
0320e25227 Mart*0098       _RL  c0         (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
31f96e9372 Jean*0099 C     tau_d :: dissipation parameter (see Olbers and Eden 2013, eq.12)
                0100       _RL  tau_d      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0320e25227 Mart*0101       _RL  forc       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0102       _RL  gm_forc    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31f96e9372 Jean*0103 #ifdef ALLOW_DIAGNOSTICS
                0104       _RL  osborn_diff(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0105 #endif
31a3206180 Mart*0106 CEOP
                0107 
                0108       iMin = 2-OLx
                0109       iMax = sNx+OLx-1
                0110       jMin = 2-OLy
                0111       jMax = sNy+OLy-1
0320e25227 Mart*0112 
                0113       IF ( usingPCoords ) THEN
                0114        kSrf = Nr
                0115        kTop = Nr
                0116       ELSE
                0117        kSrf =  1
                0118        kTop =  2
                0119       ENDIF
                0120       deltaTloc = dTtracerLev(kSrf)
                0121 
                0122       coordFac = 1. _d 0
                0123       IF ( usingPCoords) coordFac = gravity * rhoConst
                0124       recip_coordFac = 1./coordFac
31a3206180 Mart*0125 
5b0716a6b3 Mart*0126       twoOverPi     = 2. _d 0/PI
                0127       pijstar       = PI*IDEMIX_jstar
                0128       recip_pijstar = 1. _d 0 / pijstar
31f96e9372 Jean*0129 
31a3206180 Mart*0130 C     Initialize local fields
                0131       DO k = 1, Nr
                0132        DO j=1-OLy,sNy+OLy
                0133         DO i=1-OLx,sNx+OLx
31f96e9372 Jean*0134          gTKE(i,j,k)    = 0. _d 0
31a3206180 Mart*0135          Nsquare(i,j,k) = 0. _d 0
                0136          delta(i,j,k)   = 0. _d 0
                0137          a3d(i,j,k)     = 0. _d 0
                0138          b3d(i,j,k)     = 1. _d 0
                0139          c3d(i,j,k)     = 0. _d 0
                0140          c0(i,j,k)      = 0. _d 0
31f96e9372 Jean*0141          v0(i,j,k)      = 0. _d 0
                0142          tau_d(i,j,k)   = 0. _d 0
31a3206180 Mart*0143          forc(i,j,k)    = 0. _d 0
31f96e9372 Jean*0144 #ifdef ALLOW_DIAGNOSTICS
                0145          osborn_diff(i,j,k) = 0. _d 0
                0146 #endif
31a3206180 Mart*0147         ENDDO
                0148        ENDDO
                0149       ENDDO
                0150       DO j=1-OLy,sNy+OLy
                0151        DO i=1-OLx,sNx+OLx
                0152          dfx(i,j) = 0. _d 0
                0153          dfy(i,j) = 0. _d 0
                0154          bN0(i,j) = 0. _d 0
                0155          gm_forc(i,j) = 0. _d 0
                0156        ENDDO
                0157       ENDDO
                0158 c-----------------------------------------------------------------------
                0159 c     allow for IW everywhere by limiting buoyancy freq.
                0160 c-----------------------------------------------------------------------
1161225441 Jean*0161       DO k=2,Nr
31a3206180 Mart*0162        DO j=1-OLy,sNy+OLy
                0163         DO i=1-OLx,sNx+OLx
                0164          Nsquare(i,j,k) = gravity*gravitySign*recip_rhoConst
0320e25227 Mart*0165      &                  * sigmaR(i,j,k) * coordFac
5b0716a6b3 Mart*0166 #ifdef GGL90_IDEMIX_CVMIX_VERSION
                0167          Nsquare(i,j,k)= MAX( 0. _d 0, Nsquare(i,j,k) )
                0168 #else
0320e25227 Mart*0169          fxb = MAX( 1. _d -6, ABS( fCori(i,j,bi,bj) ))
                0170          Nsquare(i,j,k)= MAX( 100.*fxb*fxb, Nsquare(i,j,k) )
1161225441 Jean*0171      &                 *maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
5b0716a6b3 Mart*0172 #endif
1161225441 Jean*0173         ENDDO
                0174        ENDDO
                0175       ENDDO
31a3206180 Mart*0176 c-----------------------------------------------------------------------
                0177 c     vertically integrated N
                0178 c-----------------------------------------------------------------------
1161225441 Jean*0179       DO k=2,Nr
31a3206180 Mart*0180        DO j=1-OLy,sNy+OLy
                0181         DO i=1-OLx,sNx+OLx
                0182            bN0(i,j)=bN0(i,j)
0320e25227 Mart*0183      &       +SQRT(Nsquare(i,j,k))*drC(k)*recip_coordFac*hFacI(i,j,k)
1161225441 Jean*0184         ENDDO
                0185        ENDDO
                0186       ENDDO
31a3206180 Mart*0187 c-----------------------------------------------------------------------
                0188 c     vertical and horizontal group velocities
                0189 c     and constant for dissipation
                0190 c-----------------------------------------------------------------------
1161225441 Jean*0191       DO k=2,Nr
31a3206180 Mart*0192        DO j=1-OLy,sNy+OLy
                0193         DO i=1-OLx,sNx+OLx
5b0716a6b3 Mart*0194 #ifdef GGL90_IDEMIX_CVMIX_VERSION
                0195           fxb   = ABS( fCori(i,j,bi,bj) )
                0196           fxa   = SQRT(Nsquare(i,j,k))/(1. _d -22 + fxb)
                0197           cstar = MAX(1. _d -2, bN0(i,j)*recip_pijstar)
                0198 #else
                0199 C     cstar is not limited from below (to say 1e-2), instead Nsquare is
                0200 C     limited from below to (10 * max(1e-6, fCori))**2
                0201           fxb   = MAX( 1. _d -6, ABS( fCori(i,j,bi,bj) ))
                0202           fxa   = SQRT(Nsquare(i,j,k))/fxb
                0203           cstar = bN0(i,j)*recip_pijstar
                0204 #endif
0320e25227 Mart*0205           c0(i,j,k)=MAX(0. _d 0,
31f96e9372 Jean*0206      &             cstar*IDEMIX_gamma*IDEMIX_gofx2(fxa,twoOverPI))
                0207           v0(i,j,k)=MAX(0. _d 0,
                0208      &             cstar*IDEMIX_gamma*IDEMIX_hofx1(fxa,twoOverPI))
5b0716a6b3 Mart*0209 C     next two lines: fxc = ACOSH( MAX(1,fxa) )
0320e25227 Mart*0210           fxc = MAX( 1. _d 0 , fxa )
                0211           fxc = LOG( fxc + SQRT( fxc*fxc -1.))
5b0716a6b3 Mart*0212 #ifdef GGL90_IDEMIX_CVMIX_VERSION
                0213           tau_d(i,j,k) = MAX( 1. _d -4, IDEMIX_mu0*fxb*fxc/cstar**2 )
                0214 #else
                0215           tau_d(i,j,k) = IDEMIX_mu0*fxb*fxc
                0216      &         * ( pijstar/(GGL90eps+bN0(i,j)) )**2
                0217 #endif
1161225441 Jean*0218         ENDDO
                0219        ENDDO
                0220       ENDDO
555f1df2de Mart*0221       IF ( IDEMIX_tau_h .GT. 0. _d 0 ) THEN
fa8d87e2db Jean*0222 C     horizontal diffusion of IW energy can become unstable for long
555f1df2de Mart*0223 C     time steps, so limit horizontal group velocity to satisfy simple
                0224 C     CFL-like criterion:
                0225 C     tau_h V0**2 *dt/dx**2 < 0.25 <=> V0 < sqrt( 0.25 * dx**2/(dt*tau_h) )
                0226        DO k=2,Nr
                0227         DO j=1-OLy,sNy+OLy
                0228          DO i=1-OLx,sNx+OLx
0320e25227 Mart*0229           fxa = SQRT( 1. _d 0/( deltaTloc * IDEMIX_tau_h ) )
                0230           fxb = 0.5*MIN( _dxF(i,j,bi,bj), _dyF(i,j,bi,bj) )*fxa
31f96e9372 Jean*0231           v0(i,j,k) = MIN( v0(i,j,k), fxb )
555f1df2de Mart*0232          ENDDO
                0233         ENDDO
                0234        ENDDO
                0235       ENDIF
31a3206180 Mart*0236 c-----------------------------------------------------------------------
                0237 c     forcing by mesoscale GM
                0238 c-----------------------------------------------------------------------
                0239 
cdaa8119a1 Jean*0240 c     vertically integrated forcing
31a3206180 Mart*0241 #ifdef ALLOW_GMREDI
0320e25227 Mart*0242       IF (useGmredi) THEN
31a3206180 Mart*0243 #ifdef GM_EG_PROGNOSTIC
1161225441 Jean*0244        DO k=1,Nr
31a3206180 Mart*0245         DO j=1-OLy,sNy+OLy
                0246          DO i=1-OLx,sNx+OLx
                0247            gm_forc(i,j) = gm_forc(i,j)
0320e25227 Mart*0248      &         +GM_EG_diss(i,j,k,bi,bj)
                0249      &         *drF(k)*recip_coordFac*hFacC(i,j,k,bi,bj)
1161225441 Jean*0250          ENDDO
                0251         ENDDO
                0252        ENDDO
0320e25227 Mart*0253 #else /* GM_EG_PROGNOSTIC */
1161225441 Jean*0254        DO k=2,Nr
31a3206180 Mart*0255         DO j=1-OLy,sNy+OLy
                0256          DO i=1-OLx,sNx+OLx
                0257            gm_forc(i,j) = gm_forc(i,j)
0320e25227 Mart*0258      &              +MAX( 0. _d 0,Kwz(i,j,k,bi,bj)*Nsquare(i,j,k) )
                0259      &               *drC(k)*recip_coordFac*hFacI(i,j,k)
1161225441 Jean*0260          ENDDO
                0261         ENDDO
                0262        ENDDO
0320e25227 Mart*0263 #endif /* GM_EG_PROGNOSTIC */
                0264       ENDIF
31a3206180 Mart*0265 
0320e25227 Mart*0266       IF (IDEMIX_include_GM .AND. useGmredi) THEN
31a3206180 Mart*0267 c      inject locally
                0268 #ifdef GM_EG_PROGNOSTIC
1161225441 Jean*0269        DO k=2,Nr
31a3206180 Mart*0270         DO j=1-OLy,sNy+OLy
                0271          DO i=1-OLx,sNx+OLx
                0272           forc(i,j,k) = forc(i,j,k)
                0273      &              +.5 _d 0*(GM_EG_diss(i,j,k,bi,bj)+
                0274      &                        GM_EG_diss(i,j,k-1,bi,bj))
1161225441 Jean*0275          ENDDO
                0276         ENDDO
                0277        ENDDO
0320e25227 Mart*0278 #else /* GM_EG_PROGNOSTIC */
1161225441 Jean*0279        DO k=2,Nr
31a3206180 Mart*0280         DO j=1-OLy,sNy+OLy
                0281          DO i=1-OLx,sNx+OLx
                0282           forc(i,j,k) = forc(i,j,k)
0320e25227 Mart*0283      &              +MAX( 0. _d 0,Kwz(i,j,k,bi,bj)*Nsquare(i,j,k) )
1161225441 Jean*0284          ENDDO
                0285         ENDDO
                0286        ENDDO
0320e25227 Mart*0287 #endif /* GM_EG_PROGNOSTIC */
                0288       ENDIF
31a3206180 Mart*0289 
0320e25227 Mart*0290       IF (IDEMIX_include_GM_bottom .AND. useGmredi) THEN
31a3206180 Mart*0291 c      inject at bottom box only
                0292        DO j=1-OLy,sNy+OLy
                0293         DO i=1-OLx,sNx+OLx
0320e25227 Mart*0294          IF ( usingPCoords ) THEN
                0295           kBot=MIN(kSurfC(i,j,bi,bj)+1,Nr)
                0296          ELSE
                0297           kBot   = MAX(kLowC(i,j,bi,bj),1)
                0298          ENDIF
                0299          forc(i,j,kBot)=forc(i,j,kBot)
                0300      &     + gm_forc(i,j)*recip_drC(kBot)*coordFac
                0301      &                  *recip_hFacI(i,j,kBot)
1161225441 Jean*0302         ENDDO
                0303        ENDDO
0320e25227 Mart*0304       ENDIF
                0305 #endif /* ALLOW_GMREDI */
31a3206180 Mart*0306 
                0307 c-----------------------------------------------------------------------
                0308 c     horizontal diffusion of IW energy
                0309 c-----------------------------------------------------------------------
1161225441 Jean*0310        DO k=2,Nr
0320e25227 Mart*0311         kl = k
                0312         IF ( usingPCoords ) kl = k - 1
31a3206180 Mart*0313         DO j=1-OLy,sNy+OLy
                0314          dfx(1-OLx,j)=0. _d 0
                0315          DO i=1-OLx+1,sNx+OLx
                0316           fxa = IDEMIX_tau_h*0.5 _d 0*(
31f96e9372 Jean*0317      &        v0(i-1,j,k)*maskC(i-1,j,kl,bi,bj)
                0318      &       +v0(i  ,j,k)*maskC(i  ,j,kl,bi,bj))
31a3206180 Mart*0319           dfx(i,j) = -fxa*_dyG(i,j,bi,bj)*drC(k)
0320e25227 Mart*0320      &                *(MIN(.5 _d 0,_hFacW(i,j,k-1,bi,bj) ) +
                0321      &                  MIN(.5 _d 0,_hFacW(i,j,k  ,bi,bj) ) )
31a3206180 Mart*0322      &      *_recip_dxC(i,j,bi,bj)
31f96e9372 Jean*0323      &      *(v0(i  ,j,k)*IDEMIX_E(i  ,j,k,bi,bj)
                0324      &       -v0(i-1,j,k)*IDEMIX_E(i-1,j,k,bi,bj))
0320e25227 Mart*0325      &         *maskW(i,j,kl,bi,bj) ! paranoia setting
31a3206180 Mart*0326          ENDDO
                0327         ENDDO
                0328         DO i=1-OLx,sNx+OLx
                0329          dfy(i,1-OLy)=0. _d 0
                0330         ENDDO
                0331         DO j=1-OLy+1,sNy+OLy
                0332          DO i=1-OLx,sNx+OLx
                0333           fxa = IDEMIX_tau_h*0.5 _d 0*(
31f96e9372 Jean*0334      &        v0(i,j  ,k)*maskC(i,j  ,kl,bi,bj)
                0335      &       +v0(i,j-1,k)*maskC(i,j-1,kl,bi,bj) )
31a3206180 Mart*0336           dfy(i,j) = -fxa*_dxG(i,j,bi,bj)*drC(k)
0320e25227 Mart*0337      &                *(MIN(.5 _d 0,_hFacS(i,j,k-1,bi,bj) ) +
                0338      &                  MIN(.5 _d 0,_hFacS(i,j,k  ,bi,bj) ) )
31a3206180 Mart*0339      &      *_recip_dyC(i,j,bi,bj)
31f96e9372 Jean*0340      &      *(v0(i,j  ,k)*IDEMIX_E(i,j  ,k,bi,bj)
                0341      &       -v0(i,j-1,k)*IDEMIX_E(i,j-1,k,bi,bj))
0320e25227 Mart*0342      &         *maskS(i,j,kl,bi,bj) ! paranoia setting
31a3206180 Mart*0343          ENDDO
                0344         ENDDO
                0345 c-----------------------------------------------------------------------
                0346 C     Compute divergence of fluxes, add time tendency
                0347 c-----------------------------------------------------------------------
                0348         DO j=jMin,jMax
                0349          DO i=iMin,iMax
                0350           IDEMIX_E(i,j,k,bi,bj) = IDEMIX_E(i,j,k,bi,bj)
0320e25227 Mart*0351      &       + deltaTloc*(-recip_drC(k)*recip_rA(i,j,bi,bj)
cdafb98dea Mart*0352      &                   *recip_hFacI(i,j,k)
31a3206180 Mart*0353      &         *((dfx(i+1,j)-dfx(i,j))+(dfy(i,j+1)-dfy(i,j)) )  )
0320e25227 Mart*0354      &         *maskC(i,j,kl,bi,bj) ! paranoia setting
31a3206180 Mart*0355          ENDDO
                0356         ENDDO
1161225441 Jean*0357        ENDDO ! k loop
31a3206180 Mart*0358 c-----------------------------------------------------------------------
                0359 c      add interior forcing e.g. by mesoscale GM
                0360 c-----------------------------------------------------------------------
1161225441 Jean*0361       DO k=2,Nr
31a3206180 Mart*0362        DO j=jMin,jMax
                0363         DO i=iMin,iMax
                0364           IDEMIX_E(i,j,k,bi,bj) = IDEMIX_E(i,j,k,bi,bj)
0320e25227 Mart*0365      &                      + forc(i,j,k)*deltaTloc
1161225441 Jean*0366         ENDDO
                0367        ENDDO
                0368       ENDDO
31a3206180 Mart*0369 c-----------------------------------------------------------------------
                0370 c      solve vertical diffusion implicitly
                0371 c-----------------------------------------------------------------------
                0372 
                0373 C     delta_k = dt tau_v /drF_k (c_k+c_k+1)/2
0320e25227 Mart*0374 C     delta(1) and delta(Nr) are zero by initialisation
31a3206180 Mart*0375       DO k=2,Nr-1
                0376        DO j=jMin,jMax
                0377         DO i=iMin,iMax
0320e25227 Mart*0378          delta(i,j,k)  = deltaTloc*IDEMIX_tau_v
                0379      &        *recip_drF(k)*coordFac*recip_hFacC(i,j,k,bi,bj)
                0380      &        *.5 _d 0*(c0(i,j,k)+c0(i,j,k+1))
31a3206180 Mart*0381         ENDDO
                0382        ENDDO
                0383       ENDDO
0320e25227 Mart*0384       IF ( usingPCoords ) THEN
                0385        DO j=jMin,jMax
                0386         DO i=iMin,iMax
                0387          kBot = MIN(kSurfC(i,j,bi,bj),Nr)
                0388          delta(i,j,kBot) = 0. _d 0
                0389         ENDDO
31a3206180 Mart*0390        ENDDO
0320e25227 Mart*0391       ELSE
                0392        DO j=jMin,jMax
                0393         DO i=iMin,iMax
                0394          kBot = MAX(kLowC(i,j,bi,bj),1)
                0395          delta(i,j,kBot) = 0. _d 0
                0396         ENDDO
                0397        ENDDO
                0398       ENDIF
31a3206180 Mart*0399 
                0400 C--   Lower diagonal  for E_(k-1) : -delta_k-1 c_k-1/drC_k
0320e25227 Mart*0401 C     but leaving the contribution of c0_k-1 for later
                0402       DO k=2,Nr
31a3206180 Mart*0403        DO j=jMin,jMax
                0404         DO i=iMin,iMax
                0405 C-       No need for maskC(k-1) with recip_hFacC(k-1) in delta(k-1)
0320e25227 Mart*0406          a3d(i,j,k) = -delta(i,j,k-1)
                0407      &        *recip_drC(k)*coordFac*recip_hFacI(i,j,k)
                0408      &        *maskC(i,j,k,bi,bj)
31a3206180 Mart*0409         ENDDO
                0410        ENDDO
                0411       ENDDO
                0412 
cdaa8119a1 Jean*0413 C--   Upper diagonal for E_(k+1):  delta_k c_k+1/drC_k
0320e25227 Mart*0414 C     but leaving the contribution of c0_k+1 for later
31a3206180 Mart*0415       DO k=2,Nr
                0416        DO j=jMin,jMax
                0417         DO i=iMin,iMax
                0418 C-       No need for maskC(k) with recip_hFacC(k) in delta(k)
0320e25227 Mart*0419          c3d(i,j,k) = -delta(i,j,k)
                0420      &        *recip_drC(k)*coordFac*recip_hFacI(i,j,k)
31a3206180 Mart*0421      &        *maskC(i,j,k-1,bi,bj)
                0422         ENDDO
                0423        ENDDO
                0424       ENDDO
                0425 
0320e25227 Mart*0426 C     treat bottom and surface boundaries for coeffients
                0427 C     of upper and lower diagonal by masking
                0428       IF ( usingPCoords ) THEN
                0429        DO j=jMin,jMax
                0430         DO i=iMin,iMax
                0431 C     a3d at bottom is zero
                0432          kBot=MIN(kSurfC(i,j,bi,bj)+1,Nr)
                0433          a3d(i,j,kBot) = 0. _d 0
                0434 C     for p-coords, c3d is zero at the surface, too
                0435          c3d(i,j,Nr)   = 0. _d 0
                0436         ENDDO
                0437        ENDDO
                0438       ELSEIF ( usingZCoords ) THEN
                0439        DO j=jMin,jMax
                0440         DO i=iMin,iMax
                0441 C     c3d at bottom is zero
                0442          kBot = MAX(kLowC(i,j,bi,bj),1)
                0443          c3d(i,j,kBot) = 0. _d 0
                0444 C     for z-coords, a3d is zero at the surface (level 2), too
                0445          a3d(i,j,kTop) = 0. _d 0
                0446         ENDDO
                0447        ENDDO
                0448       ENDIF
                0449 C     For k=1 there is nothing to solve for
cdaa8119a1 Jean*0450       DO j=jMin,jMax
31a3206180 Mart*0451        DO i=iMin,iMax
0320e25227 Mart*0452 C     so that both off-diagonal coefficients are zero
                0453         a3d(i,j,1) = 0. _d 0
                0454         c3d(i,j,1) = 0. _d 0
                0455 C     and the main diagonal is one (for stability)
                0456         b3d(i,j,1) = 1. _d 0
1161225441 Jean*0457        ENDDO
                0458       ENDDO
0320e25227 Mart*0459 
                0460 C--   Center diagonal
31a3206180 Mart*0461       DO k=2,Nr
                0462        DO j=jMin,jMax
                0463         DO i=iMin,iMax
31f96e9372 Jean*0464          b3d(i,j,k) = 1. _d 0+deltaTloc*tau_d(i,j,k)
31a3206180 Mart*0465      &         *IDEMIX_E(i,j,k,bi,bj)
0320e25227 Mart*0466      &         *maskC(i,j,k,bi,bj)*maskC(i,j,k-1,bi,bj)
                0467      &        - ( a3d(i,j,k) + c3d(i,j,k) ) * c0(i,j,k)
                0468         ENDDO
31a3206180 Mart*0469        ENDDO
                0470       ENDDO
                0471 
0320e25227 Mart*0472 C--   Complete computation of lower and upper diagonal after they have
                0473 C     been used for the center diagonal: add the contribution of c0_k+/-1
                0474       DO k=2,Nr
                0475        kp1=MIN(k+1,Nr)
                0476        km1=MAX(k-1,2)
                0477        DO j=jMin,jMax
                0478         DO i=iMin,iMax
                0479          a3d(i,j,k) = a3d(i,j,k)*c0(i,j,km1)
                0480          c3d(i,j,k) = c3d(i,j,k)*c0(i,j,kp1)
                0481         ENDDO
31a3206180 Mart*0482        ENDDO
                0483       ENDDO
                0484 
0320e25227 Mart*0485 C--   Apply flux boundary condition
                0486       kl = kTop
                0487       IF ( usingPCoords ) kl = kTop - 1
31a3206180 Mart*0488       DO j=jMin,jMax
                0489        DO i=iMin,iMax
0320e25227 Mart*0490         IDEMIX_E(i,j,kTop,bi,bj) = IDEMIX_E(i,j,kTop,bi,bj)
                0491      &       +deltaTloc*IDEMIX_F_s(i,j,bi,bj)
                0492      &       *recip_drC(kTop)*coordFac*recip_hFacI(i,j,kTop)
                0493      &       *maskC(i,j,kl,bi,bj)
31a3206180 Mart*0494        ENDDO
                0495       ENDDO
0320e25227 Mart*0496       IF ( usingZCoords) THEN
                0497        DO j=jMin,jMax
                0498         DO i=iMin,iMax
                0499          kBot = MAX(kLowC(i,j,bi,bj),1)
                0500          IDEMIX_E(i,j,kBot,bi,bj) = IDEMIX_E(i,j,kBot,bi,bj)
                0501      &        -deltaTloc*IDEMIX_F_b(i,j,bi,bj)
                0502      &        *recip_drC(kBot)*coordFac*recip_hFacI(i,j,kBot)
                0503      &        *maskC(i,j,kBot,bi,bj)
                0504         ENDDO
                0505        ENDDO
                0506       ELSEIF ( usingPCoords ) THEN
                0507        DO j=jMin,jMax
                0508         DO i=iMin,iMax
                0509          kBot = MIN(kSurfC(i,j,bi,bj)+1,Nr)
                0510          IDEMIX_E(i,j,kBot,bi,bj) = IDEMIX_E(i,j,kBot,bi,bj)
                0511      &        -deltaTloc*IDEMIX_F_b(i,j,bi,bj)
                0512      &        *recip_drC(kBot)*coordFac*recip_hFacI(i,j,kBot)
                0513      &        *maskC(i,j,kBot-1,bi,bj)
                0514         ENDDO
                0515        ENDDO
                0516       ENDIF
31a3206180 Mart*0517 
                0518 C     solve tri-diagonal system
fa8d87e2db Jean*0519       errCode = -1
31a3206180 Mart*0520       CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
                0521      I                        a3d, b3d, c3d,
                0522      U                        IDEMIX_E(1-OLx,1-OLy,1,bi,bj),
                0523      O                        errCode,
                0524      I                        bi, bj, myThid )
                0525 
5b0716a6b3 Mart*0526 C     generate TKE tendency due to dissipation of IW energy (output)
31f96e9372 Jean*0527       DO k=2,Nr
                0528        DO j=jMin,jMax
                0529         DO i=iMin,iMax
                0530          gTKE(i,j,k) =
                0531      &        tau_d(i,j,k)*IDEMIX_E(i,j,k,bi,bj)*IDEMIX_E(i,j,k,bi,bj)
                0532 C-    to reproduce older results (before adding output arguement gTKE)
                0533 c    &        tau_d(i,j,k)*IDEMIX_E(i,j,k,bi,bj)**2
                0534         ENDDO
                0535        ENDDO
                0536       ENDDO
                0537 
31a3206180 Mart*0538 #ifdef ALLOW_DIAGNOSTICS
                0539       IF ( useDiagnostics ) THEN
                0540 c-----------------------------------------------------------------------
                0541 c     compute diffusivity due to internal wave breaking
                0542 c     assuming local Osborn-Cox balance model
                0543 c     kept for diagnostics only
                0544 c-----------------------------------------------------------------------
31f96e9372 Jean*0545        IF ( DIAGNOSTICS_IS_ON('IDEMIX_K',myThid) ) THEN
                0546         DO k=2,Nr
                0547          DO j=jMin,jMax
                0548           DO i=iMin,iMax
                0549            osborn_diff(i,j,k) = IDEMIX_mixing_efficiency * gTKE(i,j,k)
                0550      &          /MAX(1. _d -12,Nsquare(i,j,k))*maskC(i,j,k,bi,bj)
                0551            osborn_diff(i,j,k) = MIN(IDEMIX_diff_max,osborn_diff(i,j,k))
                0552           ENDDO
31a3206180 Mart*0553          ENDDO
                0554         ENDDO
31f96e9372 Jean*0555         CALL DIAGNOSTICS_FILL( osborn_diff ,'IDEMIX_K',
                0556      &                          0, Nr, 2, bi, bj, myThid )
                0557        ENDIF
                0558        CALL DIAGNOSTICS_FILL( IDEMIX_E ,'IDEMIX_E',0,Nr,1,bi,bj,myThid)
5b0716a6b3 Mart*0559        CALL DIAGNOSTICS_FILL( gTKE ,    'IDEMgTKE',0,Nr,1,bi,bj,myThid)
31f96e9372 Jean*0560        CALL DIAGNOSTICS_FILL( tau_d,    'IDEMIX_t',0,Nr,2,bi,bj,myThid)
                0561        CALL DIAGNOSTICS_FILL( v0,       'IDEMIX_v',0,Nr,2,bi,bj,myThid)
                0562        CALL DIAGNOSTICS_FILL( c0 ,      'IDEMIX_c',0,Nr,2,bi,bj,myThid)
                0563        CALL DIAGNOSTICS_FILL( forc,     'IDEMIX_F',0,Nr,2,bi,bj,myThid)
                0564        CALL DIAGNOSTICS_FILL(IDEMIX_F_b,'IDEM_F_b',0, 1,1,bi,bj,myThid)
                0565        CALL DIAGNOSTICS_FILL(IDEMIX_F_s,'IDEM_F_s',0, 1,1,bi,bj,myThid)
                0566 # ifdef ALLOW_GMREDI
                0567        IF (useGmredi) THEN
                0568         CALL DIAGNOSTICS_FILL( gm_forc, 'IDEM_F_g',0, 1,2,bi,bj,myThid)
                0569        ENDIF
                0570 # endif
31a3206180 Mart*0571       ENDIF
                0572 #endif /* ALLOW_DIAGNOSTICS */
                0573 
                0574 #endif /* ALLOW_GGL90_IDEMIX */
                0575       RETURN
                0576       END
                0577 
                0578 #ifdef ALLOW_GGL90_IDEMIX
                0579 C     helper functions
cdaa8119a1 Jean*0580 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
31f96e9372 Jean*0581       _RL FUNCTION IDEMIX_gofx2(xx,toPI)
cdaa8119a1 Jean*0582       IMPLICIT NONE
31f96e9372 Jean*0583       _RL xx
                0584       _RL toPI                  ! 2.d0/PI
                0585       _RL x,c
cdaa8119a1 Jean*0586       x=MAX(3.d0,xx)
31f96e9372 Jean*0587       c= 1.d0-toPI*ASIN(1.d0/x)
                0588       IDEMIX_gofx2 = toPI/c*0.9d0*x**(-2.d0/3.d0)*(1.-EXP(-x/4.3d0))
cdaa8119a1 Jean*0589       END
                0590 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
31f96e9372 Jean*0591       _RL FUNCTION IDEMIX_hofx1(x,toPI)
cdaa8119a1 Jean*0592       IMPLICIT NONE
                0593       _RL x
31f96e9372 Jean*0594       _RL toPI                  ! 2.d0/PI
                0595       IDEMIX_hofx1 = toPI/(1.d0-toPI*ASIN(1.d0/MAX(1.01d0,x)))
                0596      &     *(x-1.d0)/(x+1.d0)
cdaa8119a1 Jean*0597       END
31a3206180 Mart*0598 #endif /* ALLOW_GGL90_IDEMIX */
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!