# Changeset 53 for palm/trunk/SOURCE/wall_fluxes.f90

Ignore:
Timestamp:
Mar 7, 2007 12:33:47 PM (15 years ago)
Message:

preliminary version, several changes to be explained later

File:
1 edited

Unmodified
Removed
• ## palm/trunk/SOURCE/wall_fluxes.f90

 r52 IMPLICIT NONE INTEGER ::  i, j, k, nzb_w, nzt_w REAL    ::  a, b, c1, c2, h1, h2, delta_p INTEGER ::  i, j, k, nzb_w, nzt_w, wall_index REAL    ::  a, b, c1, c2, h1, h2, zp REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts delta_p   = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) wall_flux = 0.0 zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) wall_flux  = 0.0 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) ! !-- the relevant variable is defined DO  k = nzb_w, nzt_w ! !--    (1) Compute rifs, u_i, v_i, ws, pt' and w'pt' rifs    = rif_wall(k,j,i,NINT(a+2*b+3*c1+4*c2)) u_i     = a * u(k,j,i) +  & c1 * 0.25 * ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) v_i     = b * v(k,j,i) +  & c2 * 0.25 * ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) ws      = ( c1 + c2 ) * w(k,j,i) +  & a * 0.25 * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) + & b * 0.25 * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) pt_i    = 0.5 * ( pt(k,j,i) +  & a *  pt(k,j,i-1) + b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) ) pts     = pt_i - hom(k,1,4,0) wspts   = ws * pts rifs  = rif_wall(k,j,i,wall_index) u_i   = a * u(k,j,i) + c1 * 0.25 * & ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) v_i   = b * v(k,j,i) + c2 * 0.25 * & ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                             & a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) & + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) & ) pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + b * pt(k,j-1,i)  & + ( c1 + c2 ) * pt(k+1,j,i) ) pts   = pt_i - hom(k,1,4,0) wspts = ws * pts ! !--    (2) Compute wall-parallel absolute velocity vel_total vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 ) ! !--    (3) Compute wall friction velocity us_wall IF ( rifs >= 0.0 )  THEN ! !--       Stable stratification (and neutral) us_wall = kappa * vel_total / (                         & LOG( delta_p / z0(j,i) ) +                    & 5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p  & us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +               & 5.0 * rifs * ( zp - z0(j,i) ) / zp  & ) ELSE ! !--       Unstable stratification h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / delta_p * z0(j,i) ) ) h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) ! !--       If a borderline case occurs, the formula for stable stratification !--       argument of the logarithm. IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN us_wall = kappa * vel_total / (                           & LOG( delta_p / z0(j,i) ) +                         & 5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p & ) us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             & 5.0 * rifs * ( zp - z0(j,i) ) / zp & ) ELSE us_wall = kappa * vel_total / (                           & us_wall = kappa * vel_total / (                                   & LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + & 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                  & ) ) ENDIF ENDIF ! !--    (4) Compute delta_p/L (corresponds to neutral Richardson flux number ENDIF ! !--    (4) Compute zp/L (corresponds to neutral Richardson flux number !--        rifs) rifs = -1.0 * delta_p * kappa * g * wspts / & ( pt_i * ( us_wall**3 + 1E-30 ) ) rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * ( us_wall**3 + 1E-30 ) ) ! !--    Limit the value range of the Richardson numbers. IF ( rifs < rif_min )  rifs = rif_min IF ( rifs > rif_max )  rifs = rif_max ! !--    (5) Compute wall_flux (u'v', v'u', w'v', or w'u') IF ( rifs >= 0.0 )  THEN ! !--       Stable stratification (and neutral) wall_flux(k) = kappa *  & ( a * u(k,j,i) + b * v(k,j,i) + (c1 + c2 ) * w(k,j,i) ) / (  & LOG( delta_p / z0(j,i) ) +                                & 5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p        & ) ELSE wall_flux(k) = kappa *                                          & ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & (  LOG( zp / z0(j,i) ) +                         & 5.0 * rifs * ( zp - z0(j,i) ) / zp            & ) ELSE ! !--       Unstable stratification h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / delta_p * z0(j,i) ) ) h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) ! !--       If a borderline case occurs, the formula for stable stratification !--       argument of the logarithm. IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN wall_flux(k) = kappa *  & ( a * u(k,j,i) + b * v(k,j,i) + (c1 + c2 ) * w(k,j,i) ) / ( & LOG( delta_p / z0(j,i) ) +                                & 5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p        & ) wall_flux(k) = kappa *                                          & ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & ( LOG( zp / z0(j,i) ) +                          & 5.0 * rifs * ( zp - z0(j,i) ) / zp             & ) ELSE wall_flux(k) = kappa *  & ( a * u(k,j,i) + b * v(k,j,i) + (c1 + c2 ) * w(k,j,i) ) / (  & LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) +          & 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                            & ) wall_flux(k) = kappa *                                            & ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /   & ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) & + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              & ) ENDIF ENDIF wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) ) ! !--    store rifs for next time step rif_wall(k,j,i,NINT(a+2*b+3*c1+4*c2)) = rifs rif_wall(k,j,i,wall_index) = rifs ENDDO END SUBROUTINE wall_fluxes SUBROUTINE wall_fluxes_e( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 ) !------------------------------------------------------------------------------! ! Actual revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! Initial version (2007/03/07) ! ! Description: ! ------------ ! Calculates momentum fluxes at vertical walls for routine production_e ! assuming Monin-Obukhov similarity. ! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0). !------------------------------------------------------------------------------! USE arrays_3d USE control_parameters USE grid_variables USE indices USE statistics USE user IMPLICIT NONE INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index REAL    ::  a, b, c1, c2, h1, h2, zp REAL ::  rifs REAL, DIMENSION(nzb:nzt+1) ::  wall_flux zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) wall_flux  = 0.0 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) ! !-- All subsequent variables are computed for the respective location where !-- the relevant variable is defined DO  k = nzb_w, nzt_w ! !--    (1) Compute rifs IF ( k == nzb_w )  THEN kk = nzb_w ELSE kk = k-1 ENDIF rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                         & a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + & c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    & ) ! !--    Skip (2) to (4) of wall_fluxes, because here rifs is already available !--    from (1) ! !--    (5) Compute wall_flux (u'v', v'u', w'v', or w'u') IF ( rifs >= 0.0 )  THEN ! !--       Stable stratification (and neutral) wall_flux(k) = kappa *                                          & ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & (  LOG( zp / z0(j,i) ) +                         & 5.0 * rifs * ( zp - z0(j,i) ) / zp            & ) ELSE ! !--       Unstable stratification h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) ! !--       If a borderline case occurs, the formula for stable stratification !--       must be used anyway, or else a zero division would occur in the !--       argument of the logarithm. IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN wall_flux(k) = kappa *                                          & ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & ( LOG( zp / z0(j,i) ) +                          & 5.0 * rifs * ( zp - z0(j,i) ) / zp             & ) ELSE wall_flux(k) = kappa *                                            & ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /   & ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) & + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) )              & ) ENDIF ENDIF wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) ) ! !--    store rifs for next time step rif_wall(k,j,i,wall_index) = rifs ENDDO END SUBROUTINE wall_fluxes_e
Note: See TracChangeset for help on using the changeset viewer.