/[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.2 by m_bates, Fri Jun 21 21:56:18 2013 UTC revision 1.3 by m_bates, Thu Jun 27 14:51:40 2013 UTC
# Line 70  C      _RL dSigmaLoX,dSigmaLoY Line 70  C      _RL dSigmaLoX,dSigmaLoY
70        _RL sigx, sigy, sigz        _RL sigx, sigy, sigz
71        _RL hsurf        _RL hsurf
72        _RL small        _RL small
73          _RL KPV(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
74        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL SlopeX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
75        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL SlopeY(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
76        _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL dSigmaDr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
# Line 445  C            Enforce lower and upper bou Line 446  C            Enforce lower and upper bou
446  C     ======================================  C     ======================================
447  C     Find the PV gradient  C     Find the PV gradient
448  C     ======================================  C     ======================================
449  C     Find the mixed layer depth as per Large et al. (1997)  C     Calculate the surface layer thickness.
450  C     Start with the surface density  C     Use hMixLayer (calculated in model/src/calc_oce_mxlayer)
451        CALL FIND_RHO_2D(  C     for the mixed layer depth.
      I                 iMin, iMax, jMin, jMax, 1,  
      I                 theta(1-OLx,1-OLy,1,bi,bj),  
      I                 salt (1-OLx,1-OLy,1,bi,bj),  
      O                 rhosurf,  
      I                 1, bi, bj, myThid )  
452    
453  C     Set the minimum surface layer depth to GM_K3D_surfMinDepth  C     Enforce a minimum surface layer depth
       DO k=1,Nr  
        IF (rF(k+1).LE.-GM_K3D_surfMinDepth) THEN  
          minsurfk=k  
          EXIT  
        ENDIF  
       ENDDO  
454        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
455         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
456          surfk(i,j) = min( minsurfk, kLow_C(i,j) )          surfkz(i,j) = MIN(-GM_K3D_surfMinDepth,-hMixLayer(i,j,bi,bj))
457          surfkz(i,j) = rF(minsurfk+1)          surfkz(i,j) = MAX(surfkz(i,j),R_low(i,j,bi,bj))
458          Dbz(i,j) = zeroRL          IF(maskC(i,j,1,bi,bj).EQ.0.0) surfkz(i,j)=0.0
459            surfk(i,j) = 0
460         ENDDO         ENDDO
461        ENDDO        ENDDO
462          DO k=0,Nr
       DO k=2,Nr  
 C     Calculate the surface layer depth.  
 C     If the mixed layer is deeper than GM_K3D_surfMinDepth then  
 C     set the surface layer depth to the mixed layer depth  
        CALL FIND_RHO_2D(  
      I                  iMin, iMax, jMin, jMax, 1,  
      I                  theta(1-OLx,1-OLy,k,bi,bj),  
      I                  salt (1-OLx,1-OLy,k,bi,bj),  
      O                  rhodeep,  
      I                  k, bi, bj, myThid )  
463         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
464          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
465           IF (maskC(i,j,k,bi,bj).NE.zeroRL) THEN           IF (rF(k).GT.surfkz(i,j) .AND. surfkz(i,j).GE.rF(k+1))
466             newDbz = ( rhosurf(i,j)-rhodeep(i,j) )/rC(k)       &        surfk(i,j) = k
            IF (newDbz.GT.Dbz(i,j)) THEN  
              Dbz(i,j) = newDbz  
              IF (surfk(i,j).LT.k) THEN  
                surfk(i,j) = k  
                surfkz(i,j) = rF(k+1)  
              ENDIF  
            ENDIF  
          ENDIF  
467          ENDDO          ENDDO
468         ENDDO         ENDDO
469        ENDDO        ENDDO
470          
471  C     Calculate the isopycnal slope  C     Calculate the isopycnal slope
472        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
473         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
# Line 536  C     &          slope = SIGN(GM_K3D_max Line 509  C     &          slope = SIGN(GM_K3D_max
509         ENDDO         ENDDO
510        ENDDO        ENDDO
511    
512    CC     Try something different
513    C      DO k=1,Nr
514    CC        Stuff for slope limiting
515    C       DO j=1-Oly+1,sNy+Oly
516    C        DO i=1-Olx+1,sNx+Olx
517    C         slpX(i,j) = sigmaX(i,j,k)
518    C         slpY(i,j) = sigmaY(i,j,k)
519    C         dSigmaDrW(i,j)=op5*( sigmaR(i-1,j,k)+sigmaR(i,j,k) )
520    C     &        *maskW(i,j,k,bi,bj)
521    C         dSigmaDrS(i,j)=op5*( sigmaR(i,j-1,k)+sigmaR(i,j,k) )
522    C     &        *maskS(i,j,k,bi,bj)
523    C         dummy(i,j) = zeroRL
524    C        ENDDO
525    C       ENDDO
526    C       CALL GMREDI_SLOPE_PSI(
527    C     O      tprX, tprY,
528    C     U      slpX, slpY,
529    C     U      dSigmaDrW, dSigmaDrS,
530    C     I      dummy, dummy, rF(k), k,
531    C     I      bi, bj, myThid)
532    C       DO j=1-Oly,sNy+Oly
533    C        DO i=1-Olx,sNx+Olx
534    C         SlopeX(i,j,k) = slpX(i,j)
535    C         SlopeY(i,j,k) = slpY(i,j)
536    C        ENDDO
537    C       ENDDO
538    C      ENDDO
539    
540  C     Calculate the thickness flux  C     Calculate the thickness flux
541  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)
542        k=Nr        k=Nr
# Line 557  C     We still need to give special trea Line 558  C     We still need to give special trea
558        DO k=Nr-1,1,-1        DO k=Nr-1,1,-1
559         DO j=1-Oly,sNy+Oly-1         DO j=1-Oly,sNy+Oly-1
560          DO i=1-Olx,sNx+Olx-1          DO i=1-Olx,sNx+Olx-1
561           IF(k.LE.surfk(i,j)) THEN           IF(k.LE.surfk(i,j) .AND. .NOT. GM_K3D_likeGM) THEN
562  C          We are in the surface layer, so set the thickness flux  C          We are in the surface layer, so set the thickness flux
563  C          based on the average slope over the surface layer  C          based on the average slope over the surface layer
564  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 664  C     It does not affect the end result Line 665  C     It does not affect the end result
665        ENDDO        ENDDO
666    
667        IF(GM_K3D_likeGM) THEN        IF(GM_K3D_likeGM) THEN
 C       To make it similar to GM, we do not do any smoothing  
         m=1  
         DO j=1-Oly,sNy+Oly  
          DO i=1-Olx,sNx+Olx  
           XimX(1,i,j) = zeroRL  
           XimY(1,i,j) = zeroRL  
          ENDDO  
         ENDDO  
668          DO k=1,Nr          DO k=1,Nr
669           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
670            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
671             IF (GM_K3D_likeGM) THEN             KPV(i,j,k) = GM_K3D_constK
              Xix(i,j,k)=-GM_K3D_constK*gradqx(i,j,k)  
              Xiy(i,j,k)=-GM_K3D_constK*gradqy(i,j,k)  
            ELSE  
              Xix(i,j,k)=-K3D(i,j,k,bi,bj)*gradqx(i,j,k)  
              Xiy(i,j,k)=-K3D(i,j,k,bi,bj)*gradqy(i,j,k)  
            ENDIF  
            XimX(1,i,j) = XimX(m,i,j)  
      &          + Xix(i,j,k)*drF(k)*hfacW(i,j,k,bi,bj)  
            XimY(1,i,j) = XimY(m,i,j)  
      &          + Xiy(i,j,k)*drF(k)*hfacS(i,j,k,bi,bj)  
672            ENDDO            ENDDO
673           ENDDO           ENDDO
674          ENDDO          ENDDO
675          DO j=1-Oly,sNy+Oly        ELSE
676           DO i=1-Olx,sNx+Olx          DO k=1,Nr
677  C         minus comes from rlow being negative.           DO j=1-Oly,sNy+Oly
678            IF(rLowW(i,j,bi,bj).LT.0.0)            DO i=1-Olx,sNx+Olx
679       &         XimX(1,i,j) = -XimX(1,i,j)/rLowW(i,j,bi,bj)             KPV(i,j,k) = K3D(i,j,k,bi,bj)
680            IF(rLowS(i,j,bi,bj).LT.0.0)            ENDDO
      &         XimY(1,i,j) = -XimY(1,i,j)/rLowS(i,j,bi,bj)  
681           ENDDO           ENDDO
682          ENDDO          ENDDO
683          ENDIF
684    
685          IF (.NOT. GM_K3D_smooth) THEN
686    C       Do not expand K grad(q) => no smoothing
687    C       May only be done with a constant K, otherwise the
688    C       integral constraint is violated.
689          DO k=1,Nr          DO k=1,Nr
690           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
691            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
692             Xix(i,j,k) = Xix(i,j,k) - maskW(i,j,k,bi,bj)*XimX(1,i,j)             Xix(i,j,k) = -maskW(i,j,k,bi,bj)*KPV(i,j,k)*gradqx(i,j,k)
693             Xiy(i,j,k) = Xiy(i,j,k) - maskS(i,j,k,bi,bj)*XimY(1,i,j)             Xiy(i,j,k) = -maskS(i,j,k,bi,bj)*KPV(i,j,k)*gradqy(i,j,k)
694            ENDDO            ENDDO
695           ENDDO           ENDDO
696          ENDDO          ENDDO
697    
698        ELSE        ELSE
699  C       Expand K grad(q) in terms of baroclinic modes  C       Expand K grad(q) in terms of baroclinic modes to smooth
700    C       and satisfy the integral constraint
701    
702  C       Start with the X direction  C       Start with the X direction
703  C       ------------------------------  C       ------------------------------
# Line 733  C       Calculate Xi_m at the west face Line 721  C       Calculate Xi_m at the west face
721           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
722            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
723             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
             IF (GM_K3D_likeGM) THEN  
               XimX(m,i,j) = XimX(m,i,j)  
      &             - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)  
      &             *GM_K3D_constK*gradqx(i,j,k)*vec(m,i,j,k)  
             ELSE  
724                XimX(m,i,j) = XimX(m,i,j)                XimX(m,i,j) = XimX(m,i,j)
725       &             - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)       &             - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)
726       &             *K3D(i,j,k,bi,bj)*gradqx(i,j,k)*vec(m,i,j,k)       &             *KPV(i,j,k)*gradqx(i,j,k)*vec(m,i,j,k)
             ENDIF  
