C $Header: /u/gcmpack/MITgcm/pkg/my82/my82_calc.F,v 1.9 2015/02/23 21:20:15 jmc Exp $ C $Name: $ #include "MY82_OPTIONS.h" CBOP C !ROUTINE: MY82_CALC C !INTERFACE: ====================================================== subroutine MY82_CALC( I bi, bj, sigmaR, myTime, myIter, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE MY82_CALC | C | o Compute all MY82 fields defined in MY82.h | C *==========================================================* C | This subroutine is based on SPEM code | C *==========================================================* IMPLICIT NONE C 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: ============================================================ #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "MY82.h" #ifdef ALLOW_3D_DIFFKR # include "DYNVARS.h" #endif 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_MY82 C !LOCAL VARIABLES: ==================================================== c Local constants C imin, imax, jmin, jmax - array computation indices C RiNumber - Richardson Number = -GH/GM C GH - buoyancy freqency after call to Ri_number, later C GH of M. Satoh, Eq. (11.3.45) C GM - vertical shear of velocity after call Ri_number, C later GM of M. Satoh, Eq. (11.3.44) INTEGER I, J, K INTEGER iMin ,iMax ,jMin ,jMax _RL RiFlux, RiTmp _RL SHtmp, bTmp, tkesquare, tkel _RL RiNumber(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL GH(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL GM(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL SH(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL SM(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) _RL tke(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr) CEOP iMin = 2-OLx iMax = sNx+OLx-1 jMin = 2-OLy jMax = sNy+OLy-1 C Initialize local fields DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx GH(I,J) = 0. _d 0 GM(I,J) = 0. _d 0 ENDDO ENDDO DO K = 1, Nr DO J=1-OLy,sNy+OLy DO I=1-OLx,sNx+OLx SH(I,J,K) = 0. _d 0 SM(I,J,K) = 0. _d 0 tke(I,J,K) = 0. _d 0 ENDDO ENDDO ENDDO C first k-loop C compute turbulent kinetic energy from richardson number DO K = 2, Nr CALL MY82_RI_NUMBER( I bi, bj, K, iMin, iMax, jMin, jMax, O RiNumber, GH, GM, I myTime, myThid ) DO J=jMin,jMax DO I=iMin,iMax RiTmp = MIN(RiNumber(I,J),RiMax) btmp = beta1+beta4*RiTmp C M. Satoh, Atmospheric Circulation Dynamics and General C Circulation models, Springer, 2004: Eq. (11.3.60) RiFlux =( btmp - SQRT(btmp*btmp - 4. _d 0 *beta2*beta3*RiTmp) ) & /(2. _d 0*beta2) C M. Satoh: Eq. (11.3.58) SHtmp = (alpha1-alpha2*RiFlux)/(1. _d 0-RiFlux) SH(I,J,K) = SHtmp SM(I,J,K) = SHtmp*(beta1-beta2*RiFlux)/(beta3-beta4*RiFlux) C M. Satoh: Eq. (11.3.53/55) tkesquare = MAX( 0. _d 0, & b1*(SH(I,J,K)*GH(I,J) + SM(I,J,K)*GM(I,J)) ) tke(I,J,K) = sqrt(tkesquare) CML if ( k.eq.2 .and. i.ge.1 .and. i.le.sNx .and. j.eq.1) CML & print '(A,3I3,8E13.5)', 'ml-my82', I,J,K, RiNumber(I,J), CML & RiTmp,RiFlux, CML & SH(I,J,K), SM(I,J,K), GH(I,J), GM(I,J), tke(I,J,K) ENDDO ENDDO C end of first k-loop ENDDO C re-initilialize GM and GH for abuse, they no longer have C the meaning of shear and negative buoyancy frequency DO J=jMin,jMax DO I=iMin,iMax GH(I,J) = 0. _d 0 GM(I,J) = 0. _d 0 ENDDO ENDDO C Find boundary length scale from energy weighted mean. C This is something like "center of mass" of the vertical C tke-distribution C begin second k-loop DO K = 2, Nr DO J=jMin,jMax DO I=iMin,iMax GM(I,J) = GM(I,J) + tke(I,J,K)*rF(K) GH(I,J) = GH(I,J) + tke(I,J,K) ENDDO ENDDO C end of second k-loop END
DO C compute boundary length scale MYhbl DO J=jMin,jMax DO I=iMin,iMax IF ( GH(I,J) .EQ. 0. _d 0 ) THEN MYhbl(I,J,bi,bj) = 0. _d 0 ELSE MYhbl(I,J,bi,bj) = -GM(I,J)/GH(I,J)*MYhblScale ENDIF ENDDO ENDDO C begin third k-loop DO K = 1, Nr C integrate tke to find integral length scale DO J=jMin,jMax DO I=iMin,iMax tkel = MYhbl(I,J,bi,bj)*tke(I,J,K) C M. Satoh: Eq. (11.3.43) MYviscAr(I,J,K,bi,bj) = MYhbl(I,J,bi,bj)*tkel*SM(I,J,K) MYdiffKr(I,J,K,bi,bj) = MYhbl(I,J,bi,bj)*tkel*SH(I,J,K) C Set a minium (= background) value MYviscAr(I,J,K,bi,bj) = MAX(MYviscAr(I,J,K,bi,bj), & viscArnr(k) ) MYdiffKr(I,J,K,bi,bj) = MAX(MYdiffKr(I,J,K,bi,bj), #ifdef ALLOW_3D_DIFFKR & diffKr(i,j,k,bi,bj) ) #else & diffKrNrS(k) ) #endif C Set a maximum and mask land point MYviscAr(I,J,K,bi,bj) = MIN(MYviscAr(I,J,K,bi,bj),MYviscMax) & * maskC(I,J,K,bi,bj) MYdiffKr(I,J,K,bi,bj) = MIN(MYdiffKr(I,J,K,bi,bj),MYdiffMax) & * maskC(I,J,K,bi,bj) ENDDO ENDDO C end third k-loop ENDDO #endif /* ALLOW_MY82 */ RETURN END