pkg/pp81/pp81_ri_number.F

08be60903a Mart*0001 #include "PP81_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: PP81_RI_NUMBER
                0005
                0006 C !INTERFACE: ==========================================================
4e4ad91a39 Jean*0007       SUBROUTINE PP81_RI_NUMBER(
94c8eb5701 Jean*0008      I     bi, bj, K, iMin, iMax, jMin, jMax,
08be60903a Mart*0009      O     RiNumber,
                0010      I     myTime, myThid )
                0011
                0012 C !DESCRIPTION: \bv
4e4ad91a39 Jean*0013 C     *==========================================================*
08be60903a Mart*0014 C     | SUBROUTINE PP81_RI_NUMBER                                |
                0015 C     | o Compute gradient Richardson number for Pacanowski and  |
                0016 C     |   Philander (1981) mixing scheme                         |
4e4ad91a39 Jean*0017 C     *==========================================================*
                0018 C \ev
08be60903a Mart*0019
                0020 C !USES: ===============================================================
4e4ad91a39 Jean*0021       IMPLICIT NONE
08be60903a Mart*0022 #include "SIZE.h"
                0023 #include "EEPARAMS.h"
                0024 #include "PARAMS.h"
4e4ad91a39 Jean*0025 #include "GRID.h"
08be60903a Mart*0026 #include "DYNVARS.h"
                0027 #include "PP81.h"
                0028
                0029 C !INPUT PARAMETERS: ===================================================
                0030 C Routine arguments
                0031 C     bi, bj - array indices on which to apply calculations
