/[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.16 by m_bates, Fri Mar 28 04:22:11 2014 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 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 117  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 131  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 144  C     surfkz   :: Depth of surface layer Line 150  C     surfkz   :: Depth of surface layer
150        _RL fCoriV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL fCoriV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
151        _RL surfkz(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL surfkz(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
152    
153    C     centreX,centreY :: used for calculating averages at centre of cell
154    C     numerator,denominator :: of the renormalisation factor
155    C     uInt      :: column integral of u velocity (sum u*dz)
156    C     vInt      :: column integral of v velocity (sum v*dz)
157    C     KdqdxInt  :: column integral of K*dqdx (sum K*dqdx*dz)
158    C     KdqdyInt  :: column integral of K*dqdy (sum K*dqdy*dz)
159    C     uKdqdyInt :: column integral of u*K*dqdy (sum u*K*dqdy*dz)
160    C     vKdqdxInt :: column integral of v*K*dqdx (sum v*K*dqdx*dz)
161    C     uXiyInt   :: column integral of u*Xiy (sum u*Xiy*dz)
162    C     vXixInt   :: column integral of v*Xix (sum v*Xix*dz)
163    C     Renorm    :: renormalisation factor at the centre of a cell
164    C     RenormU   :: renormalisation factor at the western face of a cell
165    C     RenormV   :: renormalisation factor at the southern face of a cell
166          _RL centreX, centreY
167          _RL numerator, denominator
168          _RL uInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
169          _RL vInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
170          _RL KdqdxInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
171          _RL KdqdyInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
172          _RL uKdqdyInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
173          _RL vKdqdxInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
174          _RL uXiyInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
175          _RL vXixInt(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
176          _RL Renorm(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
177          _RL RenormU(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
178          _RL RenormV(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
179    
180  C     gradqx   :: Potential vorticity gradient in x direction  C     gradqx   :: Potential vorticity gradient in x direction
181  C     gradqy   :: Potential vorticity gradient in y direction  C     gradqy   :: Potential vorticity gradient in y direction
182  C     XimX     :: Vertical integral of phi_m*K*gradqx  C     XimX     :: Vertical integral of phi_m*K*gradqx
# Line 258  C       Not limited Line 291  C       Not limited
291  C     Zeroing some cumulative fields  C     Zeroing some cumulative fields
292        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
293         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
294          eady(i,j)  = zeroRL          eady(i,j)   = zeroRL
295          BVint(i,j) = zeroRL          BVint(i,j)  = zeroRL
296          Ubaro(i,j) = zeroRL          Ubaro(i,j)  = zeroRL
297            deltaH(i,j) = zeroRL
298           ENDDO
299          ENDDO
300          DO k=1,Nr
301           DO j=1-Oly,sNy+Oly
302            DO i=1-Olx,sNx+Olx
303             slopeC(i,j,k)=zeroRL
304            ENDDO
305         ENDDO         ENDDO
306        ENDDO        ENDDO
307    
# Line 382  C     Calculate the Eady growth rate Line 423  C     Calculate the Eady growth rate
423  C     ===============================  C     ===============================
424        DO k=1,Nr        DO k=1,Nr
425    
426           DO j=1-Oly,sNy+Oly-1
427            DO i=1-Olx,sNx+Olx-1
428             M4loc(i,j,k) = g_reciprho_sq*( dSigmaDx(i,j,k)**2
429         &                                 +dSigmaDy(i,j,k)**2 )
430             IF (k.NE.kLow_C(i,j)) THEN
431               N2loc(i,j,k) = op5*(N2(i,j,k)+N2(i,j,k+1))
432             ELSE
433               N2loc(i,j,k) = op5*N2(i,j,k)
434             ENDIF
435            ENDDO
436           ENDDO
437  C      The bottom of the grid cell is shallower than the top  C      The bottom of the grid cell is shallower than the top
438  C      integration level, so, advance the depth.  C      integration level, so, advance the depth.
439         IF (-rF(k+1).LE. GM_K3D_EadyMinDepth) CYCLE         IF (-rF(k+1) .LE. GM_K3D_EadyMinDepth) CYCLE
440    
441  C      Do not bother going any deeper since the top of the  C      Do not bother going any deeper since the top of the
442  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 445  C      cell is deeper than the bottom in
445  C      We are in the integration depth range  C      We are in the integration depth range
446         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly-1
447          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx-1
448           IF (kLow_C(i,j).GE.k) THEN           IF ( (kLow_C(i,j).GE.k) .AND.
449             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)  
450    
451               slopeC(i,j,k) = SQRT(M4loc(i,j,k))/N2loc(i,j,k)
452  C          Limit the slope.  Note, this is not all the Eady calculations.  C          Limit the slope.  Note, this is not all the Eady calculations.
453             IF (slope.LE.GM_K3D_maxSlope) THEN             IF (slopeC(i,j,k).LE.GM_maxSlope) THEN
454               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/(N2loc)  
455             ELSE             ELSE
456               eady(i,j) = eady(i,j)               slopeC(i,j,k) = GM_maxslope
457       &            + hfacC(i,j,k,bi,bj)*drF(k)*SQRT(M4loc)               M4onN2(i,j,k) = SQRT(M4loc(i,j,k))*GM_maxslope
      &              *GM_K3D_maxSlope*GM_K3D_maxSlope  
458             ENDIF             ENDIF
459               eady(i,j)   = eady(i,j)  
460         &                   + hfacC(i,j,k,bi,bj)*drF(k)*M4onN2(i,j,k)
461               deltaH(i,j) = deltaH(i,j) + drF(k)
462           ENDIF           ENDIF
463          ENDDO          ENDDO
464         ENDDO         ENDDO
# Line 419  C          Limit the slope.  Note, this Line 466  C          Limit the slope.  Note, this
466    
467        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
468         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
469  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
470  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
471  C       to avoid floating point exceptions.  These points are taken care  C       integration, we set the Eady growth rate to something small
472  C       of by setting K3D=GM_K3D_smallK later.  C       to avoid floating point exceptions.
473          IF (-r_Low(i,j,bi,bj).LE.GM_K3D_EadyMinDepth) THEN  C       Later, these areas will be given a small diffusivity.
474            IF (deltaH(i,j).EQ.zeroRL) THEN
475            eady(i,j) = small            eady(i,j) = small
476    
477  C       Otherwise, multiply eady by the various constants to get the  C       Otherwise, divide over the integration and take the square root
478  C       growth rate.  C       to actually find the Eady growth rate.
479          ELSE          ELSE
480            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)  
481                        
482          ENDIF          ENDIF
483    
# Line 444  C     ================================== Line 490  C     ==================================
490        DO j=1-Oly+1,sNy+Oly        DO j=1-Oly+1,sNy+Oly
491         DO i=1-Olx+1,sNx+Olx-1         DO i=1-Olx+1,sNx+Olx-1
492  C       Calculate the Visbeck velocity  C       Calculate the Visbeck velocity
493          Rurms(i,j) = MIN(Rdef(i,j,bi,bj),GM_K3D_maxLurms)          Rurms(i,j) = MIN(Rdef(i,j,bi,bj),GM_K3D_Rmax)
494          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)
495  C       Set the bottom urms to zero  C       Set the bottom urms to zero
496          k=kLow_C(i,j)          k=kLow_C(i,j)
# Line 455  C       Calculate the Rhines scale Line 501  C       Calculate the Rhines scale
501    
502  C       Calculate the estimated length scale  C       Calculate the estimated length scale
503          Rmix(i,j) = MIN(Rdef(i,j,bi,bj), RRhines(i,j))          Rmix(i,j) = MIN(Rdef(i,j,bi,bj), RRhines(i,j))
504            Rmix(i,j) = MAX(Rmix(i,j),GM_K3D_Rmin)
505    
506  C       Calculate the Doppler shifted long Rossby wave speed  C       Calculate the Doppler shifted long Rossby wave speed
507  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 478  C     Calculate the diffusivity (on the Line 525  C     Calculate the diffusivity (on the
525         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
526          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
527           IF (k.LE.kLow_C(i,j)) THEN           IF (k.LE.kLow_C(i,j)) THEN
528             IF (-r_Low(i,j,bi,bj).LE.GM_K3D_EadyMinDepth) THEN             IF (deltaH(i,j).EQ.zeroRL) THEN
529               K3D(i,j,k,bi,bj) = GM_K3D_smallK               K3D(i,j,k,bi,bj) = GM_K3D_smallK
530             ELSE             ELSE
531               IF (urms(i,j,k).EQ.0.0) THEN               IF (urms(i,j,k).EQ.0.0) THEN
532                 K3D(i,j,k,bi,bj) = GM_K3D_smallK                 K3D(i,j,k,bi,bj) = GM_K3D_smallK
533               ELSE               ELSE
534                 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)
535                 supp(i,j,k) = 1/( 1 + 4*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)
536    C              2*Rmix gives the diameter
537                 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)
538       &                            *Rmix(i,j)*supp(i,j,k)       &                            *2*Rmix(i,j)*supp(i,j,k)
539               ENDIF               ENDIF
540    
541  C            Enforce lower and upper bounds on the diffusivity  C            Enforce lower and upper bounds on the diffusivity
# Line 547  C          Calculate the zonal slope at Line 595  C          Calculate the zonal slope at
595             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 )
596             sigx  = op5*( sigmaX(i,j,k)+sigmaX(i,j,k-1) )             sigx  = op5*( sigmaX(i,j,k)+sigmaX(i,j,k-1) )
597             slope = sigx/sigz             slope = sigx/sigz
 C           IF(ABS(slope).GT.GM_K3D_maxSlope)  
 C     &          slope = SIGN(GM_K3D_maxSlope,slope)  
598             IF(ABS(slope).GT.GM_maxSlope)             IF(ABS(slope).GT.GM_maxSlope)
599       &          slope = SIGN(GM_maxSlope,slope)       &          slope = SIGN(GM_maxSlope,slope)
600             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 557  C          Calculate the meridional slop Line 603  C          Calculate the meridional slop
603             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 )
604             sigy  = op5*( sigmaY(i,j,k) + sigmaY(i,j,k-1) )             sigy  = op5*( sigmaY(i,j,k) + sigmaY(i,j,k-1) )
605             slope = sigy/sigz             slope = sigy/sigz
 C           IF(ABS(slope).GT.GM_K3D_maxSlope)  
 C     &          slope = SIGN(GM_K3D_maxSlope,slope)  
606             IF(ABS(slope).GT.GM_maxSlope)             IF(ABS(slope).GT.GM_maxSlope)
607       &          slope = SIGN(GM_maxSlope,slope)       &          slope = SIGN(GM_maxSlope,slope)
608             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 567  C     &          slope = SIGN(GM_K3D_max Line 611  C     &          slope = SIGN(GM_K3D_max
611         ENDDO         ENDDO
612        ENDDO        ENDDO
613    
614  C     Calculate the thickness flux  C     Calculate the thickness flux and diffusivity.  These may be altered later
615    C     depending on namelist options.
616  C     Enforce a zero slope bottom boundary condition for the bottom most cells (k=Nr)  C     Enforce a zero slope bottom boundary condition for the bottom most cells (k=Nr)
617        k=Nr        k=Nr
618        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
# Line 578  C          Zonal thickness flux at the w Line 623  C          Zonal thickness flux at the w
623  C          Meridional thickness flux at the southern cell face  C          Meridional thickness flux at the southern cell face
624             tfluxY(i,j,k) = -fCoriV(i,j)*SlopeY(i,j,k)             tfluxY(i,j,k) = -fCoriV(i,j)*SlopeY(i,j,k)
625       &                     *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)       &                     *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
626    
627    C          Use the interior diffusivity. Note: if we are using a
628    C          constant diffusivity KPV is overwritten later
629               KPV(i,j,k) = K3D(i,j,k,bi,bj)
630    
631         ENDDO         ENDDO
632        ENDDO        ENDDO
633    
634  C     Calculate the thickness flux for other cells (k<Nr)  C     Calculate the thickness flux and diffusivity and for other cells (k<Nr)
 C     SlopeX and SlopeY are zero in dry cells, so, the bottom boundary  
 C     condition that the slope is zero is taken care of.  
 C     We still need to give special treatment for the surface layer however.  
635        DO k=Nr-1,1,-1        DO k=Nr-1,1,-1
636         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly
637          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx
638           IF(k.LE.surfk(i,j) .AND. .NOT. GM_K3D_likeGM) THEN  C        Zonal thickness flux at the western cell face
639  C          We are in the surface layer, so set the thickness flux           tfluxX(i,j,k)=-fCoriU(i,j)*(SlopeX(i,j,k)-SlopeX(i,j,k+1))
640  C          based on the average slope over the surface layer       &        *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)
641  C          If we are on the edge of a "cliff" the surface layer at the       &        *maskW(i,j,k,bi,bj)
642  C          centre of the grid point could be deeper than the U or V point.            
643  C          So, we ensure that we always take a sensible slope.  C        Meridional thickness flux at the southern cell face
644             IF(kLow_U(i,j).LT.surfk(i,j)) THEN           tfluxY(i,j,k)=-fCoriV(i,j)*(SlopeY(i,j,k)-SlopeY(i,j,k+1))
645               kk=kLow_U(i,j)       &        *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
646               hsurf = -rLowW(i,j,bi,bj)       &        *maskS(i,j,k,bi,bj)
647             ELSE          
648               kk=surfk(i,j)  C        Use the interior diffusivity. Note: if we are using a
649               hsurf = -surfkz(i,j)  C        constant diffusivity KPV is overwritten later
650             ENDIF           KPV(i,j,k) = K3D(i,j,k,bi,bj)
651             IF(kk.GT.0) THEN          ENDDO
652               IF(kk.EQ.Nr) THEN         ENDDO
653                 tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)        ENDDO
654       &              *SlopeX(i,j,kk)/hsurf  
655               ELSE  
656                 tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)  C     Special treatment for the surface layer (if necessary), which overwrites
657       &              *( SlopeX(i,j,kk)-SlopeX(i,j,kk+1) )/hsurf  C     values in the previous loops.
658          IF (GM_K3D_ThickSheet .OR. GM_K3D_surfK) THEN
659            DO k=Nr-1,1,-1
660             DO j=1-Oly,sNy+Oly
661              DO i=1-Olx,sNx+Olx
662               IF(k.LE.surfk(i,j)) THEN
663    C            We are in the surface layer.  Change the thickness flux
664    C            and diffusivity as necessary.
665    
666                 IF (GM_K3D_ThickSheet) THEN
667    C              We are in the surface layer, so set the thickness flux
668    C              based on the average slope over the surface layer
669    C              If we are on the edge of a "cliff" the surface layer at the
670    C              centre of the grid point could be deeper than the U or V point.  
671    C              So, we ensure that we always take a sensible slope.
672                   IF(kLow_U(i,j).LT.surfk(i,j)) THEN
673                     kk=kLow_U(i,j)
674                     hsurf = -rLowW(i,j,bi,bj)
675                   ELSE
676                     kk=surfk(i,j)
677                     hsurf = -surfkz(i,j)
678                   ENDIF
679                   IF(kk.GT.0) THEN
680                     IF(kk.EQ.Nr) THEN
681                       tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)
682         &                  *SlopeX(i,j,kk)/hsurf
683                     ELSE
684                       tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)
685         &                  *( SlopeX(i,j,kk)-SlopeX(i,j,kk+1) )/hsurf
686                     ENDIF
687                   ELSE
688                     tfluxX(i,j,k) = zeroRL
689                   ENDIF
690                  
691                   IF(kLow_V(i,j).LT.surfk(i,j)) THEN
692                     kk=kLow_V(i,j)
693                     hsurf = -rLowS(i,j,bi,bj)
694                   ELSE
695                     kk=surfk(i,j)
696                     hsurf = -surfkz(i,j)
697                   ENDIF
698                   IF(kk.GT.0) THEN
699                     IF(kk.EQ.Nr) THEN
700                       tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)
701         &                  *SlopeY(i,j,kk)/hsurf
702                     ELSE
703                       tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)
704         &                  *( SlopeY(i,j,kk)-SlopeY(i,j,kk+1) )/hsurf
705                     ENDIF
706                   ELSE
707                     tfluxY(i,j,k) = zeroRL
708                   ENDIF
709               ENDIF               ENDIF
            ELSE  
              tfluxX(i,j,k) = zeroRL  
            ENDIF  
