/[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.40 by mlosch, Tue May 1 04:09:25 2007 UTC revision 1.41 by mlosch, Thu May 3 14:51:05 2007 UTC
# Line 343  CADJ store TTALPHA, SSBETA    = comlev1_ Line 343  CADJ store TTALPHA, SSBETA    = comlev1_
343  #endif  #endif
344  cph)  cph)
345    
 c------------------------------------------------------------------------  
 c     friction velocity, turbulent and radiative surface buoyancy forcing  
 c     -------------------------------------------------------------------  
 c     taux / rho = surfaceForcingU                               (N/m^2)  
 c     tauy / rho = surfaceForcingV                               (N/m^2)  
 c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)  
 c     bo    = - g * ( alpha*surfaceForcingT +  
 c                     beta *surfaceForcingS ) / rho            (m^2/s^3)  
 c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)  
 c------------------------------------------------------------------------  
   
 c initialize arrays to zero  
       DO j = 1-OLy, sNy+OLy  
          DO i = 1-OLx, sNx+OLx  
             ustar(i,j) = p0  
             bo   (I,J) = p0  
             bosol(I,J) = p0  
          END DO  
       END DO  
   
       DO j = jmin, jmax  
        jp1 = j + 1  
        DO i = imin, imax  
         ip1 = i+1  
         work3(i,j) =  
      &   (surfaceForcingU(i,j,bi,bj) + surfaceForcingU(ip1,j,bi,bj)) *  
      &   (surfaceForcingU(i,j,bi,bj) + surfaceForcingU(ip1,j,bi,bj)) +  
      &   (surfaceForcingV(i,j,bi,bj) + surfaceForcingV(i,jp1,bi,bj)) *  
      &   (surfaceForcingV(i,j,bi,bj) + surfaceForcingV(i,jp1,bi,bj))  
        END DO  
       END DO  
 cph(  
 CADJ store work3 = comlev1_kpp, key = ikppkey  
 cph)  
       DO j = jmin, jmax  
        jp1 = j + 1  
        DO i = imin, imax  
         ip1 = i+1  
   
         if ( work3(i,j) .lt. (phepsi*phepsi*drF(1)*drF(1)) ) then  
            ustar(i,j) = SQRT( phepsi * p5 * drF(1) )  
         else  
            tempVar2 =  SQRT( work3(i,j) ) * p5  
            ustar(i,j) = SQRT( tempVar2 )  
         endif  
   
        END DO  
       END DO  
   
346  CML#ifdef ALLOW_SHELFICE  CML#ifdef ALLOW_SHELFICE
347  CMLC     For the pbl parameterisation to work underneath the ice shelves  CMLC     For the pbl parameterisation to work underneath the ice shelves
348  CMLC     it needs to know the surface (ice-ocean) fluxes. However, masking  CMLC     it needs to know the surface (ice-ocean) fluxes. However, masking
349  CMLC     and indexing problems make this part of the code not work  CMLC     and indexing problems make this part of the code not work
350  CMLC     underneath the ice shelves and the following lines are only here  CMLC     underneath the ice shelves and the following lines are only here
351  CMLC     to remind me that this still needs to be sorted out.  CMLC     to remind me that this still needs to be sorted out.
352  CML      IF ( useShelfIce ) THEN  CML      shelfIceFac = 0. _d 0
353    CML      IF ( useShelfIce ) selfIceFac = 1. _d 0
354  CML      DO j = jmin, jmax  CML      DO j = jmin, jmax
 CML       jp1 = j + 1  
355  CML       DO i = imin, imax  CML       DO i = imin, imax
356  CML        ip1 = i+1  CML        surfForcT = surfaceForcingT(i,j,bi,bj)
357  CML        k   = kTopC(I,J,bi,bj)  CML     &       + shelficeForcingT(i,j,bi,bj) * shelfIceFac
358  CML        bo(I,J) = - gravity *  CML        surfForcS = surfaceForcingS(i,j,bi,bj)
359  CML     &       ( TTALPHA(I,J,K) * ( maskC(I,J,1,bi,bj)  CML     &       + shelficeForcingS(i,j,bi,bj) * shelfIceFac
 CML     &       * ( surfaceForcingT(i,j,bi,bj)  
 CML     &       + surfaceForcingTice(i,j,bi,bj) )  
 CML     &       + shelficeForcingT(i,j,bi,bj) )  
 CML     &       + SSBETA(I,J,K) * ( maskC(I,J,1,bi,bj)  
 CML     &       * surfaceForcingS(i,j,bi,bj)  
 CML     &       + shelficeForcingS(i,j,bi,bj) ) )  
 CML     &       / work2(I,J)  
 CML  
 CML        bosol(I,J) = gravity * TTALPHA(I,J,1) * Qsw(i,j,bi,bj)  
 CML     &       * maskC(I,J,1,bi,bj) * recip_Cp*recip_rhoConst  
 CML     &       / work2(I,J)  
 CML  
