MODULE wall_fluxes_mod !------------------------------------------------------------------------------! ! Actual revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: wall_fluxes.f90 75 2007-03-22 09:54:05Z raasch $ ! Initial version (2007/03/07) ! ! Description: ! ------------ ! Calculates momentum fluxes at vertical walls assuming Monin-Obukhov ! similarity. ! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0). ! The all-gridpoint version of wall_fluxes_e is not used so far, because ! it gives slightly different results from the ij-version for some unknown ! reason. !------------------------------------------------------------------------------! PRIVATE PUBLIC wall_fluxes, wall_fluxes_e INTERFACE wall_fluxes MODULE PROCEDURE wall_fluxes MODULE PROCEDURE wall_fluxes_ij END INTERFACE wall_fluxes INTERFACE wall_fluxes_e MODULE PROCEDURE wall_fluxes_e MODULE PROCEDURE wall_fluxes_e_ij END INTERFACE wall_fluxes_e CONTAINS !------------------------------------------------------------------------------! ! Call for all grid points !------------------------------------------------------------------------------! SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner, & nzb_uvw_outer, wall ) USE arrays_3d USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, j, k, wall_index INTEGER, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1) :: nzb_uvw_inner, & nzb_uvw_outer REAL :: a, b, c1, c2, h1, h2, zp REAL :: pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts REAL, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1) :: wall REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: 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 ) DO i = nxl, nxr DO j = nys, nyn IF ( wall(j,i) /= 0.0 ) THEN ! !-- All subsequent variables are computed for the respective !-- location where the relevant variable is defined DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i) ! !-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt' 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( 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 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 5.0 * rifs * ( zp - z0(j,i) ) / zp & ) ELSE 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 zp/L (corresponds to neutral Richardson flux !-- number rifs) rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * & ( us_wall**3 + 1E-30 ) ) ! !-- Limit the value range of the Richardson numbers. !-- This is necessary for very small velocities (u,w --> 0), !-- because the absolute value of rif can then become very !-- large, which in consequence would result in very large !-- shear stresses and very small momentum fluxes (both are !-- generally unrealistic). 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,j,i) = 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,j,i) = 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,j,i) = 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,j,i) = -wall_flux(k,j,i) * ABS(wall_flux(k,j,i)) ! !-- store rifs for next time step rif_wall(k,j,i,wall_index) = rifs ENDDO ENDIF ENDDO ENDDO END SUBROUTINE wall_fluxes !------------------------------------------------------------------------------! ! Call for all grid point i,j !------------------------------------------------------------------------------! SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 ) USE arrays_3d USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE 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 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, u_i, v_i, ws, pt' and w'pt' 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( 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 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 5.0 * rifs * ( zp - z0(j,i) ) / zp & ) ELSE 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 zp/L (corresponds to neutral Richardson flux number !-- rifs) rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) ) ! !-- Limit the value range of the Richardson numbers. !-- This is necessary for very small velocities (u,w --> 0), because !-- the absolute value of rif can then become very large, which in !-- consequence would result in very large shear stresses and very !-- small momentum fluxes (both are generally unrealistic). 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( 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_ij !------------------------------------------------------------------------------! ! Call for all grid points !------------------------------------------------------------------------------! SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall ) !------------------------------------------------------------------------------! ! 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 IMPLICIT NONE INTEGER :: i, j, k, kk, wall_index REAL :: a, b, c1, c2, h1, h2, vel_zp, zp REAL :: rifs REAL, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1) :: wall REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: 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 ) DO i = nxl, nxr DO j = nys, nyn IF ( wall(j,i) /= 0.0 ) THEN ! !-- All subsequent variables are computed for the respective !-- location where the relevant variable is defined DO k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2 ! !-- (1) Compute rifs IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN kk = nzb_diff_s_inner(j,i)-1 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') vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + & b * ( v(k,j,i) + v(k,j+1,i) ) + & (c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) & ) IF ( rifs >= 0.0 ) THEN ! !-- Stable stratification (and neutral) wall_flux(k,j,i) = kappa * vel_zp / & ( 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,j,i) = kappa * vel_zp / & ( LOG( zp / z0(j,i) ) + & 5.0 * rifs * ( zp - z0(j,i) ) / zp & ) ELSE wall_flux(k,j,i) = kappa * vel_zp / & ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) & + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & ) ENDIF ENDIF wall_flux(k,j,i) = wall_flux(k,j,i) * ABS( wall_flux(k,j,i) ) ! !-- Store rifs for next time step rif_wall(k,j,i,wall_index) = rifs ENDDO ENDIF ENDDO ENDDO END SUBROUTINE wall_fluxes_e !------------------------------------------------------------------------------! ! Call for grid point i,j !------------------------------------------------------------------------------! SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 ) USE arrays_3d USE control_parameters USE grid_variables USE indices USE statistics IMPLICIT NONE INTEGER :: i, j, k, kk, nzb_w, nzt_w, wall_index REAL :: a, b, c1, c2, h1, h2, vel_zp, 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') vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + & b * ( v(k,j,i) + v(k,j+1,i) ) + & (c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) & ) IF ( rifs >= 0.0 ) THEN ! !-- Stable stratification (and neutral) wall_flux(k) = kappa * vel_zp / & ( 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 * vel_zp / & ( LOG( zp / z0(j,i) ) + & 5.0 * rifs * ( zp - z0(j,i) ) / zp & ) ELSE wall_flux(k) = kappa * vel_zp / & ( 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_ij END MODULE wall_fluxes_mod