/[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.1 by adcroft, Wed Jun 21 19:45:47 2000 UTC revision 1.2 by heimbach, Tue Sep 12 18:14:32 2000 UTC
# Line 83  c     KPPfrac     - Fraction of short-wa Line 83  c     KPPfrac     - Fraction of short-wa
83    
84  c--   KPP_CALC computes vertical viscosity and diffusivity for region  c--   KPP_CALC computes vertical viscosity and diffusivity for region
85  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
86  c     values of uVel, vVel, fu, fv in the region (-2:sNx+4,-2:sNy+4).  c     values of uVel, vVel, SurfaceTendencyU, SurfaceTendencyV in the
87    c     region (-2:sNx+4,-2:sNy+4).
88  c     Hence overlap region needs to be set OLx=4, OLy=4.  c     Hence overlap region needs to be set OLx=4, OLy=4.
89  c     When option FRUGAL_KPP is used, computation in overlap regions  c     When option FRUGAL_KPP is used, computation in overlap regions
90  c     is replaced with exchange calls hence reducing overlap requirements  c     is replaced with exchange calls hence reducing overlap requirements
# Line 101  c     to OLx=1, OLy=1. Line 102  c     to OLx=1, OLy=1.
102  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
103  #include "tamc.h"  #include "tamc.h"
104  #include "tamc_keys.h"  #include "tamc_keys.h"
105          INTEGER    isbyte
106          PARAMETER( isbyte = 4 )
107  #else /* ALLOW_AUTODIFF_TAMC */  #else /* ALLOW_AUTODIFF_TAMC */
108        integer ikey        integer ikey
109  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 192  c     mixing process switches Line 195  c     mixing process switches
195        logical    lri        logical    lri
196        parameter( lri = .true. )        parameter( lri = .true. )
197    
198          _RS        m1
199          parameter( m1=-1.0)
200        _RS        p0    , p5    , p25     , p125      , p0625        _RS        p0    , p5    , p25     , p125      , p0625
201        parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )        parameter( p0=0.0, p5=0.5, p25=0.25, p125=0.125, p0625=0.0625 )
202                
203        _RS     tempVar        _RL     tempVar
204        integer i, j, k, kp1, im1, ip1, jm1, jp1        integer i, j, k, kp1, im1, ip1, jm1, jp1
205    
206  #ifdef KPP_ESTIMATE_UREF  #ifdef KPP_ESTIMATE_UREF
207        _RS     dBdz1, dBdz2, ustarX, ustarY        _RS     dBdz1, dBdz2, ustarX, ustarY
208  #endif  #endif
209    
       IF (use_KPPmixing) THEN  
   
 CADJ STORE fu     (:,:  ,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE fv     (:,:  ,bi,bj)  = uvtape, key = ikey, byte = isbyte  
   
210  c     Check to see if new vertical mixing coefficient should be computed now?  c     Check to see if new vertical mixing coefficient should be computed now?
211        IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,myTime-deltaTClock) .OR.        IF ( DIFFERENT_MULTIPLE(kpp_freq,myTime,myTime-deltaTClock) .OR.
212       1     myTime .EQ. startTime ) THEN       1     myTime .EQ. startTime ) THEN
# Line 343  c     so that the subroutine "bldepth" w Line 343  c     so that the subroutine "bldepth" w
343  c------------------------------------------------------------------------  c------------------------------------------------------------------------
344  c     friction velocity, turbulent and radiative surface buoyancy forcing  c     friction velocity, turbulent and radiative surface buoyancy forcing
345  c     -------------------------------------------------------------------  c     -------------------------------------------------------------------
346  c     taux / rho = fu * delZ(1)                                   (N/m^2)  c     taux / rho = SurfaceTendencyU * delZ(1)                     (N/m^2)
347  c     tauy / rho = fv * delZ(1)                                   (N/m^2)  c     tauy / rho = SurfaceTendencyV * delZ(1)                     (N/m^2)
348  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                 (m/s)  c     ustar = sqrt( sqrt( taux^2 + tauy^2 ) / rho )                 (m/s)
349  c     bo    = - g * (alpha*Qnet + beta*EmPmR) * delZ(1) / rho   (m^2/s^3)  c     bo    = - g * ( alpha*SurfaceTendencyT +
350    c                     beta *SurfaceTendencyS ) * delZ(1) / rho  (m^2/s^3)
351  c     bosol = - g * alpha * Qsw * delZ(1) / rho                 (m^2/s^3)  c     bosol = - g * alpha * Qsw * delZ(1) / rho                 (m^2/s^3)
352  c------------------------------------------------------------------------  c------------------------------------------------------------------------
353    
# Line 356  c--------------------------------------- Line 357  c---------------------------------------
357           DO i = 1, sNx           DO i = 1, sNx
358  #else  #else
359        DO j = jmin, jmax        DO j = jmin, jmax
360           jp1 = j + 1         jp1 = j + 1
361           DO i = imin, imax         DO i = imin, imax
362  #endif  #endif
363              ip1 = i+1          ip1 = i+1
364              ustar(i,j) = SQRT( SQRT(          tempVar =
365       &           (fu(i,j,bi,bj) + fu(ip1,j,bi,bj)) * p5 * delZ(1) *       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) *
366       &           (fu(i,j,bi,bj) + fu(ip1,j,bi,bj)) * p5 * delZ(1) +       &   (SurfaceTendencyU(i,j,bi,bj) + SurfaceTendencyU(ip1,j,bi,bj)) +
367       &           (fv(i,j,bi,bj) + fv(i,jp1,bi,bj)) * p5 * delZ(1) *       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj)) *
368       &           (fv(i,j,bi,bj) + fv(i,jp1,bi,bj)) * p5 * delZ(1) ))       &   (SurfaceTendencyV(i,j,bi,bj) + SurfaceTendencyV(i,jp1,bi,bj))
369              bo(I,J) = - gravity *          if ( tempVar .lt. (epsln*epsln) ) then
370       &           ( vddiff(I,J,1,1) * Qnet(i,j,bi,bj) +             ustar(i,j) = SQRT( epsln * p5 * delZ(1) )
371       &             vddiff(I,J,1,2) * EmPmR(i,j,bi,bj)          else
372       &           ) *             ustar(i,j) = SQRT( SQRT( tempVar ) * p5 * delZ(1) )
373       &           delZ(1) / work2(I,J)          endif
374              bosol(I,J) = - gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *          bo(I,J) = - gravity *
375       &           delZ(1) / work2(I,J)       &       ( vddiff(I,J,1,1) * SurfaceTendencyT(i,j,bi,bj) +
376           END DO       &         vddiff(I,J,1,2) * SurfaceTendencyS(i,j,bi,bj)
377         &       ) *
378         &       delZ(1) / work2(I,J)
379            bosol(I,J) = - gravity * vddiff(I,J,1,1) * Qsw(i,j,bi,bj) *
380         &       delZ(1) / work2(I,J)
381           END DO
382        END DO        END DO
383    
384  #ifndef FRUGAL_KPP  #ifndef FRUGAL_KPP
# Line 459  c     Linearly interpolate to find hMix. Line 465  c     Linearly interpolate to find hMix.
465              ENDIF              ENDIF
466    
467  c     Compute roughness length scale z0 subject to 0 < z0  c     Compute roughness length scale z0 subject to 0 < z0
468                 tempVar = SQRT ( p5 * (                 tempVar = p5 * (
469       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) *
470       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +       &              (uVel(i,  j,  1,bi,bj)-uVel(i,  j,  2,bi,bj)) +
471       &              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) *       &              (uVel(ip1,j,  1,bi,bj)-uVel(ip1,j,  2,bi,bj)) *
# Line 467  c     Compute roughness length scale z0 Line 473  c     Compute roughness length scale z0
473       &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) *       &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) *
474       &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +       &              (vVel(i,  j,  1,bi,bj)-vVel(i,  j,  2,bi,bj)) +
475       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) *
476       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) ) )       &              (vVel(i,  jp1,1,bi,bj)-vVel(i,  jp1,2,bi,bj)) )
477                   if ( tempVar .lt. (epsln*epsln) ) then
478                      tempVar = epsln
479                   else
480                      tempVar = SQRT ( tempVar )
481                   endif
482                 z0(i,j) = rF(2) *                 z0(i,j) = rF(2) *
483       &                   ( rF(3) * LOG ( rF(3) / rF(2) ) /       &                   ( rF(3) * LOG ( rF(3) / rF(2) ) /
484       &                     ( rF(3) - rF(2) ) -       &                     ( rF(3) - rF(2) ) -
# Line 485  c     Estimate reference velocity uRef a Line 496  c     Estimate reference velocity uRef a
496                 vRef(i,j) = p5 *                 vRef(i,j) = p5 *
497       &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )       &                     ( vVel(i,j,1,bi,bj) + vVel(i,jp1,1,bi,bj) )
498                 IF ( zRef(i,j) .LT. drF(1) ) THEN                 IF ( zRef(i,j) .LT. drF(1) ) THEN
499                    ustarX = ( fu(i,j,bi,bj) + fu(ip1,j,bi,bj) ) * p5                    ustarX = ( SurfaceTendencyU(i,  j,bi,bj) +
500                    ustarY = ( fv(i,j,bi,bj) + fu(i,jp1,bi,bj) ) * p5       &                       SurfaceTendencyU(ip1,j,bi,bj) ) * p5
501                    tempVar = MAX ( phepsi,                    ustarY = ( SurfaceTendencyV(i,j,  bi,bj) +
502       &                 SQRT ( ustarX * ustarX + ustarY * ustarY ) )       &                       SurfaceTendencyU(i,jp1,bi,bj) ) * p5
503                      tempVar = ustarX * ustarX + ustarY * ustarY
504                      if ( tempVar .lt. (epsln*epsln) ) then
505                         tempVar = epsln
506                      else
507                         tempVar = SQRT ( tempVar )
508                      endif
509                    tempVar = ustar(i,j) *                    tempVar = ustar(i,j) *
510       &                 ( LOG ( zRef(i,j) / rF(2) ) +       &                 ( LOG ( zRef(i,j) / rF(2) ) +
511       &                 z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) /       &                 z0(i,j) / zRef(i,j) - z0(i,j) / rF(2) ) /
# Line 761  c--------------------------------------- Line 778  c---------------------------------------
778  #endif  #endif
779         ENDIF         ENDIF
780    
781  CADJ STORE vddiff, ghat  = uvtape, key = ikey  #ifdef ALLOW_AUTODIFF_TAMC
782    CADJ STORE vddiff, ghat  = comlev1_kpp, key = ikey
783    #endif /* ALLOW_AUTODIFF_TAMC */
784    
785  c-----------------------------------------------------------------------  c-----------------------------------------------------------------------
786  c     zero out land values,  c     zero out land values,
# Line 847  C     the bottom of the mixing layer. Line 866  C     the bottom of the mixing layer.
866           ENDDO           ENDDO
867        ENDDO        ENDDO
868        CALL SWFRAC(        CALL SWFRAC(
869       I     (sNx+2*OLx)*(sNy+2*OLy), -1., worka,       I     (sNx+2*OLx)*(sNy+2*OLy), m1, worka,
870       O     workb )       O     workb )
871        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
872           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
873              KPPfrac(i,j,bi,bj) = workb(i,j)              KPPfrac(i,j,bi,bj) = workb(i,j)
874           ENDDO           ENDDO
875        ENDDO        ENDDO
           
 CADJ STORE KPPhbl    (:,:  ,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE KPPghat   (:,:,:,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE KPPviscAz (:,:,:,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE KPPdiffKzT(:,:,:,bi,bj)  = uvtape, key = ikey, byte = isbyte  
 CADJ STORE KPPdiffKzS(:,:,:,bi,bj)  = uvtape, key = ikey, byte = isbyte  
876    
877        ENDIF        ENDIF
       ENDIF  
878    
879  #endif ALLOW_KPP  #endif ALLOW_KPP
880    
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22