pkg/my82/my82_calc.F

08be60903a Mart*0001 #include "MY82_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: MY82_CALC
                0005
                0006 C !INTERFACE: ======================================================
                0007       subroutine MY82_CALC(
dc6107c029 Jean*0008      I                bi, bj, sigmaR, myTime, myIter, myThid )
08be60903a Mart*0009
                0010 C !DESCRIPTION: \bv
5e48dccc42 Jean*0011 C     *==========================================================*
08be60903a Mart*0012 C     | SUBROUTINE MY82_CALC                                     |
                0013 C     | o Compute all MY82 fields defined in MY82.h              |
5e48dccc42 Jean*0014 C     *==========================================================*
08be60903a Mart*0015 C     | This subroutine is based on SPEM code                    |
5e48dccc42 Jean*0016 C     *==========================================================*
08be60903a Mart*0017       IMPLICIT NONE
                0018 C
                0019 C--------------------------------------------------------------------
                0020
                0021 C global parameters updated by pp_calc
016b84c482 Mart*0022 C     MYviscAz  :: MY82 eddy viscosity coefficient              (m^2/s)
                0023 C     MYdiffKzT :: MY82 diffusion coefficient for temperature   (m^2/s)
08be60903a Mart*0024 C
                0025 C \ev
                0026
                0027 C !USES: ============================================================
                0028 #include "SIZE.h"
                0029 #include "EEPARAMS.h"
                0030 #include "PARAMS.h"
                0031 #include "GRID.h"
dbbea2be4e Jean*0032 #include "MY82.h"
                0033 #ifdef ALLOW_3D_DIFFKR
                0034 # include "DYNVARS.h"
                0035 #endif
08be60903a Mart*0036
                0037 C !INPUT PARAMETERS: ===================================================
dc6107c029 Jean*0038 C Routine arguments
                0039 C     bi, bj :: Current tile indices
                0040 C     sigmaR :: Vertical gradient of iso-neutral density
                0041 C     myTime :: Current time in simulation
                0042 C     myIter :: Current time-step number
                0043 C     myThid :: My Thread Id number
