/[MITgcm]/MITgcm/pkg/kpp/kpp_calc.F
ViewVC logotype

Diff of /MITgcm/pkg/kpp/kpp_calc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.27 by heimbach, Tue May 11 20:57:08 2004 UTC revision 1.28 by jmc, Sun Jul 18 01:20:22 2004 UTC
# Line 90  c     KPPfrac     - Fraction of short-wa Line 90  c     KPPfrac     - Fraction of short-wa
90    
91  c--   KPP_CALC computes vertical viscosity and diffusivity for region  c--   KPP_CALC computes vertical viscosity and diffusivity for region
92  c     (-2:sNx+3,-2:sNy+3) as required by CALC_DIFFUSIVITY and requires  c     (-2:sNx+3,-2:sNy+3) as required by CALC_DIFFUSIVITY and requires
93  c     values of uVel, vVel, SurfaceTendencyU, SurfaceTendencyV in the  c     values of uVel, vVel, surfaceForcingU, surfaceForcingV in the
94  c     region (-2:sNx+4,-2:sNy+4).  c     region (-2:sNx+4,-2:sNy+4).
95  c     Hence overlap region needs to be set OLx=4, OLy=4.  c     Hence overlap region needs to be set OLx=4, OLy=4.
96  c     When option FRUGAL_KPP is used, computation in overlap regions  c     When option FRUGAL_KPP is used, computation in overlap regions
# Line 305  c     horizontally smooth density relate Line 305  c     horizontally smooth density relate
305    
306  c     zero out dbloc over land points (so that the convective  c     zero out dbloc over land points (so that the convective
307  c     part of the interior mixing can be diagnosed)  c     part of the interior mixing can be diagnosed)
308                 dbloc(i,j,k) = dbloc(i,j,k) * pMask(i,j,k,bi,bj)                 dbloc(i,j,k) = dbloc(i,j,k) * maskC(i,j,k,bi,bj)
309                 ghat(i,j,k)  = ghat(i,j,k)  * pMask(i,j,k,bi,bj)                 ghat(i,j,k)  = ghat(i,j,k)  * maskC(i,j,k,bi,bj)
310                 Ritop(i,j,k) = Ritop(i,j,k) * pMask(i,j,k,bi,bj)                 Ritop(i,j,k) = Ritop(i,j,k) * maskC(i,j,k,bi,bj)
311                 if(k.eq.nzmax(i,j,bi,bj)) then                 if(k.eq.nzmax(i,j,bi,bj)) then
312                    dbloc(i,j,k) = p0                    dbloc(i,j,k) = p0
313                    ghat(i,j,k)  = p0                    ghat(i,j,k)  = p0
# Line 335  cph) Line 335  cph)
335  c------------------------------------------------------------------------  c------------------------------------------------------------------------
336  c     friction velocity, turbulent and radiative surface buoyancy forcing  c     friction velocity, turbulent and radiative surface buoyancy forcing
337  c     -------------------------------------------------------------------  c     -------------------------------------------------------------------
338  c     taux / rho = SurfaceTendencyU * drF(1)                     (N/m^2)  c     taux / rho = surfaceForcingU                               (N/m^2)
339  c     tauy / rho = SurfaceTendencyV * drF(1)                     (N/m^2)  c     tauy / rho = surfaceForcingV                               (N/m^2)
340  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)
341  c     bo    = - g * ( alpha*SurfaceTendencyT +  c     bo    = - g * ( alpha*surfaceForcingT +
342  c                     beta *SurfaceTendencyS ) * drF(1) / rho  (m^2/s^3)  c                     beta *surfaceForcingS ) / rho            (m^2/s^3)
343  c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)  c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)
344  c------------------------------------------------------------------------  c------------------------------------------------------------------------
345    
# Line 357  c initialize arrays to zero Line 357  c initialize arrays to zero
357         DO i = imin, imax         DO i = imin, imax
358          ip1 = i+1          ip1 = i+1
359          work3(i,j) =          work3(i,j) =
360       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) *       &   (surfaceForcingU(i,j,bi,bj) + surfaceForcingU(ip1,j,bi,bj)) *
361       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) +       &   (surfaceForcingU(i,j,bi,bj) + surfaceForcingU(ip1,j,bi,bj)) +
362       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) *       &   (surfaceForcingV(i,j,bi,bj) + surfaceForcingV(i,jp1,bi,bj)) *
363       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj))       &   (surfaceForcingV(i,j,bi,bj) + surfaceForcingV(i,jp1,bi,bj))
364         END DO         END DO
365        END DO        END DO
366  cph(  cph(
# Line 371  cph) Line 371  cph)
371         DO i = imin, imax         DO i = imin, imax
372          ip1 = i+1          ip1 = i+1
373    
374          if ( work3(i,j) .lt. (phepsi*phepsi) ) then          if ( work3(i,j) .lt. (phepsi*phepsi*drF(1)*drF(1)) ) then
375             ustar(i,j) = SQRT( phepsi * p5 * drF(1) )             ustar(i,j) = SQRT( phepsi * p5 * drF(1) )
376          else          else
377             tempVar2 =  SQRT( work3(i,j) ) * p5 * drF(1)             tempVar2 =  SQRT( work3(i,j) ) * p5
378             ustar(i,j) = SQRT( tempVar2 )             ustar(i,j) = SQRT( tempVar2 )
379          endif          endif
380    
381          bo(I,J) = - gravity *          bo(I,J) = - gravity *
382       &       ( vddiff(I,J,1,1) * (SurfaceTendencyT(i,j,bi,bj)+       &       ( vddiff(I,J,1,1) * (surfaceForcingT(i,j,bi,bj)+
383       &       SurfaceTendencyTice(i,j,bi,bj)) +       &       surfaceForcingTice(i,j,bi,bj)) +
384       &       vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj) ) *       &       vddiff(I,J,1,2) * surfaceForcingS(i,j,bi,bj) )
385       &       drF(1) / work2(I,J)       &       / work2(I,J)
386    
387          bosol(I,J) = gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *          bosol(I,J) = gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *
388       &       recip_Cp*recip_rhoConst*recip_dRf(1) *       &       recip_Cp*recip_rhoConst
389       &       drF(1) / work2(I,J)       &       / work2(I,J)
390    
391         END DO         END DO
392        END DO        END DO
# Line 486  c     Estimate reference velocity uRef a Line 486  c     Estimate reference velocity uRef a
486                 vRef(i,j) = p5 *                 vRef(i,j) = p5 *
487       &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )       &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )
488                 IF ( zRef(i,j) .LT. drF(1) ) THEN                 IF ( zRef(i,j) .LT. drF(1) ) THEN
489                    ustarX = ( SurfaceTendencyU(i,  j,bi,bj) +                    ustarX = ( surfaceForcingU(i,  j,bi,bj) +
490       &                       SurfaceTendencyU(ip1,j,bi,bj) ) * p5       &                       surfaceForcingU(ip1,j,bi,bj) ) * p5
491                    ustarY = ( SurfaceTendencyV(i,j,  bi,bj) +       &                   *recip_drF(1)
492       &                       SurfaceTendencyU(i,jp1,bi,bj) ) * p5                    ustarY = ( surfaceForcingV(i,j,  bi,bj) +
493         &                       surfaceForcingU(i,jp1,bi,bj) ) * p5
494         &                   *recip_drF(1)
495                    tempVar1 = ustarX * ustarX + ustarY * ustarY                    tempVar1 = ustarX * ustarX + ustarY * ustarY
496                    if ( tempVar1 .lt. (epsln*epsln) ) then                    if ( tempVar1 .lt. (epsln*epsln) ) then
497                       tempVar2 = epsln                       tempVar2 = epsln
# Line 665  c--------------------------------------- Line 667  c---------------------------------------
667        DO j = jmin, jmax        DO j = jmin, jmax
668         DO i = imin, imax         DO i = imin, imax
669          DO k = 1, Nr          DO k = 1, Nr
670           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * pMask(i,j,k,bi,bj)           KPPviscAz(i,j,k,bi,bj) = vddiff(i,j,k-1,1) * maskC(i,j,k,bi,bj)
671           KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * pMask(i,j,k,bi,bj)           KPPdiffKzS(i,j,k,bi,bj)= vddiff(i,j,k-1,2) * maskC(i,j,k,bi,bj)
672           KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * pMask(i,j,k,bi,bj)           KPPdiffKzT(i,j,k,bi,bj)= vddiff(i,j,k-1,3) * maskC(i,j,k,bi,bj)
673           KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * pMask(i,j,k,bi,bj)           KPPghat(i,j,k,bi,bj)   = ghat(i,j,k)       * maskC(i,j,k,bi,bj)
674          END DO          END DO
675          KPPhbl(i,j,bi,bj) = hbl(i,j) * pMask(i,j,1,bi,bj)          KPPhbl(i,j,bi,bj) = hbl(i,j) * maskC(i,j,1,bi,bj)
676         END DO         END DO
677        END DO        END DO
678  #ifdef FRUGAL_KPP  #ifdef FRUGAL_KPP
Legend:
Removed from v.1.27  
changed lines
  Added in v.1.28
  ViewVC Help
Powered by ViewVC 1.1.22