/[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.13 by m_bates, Mon Oct 21 18:46:05 2013 UTC revision 1.14 by m_bates, Wed Jan 1 23:20:48 2014 UTC
# Line 584  C          Calculate the meridional slop Line 584  C          Calculate the meridional slop
584         ENDDO         ENDDO
585        ENDDO        ENDDO
586    
587  C     Calculate the thickness flux  C     Calculate the thickness flux and diffusivity.  These may be altered later
588    C     depending on namelist options.
589  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)
590        k=Nr        k=Nr
591        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
# Line 595  C          Zonal thickness flux at the w Line 596  C          Zonal thickness flux at the w
596  C          Meridional thickness flux at the southern cell face  C          Meridional thickness flux at the southern cell face
597             tfluxY(i,j,k) = -fCoriV(i,j)*SlopeY(i,j,k)             tfluxY(i,j,k) = -fCoriV(i,j)*SlopeY(i,j,k)
598       &                     *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)       &                     *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
599    
600    C          Use the interior diffusivity. Note: if we are using a
601    C          constant diffusivity KPV is overwritten later
602               KPV(i,j,k) = K3D(i,j,k,bi,bj)
603    
604         ENDDO         ENDDO
605        ENDDO        ENDDO
606    
607  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.  
608        DO k=Nr-1,1,-1        DO k=Nr-1,1,-1
609         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly
610          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx
611           IF(k.LE.surfk(i,j) .AND. GM_K3D_ThickSheet) THEN  C        Zonal thickness flux at the western cell face
612  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))
613  C          based on the average slope over the surface layer       &        *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)
614  C          If we are on the edge of a "cliff" the surface layer at the       &        *maskW(i,j,k,bi,bj)
615  C          centre of the grid point could be deeper than the U or V point.            
616  C          So, we ensure that we always take a sensible slope.  C        Meridional thickness flux at the southern cell face
617             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))
618               kk=kLow_U(i,j)       &        *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
619               hsurf = -rLowW(i,j,bi,bj)       &        *maskS(i,j,k,bi,bj)
620             ELSE          
621               kk=surfk(i,j)  C        Use the interior diffusivity. Note: if we are using a
622               hsurf = -surfkz(i,j)  C        constant diffusivity KPV is overwritten later
623             ENDIF           KPV(i,j,k) = K3D(i,j,k,bi,bj)
624             IF(kk.GT.0) THEN          ENDDO
625               IF(kk.EQ.Nr) THEN         ENDDO
626                 tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)        ENDDO
627       &              *SlopeX(i,j,kk)/hsurf  
628               ELSE  
629                 tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)  C     Special treatment for the surface layer (if necessary), which overwrites
630       &              *( SlopeX(i,j,kk)-SlopeX(i,j,kk+1) )/hsurf  C     values in the previous loops.
631          IF (GM_K3D_ThickSheet .OR. GM_K3D_surfK) THEN
632            DO k=Nr-1,1,-1
633             DO j=1-Oly,sNy+Oly
634              DO i=1-Olx,sNx+Olx
635               IF(k.LE.surfk(i,j)) THEN
636    C            We are in the surface layer.  Change the thickness flux
637    C            and diffusivity as necessary.
638    
639                 IF (GM_K3D_ThickSheet) THEN
640    C              We are in the surface layer, so set the thickness flux
641    C              based on the average slope over the surface layer
642    C              If we are on the edge of a "cliff" the surface layer at the
643    C              centre of the grid point could be deeper than the U or V point.  
644    C              So, we ensure that we always take a sensible slope.
645                   IF(kLow_U(i,j).LT.surfk(i,j)) THEN
646                     kk=kLow_U(i,j)
647                     hsurf = -rLowW(i,j,bi,bj)
648                   ELSE
649                     kk=surfk(i,j)
650                     hsurf = -surfkz(i,j)
651                   ENDIF
652                   IF(kk.GT.0) THEN
653                     IF(kk.EQ.Nr) THEN
654                       tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)
655         &                  *SlopeX(i,j,kk)/hsurf
656                     ELSE
657                       tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)
658         &                  *( SlopeX(i,j,kk)-SlopeX(i,j,kk+1) )/hsurf
659                     ENDIF
660                   ELSE
661                     tfluxX(i,j,k) = zeroRL
662                   ENDIF
663                  
664                   IF(kLow_V(i,j).LT.surfk(i,j)) THEN
665                     kk=kLow_V(i,j)
666                     hsurf = -rLowS(i,j,bi,bj)
667                   ELSE
668                     kk=surfk(i,j)
669                     hsurf = -surfkz(i,j)
670                   ENDIF
671                   IF(kk.GT.0) THEN
672                     IF(kk.EQ.Nr) THEN
673                       tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)
674         &                  *SlopeY(i,j,kk)/hsurf
675                     ELSE
676                       tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)
677         &                  *( SlopeY(i,j,kk)-SlopeY(i,j,kk+1) )/hsurf
678                     ENDIF
679                   ELSE
680                     tfluxY(i,j,k) = zeroRL
681                   ENDIF
682               ENDIF               ENDIF
            ELSE  
              tfluxX(i,j,k) = zeroRL  
            ENDIF  