94c8eb5701 Jean*0032 C     iMin, iMax, jMin, jMax
08be60903a Mart*0033 C            - array boundaries
                0034 C     k      - depth level
                0035 C     myTime - Current time in simulation
                0036 C     RiNumber - (output) Richardson number
                0037       INTEGER bi, bj, iMin, iMax, jMin, jMax, k
                0038       INTEGER myThid
                0039       _RL     myTime
                0040       _RL     RiNumber(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0041
                0042 #ifdef ALLOW_PP81
                0043
                0044 C !LOCAL VARIABLES: ====================================================
                0045 C     I,J,kUp      - loop indices
                0046 C     p0-125       - averaging coefficients
                0047 C     tempu, tempv - temporary variables
                0048 C     rhoK, rhoKm1 - Density below and above current interface
                0049 C     epsilon      - small number
94c8eb5701 Jean*0050 C     RiFlux       - denominator of Richardson number
08be60903a Mart*0051 C     BuoyFreq     - buoyancy frequency
                0052       INTEGER I,J,Km1
94c8eb5701 Jean*0053       _RL        p5    , p125
08be60903a Mart*0054       PARAMETER( p5=0.5, p125=0.125 )
                0055       _RL tempu, tempv
                0056       _RL epsilon
                0057       PARAMETER    (  epsilon = 1.D-10 )
94c8eb5701 Jean*0058       _RL RiFlux
08be60903a Mart*0059       _RL buoyFreq
                0060       _RL rhoKm1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0061       _RL rhoK     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0062 #ifdef PP81_SMOOTH_RI
                0063       _RL RiTmp   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0064 #endif /* PP81_SMOOTH_RI */
94c8eb5701 Jean*0065 CEOP
08be60903a Mart*0066
                0067       Km1   = MAX(1,K-1)
                0068 C     Not sure what is correct for pressure coordinates:
                0069 C     RiFlux is always correct because it quadratic
                0070 C     buoyFreq should also be correct in pressure coordinates:
                0071 C     N^2=g^2drho/dp and K=1 at the bottom leads to the correct sign
                0072 C     So the following is wrong:
                0073 CML      IF ( buoyancyRelation .EQ. 'OCEANIC') THEN
                0074 CML       kUp   = MAX(1,K-1)
                0075 CML       kDown = K
                0076 CML      ELSEIF  ( buoyancyRelation .EQ. 'OCEANIP') THEN
                0077 CML       kUp   = K
                0078 CML       kDown = MAX(1,K-1)
                0079 CML      ELSE
                0080 CML       STOP 'PP81_RI_NUMBER: We should never reach this point'
                0081 CML      ENDIF
                0082 C     preparation: find density a kth and k-1st level
94c8eb5701 Jean*0083       CALL FIND_RHO_2D(
                0084      I     iMin, iMax, jMin, jMax, K,
                0085      I     theta(1-OLx,1-OLy,Km1,bi,bj), salt(1-OLx,1-OLy,Km1,bi,bj),
08be60903a Mart*0086      O     rhoKm1,
94c8eb5701 Jean*0087      I     Km1, bi, bj, myThid )
                0088       CALL FIND_RHO_2D(
                0089      I     iMin, iMax, jMin, jMax, K,
                0090      I     theta(1-OLx,1-OLy,K,bi,bj), salt(1-OLx,1-OLy,K,bi,bj),
08be60903a Mart*0091      O     rhoK,
94c8eb5701 Jean*0092      I     K, bi, bj, myThid )
08be60903a Mart*0093
                0094 C     first step:  calculate flux Richardson number.
                0095 C     calculate (du/dz)^2 + (dv/dz)^2 on W (between T) points.
                0096       DO J= jMin, jMax
                0097        DO I = iMin, iMax
                0098         tempu= .5*(  uVel(I,J,Km1,bi,bj)+uVel(I+1,J,Km1,bi,bj)
                0099      &            - (uVel(I,J,K  ,bi,bj)+uVel(I+1,J,K  ,bi,bj)) )
                0100      &       *recip_drC(K)
                0101         tempv= .5*(  vVel(I,J,Km1,bi,bj)+vVel(I,J+1,Km1,bi,bj)
                0102      &            - (vVel(I,J,K  ,bi,bj)+vVel(I,J+1,K  ,bi,bj)) )
                0103      &       *recip_drC(K)
                0104         RiFlux = tempu*tempu+tempv*tempv
                0105
                0106 C
                0107 C     calculate - g*(drho/dz)/rho0= N^2  .
                0108 C
0b1017b546 Jean*0109         buoyFreq = - gravity*mass2rUnit *
08be60903a Mart*0110      &       (rhoKm1(I,J) - rhoK(I,J))*recip_drC(K)
94c8eb5701 Jean*0111 C
08be60903a Mart*0112 C     calculates gradient Richardson number, bounded by
                0113 C     a very large negative value.
94c8eb5701 Jean*0114 C
08be60903a Mart*0115         RiNumber(I,J) = buoyFreq/max(RiFlux,epsilon)
                0116        ENDDO
                0117       ENDDO
94c8eb5701 Jean*0118
08be60903a Mart*0119 #ifdef PP81_SMOOTH_RI
                0120 C     average Richardson number horizontally
                0121       DO J=jMin,jMax
                0122        DO I=iMin,iMax
                0123         RiTmp(I,J) = RiNumber(I,J)
                0124        ENDDO
94c8eb5701 Jean*0125       ENDDO
4e4ad91a39 Jean*0126       DO J=1-OLy+1,sNy+OLy-1
                0127        DO I=1-OLx+1,sNx+OLx-1
94c8eb5701 Jean*0128         RiNumber(I,J) = p5*RiTmp(I,J)
                0129      &       + p125*RiTmp(I-1,J) + p125*RiTmp(I+1,J)
                0130      &       + p125*RiTmp(I,J-1) + p125*RiTmp(I,J+1)
08be60903a Mart*0131        ENDDO
94c8eb5701 Jean*0132       ENDDO
08be60903a Mart*0133 #endif /* PP81_SMOOTH_RI */
                0134
                0135 #endif /* ALLOW_PP81 */
                0136
                0137       RETURN
                0138       END