727             ENDDO             ENDDO
728            ENDDO            ENDDO
729           ENDDO           ENDDO
# Line 786  C     Calculate the eigenvectors at the Line 768  C     Calculate the eigenvectors at the
768           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
769            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
770             DO m=1,GM_K3D_NModes             DO m=1,GM_K3D_NModes
771              IF (GM_K3D_likeGM) THEN              XimY(m,i,j) = XimY(m,i,j)
772                XimY(m,i,j) = XimY(m,i,j)       &           - drF(k)*hfacS(i,j,k,bi,bj)
773       &             - maskS(i,j,k,bi,bj)*drF(k)*hfacS(i,j,k,bi,bj)       &           *KPV(i,j,k)*gradqy(i,j,k)*vec(m,i,j,k)
      &             *GM_K3D_constK*gradqy(i,j,k)*vec(m,i,j,k)  
             ELSE  
               XimY(m,i,j) = XimY(m,i,j)  
      &             - drF(k)*hfacS(i,j,k,bi,bj)  
      &             *K3D(i,j,k,bi,bj)*gradqy(i,j,k)*vec(m,i,j,k)  
             ENDIF  
774             ENDDO             ENDDO
775            ENDDO            ENDDO
776           ENDDO           ENDDO
# Line 819  C     Calculate Xi for Y direction at th Line 795  C     Calculate Xi for Y direction at th
795           ENDDO           ENDDO
796          ENDDO          ENDDO
797    
798    C     ENDIF GM_K3D_likeGM
799        ENDIF        ENDIF
800    
801    
# Line 925  C     Diagnostics Line 902  C     Diagnostics
902          CALL DIAGNOSTICS_FILL(Kdef,   'GM_KDEF ',0, 1,0,1,1,myThid)          CALL DIAGNOSTICS_FILL(Kdef,   'GM_KDEF ',0, 1,0,1,1,myThid)
903        ENDIF        ENDIF
904  #endif  #endif
         
 C     Even when using a constant K, we diagnose the K that we would  
 C     use.  So, after the diagnostics are written, we change the  
 C     array K3D to be equal to GM_K3D_constK for use in the Redi tensor.  
       IF (GM_K3D_likeGM) THEN  
         DO k=1,Nr  
          DO j=1-Oly,sNy+Oly  
           DO i=1-Olx,sNx+Olx  
            K3D(i,j,k,bi,bj) = GM_K3D_constK  
           ENDDO  
          ENDDO  
         ENDDO  
       ENDIF  
905    
906  #endif /* GM_K3D */  #endif /* GM_K3D */
907        RETURN        RETURN
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22