pkg/pp81/pp81_calc.F

08be60903a Mart*0001 #include "PP81_OPTIONS.h"
                0002
                0003 CBOP
                0004 C !ROUTINE: PP81_CALC
                0005
                0006 C !INTERFACE: =======================================================
                0007       subroutine PP81_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 PP81_CALC                                     |
                0013 C     | o Compute all PP81 fields defined in PP81.h              |
5e48dccc42 Jean*0014 C     *==========================================================*
08be60903a Mart*0015 C     | This subroutine is based on SPEM code                    |
5e48dccc42 Jean*0016 C     *==========================================================*
08be60903a Mart*0017
                0018 C global parameters updated by pp_calc
5e48dccc42 Jean*0019 C     PPviscAz  :: PP eddy viscosity coefficient              (m^2/s)
                0020 C     PPdiffKzT :: PP diffusion coefficient for temperature   (m^2/s)
08be60903a Mart*0021 C
                0022 C \ev
                0023
                0024 C !USES: ============================================================
5e48dccc42 Jean*0025       IMPLICIT NONE
08be60903a Mart*0026 #include "SIZE.h"
                0027 #include "EEPARAMS.h"
                0028 #include "PARAMS.h"
                0029 #include "GRID.h"
40cf87ade1 Jean*0030 #ifdef ALLOW_3D_DIFFKR
                0031 # include "DYNVARS.h"
                0032 #endif
                0033 #include "PP81.h"
                0034 c#ifdef ALLOW_AUTODIFF_TAMC
                0035 c#include "tamc.h"
                0036 c#endif /* ALLOW_AUTODIFF_TAMC */
08be60903a Mart*0037
                0038 C !INPUT PARAMETERS: ===================================================
dc6107c029 Jean*0039 C Routine arguments
                0040 C     bi, bj :: Current tile indices
                0041 C     sigmaR :: Vertical gradient of iso-neutral density
                0042 C     myTime :: Current time in simulation
                0043 C     myIter :: Current time-step number
                0044 C     myThid :: My Thread Id number
08be60903a Mart*0045       INTEGER bi, bj
dc6107c029 Jean*0046       _RL     sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
08be60903a Mart*0047       _RL     myTime
dc6107c029 Jean*0048       INTEGER myIter
5e48dccc42 Jean*0049       INTEGER myThid
08be60903a Mart*0050
                0051 #ifdef ALLOW_PP81
                0052 C !LOCAL VARIABLES: ====================================================
                0053 c Local constants
                0054 C     imin, imax, jmin, jmax  - array computation indices
                0055 C     RiNumber                - Richardson Number
                0056       INTEGER I, J, K
                0057       INTEGER   iMin ,iMax ,jMin ,jMax
                0058       _RL     denom, PPviscTmp
                0059       _RL     RiNumber(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0060 CEOP
                0061       iMin = 2-OLx
                0062       iMax = sNx+OLx-1
                0063       jMin = 2-OLy
                0064       jMax = sNy+OLy-1
                0065
                0066       DO K = 2, Nr
                0067        CALL PP81_RI_NUMBER(
5e48dccc42 Jean*0068      I      bi, bj, K, iMin, iMax, jMin, jMax,
                0069      O      RiNumber,
08be60903a Mart*0070      I      myTime, myThid )
                0071        DO J=jMin,jMax
                0072         DO I=iMin,iMax
                0073          IF ( RiNumber(I,J) .LT. RiLimit ) THEN
                0074           denom = 1.0 + PPalpha*RiLimit
                0075           PPviscTmp = PPviscMax
                0076          ELSE
                0077           denom = 1.0 + PPalpha*RiNumber(I,J)
                0078           PPviscTmp = PPnu0/(denom**PPnRi)
                0079          ENDIF
                0080 C     assign a minium ( = background ) value
5e48dccc42 Jean*0081          PPviscAr(I,J,K,bi,bj) = MAX(PPviscTmp,viscArNr(k))
08be60903a Mart*0082          PPdiffKr(I,J,K,bi,bj) = MAX(PPviscAr(I,J,K,bi,bj)/denom,
40cf87ade1 Jean*0083 #ifdef ALLOW_3D_DIFFKR
                0084      &                               diffKr(i,j,k,bi,bj) )
                0085 #else
78524d1402 Jean*0086      &                               diffKrNrS(k) )
40cf87ade1 Jean*0087 #endif
08be60903a Mart*0088 CML         if ( k.eq.2 .and. i.ge.1 .and. i.le.sNx .and. j.eq.1)
                0089 CML     &        print '(A,3I3,5E14.5)', 'ml-pp81', I,J,K, RiLimit,
                0090 CML     &        RiNumber(I,J),denom,
                0091 CML     &        PPviscAr(I,J,K,bi,bj), PPdiffKr(I,J,K,bi,bj)
                0092         ENDDO
                0093        ENDDO
                0094 #ifdef ALLOW_PP81_LOWERBOUND
5e48dccc42 Jean*0095 CRT   This is where the lower limit for subsurface layers
                0096 CRT   (BRIOS special) is set.
6709240de5 Jean*0097        IF ( (usingZCoords .AND. K .EQ. 2) .OR.
                0098      &      (usingPCoords .AND. K .EQ. Nr) ) THEN
08be60903a Mart*0099         DO J=jMin,jMax
                0100          DO I=iMin,iMax
                0101           PPviscAr(I,J,K,bi,bj) = MAX(PPviscMin,PPviscAr(I,J,K,bi,bj))
                0102           PPdiffKr(I,J,K,bi,bj) = MAX(PPdiffMin,PPdiffKr(I,J,K,bi,bj))
5e48dccc42 Jean*0103          ENDDO
08be60903a Mart*0104         ENDDO
                0105        ENDIF
                0106 #endif /* ALLOW_PP81_LOWERBOUND */
                0107 C     Mask land points
                0108        DO J=jMin,jMax
                0109         DO I=iMin,iMax
                0110          PPviscAr(I,J,K,bi,bj) = PPviscAr(I,J,K,bi,bj)
                0111      &        * maskC(I,J,K,bi,bj)
                0112          PPdiffKr(I,J,K,bi,bj) = PPdiffKr(I,J,K,bi,bj)
                0113      &        * maskC(I,J,K,bi,bj)
5e48dccc42 Jean*0114         ENDDO
08be60903a Mart*0115        ENDDO
                0116 C     end K-loop
                0117       ENDDO
                0118
016b84c482 Mart*0119 #ifdef ALLOW_DIAGNOSTICS
                0120       IF ( useDiagnostics ) THEN
                0121        CALL DIAGNOSTICS_FILL(PPviscAr,'PPviscAr',0,Nr,1,bi,bj,myThid)
                0122        CALL DIAGNOSTICS_FILL(PPdiffKr,'PPdiffKr',0,Nr,1,bi,bj,myThid)
                0123       ENDIF
                0124 #endif /* ALLOW_DIAGNOSTICS */
                0125
08be60903a Mart*0126 #endif /* ALLOW_PP81 */
                0127
                0128       RETURN
                0129       END