pkg/my82/my82_ri_number.F

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