/[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.11 by m_bates, Sat Sep 28 17:59:30 2013 UTC revision 1.13 by m_bates, Mon Oct 21 18:46:05 2013 UTC
# Line 81  C     small  :: a small number (to avoid Line 81  C     small  :: a small number (to avoid
81        _RL deltaH(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL deltaH(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
82        _RL g_reciprho_sq        _RL g_reciprho_sq
83        _RL M4loc(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL M4loc(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
84          _RL M4onN2(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
85        _RL maxDRhoDz        _RL maxDRhoDz
86        _RL sigx, sigy, sigz        _RL sigx, sigy, sigz
87        _RL hsurf        _RL hsurf
# Line 118  C     Ubaro   :: Barotropic velocity (m/ Line 119  C     Ubaro   :: Barotropic velocity (m/
119    
120  C     Rmid     :: Rossby radius (m)  C     Rmid     :: Rossby radius (m)
121  C     KPV      :: Diffusivity (m**2/s)  C     KPV      :: Diffusivity (m**2/s)
122    C     Kdqdx    :: diffusivity multiplied by zonal PV gradient
123    C     Kdqdy    :: diffusivity multiplied by meridional PV gradient
124  C     SlopeX   :: isopycnal slope in x direction  C     SlopeX   :: isopycnal slope in x direction
125  C     SlopeY   :: isopycnal slope in y direction  C     SlopeY   :: isopycnal slope in y direction
126  C     dSigmaDx :: sigmaX averaged onto tracer grid  C     dSigmaDx :: sigmaX averaged onto tracer grid
# Line 132  C     coriV    :: As for cori, but, at V Line 135  C     coriV    :: As for cori, but, at V
135  C     surfkz   :: Depth of surface layer (in r units)  C     surfkz   :: Depth of surface layer (in r units)
136        _RL Rmid(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL Rmid(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
137        _RL KPV(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KPV(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
138          _RL Kdqdy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
139          _RL Kdqdx(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
140        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
141        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
142        _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL dSigmaDx(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
# Line 416  C      We are in the integration depth r Line 421  C      We are in the integration depth r
421           IF ( (kLow_C(i,j).GE.k) .AND.           IF ( (kLow_C(i,j).GE.k) .AND.
422       &        (-hMixLayer(i,j,bi,bj).LE.-rC(k)) ) THEN       &        (-hMixLayer(i,j,bi,bj).LE.-rC(k)) ) THEN
423    
424             slopeC(i,j,k) = SQRT(SQRT(M4loc(i,j,k))/N2loc(i,j,k))             slopeC(i,j,k) = SQRT(M4loc(i,j,k))/N2loc(i,j,k)
425  C          Limit the slope.  Note, this is not all the Eady calculations.  C          Limit the slope.  Note, this is not all the Eady calculations.
426             IF (slopeC(i,j,k).LE.GM_K3D_maxSlope) THEN             IF (slopeC(i,j,k).LE.GM_maxSlope) THEN
427               eady(i,j) = eady(i,j)               M4onN2(i,j,k) = M4loc(i,j,k)/N2loc(i,j,k)
      &           + hfacC(i,j,k,bi,bj)*drF(k)*M4loc(i,j,k)/(N2loc(i,j,k))  
428             ELSE             ELSE
429               slopeC(i,j,k) = GM_K3D_maxSlope               slopeC(i,j,k) = GM_maxslope
430               eady(i,j) = eady(i,j)               M4onN2(i,j,k) = SQRT(M4loc(i,j,k))*GM_maxslope
      &            + hfacC(i,j,k,bi,bj)*drF(k)*SQRT(M4loc(i,j,k))  
      &              *GM_K3D_maxSlope*GM_K3D_maxSlope  
431             ENDIF             ENDIF
432               eady(i,j)   = eady(i,j)  
433         &                   + hfacC(i,j,k,bi,bj)*drF(k)*M4onN2(i,j,k)
434             deltaH(i,j) = deltaH(i,j) + drF(k)             deltaH(i,j) = deltaH(i,j) + drF(k)
435           ENDIF           ENDIF
436          ENDDO          ENDDO
# Line 459  C     ================================== Line 463  C     ==================================
463        DO j=1-Oly+1,sNy+Oly        DO j=1-Oly+1,sNy+Oly
464         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
465  C       Calculate the Visbeck velocity  C       Calculate the Visbeck velocity
466          Rurms(i,j) = MIN(Rdef(i,j,bi,bj),GM_K3D_maxLurms)          Rurms(i,j) = MIN(Rdef(i,j,bi,bj),GM_K3D_Rmax)
         Rurms(i,j) = MAX(Rurms(i,j),GM_K3D_minLurms)  
467          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)
468  C       Set the bottom urms to zero  C       Set the bottom urms to zero
469          k=kLow_C(i,j)          k=kLow_C(i,j)
# Line 471  C       Calculate the Rhines scale Line 474  C       Calculate the Rhines scale
474    
475  C       Calculate the estimated length scale  C       Calculate the estimated length scale
476          Rmix(i,j) = MIN(Rdef(i,j,bi,bj), RRhines(i,j))          Rmix(i,j) = MIN(Rdef(i,j,bi,bj), RRhines(i,j))
477            Rmix(i,j) = MAX(Rmix(i,j),GM_K3D_Rmin)
478    
479  C       Calculate the Doppler shifted long Rossby wave speed  C       Calculate the Doppler shifted long Rossby wave speed
480  C       Ubaro is on the U grid so we must average onto the M grid.  C       Ubaro is on the U grid so we must average onto the M grid.
# Line 500  C     Calculate the diffusivity (on the Line 504  C     Calculate the diffusivity (on the
504               IF (urms(i,j,k).EQ.0.0) THEN               IF (urms(i,j,k).EQ.0.0) THEN
505                 K3D(i,j,k,bi,bj) = GM_K3D_smallK                 K3D(i,j,k,bi,bj) = GM_K3D_smallK
506               ELSE               ELSE
507                 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)
508                 supp(i,j,k) = 1/( 1 + 10*umc(i,j,k)**2/urms(i,j,1)**2 )                 supp(i,j,k)=1/(1+GM_K3D_b1*umc(i,j,k)**2/urms(i,j,1)**2)
509    C              2*Rmix gives the diameter
510                 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)
511       &                            *Rmix(i,j)*supp(i,j,k)       &                            *2*Rmix(i,j)*supp(i,j,k)
512               ENDIF               ENDIF
513    
514  C            Enforce lower and upper bounds on the diffusivity  C            Enforce lower and upper bounds on the diffusivity
# Line 563  C          Calculate the zonal slope at Line 568  C          Calculate the zonal slope at
568             sigz  = MIN( op5*(sigmaR(i,j,k)+sigmaR(i-1,j,k)), maxDRhoDz )             sigz  = MIN( op5*(sigmaR(i,j,k)+sigmaR(i-1,j,k)), maxDRhoDz )
569             sigx  = op5*( sigmaX(i,j,k)+sigmaX(i,j,k-1) )             sigx  = op5*( sigmaX(i,j,k)+sigmaX(i,j,k-1) )
570             slope = sigx/sigz             slope = sigx/sigz
 C           IF(ABS(slope).GT.GM_K3D_maxSlope)  
 C     &          slope = SIGN(GM_K3D_maxSlope,slope)  
571             IF(ABS(slope).GT.GM_maxSlope)             IF(ABS(slope).GT.GM_maxSlope)
572       &          slope = SIGN(GM_maxSlope,slope)       &          slope = SIGN(GM_maxSlope,slope)
573             SlopeX(i,j,k)=-maskW(i,j,k-1,bi,bj)*maskW(i,j,k,bi,bj)*slope             SlopeX(i,j,k)=-maskW(i,j,k-1,bi,bj)*maskW(i,j,k,bi,bj)*slope
# Line 573  C          Calculate the meridional slop Line 576  C          Calculate the meridional slop
576             sigz  = MIN( op5*(sigmaR(i,j,k)+sigmaR(i,j-1,k)), maxDRhoDz )             sigz  = MIN( op5*(sigmaR(i,j,k)+sigmaR(i,j-1,k)), maxDRhoDz )
577             sigy  = op5*( sigmaY(i,j,k) + sigmaY(i,j,k-1) )             sigy  = op5*( sigmaY(i,j,k) + sigmaY(i,j,k-1) )
578             slope = sigy/sigz             slope = sigy/sigz
 C           IF(ABS(slope).GT.GM_K3D_maxSlope)  
 C     &          slope = SIGN(GM_K3D_maxSlope,slope)  
579             IF(ABS(slope).GT.GM_maxSlope)             IF(ABS(slope).GT.GM_maxSlope)
580       &          slope = SIGN(GM_maxSlope,slope)       &          slope = SIGN(GM_maxSlope,slope)
581             SlopeY(i,j,k)=-maskS(i,j,k-1,bi,bj)*maskS(i,j,k,bi,bj)*slope             SlopeY(i,j,k)=-maskS(i,j,k-1,bi,bj)*maskS(i,j,k,bi,bj)*slope
# Line 604  C     We still need to give special trea Line 605  C     We still need to give special trea
605        DO k=Nr-1,1,-1        DO k=Nr-1,1,-1
606         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly-1
607          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx-1
608           IF(k.LE.surfk(i,j) .AND. GM_K3D_PVsheet) THEN           IF(k.LE.surfk(i,j) .AND. GM_K3D_ThickSheet) THEN
609  C          We are in the surface layer, so set the thickness flux  C          We are in the surface layer, so set the thickness flux
610  C          based on the average slope over the surface layer  C          based on the average slope over the surface layer
611  C          If we are on the edge of a "cliff" the surface layer at the  C          If we are on the edge of a "cliff" the surface layer at the
# Line 782  C       Calculate Xi_m at the west face Line 783  C       Calculate Xi_m at the west face
783           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
784            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
785             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
786                XimX(m,i,j) = XimX(m,i,j)              Kdqdx(i,j,k) = KPV(i,j,k)*gradqx(i,j,k)
787       &             - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)              XimX(m,i,j) = XimX(m,i,j)
788       &             *KPV(i,j,k)*gradqx(i,j,k)*modesW(m,i,j,k,bi,bj)       &           - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)
789         &           *Kdqdx(i,j,k)*modesW(m,i,j,k,bi,bj)
790             ENDDO             ENDDO
791            ENDDO            ENDDO
792           ENDDO           ENDDO
# Line 832  C     Calculate the eigenvectors at the Line 834  C     Calculate the eigenvectors at the
834           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
835            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
836             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
837                Kdqdy(i,j,k) = KPV(i,j,k)*gradqy(i,j,k)
838              XimY(m,i,j) = XimY(m,i,j)              XimY(m,i,j) = XimY(m,i,j)
839       &           - drF(k)*hfacS(i,j,k,bi,bj)       &           - drF(k)*hfacS(i,j,k,bi,bj)
840       &           *KPV(i,j,k)*gradqy(i,j,k)*modesS(m,i,j,k,bi,bj)       &           *Kdqdy(i,j,k)*modesS(m,i,j,k,bi,bj)
841             ENDDO             ENDDO
842            ENDDO            ENDDO
843           ENDDO           ENDDO
# Line 945  C     Diagnostics Line 948  C     Diagnostics
948          CALL DIAGNOSTICS_FILL(tfluxY, 'GM_TFLXY',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(tfluxY, 'GM_TFLXY',0,Nr,0,1,1,myThid)
949          CALL DIAGNOSTICS_FILL(gradqx, 'GM_dqdx ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(gradqx, 'GM_dqdx ',0,Nr,0,1,1,myThid)
950          CALL DIAGNOSTICS_FILL(gradqy, 'GM_dqdy ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(gradqy, 'GM_dqdy ',0,Nr,0,1,1,myThid)
951            CALL DIAGNOSTICS_FILL(Kdqdy,  'GM_Kdqdy',0,Nr,0,1,1,myThid)
952            CALL DIAGNOSTICS_FILL(Kdqdx,  'GM_Kdqdx',0,Nr,0,1,1,myThid)
953          CALL DIAGNOSTICS_FILL(surfkz, 'GM_SFLYR',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(surfkz, 'GM_SFLYR',0, 1,0,1,1,myThid)
954          CALL DIAGNOSTICS_FILL(ustar,  'GM_USTAR',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(ustar,  'GM_USTAR',0,Nr,0,1,1,myThid)
955          CALL DIAGNOSTICS_FILL(vstar,  'GM_VSTAR',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(vstar,  'GM_VSTAR',0,Nr,0,1,1,myThid)
# Line 954  C     Diagnostics Line 959  C     Diagnostics
959       &                                'GM_MODEC',0,Nr,0,1,1,myThid)       &                                'GM_MODEC',0,Nr,0,1,1,myThid)
960          CALL DIAGNOSTICS_FILL(M4loc,  'GM_M4   ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(M4loc,  'GM_M4   ',0,Nr,0,1,1,myThid)
961          CALL DIAGNOSTICS_FILL(N2loc,  'GM_N2   ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(N2loc,  'GM_N2   ',0,Nr,0,1,1,myThid)
962            CALL DIAGNOSTICS_FILL(M4onN2, 'GM_M4_N2',0,Nr,0,1,1,myThid)
963          CALL DIAGNOSTICS_FILL(slopeC, 'GM_SLOPE',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(slopeC, 'GM_SLOPE',0,Nr,0,1,1,myThid)
964    
965        ENDIF        ENDIF
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22