360  CML       END DO  CML       END DO
361  CML      END DO  CML      END DO
 CML      ELSE  
 CML#endif /* ALLOW_SHELFICE */  
       DO j = jmin, jmax  
        jp1 = j + 1  
        DO i = imin, imax  
         ip1 = i+1  
   
         bo(I,J) = - gravity *  
      &       ( TTALPHA(I,J,1) * (surfaceForcingT(i,j,bi,bj)+  
      &       surfaceForcingTice(i,j,bi,bj)) +  
      &       SSBETA(I,J,1) * surfaceForcingS(i,j,bi,bj) )  
      &       / work2(I,J)  
   
         bosol(I,J) = gravity * TTALPHA(I,J,1) * Qsw(i,j,bi,bj) *  
      &       recip_Cp*recip_rhoConst  
      &       / work2(I,J)  
   
        END DO  
       END DO  
 CML#ifdef ALLOW_SHELFICE  
 CML      ENDIF  
362  CML#endif /* ALLOW_SHELFICE */  CML#endif /* ALLOW_SHELFICE */
         
 cph(  
 CADJ store ustar = comlev1_kpp, key = ikppkey  
 cph)  
363    
364  c------------------------------------------------------------------------  c------------------------------------------------------------------------
365    c     friction velocity, turbulent and radiative surface buoyancy forcing
366    c     -------------------------------------------------------------------
367    c     taux / rho = surfaceForcingU                               (N/m^2)
368    c     tauy / rho = surfaceForcingV                               (N/m^2)
369    c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                (m/s)
370    c     bo    = - g * ( alpha*surfaceForcingT +
371    c                     beta *surfaceForcingS ) / rho            (m^2/s^3)
372    c     bosol = - g * alpha * Qsw * drF(1) / rho                 (m^2/s^3)
373    c------------------------------------------------------------------------
374  c     velocity shear  c     velocity shear
375  c     --------------  c     --------------
376  c     Get velocity shear squared, averaged from "u,v-grid"  c     Get velocity shear squared, averaged from "u,v-grid"
377  c     onto "t-grid" (in (m/s)**2):  c     onto "t-grid" (in (m/s)**2):
378  c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2      at grid levels  c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2      at grid levels
379  c     shsq(k)=(U(k)-U(k+1))**2+(V(k)-V(k+1))**2  at interfaces  c     shsq(k)=(U(k)-U(k+1))**2+(V(k)-V(k+1))**2  at interfaces
380    c    
381    c     note: Vref can depend on the surface fluxes that is why we compute
382    c     dVsq in the subroutine that does the surface related stuff
383    c     (admittedly this is a bit messy)
384  c------------------------------------------------------------------------  c------------------------------------------------------------------------
385          
386          CALL KPP_FORCING_SURF(
387         I     work2, surfaceForcingU, surfaceForcingV,
388         I     surfaceForcingT, surfaceForcingS, surfaceForcingTice,
389         I     Qsw, ttalpha, ssbeta,  
390         O     ustar, bo, bosol, dVsq,
391         I     ikppkey, iMin, iMax, jMin, jMax, bi, bj, myTime, myThid )
392    
393    CMLcph(
394    CMLCADJ store ustar = comlev1_kpp, key = ikppkey
395    CMLcph)
396    
397  c initialize arrays to zero  c initialize arrays to zero
398        DO k = 1, Nr        DO k = 1, Nr
399           DO j = 1-OLy, sNy+OLy         DO j = 1-OLy, sNy+OLy
400              DO i = 1-OLx, sNx+OLx          DO i = 1-OLx, sNx+OLx
401                 shsq(i,j,k) = p0           shsq(i,j,k) = p0
402                 dVsq(i,j,k) = p0          END DO
403              END DO         END DO
          END DO  
       END DO  
   
 c     dVsq computation  
   
 #ifdef KPP_ESTIMATE_UREF  
   
 c     Get rid of vertical resolution dependence of dVsq term by  
 c     estimating a surface velocity that is independent of first level  
 c     thickness in the model.  First determine mixed layer depth hMix.  
 c     Second zRef = espilon * hMix.  Third determine roughness length  
 c     scale z0.  Third estimate reference velocity.  
   
       DO j = jmin, jmax  
          jp1 = j + 1  
          DO i = imin, imax  
             ip1 = i + 1  
   
 c     Determine mixed layer depth hMix as the shallowest depth at which  
 c     dB/dz exceeds 5.2e-5 s^-2.  
             work1(i,j) = nzmax(i,j,bi,bj)  
             DO k = 1, Nr  
                IF ( k .LT. nzmax(i,j,bi,bj) .AND.  
      &              maskC(I,J,k,bi,bj) .GT. 0. .AND.  
      &              dbloc(i,j,k) / drC(k+1) .GT. dB_dz )  
      &              work1(i,j) = k  
             END DO  
   
 c     Linearly interpolate to find hMix.  
             k = work1(i,j)  
             IF ( k .EQ. 0 .OR. nzmax(i,j,bi,bj) .EQ. 1 ) THEN  
                zRef(i,j) = p0  
             ELSEIF ( k .EQ. 1) THEN  
                dBdz2 = dbloc(i,j,1) / drC(2)  
                zRef(i,j) = drF(1) * dB_dz / dBdz2  
             ELSEIF ( k .LT. nzmax(i,j,bi,bj) ) THEN  
                dBdz1 = dbloc(i,j,k-1) / drC(k  )  
                dBdz2 = dbloc(i,j,k  ) / drC(k+1)  
                zRef(i,j) = rF(k) + drF(k) * (dB_dz - dBdz1) /  
      &                     MAX ( phepsi, dBdz2 - dBdz1 )  
             ELSE  
                zRef(i,j) = rF(k+1)  
             ENDIF  
   
 c     Compute roughness length scale z0 subject to 0 < z0  
                tempVar1 = p5 * (  
      &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *  
      &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +  
      &              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) *  
      &              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) +  
      &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) *  
      &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +  
      &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *  
      &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) )  
                if ( tempVar1 .lt. (epsln*epsln) ) then  
                   tempVar2 = epsln  
                else  
                   tempVar2 = SQRT ( tempVar1 )  
                endif  
                z0(i,j) = rF(2) *  
      &                   ( rF(3) * LOG ( rF(3) / rF(2) ) /  
      &                     ( rF(3) - rF(2) ) -  
      &                     tempVar2 * vonK /  
      &                     MAX ( ustar(i,j), phepsi ) )  
                z0(i,j) = MAX ( z0(i,j), phepsi )  
   
 c     zRef is set to 0.1 * hMix subject to z0 <= zRef <= drF(1)  
                zRef(i,j) = MAX ( epsilon * zRef(i,j), z0(i,j) )  
                zRef(i,j) = MIN ( zRef(i,j), drF(1) )  
   
 c     Estimate reference velocity uRef and vRef.  
                uRef(i,j) = p5 *  
      &                     ( uVel(i,j,1,bi,bj) + uVel(ip1,j,1,bi,bj) )  
                vRef(i,j) = p5 *  
      &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )  
                IF ( zRef(i,j) .LT. drF(1) ) THEN  
                   ustarX = ( surfaceForcingU(i,  j,bi,bj) +  
      &                       surfaceForcingU(ip1,j,bi,bj) ) * p5  
      &                   *recip_drF(1)  
                   ustarY = ( surfaceForcingV(i,j,  bi,bj) +  
      &                       surfaceForcingV(i,jp1,bi,bj) ) * p5  
      &                   *recip_drF(1)  
                   tempVar1 = ustarX * ustarX + ustarY * ustarY  
                   if ( tempVar1 .lt. (epsln*epsln) ) then  
                      tempVar2 = epsln  
                   else  
                      tempVar2 = SQRT ( tempVar1 )  
                   endif  
                   tempVar2 = ustar(i,j) *  
      &                 ( LOG ( zRef(i,j) / rF(2) ) +  
      &                 z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) /  
      &                 vonK / tempVar2  
                   uRef(i,j) = uRef(i,j) + ustarX * tempVar2  
                   vRef(i,j) = vRef(i,j) + ustarY * tempVar2  
                ENDIF  
   
          END DO  
       END DO  
   
       DO k = 1, Nr  
          DO j = jmin, jmax  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = imin, imax  
                im1 = i - 1  
                ip1 = i + 1  
                dVsq(i,j,k) = p5 * (  
      $              (uRef(i,j) - uVel(i,  j,  k,bi,bj)) *  
      $              (uRef(i,j) - uVel(i,  j,  k,bi,bj)) +  
      $              (uRef(i,j) - uVel(ip1,j,  k,bi,bj)) *  
      $              (uRef(i,j) - uVel(ip1,j,  k,bi,bj)) +  
      $              (vRef(i,j) - vVel(i,  j,  k,bi,bj)) *  
      $              (vRef(i,j) - vVel(i,  j,  k,bi,bj)) +  
      $              (vRef(i,j) - vVel(i,  jp1,k,bi,bj)) *  
      $              (vRef(i,j) - vVel(i,  jp1,k,bi,bj)) )  
 #ifdef KPP_SMOOTH_DVSQ  
                dVsq(i,j,k) = p5 * dVsq(i,j,k) + p125 * (  
      $              (uRef(i,j) - uVel(i,  jm1,k,bi,bj)) *  
      $              (uRef(i,j) - uVel(i,  jm1,k,bi,bj)) +  
      $              (uRef(i,j) - uVel(ip1,jm1,k,bi,bj)) *  
      $              (uRef(i,j) - uVel(ip1,jm1,k,bi,bj)) +  
      $              (uRef(i,j) - uVel(i,  jp1,k,bi,bj)) *  
      $              (uRef(i,j) - uVel(i,  jp1,k,bi,bj)) +  
      $              (uRef(i,j) - uVel(ip1,jp1,k,bi,bj)) *  
      $              (uRef(i,j) - uVel(ip1,jp1,k,bi,bj)) +  
      $              (vRef(i,j) - vVel(im1,j,  k,bi,bj)) *  
      $              (vRef(i,j) - vVel(im1,j,  k,bi,bj)) +  
      $              (vRef(i,j) - vVel(im1,jp1,k,bi,bj)) *  
      $              (vRef(i,j) - vVel(im1,jp1,k,bi,bj)) +  
      $              (vRef(i,j) - vVel(ip1,j,  k,bi,bj)) *  
      $              (vRef(i,j) - vVel(ip1,j,  k,bi,bj)) +  
      $              (vRef(i,j) - vVel(ip1,jp1,k,bi,bj)) *  
      $              (vRef(i,j) - vVel(ip1,jp1,k,bi,bj)) )  
 #endif /* KPP_SMOOTH_DVSQ */  
             END DO  
          END DO  
       END DO  
   
 #else /* KPP_ESTIMATE_UREF */  
   
       DO k = 1, Nr  
          DO j = jmin, jmax  
             jm1 = j - 1  
             jp1 = j + 1  
             DO i = imin, imax  
                im1 = i - 1  
                ip1 = i + 1  
                dVsq(i,j,k) = p5 * (  
      $              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  k,bi,bj)) *  
      $              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  k,bi,bj)) +  
      $              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  k,bi,bj)) *  
      $              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  k,bi,bj)) +  
      $              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  k,bi,bj)) *  
      $              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  k,bi,bj)) +  
      $              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,k,bi,bj)) *  
      $              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,k,bi,bj)) )  
 #ifdef KPP_SMOOTH_DVSQ  
                dVsq(i,j,k) = p5 * dVsq(i,j,k) + p125 * (  
      $              (uVel(i,  jm1,1,bi,bj)-uVel(i,  jm1,k,bi,bj)) *  
      $              (uVel(i,  jm1,1,bi,bj)-uVel(i,  jm1,k,bi,bj)) +  
      $              (uVel(ip1,jm1,1,bi,bj)-uVel(ip1,jm1,k,bi,bj)) *  
      $              (uVel(ip1,jm1,1,bi,bj)-uVel(ip1,jm1,k,bi,bj)) +  
      $              (uVel(i,  jp1,1,bi,bj)-uVel(i,  jp1,k,bi,bj)) *  
      $              (uVel(i,  jp1,1,bi,bj)-uVel(i,  jp1,k,bi,bj)) +  
      $              (uVel(ip1,jp1,1,bi,bj)-uVel(ip1,jp1,k,bi,bj)) *  
      $              (uVel(ip1,jp1,1,bi,bj)-uVel(ip1,jp1,k,bi,bj)) +  
      $              (vVel(im1,j,  1,bi,bj)-vVel(im1,j,  k,bi,bj)) *  
      $              (vVel(im1,j,  1,bi,bj)-vVel(im1,j,  k,bi,bj)) +  
      $              (vVel(im1,jp1,1,bi,bj)-vVel(im1,jp1,k,bi,bj)) *  
      $              (vVel(im1,jp1,1,bi,bj)-vVel(im1,jp1,k,bi,bj)) +  
      $              (vVel(ip1,j,  1,bi,bj)-vVel(ip1,j,  k,bi,bj)) *  
      $              (vVel(ip1,j,  1,bi,bj)-vVel(ip1,j,  k,bi,bj)) +  
      $              (vVel(ip1,jp1,1,bi,bj)-vVel(ip1,jp1,k,bi,bj)) *  
      $              (vVel(ip1,jp1,1,bi,bj)-vVel(ip1,jp1,k,bi,bj)) )  
 #endif /* KPP_SMOOTH_DVSQ */  
             END DO  
          END DO  
