/[MITgcm]/MITgcm/pkg/my82/my82_ri_number.F
ViewVC logotype

Diff of /MITgcm/pkg/my82/my82_ri_number.F

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

revision 1.3 by jmc, Thu Aug 23 19:13:10 2007 UTC revision 1.4 by jmc, Mon Aug 11 22:28:06 2008 UTC
# Line 8  C !ROUTINE: MY82_RI_NUMBER Line 8  C !ROUTINE: MY82_RI_NUMBER
8    
9  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
10        subroutine MY82_RI_NUMBER(        subroutine MY82_RI_NUMBER(
11       I     bi, bj, K, iMin, iMax, jMin, jMax,       I     bi, bj, K, iMin, iMax, jMin, jMax,
12       O     RiNumber, buoyFreq, vertShear,       O     RiNumber, buoyFreq, vertShear,
13       I     myTime, myThid )       I     myTime, myThid )
14    
# Line 41  C !USES: =============================== Line 41  C !USES: ===============================
41  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
42  C Routine arguments  C Routine arguments
43  C     bi, bj - array indices on which to apply calculations  C     bi, bj - array indices on which to apply calculations
44  C     iMin, iMax, jMin, jMax  C     iMin, iMax, jMin, jMax
45  C            - array boundaries  C            - array boundaries
46  C     k      - depth level  C     k      - depth level
47  C     myTime - Current time in simulation  C     myTime - Current time in simulation
# Line 64  C     tempu, tempv - temporary variables Line 64  C     tempu, tempv - temporary variables
64  C     rhoK, rhoKm1 - Density below and above current interface  C     rhoK, rhoKm1 - Density below and above current interface
65  C     epsilon      - small number  C     epsilon      - small number
66        INTEGER I,J,Km1        INTEGER I,J,Km1
67        _RL        p5    , p125              _RL        p5    , p125
68        PARAMETER( p5=0.5D0, p125=0.125D0 )        PARAMETER( p5=0.5D0, p125=0.125D0 )
69        _RL tempu, tempv        _RL tempu, tempv
70        _RL epsilon        _RL epsilon
# Line 74  C     epsilon      - small number Line 74  C     epsilon      - small number
74  #ifdef MY82_SMOOTH_RI  #ifdef MY82_SMOOTH_RI
75        _RL RiTmp   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL RiTmp   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76  #endif /* MY82_SMOOTH_RI */  #endif /* MY82_SMOOTH_RI */
77  CEOP      CEOP
78    
79        Km1   = MAX(1,K-1)        Km1   = MAX(1,K-1)
80  C     Not sure what is correct for pressure coordinates:  C     Not sure what is correct for pressure coordinates:
# Line 92  CML      ELSE Line 92  CML      ELSE
92  CML       STOP 'MY82_RI_NUMBER: We should never reach this point'  CML       STOP 'MY82_RI_NUMBER: We should never reach this point'
93  CML      ENDIF  CML      ENDIF
94  C     preparation: find density a kth and k-1st level  C     preparation: find density a kth and k-1st level
95        CALL FIND_RHO(        CALL FIND_RHO_2D(
96       I     bi, bj, iMin, iMax, jMin, jMax, Km1, K,       I     iMin, iMax, jMin, jMax, K,
97       I     theta, salt,       I     theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),
98       O     rhoKm1,       O     rhoKm1,
99       I     myThid )       I     Km1, bi, bj, myThid )
100        CALL FIND_RHO(        CALL FIND_RHO_2D(
101       I     bi, bj, iMin, iMax, jMin, jMax, K, K,       I     iMin, iMax, jMin, jMax, K,
102       I     theta, salt,       I     theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),
103       O     rhoK,       O     rhoK,
104       I     myThid )       I     K, bi, bj, myThid )
105    
106  C     first step:  calculate flux Richardson number.  C     first step:  calculate flux Richardson number.
107  C     calculate (du/dz)^2 + (dv/dz)^2 on W (between T) points.  C     calculate (du/dz)^2 + (dv/dz)^2 on W (between T) points.
# Line 120  C     calculate g*(drho/dz)/rho0= -N^2 Line 120  C     calculate g*(drho/dz)/rho0= -N^2
120  C  C
121          buoyFreq(I,J) = gravity*mass2rUnit *          buoyFreq(I,J) = gravity*mass2rUnit *
122       &       (rhoKm1(I,J) - rhoK(I,J))*recip_drC(K)       &       (rhoKm1(I,J) - rhoK(I,J))*recip_drC(K)
123  C      C
124  C     calculates gradient Richardson number, bounded by  C     calculates gradient Richardson number, bounded by
125  C     a very large negative value.  C     a very large negative value.
126  C      C
127          RiNumber(I,J) = -buoyFreq(I,J)/max(vertShear(I,J),epsilon)          RiNumber(I,J) = -buoyFreq(I,J)/max(vertShear(I,J),epsilon)
128            
129         ENDDO         ENDDO
130        ENDDO        ENDDO
131          
132  #ifdef MY82_SMOOTH_RI  #ifdef MY82_SMOOTH_RI
133  C     average Richardson number horizontally  C     average Richardson number horizontally
134        DO J=jMin,jMax        DO J=jMin,jMax
135         DO I=iMin,iMax         DO I=iMin,iMax
136          RiTmp(I,J) = RiNumber(I,J)          RiTmp(I,J) = RiNumber(I,J)
137         ENDDO         ENDDO
138        ENDDO            ENDDO
139        DO J=1-Oly+1,sNy+Oly-1        DO J=1-Oly+1,sNy+Oly-1
140         DO I=1-Olx+1,sNx+Olx-1         DO I=1-Olx+1,sNx+Olx-1
141          RiNumber(I,J) = p5*RiTmp(I,J)          RiNumber(I,J) = p5*RiTmp(I,J)
142       &       + p125*RiTmp(I-1,J) + p125*RiTmp(I+1,J)       &       + p125*RiTmp(I-1,J) + p125*RiTmp(I+1,J)
143       &       + p125*RiTmp(I,J-1) + p125*RiTmp(I,J+1)       &       + p125*RiTmp(I,J-1) + p125*RiTmp(I,J+1)
144         ENDDO         ENDDO
145        ENDDO            ENDDO
146  #endif /* MY82_SMOOTH_RI */  #endif /* MY82_SMOOTH_RI */
147    
148  #endif /* ALLOW_MY82 */  #endif /* ALLOW_MY82 */
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22