C $Header: /u/gcmpack/MITgcm/pkg/pp81/pp81_calc.F,v 1.8 2015/02/23 21:20:15 jmc Exp $ C $Name: $ #include "PP81_OPTIONS.h" CBOP C !ROUTINE: PP81_CALC C !INTERFACE: ======================================================= subroutine PP81_CALC( I bi, bj, sigmaR, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE PP81_CALC | C | o Compute all PP81 fields defined in PP81.h | C *==========================================================* C | This subroutine is based on SPEM code | C *==========================================================* C global parameters updated by pp_calc C PPviscAz :: PP eddy viscosity coefficient (m^2/s) C PPdiffKzT :: PP diffusion coefficient for temperature (m^2/s) C C \ev C !USES: ============================================================ IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #ifdef ALLOW_3D_DIFFKR # include "DYNVARS.h" #endif #include "PP81.h" c#ifdef ALLOW_AUTODIFF_TAMC c#include "tamc.h" c#include "tamc_keys.h" c#endif /* ALLOW_AUTODIFF_TAMC */ C !INPUT PARAMETERS: =================================================== C Routine arguments C bi, bj :: Current tile indices C sigmaR :: Vertical gradient of iso-neutral density C myTime :: Current time in simulation C myIter :: Current time-step number C myThid :: My Thread Id number INTEGER bi, bj _RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL myTime INTEGER myIter INTEGER myThid #ifdef ALLOW_PP81 C !LOCAL VARIABLES: ==================================================== c Local constants C imin, imax, jmin, jmax - array computation indices C RiNumber - Richardson Number INTEGER I, J, K INTEGER iMin ,iMax ,jMin ,jMax _RL denom, PPviscTmp _RL RiNumber(1-OLx:sNx+OLx,1-OLy:sNy+OLy) CEOP iMin = 2-OLx iMax = sNx+OLx-1 jMin = 2-OLy jMax = sNy+OLy-1 DO K = 2, Nr CALL PP81_RI_NUMBER( I bi, bj, K, iMin, iMax, jMin, jMax, O RiNumber, I myTime, myThid ) DO J=jMin,jMax DO I=iMin,iMax IF ( RiNumber(I,J) .LT. RiLimit ) THEN denom = 1.0 + PPalpha*RiLimit PPviscTmp = PPviscMax ELSE denom = 1.0 + PPalpha*RiNumber(I,J) PPviscTmp = PPnu0/(denom**PPnRi) ENDIF C assign a minium ( = background ) value PPviscAr(I,J,K,bi,bj) = MAX(PPviscTmp,viscArNr(k)) PPdiffKr(I,J,K,bi,bj) = MAX(PPviscAr(I,J,K,bi,bj)/denom, #ifdef ALLOW_3D_DIFFKR & diffKr(i,j,k,bi,bj) ) #else & diffKrNrS(k) ) #endif CML if ( k.eq.2 .and. i.ge.1 .and. i.le.sNx .and. j.eq.1) CML & print '(A,3I3,5E14.5)', 'ml-pp81', I,J,K, RiLimit, CML & RiNumber(I,J),denom, CML & PPviscAr(I,J,K,bi,bj), PPdiffKr(I,J,K,bi,bj) ENDDO ENDDO #ifdef ALLOW_PP81_LOWERBOUND CRT This is where the lower limit for subsurface layers CRT (BRIOS special) is set. IF ( (usingZCoords .AND. K .EQ. 2) .OR. & (usingPCoords .AND. K .EQ. Nr) ) THEN DO J=jMin,jMax DO I=iMin,iMax PPviscAr(I,J,K,bi,bj) = MAX(PPviscMin,PPviscAr(I,J,K,bi,bj)) PPdiffKr(I,J,K,bi,bj) = MAX(PPdiffMin,PPdiffKr(I,J,K,bi,bj)) ENDDO ENDDO ENDIF #endif /* ALLOW_PP81_LOWERBOUND */ C Mask land points DO J=jMin,jMax DO I=iMin,iMax PPviscAr(I,J,K,bi,bj) = PPviscAr(I,J,K,bi,bj) & * maskC(I,J,K,bi,bj) PPdiffKr(I,J,K,bi,bj) = PPdiffKr(I,J,K,bi,bj) & * maskC(I,J,K,bi,bj) ENDDO ENDDO C end K-loop ENDDO #endif /* ALLOW_PP81 */ RETURN END