404        END DO        END DO
405    
 #endif /* KPP_ESTIMATE_UREF */  
   
406  c     shsq computation  c     shsq computation
407        DO k = 1, Nrm1        DO k = 1, Nrm1
408           kp1 = k + 1         kp1 = k + 1
409           DO j = jmin, jmax         DO j = jmin, jmax
410              jm1 = j - 1          jm1 = j - 1
411              jp1 = j + 1          jp1 = j + 1
412              DO i = imin, imax          DO i = imin, imax
413                 im1 = i - 1           im1 = i - 1
414                 ip1 = i + 1           ip1 = i + 1
415                 shsq(i,j,k) = p5 * (           shsq(i,j,k) = p5 * (
416       $              (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) *       $        (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) *
417       $              (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) +       $        (uVel(i,  j,  k,bi,bj)-uVel(i,  j,  kp1,bi,bj)) +
418       $              (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) *       $        (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) *
419       $              (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) +       $        (uVel(ip1,j,  k,bi,bj)-uVel(ip1,j,  kp1,bi,bj)) +
420       $              (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) *       $        (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) *
421       $              (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) +       $        (vVel(i,  j,  k,bi,bj)-vVel(i,  j,  kp1,bi,bj)) +
422       $              (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) *       $        (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) *
423       $              (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) )       $        (vVel(i,  jp1,k,bi,bj)-vVel(i,  jp1,kp1,bi,bj)) )
424  #ifdef KPP_SMOOTH_SHSQ  #ifdef KPP_SMOOTH_SHSQ
425                 shsq(i,j,k) = p5 * shsq(i,j,k) + p125 * (           shsq(i,j,k) = p5 * shsq(i,j,k) + p125 * (
426       $              (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) *       $        (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) *
427       $              (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) +       $        (uVel(i,  jm1,k,bi,bj)-uVel(i,  jm1,kp1,bi,bj)) +
428       $              (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) *       $        (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) *
429       $              (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) +       $        (uVel(ip1,jm1,k,bi,bj)-uVel(ip1,jm1,kp1,bi,bj)) +
430       $              (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) *       $        (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) *
431       $              (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) +       $        (uVel(i,  jp1,k,bi,bj)-uVel(i,  jp1,kp1,bi,bj)) +
432       $              (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) *       $        (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) *
433       $              (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) +       $        (uVel(ip1,jp1,k,bi,bj)-uVel(ip1,jp1,kp1,bi,bj)) +
434       $              (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) *       $        (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) *
435       $              (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) +       $        (vVel(im1,j,  k,bi,bj)-vVel(im1,j,  kp1,bi,bj)) +
436       $              (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) *       $        (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) *
437       $              (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) +       $        (vVel(im1,jp1,k,bi,bj)-vVel(im1,jp1,kp1,bi,bj)) +
438       $              (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) *       $        (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) *
439       $              (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) +       $        (vVel(ip1,j,  k,bi,bj)-vVel(ip1,j,  kp1,bi,bj)) +
440       $              (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) *       $        (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) *
441       $              (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) )       $        (vVel(ip1,jp1,k,bi,bj)-vVel(ip1,jp1,kp1,bi,bj)) )
442  #endif  #endif
443              END DO          END DO
444           END DO         END DO
445        END DO        END DO
446    
447  cph(  cph(
Legend:
Removed from v.1.40  
changed lines
  Added in v.1.41
  ViewVC Help
Powered by ViewVC 1.1.22