/[MITgcm]/MITgcm/pkg/gmredi/gmredi_k3d.F
ViewVC logotype

Diff of /MITgcm/pkg/gmredi/gmredi_k3d.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.9 by m_bates, Sun Sep 15 14:31:11 2013 UTC revision 1.10 by m_bates, Fri Sep 27 22:34:35 2013 UTC
# Line 74  C     sigy   :: local d(rho)/dy Line 74  C     sigy   :: local d(rho)/dy
74  C     sigz   :: local d(rho)/dz  C     sigz   :: local d(rho)/dz
75  C     hsurf  :: local surface layer depth  C     hsurf  :: local surface layer depth
76  C     small  :: a small number (to avoid floating point exceptions)  C     small  :: a small number (to avoid floating point exceptions)
77        _RL N2loc        _RL N2loc(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
78        _RL slope        _RL slope
79          _RL slopeC(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
80        _RL Req        _RL Req
81        _RL deltaH        _RL deltaH(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
82        _RL g_reciprho_sq        _RL g_reciprho_sq
83        _RL M4loc        _RL M4loc(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
84        _RL maxDRhoDz        _RL maxDRhoDz
85        _RL sigx, sigy, sigz        _RL sigx, sigy, sigz
86        _RL hsurf        _RL hsurf
# Line 258  C       Not limited Line 259  C       Not limited
259  C     Zeroing some cumulative fields  C     Zeroing some cumulative fields
260        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
261         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
262          eady(i,j)  = zeroRL          eady(i,j)   = zeroRL
263          BVint(i,j) = zeroRL          BVint(i,j)  = zeroRL
264          Ubaro(i,j) = zeroRL          Ubaro(i,j)  = zeroRL
265            deltaH(i,j) = zeroRL
266           ENDDO
267          ENDDO
268          DO k=1,Nr
269           DO j=1-Oly,sNy+Oly
270            DO i=1-Olx,sNx+Olx
271             slopeC(i,j,k)=zeroRL
272            ENDDO
273         ENDDO         ENDDO
274        ENDDO        ENDDO
275    
# Line 382  C     Calculate the Eady growth rate Line 391  C     Calculate the Eady growth rate
391  C     ===============================  C     ===============================
392        DO k=1,Nr        DO k=1,Nr
393    
394           DO j=1-Oly,sNy+Oly-1
395            DO i=1-Olx,sNx+Olx-1
396             M4loc(i,j,k) = g_reciprho_sq*( dSigmaDx(i,j,k)**2
397         &                                 +dSigmaDy(i,j,k)**2 )
398             IF (k.NE.kLow_C(i,j)) THEN
399               N2loc(i,j,k) = op5*(N2(i,j,k)+N2(i,j,k+1))
400             ELSE
401               N2loc(i,j,k) = op5*N2(i,j,k)
402             ENDIF
403            ENDDO
404           ENDDO
405  C      The bottom of the grid cell is shallower than the top  C      The bottom of the grid cell is shallower than the top
406  C      integration level, so, advance the depth.  C      integration level, so, advance the depth.
407         IF (-rF(k+1).LE. GM_K3D_EadyMinDepth) CYCLE         IF (-rF(k+1) .LE. GM_K3D_EadyMinDepth) CYCLE
408    
409  C      Do not bother going any deeper since the top of the  C      Do not bother going any deeper since the top of the
410  C      cell is deeper than the bottom integration level  C      cell is deeper than the bottom integration level
# Line 393  C      cell is deeper than the bottom in Line 413  C      cell is deeper than the bottom in
413  C      We are in the integration depth range  C      We are in the integration depth range
414         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly-1
415          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx-1
416           IF (kLow_C(i,j).GE.k) THEN           IF ( (kLow_C(i,j).GE.k) .AND.
417             IF (k.NE.kLow_C(i,j)) THEN       &        (-hMixLayer(i,j,bi,bj).LE.-rC(k)) ) THEN
              N2loc = op5*(N2(i,j,k)+N2(i,j,k+1))  
            ELSE  
              N2loc = op5*N2(i,j,k)  
            ENDIF  
            M4loc = g_reciprho_sq*( dSigmaDx(i,j,k)**2  
      &                            +dSigmaDy(i,j,k)**2 )  
            slope = SQRT(SQRT(M4loc)/N2loc)  
418    
419               slopeC(i,j,k) = SQRT(SQRT(M4loc(i,j,k))/N2loc(i,j,k))
420  C          Limit the slope.  Note, this is not all the Eady calculations.  C          Limit the slope.  Note, this is not all the Eady calculations.
421             IF (slope.LE.GM_K3D_maxSlope) THEN             IF (slopeC(i,j,k).LE.GM_K3D_maxSlope) THEN
422               eady(i,j) = eady(i,j)               eady(i,j) = eady(i,j)
423       &            + hfacC(i,j,k,bi,bj)*drF(k)*M4loc/(N2loc)       &           + hfacC(i,j,k,bi,bj)*drF(k)*M4loc(i,j,k)/(N2loc(i,j,k))
424             ELSE             ELSE
425                 slopeC(i,j,k) = GM_K3D_maxSlope
426               eady(i,j) = eady(i,j)               eady(i,j) = eady(i,j)
427       &            + hfacC(i,j,k,bi,bj)*drF(k)*SQRT(M4loc)       &            + hfacC(i,j,k,bi,bj)*drF(k)*SQRT(M4loc(i,j,k))
428       &              *GM_K3D_maxSlope*GM_K3D_maxSlope       &              *GM_K3D_maxSlope*GM_K3D_maxSlope
429             ENDIF             ENDIF
430               deltaH(i,j) = deltaH(i,j) + drF(k)
431           ENDIF           ENDIF
432          ENDDO          ENDDO
433         ENDDO         ENDDO
# Line 419  C          Limit the slope.  Note, this Line 435  C          Limit the slope.  Note, this
435    
436        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
437         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
438  C       If the minimum depth for the integration is deeper than ocean  C       If the minimum depth for the integration is deeper than the ocean
439  C       bottom then give the eady growth rate a dummy, non-zero value  C       bottom OR the mixed layer is deeper than the maximum depth of
440  C       to avoid floating point exceptions.  These points are taken care  C       integration, we set the Eady growth rate to something small
441  C       of by setting K3D=GM_K3D_smallK later.  C       to avoid floating point exceptions.
442          IF (-r_Low(i,j,bi,bj).LE.GM_K3D_EadyMinDepth) THEN  C       Later, these areas will be given a small diffusivity.
443            IF (deltaH(i,j).EQ.zeroRL) THEN
444            eady(i,j) = small            eady(i,j) = small
445    
446  C       Otherwise, multiply eady by the various constants to get the  C       Otherwise, divide over the integration and take the square root
447  C       growth rate.  C       to actually find the Eady growth rate.
448          ELSE          ELSE
449            deltaH = MIN(-r_low(i,j,bi,bj),GM_K3D_EadyMaxDepth)            eady(i,j) = SQRT(eady(i,j)/deltaH(i,j))
           deltaH = deltaH - GM_K3D_EadyMinDepth  
           eady(i,j) = SQRT(eady(i,j)/deltaH)  
450                        
451          ENDIF          ENDIF
452    
# Line 445  C     ================================== Line 460  C     ==================================
460         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
461  C       Calculate the Visbeck velocity  C       Calculate the Visbeck velocity
462          Rurms(i,j) = MIN(Rdef(i,j,bi,bj),GM_K3D_maxLurms)          Rurms(i,j) = MIN(Rdef(i,j,bi,bj),GM_K3D_maxLurms)
463            Rurms(i,j) = MAX(Rurms(i,j),GM_K3D_minLurms)
464          urms(i,j,1) = GM_K3D_Lambda*eady(i,j)*Rurms(i,j)          urms(i,j,1) = GM_K3D_Lambda*eady(i,j)*Rurms(i,j)
465  C       Set the bottom urms to zero  C       Set the bottom urms to zero
466          k=kLow_C(i,j)          k=kLow_C(i,j)
# Line 478  C     Calculate the diffusivity (on the Line 494  C     Calculate the diffusivity (on the
494         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
495          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
496           IF (k.LE.kLow_C(i,j)) THEN           IF (k.LE.kLow_C(i,j)) THEN
497             IF (-r_Low(i,j,bi,bj).LE.GM_K3D_EadyMinDepth) THEN             IF (deltaH(i,j).EQ.zeroRL) THEN
498               K3D(i,j,k,bi,bj) = GM_K3D_smallK               K3D(i,j,k,bi,bj) = GM_K3D_smallK
499             ELSE             ELSE
500               IF (urms(i,j,k).EQ.0.0) THEN               IF (urms(i,j,k).EQ.0.0) THEN
501                 K3D(i,j,k,bi,bj) = GM_K3D_smallK                 K3D(i,j,k,bi,bj) = GM_K3D_smallK
502               ELSE               ELSE
503                 umc(i,j,k)  = ubar(i,j,k,bi,bj) - cDopp(i,j)                 umc(i,j,k)  = ubar(i,j,k,bi,bj) - cDopp(i,j)
504                 supp(i,j,k) = 1/( 1 + 4*umc(i,j,k)**2/urms(i,j,1)**2 )                 supp(i,j,k) = 1/( 1 + 10*umc(i,j,k)**2/urms(i,j,1)**2 )
505                 K3D(i,j,k,bi,bj) = GM_K3D_gamma*urms(i,j,k)                 K3D(i,j,k,bi,bj) = GM_K3D_gamma*urms(i,j,k)
506       &                            *Rmix(i,j)*supp(i,j,k)       &                            *Rmix(i,j)*supp(i,j,k)
507               ENDIF               ENDIF
# Line 931  C     Diagnostics Line 947  C     Diagnostics
947          CALL DIAGNOSTICS_FILL(ubar,   'GM_UBAR ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(ubar,   'GM_UBAR ',0,Nr,0,1,1,myThid)
948          CALL DIAGNOSTICS_FILL(modesC(1,:,:,:,bi,bj),          CALL DIAGNOSTICS_FILL(modesC(1,:,:,:,bi,bj),
949       &                                'GM_MODEC',0,Nr,0,1,1,myThid)       &                                'GM_MODEC',0,Nr,0,1,1,myThid)
950            CALL DIAGNOSTICS_FILL(M4loc,  'GM_M4   ',0,Nr,0,1,1,myThid)
951            CALL DIAGNOSTICS_FILL(N2loc,  'GM_N2   ',0,Nr,0,1,1,myThid)
952            CALL DIAGNOSTICS_FILL(slopeC, 'GM_SLOPE',0,Nr,0,1,1,myThid)
953    
954        ENDIF        ENDIF
955  #endif  #endif
956    
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22