710    
711             IF(kLow_V(i,j).LT.surfk(i,j)) THEN               IF (GM_K3D_surfK) THEN
712               kk=kLow_V(i,j)  C              Use a constant K in the surface layer.
713               hsurf = -rLowS(i,j,bi,bj)                 KPV(i,j,k) = GM_K3D_constK
            ELSE  
              kk=surfk(i,j)  
              hsurf = -surfkz(i,j)  
            ENDIF  
            IF(kk.GT.0) THEN  
              IF(kk.EQ.Nr) THEN  
                tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)  
      &              *SlopeY(i,j,kk)/hsurf  
              ELSE  
                tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)  
      &              *( SlopeY(i,j,kk)-SlopeY(i,j,kk+1) )/hsurf  
714               ENDIF               ENDIF
            ELSE  
              tfluxY(i,j,k) = zeroRL  
715             ENDIF             ENDIF
716              ENDDO
717           ELSE           ENDDO
 C          We are not in the surface layer, so calculate the thickness  
 C          flux based on the local isopyncal slope  
   
 C          Zonal thickness flux at the western cell face  
            tfluxX(i,j,k)=-fCoriU(i,j)*(SlopeX(i,j,k)-SlopeX(i,j,k+1))  
      &                    *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)  
      &                    *maskW(i,j,k,bi,bj)  
   
 C          Meridional thickness flux at the southern cell face  
            tfluxY(i,j,k)=-fCoriV(i,j)*(SlopeY(i,j,k)-SlopeY(i,j,k+1))  
      &                    *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)  
      &                    *maskS(i,j,k,bi,bj)  
          ENDIF  