683    
684             IF(kLow_V(i,j).LT.surfk(i,j)) THEN               IF (GM_K3D_surfK) THEN
685               kk=kLow_V(i,j)  C              Use a constant K in the surface layer.
686               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  
687               ENDIF               ENDIF
            ELSE  
              tfluxY(i,j,k) = zeroRL  
688             ENDIF             ENDIF
689              ENDDO
690           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  
691          ENDDO          ENDDO
692         ENDDO        ENDIF
       ENDDO  
693    
694  C     Calculate gradq  C     Calculate gradq
695        IF (GM_K3D_likeGM .OR. GM_K3D_beta_eq_0) THEN        IF (GM_K3D_likeGM .OR. GM_K3D_beta_eq_0) THEN
# Line 719  C     It does not affect the end result Line 743  C     It does not affect the end result
743         ENDDO         ENDDO
744        ENDDO        ENDDO
745    
746  C     If GM_K3D_likeGM=.TRUE., K3D becomes a diagnostic field only  C     If GM_K3D_likeGM=.TRUE., the diffusivity for the eddy transport is
747  C     and the diffusivity is set to a constant GM_K3D_constK.  C     set to a constant equal to GM_K3D_constK.
748  C     If the diffusivity is constant the method here is the same as GM.  C     If the diffusivity is constant the method here is the same as GM.
749  C     For the Redi diffusivity K3D is set to GM_K3D_constK at the end  C     If GM_K3D_constRedi=.TRUE. K3D will be set equal to GM_K3D_constK later.
 C     of this routine (after the diagnostics are filled).  
750        IF(GM_K3D_likeGM) THEN        IF(GM_K3D_likeGM) THEN
751          DO k=1,Nr          DO k=1,Nr
752           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
# Line 732  C     of this routine (after the diagnos Line 755  C     of this routine (after the diagnos
755            ENDDO            ENDDO
756           ENDDO           ENDDO
757          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  
758        ENDIF        ENDIF
759    
760        IF (.NOT. GM_K3D_smooth) THEN        IF (.NOT. GM_K3D_smooth) THEN
# Line 931  C     ================================== Line 946  C     ==================================
946  C     Diagnostics  C     Diagnostics
947        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
948          CALL DIAGNOSTICS_FILL(K3D,    'GM_K3D  ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(K3D,    'GM_K3D  ',0,Nr,0,1,1,myThid)
949            CALL DIAGNOSTICS_FILL(KPV,    'GM_KPV  ',0,Nr,0,1,1,myThid)
950          CALL DIAGNOSTICS_FILL(urms,   'GM_URMS ',0,Nr,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(urms,   'GM_URMS ',0,Nr,0,1,1,myThid)
951          CALL DIAGNOSTICS_FILL(Rdef,   'GM_RDEF ',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Rdef,   'GM_RDEF ',0, 1,0,1,1,myThid)
952          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 966  C     Diagnostics Line 982  C     Diagnostics
982  #endif  #endif
983    
984  C     For the Redi diffusivity, we set K3D to a constant if  C     For the Redi diffusivity, we set K3D to a constant if
985  C     GM_K3D_likeGM=.TRUE. (see earlier comments)  C     GM_K3D_constRedi=.TRUE.
986        IF (GM_K3D_likeGM) THEN        IF (GM_K3D_constRedi) THEN
987          DO k=1,Nr          DO k=1,Nr
988           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
989            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
# Line 977  C     GM_K3D_likeGM=.TRUE. (see earlier Line 993  C     GM_K3D_likeGM=.TRUE. (see earlier
993          ENDDO          ENDDO
994        ENDIF        ENDIF
995    
996    #ifdef ALLOW_DIAGNOSTICS
997          IF ( useDiagnostics )
998         &     CALL DIAGNOSTICS_FILL(K3D, 'GM_K3D_T',0,Nr,0,1,1,myThid)
999    #endif
1000    
1001  #endif /* GM_K3D */  #endif /* GM_K3D */
1002        RETURN        RETURN
1003        END        END
Legend:
Removed from v.1.13  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22