08be60903a Mart*0044       INTEGER bi, bj
dc6107c029 Jean*0045       _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
08be60903a Mart*0046       _RL     myTime
dc6107c029 Jean*0047       INTEGER myIter
fb62f539dc Jean*0048       INTEGER myThid
08be60903a Mart*0049
                0050 C !LOCAL VARIABLES: ====================================================
                0051 c Local constants
                0052 C     imin, imax, jmin, jmax  - array computation indices
                0053 C     RiNumber - Richardson Number = -GH/GM
                0054 C     GH       - buoyancy freqency after call to Ri_number, later
                0055 C                GH of M. Satoh, Eq. (11.3.45)
                0056 C     GM       - vertical shear of velocity after call Ri_number,
                0057 C                later GM of M. Satoh, Eq. (11.3.44)
                0058       INTEGER I, J, K
                0059       INTEGER   iMin ,iMax ,jMin ,jMax
                0060       _RL     RiFlux, RiTmp
                0061       _RL     SHtmp, bTmp, tkesquare, tkel
                0062       _RL     RiNumber(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0063       _RL     GH(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0064       _RL     GM(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0065       _RL     SH(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0066       _RL     SM(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0067       _RL     tke(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0068 CEOP
                0069       iMin = 2-OLx
                0070       iMax = sNx+OLx-1
                0071       jMin = 2-OLy
                0072       jMax = sNy+OLy-1
                0073
                0074 C     Initialize local fields
dc6107c029 Jean*0075       DO J=1-OLy,sNy+OLy
                0076        DO I=1-OLx,sNx+OLx
fb62f539dc Jean*0077         GH(I,J) = 0. _d 0
                0078         GM(I,J) = 0. _d 0
08be60903a Mart*0079        ENDDO
                0080       ENDDO
                0081       DO K = 1, Nr
dc6107c029 Jean*0082        DO J=1-OLy,sNy+OLy
                0083         DO I=1-OLx,sNx+OLx
8fef4b1a57 Mart*0084          SH(I,J,K)  = 0. _d 0
                0085          SM(I,J,K)  = 0. _d 0
                0086          tke(I,J,K) = 0. _d 0
08be60903a Mart*0087         ENDDO
fb62f539dc Jean*0088        ENDDO
08be60903a Mart*0089       ENDDO
                0090 C     first k-loop
                0091 C     compute turbulent kinetic energy from richardson number
                0092       DO K = 2, Nr
                0093        CALL MY82_RI_NUMBER(
fb62f539dc Jean*0094      I      bi, bj, K, iMin, iMax, jMin, jMax,
08be60903a Mart*0095      O      RiNumber, GH, GM,
                0096      I      myTime, myThid )
                0097        DO J=jMin,jMax
                0098         DO I=iMin,iMax
                0099          RiTmp  = MIN(RiNumber(I,J),RiMax)
                0100          btmp   = beta1+beta4*RiTmp
                0101 C     M. Satoh, Atmospheric Circulation Dynamics and General
                0102 C     Circulation models, Springer, 2004: Eq. (11.3.60)
fb62f539dc Jean*0103          RiFlux =( btmp - SQRT(btmp*btmp - 4. _d 0 *beta2*beta3*RiTmp) )
                0104      &        /(2. _d 0*beta2)
08be60903a Mart*0105 C     M. Satoh: Eq. (11.3.58)
fb62f539dc Jean*0106          SHtmp       = (alpha1-alpha2*RiFlux)/(1. _d 0-RiFlux)
08be60903a Mart*0107          SH(I,J,K)   = SHtmp
                0108          SM(I,J,K)   = SHtmp*(beta1-beta2*RiFlux)/(beta3-beta4*RiFlux)
                0109 C     M. Satoh: Eq. (11.3.53/55)
                0110          tkesquare = MAX( 0. _d 0,
                0111      &        b1*(SH(I,J,K)*GH(I,J) + SM(I,J,K)*GM(I,J)) )
                0112          tke(I,J,K) = sqrt(tkesquare)
                0113 CML         if ( k.eq.2 .and. i.ge.1 .and. i.le.sNx .and. j.eq.1)
                0114 CML     &   print '(A,3I3,8E13.5)', 'ml-my82', I,J,K, RiNumber(I,J),
                0115 CML     &        RiTmp,RiFlux,
                0116 CML     &      SH(I,J,K), SM(I,J,K), GH(I,J), GM(I,J), tke(I,J,K)
                0117         ENDDO
                0118        ENDDO
                0119 C     end of first k-loop
fb62f539dc Jean*0120       ENDDO
08be60903a Mart*0121
fb62f539dc Jean*0122 C     re-initilialize GM and GH for abuse, they no longer have
08be60903a Mart*0123 C     the meaning of shear and negative  buoyancy frequency
                0124       DO J=jMin,jMax
                0125        DO I=iMin,iMax
8fef4b1a57 Mart*0126         GH(I,J) = 0. _d 0
                0127         GM(I,J) = 0. _d 0
08be60903a Mart*0128        ENDDO
                0129       ENDDO
                0130 C     Find boundary length scale from energy weighted mean.
                0131 C     This is something like "center of mass" of the vertical
                0132 C     tke-distribution
                0133 C     begin second k-loop
                0134       DO K = 2, Nr
                0135        DO J=jMin,jMax
                0136         DO I=iMin,iMax
                0137          GM(I,J) = GM(I,J) + tke(I,J,K)*rF(K)
fb62f539dc Jean*0138          GH(I,J) = GH(I,J) + tke(I,J,K)
08be60903a Mart*0139         ENDDO
                0140        ENDDO
                0141 C     end of second k-loop
                0142       END DO
                0143 C     compute boundary length scale MYhbl
                0144       DO J=jMin,jMax
                0145        DO I=iMin,iMax
8fef4b1a57 Mart*0146         IF ( GH(I,J) .EQ. 0. _d 0 ) THEN
                0147          MYhbl(I,J,bi,bj) = 0. _d 0
08be60903a Mart*0148         ELSE
                0149          MYhbl(I,J,bi,bj) = -GM(I,J)/GH(I,J)*MYhblScale
                0150         ENDIF
                0151        ENDDO
fb62f539dc Jean*0152       ENDDO
08be60903a Mart*0153 C     begin third k-loop
                0154       DO K = 1, Nr
                0155 C     integrate tke to find integral length scale
                0156        DO J=jMin,jMax
                0157         DO I=iMin,iMax
                0158          tkel = MYhbl(I,J,bi,bj)*tke(I,J,K)
                0159 C     M. Satoh: Eq. (11.3.43)
                0160          MYviscAr(I,J,K,bi,bj) = MYhbl(I,J,bi,bj)*tkel*SM(I,J,K)
                0161          MYdiffKr(I,J,K,bi,bj) = MYhbl(I,J,bi,bj)*tkel*SH(I,J,K)
                0162 C     Set a minium (= background) value
5e48dccc42 Jean*0163          MYviscAr(I,J,K,bi,bj) = MAX(MYviscAr(I,J,K,bi,bj),
                0164      &                               viscArnr(k) )
648edfaf30 Jean*0165          MYdiffKr(I,J,K,bi,bj) = MAX(MYdiffKr(I,J,K,bi,bj),
dbbea2be4e Jean*0166 #ifdef ALLOW_3D_DIFFKR
                0167      &                               diffKr(i,j,k,bi,bj) )
                0168 #else
78524d1402 Jean*0169      &                               diffKrNrS(k) )
dbbea2be4e Jean*0170 #endif
08be60903a Mart*0171 C     Set a maximum and mask land point
                0172          MYviscAr(I,J,K,bi,bj) = MIN(MYviscAr(I,J,K,bi,bj),MYviscMax)
                0173      &        * maskC(I,J,K,bi,bj)
                0174          MYdiffKr(I,J,K,bi,bj) = MIN(MYdiffKr(I,J,K,bi,bj),MYdiffMax)
                0175      &        * maskC(I,J,K,bi,bj)
                0176         ENDDO
fb62f539dc Jean*0177        ENDDO
08be60903a Mart*0178 C     end third k-loop
fb62f539dc Jean*0179       ENDDO
08be60903a Mart*0180
016b84c482 Mart*0181 #ifdef ALLOW_DIAGNOSTICS
                0182       IF ( useDiagnostics ) THEN
                0183        CALL DIAGNOSTICS_FILL(MYviscAr,'MYVISCAR',0,Nr,1,bi,bj,myThid)
                0184        CALL DIAGNOSTICS_FILL(MYdiffKr,'MYDIFFKR',0,Nr,1,bi,bj,myThid)
                0185        CALL DIAGNOSTICS_FILL(MYhbl   ,'MYHBL   ',0,1 ,1,bi,bj,myThid)
                0186       ENDIF
                0187 #endif /* ALLOW_DIAGNOSTICS */
08be60903a Mart*0188
                0189       RETURN
                0190       END