718          ENDDO          ENDDO
719         ENDDO        ENDIF
       ENDDO  
720    
721  C     Calculate gradq  C     Calculate gradq
722        IF (GM_K3D_likeGM .OR. GM_K3D_beta_eq_0) THEN        IF (GM_K3D_likeGM .OR. GM_K3D_beta_eq_0) THEN
# Line 702  C     It does not affect the end result Line 770  C     It does not affect the end result
770         ENDDO         ENDDO
771        ENDDO        ENDDO
772    
773    C     If GM_K3D_likeGM=.TRUE., the diffusivity for the eddy transport is
774    C     set to a constant equal to GM_K3D_constK.
775    C     If the diffusivity is constant the method here is the same as GM.
776    C     If GM_K3D_constRedi=.TRUE. K3D will be set equal to GM_K3D_constK later.
777        IF(GM_K3D_likeGM) THEN        IF(GM_K3D_likeGM) THEN
778          DO k=1,Nr          DO k=1,Nr
779           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
# Line 710  C     It does not affect the end result Line 782  C     It does not affect the end result
782            ENDDO            ENDDO
783           ENDDO           ENDDO
784          ENDDO          ENDDO
       ELSE  
         DO k=1,Nr  
          DO j=1-Oly,sNy+Oly  
           DO i=1-Olx,sNx+Olx  
            KPV(i,j,k) = K3D(i,j,k,bi,bj)  
           ENDDO  
          ENDDO  
         ENDDO  
785        ENDIF        ENDIF
786    
787        IF (.NOT. GM_K3D_smooth) THEN        IF (.NOT. GM_K3D_smooth) THEN
# Line 761  C       Calculate Xi_m at the west face Line 825  C       Calculate Xi_m at the west face
825           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
826            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
827             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
828                XimX(m,i,j) = XimX(m,i,j)              Kdqdx(i,j,k) = KPV(i,j,k)*gradqx(i,j,k)
829       &             - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)              XimX(m,i,j) = XimX(m,i,j)
830       &             *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)
831         &           *Kdqdx(i,j,k)*modesW(m,i,j,k,bi,bj)
832             ENDDO             ENDDO
833            ENDDO            ENDDO
834           ENDDO           ENDDO
# Line 811  C     Calculate the eigenvectors at the Line 876  C     Calculate the eigenvectors at the
876           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
877            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
878             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
879                Kdqdy(i,j,k) = KPV(i,j,k)*gradqy(i,j,k)
880              XimY(m,i,j) = XimY(m,i,j)              XimY(m,i,j) = XimY(m,i,j)
881       &           - drF(k)*hfacS(i,j,k,bi,bj)       &           - drF(k)*hfacS(i,j,k,bi,bj)
882       &           *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)
883             ENDDO             ENDDO
884            ENDDO            ENDDO
885           ENDDO           ENDDO
# Line 838  C     Calculate Xi for Y direction at th Line 904  C     Calculate Xi for Y direction at th
904           ENDDO           ENDDO
905          ENDDO          ENDDO
906    
907  C     ENDIF GM_K3D_likeGM  C     ENDIF (.NOT. GM_K3D_smooth)
908        ENDIF        ENDIF
909    
910    
911    C     Calculate the renormalisation factor
912          DO j=1-Oly,sNy+Oly
913           DO i=1-Olx,sNx+Olx
914            uInt(i,j)=zeroRL
915            vInt(i,j)=zeroRL
916            KdqdyInt(i,j)=zeroRL
917            KdqdxInt(i,j)=zeroRL
918            uKdqdyInt(i,j)=zeroRL
919            vKdqdxInt(i,j)=zeroRL
920            uXiyInt(i,j)=zeroRL
921            vXixInt(i,j)=zeroRL
922            Renorm(i,j)=oneRL
923            RenormU(i,j)=oneRL
924            RenormV(i,j)=oneRL
925           ENDDO
926          ENDDO
927          DO k=1,Nr
928           DO j=1-Oly,sNy+Oly-1
929            DO i=1-Olx,sNx+Olx-1
930             centreX = op5*(uVel(i,j,k,bi,bj)+uVel(i+1,j,k,bi,bj))
931             centreY = op5*(Kdqdy(i,j,k)     +Kdqdy(i,j+1,k)     )
932    C        For the numerator
933             uInt(i,j) = uInt(i,j)
934         &        + centreX*hfacC(i,j,k,bi,bj)*drF(k)
935             KdqdyInt(i,j) = KdqdyInt(i,j)
936         &        + centreY*hfacC(i,j,k,bi,bj)*drF(k)
937             uKdqdyInt(i,j) = uKdqdyInt(i,j)
938         &        + centreX*centreY*hfacC(i,j,k,bi,bj)*drF(k)
939    C        For the denominator
940             centreY = op5*(Xiy(i,j,k) + Xiy(i,j+1,k))
941             uXiyInt(i,j) = uXiyInt(i,j)
942         &        + centreX*centreY*hfacC(i,j,k,bi,bj)*drF(k)
943    
944             centreX = op5*(Kdqdx(i,j,k)     +Kdqdx(i+1,j,k))
945             centreY = op5*(vVel(i,j,k,bi,bj)+vVel(i,j+1,k,bi,bj) )
946    C        For the numerator
947             vInt(i,j) = vInt(i,j)
948         &        + centreY*hfacC(i,j,k,bi,bj)*drF(k)
949             KdqdxInt(i,j) = KdqdxInt(i,j)
950         &        + CentreX*hfacC(i,j,k,bi,bj)*drF(k)
951             vKdqdxInt(i,j) = vKdqdxInt(i,j)
952         &        + centreY*centreX*hfacC(i,j,k,bi,bj)*drF(k)
953    C        For the denominator
954             centreX = op5*(Xix(i,j,k) + Xix(i+1,j,k))
955             vXixInt(i,j) = vXixInt(i,j)
956         &        + centreY*centreX*hfacC(i,j,k,bi,bj)*drF(k)
957    
958            ENDDO
959           ENDDO
960          ENDDO
961    
962          DO j=1-Oly,sNy+Oly-1
963           DO i=1-Olx,sNx+Olx-1
964            IF (kLowC(i,j,bi,bj).GT.0) THEN
965              numerator =
966         &         (uKdqdyInt(i,j)-uInt(i,j)*KdqdyInt(i,j)/R_low(i,j,bi,bj))
967         &        -(vKdqdxInt(i,j)-vInt(i,j)*KdqdxInt(i,j)/R_low(i,j,bi,bj))
968              denominator = uXiyInt(i,j) - vXixInt(i,j)
969    C         We can have troubles with floating point exceptions if the denominator
970    C         of the renormalisation if the ocean is resting (e.g. intial conditions).
971    C         So we make the renormalisation factor one if the denominator is very small
972    C         The renormalisation factor is supposed to correct the error in the extraction of
973    C         potential energy associated with the truncation of the expansion. Thus, we
974    C         enforce a minimum value for the renormalisation factor.
975    C         We also enforce a maximum renormalisation factor.
976              IF (denominator.GT.small) THEN
977                Renorm(i,j) = ABS(numerator/denominator)
978                Renorm(i,j) = MAX(Renorm(i,j),GM_K3D_minRenorm)
979                Renorm(i,j) = MIN(Renorm(i,j),GM_K3D_maxRenorm)
980              ENDIF
981            ENDIF
982           ENDDO
983          ENDDO
984    C     Now put it back on to the velocity grids
985          DO j=1-Oly+1,sNy+Oly-1
986           DO i=1-Olx+1,sNx+Olx-1
987            RenormU(i,j) = op5*(Renorm(i-1,j)+Renorm(i,j))
988            RenormV(i,j) = op5*(Renorm(i,j-1)+Renorm(i,j))
989           ENDDO
990          ENDDO
991    
992  C     Calculate the eddy induced velocity in the X direction at the west face  C     Calculate the eddy induced velocity in the X direction at the west face
993        DO k=1,Nr        DO k=1,Nr
994         DO j=1-Oly+1,sNy+Oly         DO j=1-Oly+1,sNy+Oly
995          DO i=1-Olx+1,sNx+Olx          DO i=1-Olx+1,sNx+Olx
996           ustar(i,j,k) = -Xix(i,j,k)/coriU(i,j)           ustar(i,j,k) = -RenormU(i,j)*Xix(i,j,k)/coriU(i,j)
997          ENDDO          ENDDO
998         ENDDO         ENDDO
999        ENDDO              ENDDO      
# Line 855  C     Calculate the eddy induced velocit Line 1002  C     Calculate the eddy induced velocit
1002        DO k=1,Nr        DO k=1,Nr
1003         DO j=1-Oly+1,sNy+Oly         DO j=1-Oly+1,sNy+Oly
1004          DO i=1-Olx+1,sNx+Olx          DO i=1-Olx+1,sNx+Olx
1005           vstar(i,j,k) = -Xiy(i,j,k)/coriV(i,j)           vstar(i,j,k) = -RenormV(i,j)*Xiy(i,j,k)/coriV(i,j)
1006          ENDDO          ENDDO
1007         ENDDO         ENDDO
1008        ENDDO              ENDDO      
# Line 907  C     ================================== Line 1054  C     ==================================
1054  C     Diagnostics  C     Diagnostics
1055        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
1056          CALL DIAGNOSTICS_FILL(K3D,    'GM_K3D  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(K3D,    'GM_K3D  ',0,Nr,0,1,1,myThid)
1057            CALL DIAGNOSTICS_FILL(KPV,    'GM_KPV  ',0,Nr,0,1,1,myThid)
1058          CALL DIAGNOSTICS_FILL(urms,   'GM_URMS ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(urms,   'GM_URMS ',0,Nr,0,1,1,myThid)
1059          CALL DIAGNOSTICS_FILL(Rdef,   'GM_RDEF ',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Rdef,   'GM_RDEF ',0, 1,0,1,1,myThid)
1060          CALL DIAGNOSTICS_FILL(Rurms,  'GM_RURMS',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Rurms,  'GM_RURMS',0, 1,0,1,1,myThid)
# Line 924  C     Diagnostics Line 1072  C     Diagnostics
1072          CALL DIAGNOSTICS_FILL(tfluxY, 'GM_TFLXY',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(tfluxY, 'GM_TFLXY',0,Nr,0,1,1,myThid)
1073          CALL DIAGNOSTICS_FILL(gradqx, 'GM_dqdx ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(gradqx, 'GM_dqdx ',0,Nr,0,1,1,myThid)
1074          CALL DIAGNOSTICS_FILL(gradqy, 'GM_dqdy ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(gradqy, 'GM_dqdy ',0,Nr,0,1,1,myThid)
1075            CALL DIAGNOSTICS_FILL(Kdqdy,  'GM_Kdqdy',0,Nr,0,1,1,myThid)
1076            CALL DIAGNOSTICS_FILL(Kdqdx,  'GM_Kdqdx',0,Nr,0,1,1,myThid)
1077          CALL DIAGNOSTICS_FILL(surfkz, 'GM_SFLYR',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(surfkz, 'GM_SFLYR',0, 1,0,1,1,myThid)
1078          CALL DIAGNOSTICS_FILL(ustar,  'GM_USTAR',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(ustar,  'GM_USTAR',0,Nr,0,1,1,myThid)
1079          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 931  C     Diagnostics Line 1081  C     Diagnostics
1081          CALL DIAGNOSTICS_FILL(ubar,   'GM_UBAR ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(ubar,   'GM_UBAR ',0,Nr,0,1,1,myThid)
1082          CALL DIAGNOSTICS_FILL(modesC(1,:,:,:,bi,bj),          CALL DIAGNOSTICS_FILL(modesC(1,:,:,:,bi,bj),
1083       &                                'GM_MODEC',0,Nr,0,1,1,myThid)       &                                'GM_MODEC',0,Nr,0,1,1,myThid)
1084            CALL DIAGNOSTICS_FILL(M4loc,  'GM_M4   ',0,Nr,0,1,1,myThid)
1085            CALL DIAGNOSTICS_FILL(N2loc,  'GM_N2   ',0,Nr,0,1,1,myThid)
1086            CALL DIAGNOSTICS_FILL(M4onN2, 'GM_M4_N2',0,Nr,0,1,1,myThid)
1087            CALL DIAGNOSTICS_FILL(slopeC, 'GM_SLOPE',0,Nr,0,1,1,myThid)
1088            CALL DIAGNOSTICS_FILL(Renorm, 'GM_RENRM',0, 1,0,1,1,myThid)
1089    
1090        ENDIF        ENDIF
1091  #endif  #endif
1092    
1093    C     For the Redi diffusivity, we set K3D to a constant if
1094    C     GM_K3D_constRedi=.TRUE.
1095          IF (GM_K3D_constRedi) THEN
1096            DO k=1,Nr
1097             DO j=1-Oly,sNy+Oly
1098              DO i=1-Olx,sNx+Olx
1099               K3D(i,j,k,bi,bj) = GM_K3D_constK
1100              ENDDO
1101             ENDDO
1102            ENDDO
1103          ENDIF
1104    
1105    #ifdef ALLOW_DIAGNOSTICS
1106          IF ( useDiagnostics )
1107         &     CALL DIAGNOSTICS_FILL(K3D, 'GM_K3D_T',0,Nr,0,1,1,myThid)
1108    #endif
1109    
1110  #endif /* GM_K3D */  #endif /* GM_K3D */
1111        RETURN        RETURN
1112        END        END
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22