Changeset 56 for palm/trunk/SOURCE
- Timestamp:
- Mar 8, 2007 1:57:07 PM (18 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 6 edited
-
Makefile (modified) (2 diffs)
-
diffusion_u.f90 (modified) (4 diffs)
-
diffusion_v.f90 (modified) (4 diffs)
-
diffusion_w.f90 (modified) (5 diffs)
-
production_e.f90 (modified) (8 diffs)
-
wall_fluxes.f90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r54 r56 153 153 diffusion_e.o: modules.o 154 154 diffusion_s.o: modules.o 155 diffusion_u.o: modules.o 156 diffusion_v.o: modules.o 157 diffusion_w.o: modules.o 155 diffusion_u.o: modules.o wall_fluxes.o 156 diffusion_v.o: modules.o wall_fluxes.o 157 diffusion_w.o: modules.o wall_fluxes.o 158 158 diffusivities.o: modules.o 159 159 disturb_field.o: modules.o random_function.o … … 193 193 pres.o: modules.o poisfft.o poisfft_hybrid.o 194 194 print_1d.o: modules.o 195 production_e.o: modules.o 195 production_e.o: modules.o wall_fluxes.o 196 196 prognostic_equations.o: modules.o advec_s_pw.o advec_s_up.o advec_u_pw.o \ 197 197 advec_u_up.o advec_v_pw.o advec_v_up.o advec_w_pw.o advec_w_up.o \ -
palm/trunk/SOURCE/diffusion_u.f90
r53 r56 34 34 !------------------------------------------------------------------------------! 35 35 36 USE wall_fluxes_mod 37 36 38 PRIVATE 37 39 PUBLIC diffusion_u … … 61 63 REAL :: z0(nys-1:nyn+1,nxl-1:nxr+1) 62 64 REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) 63 REAL, DIMENSION(nzb:nzt+1) :: usvs64 65 REAL, DIMENSION(:,:), POINTER :: usws 65 66 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w 67 REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr+uxrp) :: usvs 68 69 ! 70 !-- First calculate horizontal momentum flux u'v' at vertical walls, 71 !-- if neccessary 72 IF ( topography /= 'flat' ) THEN 73 CALL wall_fluxes( usvs, 1.0, 0.0, 0.0, 0.0, uxrp, 0, nzb_u_inner, & 74 nzb_u_outer, wall_u ) 75 ENDIF 66 76 67 77 DO i = nxl, nxr+uxrp … … 103 113 !-- Wall functions at the north and south walls, respectively 104 114 IF ( wall_u(j,i) /= 0.0 ) THEN 105 106 !107 !-- Calculate the horizontal momentum flux u'v'108 CALL wall_fluxes( i, j, nzb_u_inner(j,i)+1, nzb_u_outer(j,i), &109 usvs, 1.0, 0.0, 0.0, 0.0 )110 115 111 116 DO k = nzb_u_inner(j,i)+1, nzb_u_outer(j,i) … … 139 144 + kmym_x * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 140 145 ) & 141 + wall_u(j,i) * usvs(k )&146 + wall_u(j,i) * usvs(k,j,i) & 142 147 ) * ddy 143 148 ENDDO -
palm/trunk/SOURCE/diffusion_v.f90
r53 r56 32 32 !------------------------------------------------------------------------------! 33 33 34 USE wall_fluxes_mod 35 34 36 PRIVATE 35 37 PUBLIC diffusion_v … … 59 61 REAL :: z0(nys-1:nyn+1,nxl-1:nxr+1) 60 62 REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) 61 REAL, DIMENSION(nzb:nzt+1) :: vsus62 63 REAL, DIMENSION(:,:), POINTER :: vsws 63 64 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w 65 REAL, DIMENSION(nzb:nzt+1,nys:nyn+vynp,nxl:nxr) :: vsus 66 67 ! 68 !-- First calculate horizontal momentum flux v'u' at vertical walls, 69 !-- if neccessary 70 IF ( topography /= 'flat' ) THEN 71 CALL wall_fluxes( vsus, 0.0, 1.0, 0.0, 0.0, 0, vynp, nzb_v_inner, & 72 nzb_v_outer, wall_v ) 73 ENDIF 64 74 65 75 DO i = nxl, nxr … … 101 111 !-- Wall functions at the left and right walls, respectively 102 112 IF ( wall_v(j,i) /= 0.0 ) THEN 103 104 !105 !-- Calculate the horizontal momentum flux v'u'106 CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), &107 vsus, 0.0, 1.0, 0.0, 0.0 )108 113 109 114 DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i) … … 137 142 + kmxm_y * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 138 143 ) & 139 + wall_v(j,i) * vsus(k )&144 + wall_v(j,i) * vsus(k,j,i) & 140 145 ) * ddx 141 146 ENDDO -
palm/trunk/SOURCE/diffusion_w.f90
r53 r56 28 28 ! Diffusion term of the w-component 29 29 !------------------------------------------------------------------------------! 30 31 USE wall_fluxes_mod 30 32 31 33 PRIVATE … … 59 61 REAL :: z0(nys-1:nyn+1,nxl-1:nxr+1) 60 62 REAL :: tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) 61 REAL, DIMENSION(nzb:nzt+1) :: wsus, wsvs62 63 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w 63 64 REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus, wsvs 65 66 67 ! 68 !-- First calculate horizontal momentum flux w'u' and/or w'v' at vertical 69 !-- walls, if neccessary 70 IF ( topography /= 'flat' ) THEN 71 CALL wall_fluxes( wsus, 0.0, 0.0, 0.0, 1.0, 0, 0, nzb_w_inner, & 72 nzb_w_outer, wall_w_x ) 73 CALL wall_fluxes( wsvs, 0.0, 0.0, 1.0, 0.0, 0, 0, nzb_w_inner, & 74 nzb_w_outer, wall_w_y ) 75 ENDIF 64 76 65 77 DO i = nxl, nxr … … 115 127 IF ( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0 ) THEN 116 128 117 !118 !-- Calculate the horizontal momentum fluxes w'u' and/or w'v'119 IF ( wall_w_x(j,i) /= 0.0 ) THEN120 CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, &121 nzb_w_outer(j,i), wsus, 0.0, 0.0, 0.0, &122 1.0 )123 ELSE124 wsus = 0.0125 ENDIF126 127 IF ( wall_w_y(j,i) /= 0.0 ) THEN128 CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, &129 nzb_w_outer(j,i), wsvs, 0.0, 0.0, 1.0, &130 0.0 )131 ELSE132 wsvs = 0.0133 ENDIF134 135 129 DO k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i) 136 130 ! … … 171 165 + kmxm_z * ( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) & 172 166 ) & 173 + wall_w_x(j,i) * wsus(k )&167 + wall_w_x(j,i) * wsus(k,j,i) & 174 168 ) * ddx & 175 169 + ( fwyp(j,i) * ( & … … 181 175 + kmym_z * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 182 176 ) & 183 + wall_w_y(j,i) * wsvs(k )&177 + wall_w_y(j,i) * wsvs(k,j,i) & 184 178 ) * ddy & 185 179 + 2.0 * ( & -
palm/trunk/SOURCE/production_e.f90
r55 r56 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes 6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes_e 7 7 ! 8 8 ! Former revisions: … … 31 31 !------------------------------------------------------------------------------! 32 32 33 USE wall_fluxes_mod 34 33 35 PRIVATE 34 36 PUBLIC production_e, production_e_init 35 37 36 38 LOGICAL, SAVE :: first_call = .TRUE. 37 39 … … 69 71 k1, k2, theta, temp 70 72 71 REAL, DIMENSION(nzb:nzt+1) :: usvs, vsus, wsus, wsvs 72 73 ! REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs, vsus, wsus, wsvs 74 REAL, DIMENSION(nzb:nzt+1) :: usvs, vsus, wsus, wsvs 75 76 ! 77 !-- First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at 78 !-- vertical walls, if neccessary 79 !-- So far, results are slightly different from the ij-Version. 80 !-- Therefore, ij-Version is called further below within the ij-loops. 81 ! IF ( topography /= 'flat' ) THEN 82 ! CALL wall_fluxes_e( usvs, 1.0, 0.0, 0.0, 0.0, wall_e_y ) 83 ! CALL wall_fluxes_e( wsvs, 0.0, 0.0, 1.0, 0.0, wall_e_y ) 84 ! CALL wall_fluxes_e( vsus, 0.0, 1.0, 0.0, 0.0, wall_e_x ) 85 ! CALL wall_fluxes_e( wsus, 0.0, 0.0, 0.0, 1.0, wall_e_x ) 86 ! ENDIF 73 87 74 88 ! … … 131 145 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 132 146 usvs, 1.0, 0.0, 0.0, 0.0 ) 133 134 147 dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i) 148 ! dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i) 135 149 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 136 150 wsvs, 0.0, 0.0, 1.0, 0.0 ) 137 151 dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i) 152 ! dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i) 138 153 ELSE 139 154 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & … … 147 162 vsus, 0.0, 1.0, 0.0, 0.0 ) 148 163 dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i) 164 ! dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i) 149 165 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 150 166 wsus, 0.0, 0.0, 0.0, 1.0 ) 151 167 dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i) 168 ! dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i) 152 169 ELSE 153 170 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & … … 185 202 IF ( wall_e_y(j,i) /= 0.0 ) THEN 186 203 dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i) 204 ! dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i) 187 205 dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i) 206 ! dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i) 188 207 ELSE 189 208 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & … … 195 214 IF ( wall_e_x(j,i) /= 0.0 ) THEN 196 215 dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i) 216 ! dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i) 197 217 dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i) 218 ! dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i) 198 219 ELSE 199 220 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & … … 478 499 k1, k2, theta, temp 479 500 480 REAL, DIMENSION(nzb:nzt+1) :: usvs, vsus, wsus, wsvs501 REAL, DIMENSION(nzb:nzt+1) :: usvs, vsus, wsus, wsvs 481 502 482 503 ! -
palm/trunk/SOURCE/wall_fluxes.f90
r55 r56 1 SUBROUTINE wall_fluxes( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )1 MODULE wall_fluxes_mod 2 2 !------------------------------------------------------------------------------! 3 3 ! Actual revisions: … … 15 15 ! similarity. 16 16 ! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0). 17 !------------------------------------------------------------------------------! 18 19 USE arrays_3d 20 USE control_parameters 21 USE grid_variables 22 USE indices 23 USE statistics 24 USE user 25 26 IMPLICIT NONE 27 28 INTEGER :: i, j, k, nzb_w, nzt_w, wall_index 29 REAL :: a, b, c1, c2, h1, h2, zp 30 31 REAL :: pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts 32 33 REAL, DIMENSION(nzb:nzt+1) :: wall_flux 34 35 36 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 37 wall_flux = 0.0 38 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 39 40 ! 41 !-- All subsequent variables are computed for the respective location where 42 !-- the relevant variable is defined 43 DO k = nzb_w, nzt_w 44 45 ! 46 !-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt' 47 rifs = rif_wall(k,j,i,wall_index) 48 49 u_i = a * u(k,j,i) + c1 * 0.25 * & 50 ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) 51 52 v_i = b * v(k,j,i) + c2 * 0.25 * & 53 ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) 54 55 ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( & 56 a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) & 57 + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) & 58 ) 59 pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + b * pt(k,j-1,i) & 60 + ( c1 + c2 ) * pt(k+1,j,i) ) 61 62 pts = pt_i - hom(k,1,4,0) 63 wspts = ws * pts 64 65 ! 66 !-- (2) Compute wall-parallel absolute velocity vel_total 67 vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 ) 68 69 ! 70 !-- (3) Compute wall friction velocity us_wall 71 IF ( rifs >= 0.0 ) THEN 72 73 ! 74 !-- Stable stratification (and neutral) 75 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 76 5.0 * rifs * ( zp - z0(j,i) ) / zp & 77 ) 78 ELSE 79 80 ! 81 !-- Unstable stratification 82 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 83 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 84 85 ! 86 !-- If a borderline case occurs, the formula for stable stratification 87 !-- must be used anyway, or else a zero division would occur in the 88 !-- argument of the logarithm. 89 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 17 ! The all-gridpoint version of wall_fluxes_e is not used so far, because 18 ! it gives slightly different results from the ij-version for some unknown 19 ! reason. 20 !------------------------------------------------------------------------------! 21 PRIVATE 22 PUBLIC wall_fluxes, wall_fluxes_e 23 24 INTERFACE wall_fluxes 25 MODULE PROCEDURE wall_fluxes 26 MODULE PROCEDURE wall_fluxes_ij 27 END INTERFACE wall_fluxes 28 29 INTERFACE wall_fluxes_e 30 MODULE PROCEDURE wall_fluxes_e 31 MODULE PROCEDURE wall_fluxes_e_ij 32 END INTERFACE wall_fluxes_e 33 34 CONTAINS 35 36 !------------------------------------------------------------------------------! 37 ! Call for all grid points 38 !------------------------------------------------------------------------------! 39 SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, ixp, jyp, nzb_uvw_inner, & 40 nzb_uvw_outer, wall ) 41 42 USE arrays_3d 43 USE control_parameters 44 USE grid_variables 45 USE indices 46 USE statistics 47 USE user 48 49 IMPLICIT NONE 50 51 INTEGER :: i, ixp, j, jyp, k, wall_index 52 53 INTEGER, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1) :: nzb_uvw_inner, & 54 nzb_uvw_outer 55 REAL :: a, b, c1, c2, h1, h2, zp 56 REAL :: pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts 57 58 REAL, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1) :: wall 59 REAL, DIMENSION(nzb:nzt+1,nys:nyn+jyp,nxl:nxr+ixp) :: wall_flux 60 61 62 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 63 wall_flux = 0.0 64 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 65 66 DO i = nxl, nxr+ixp 67 DO j = nys, nyn+jyp 68 69 IF ( wall(j,i) /= 0.0 ) THEN 70 ! 71 !-- All subsequent variables are computed for the respective 72 !-- location where the relevant variable is defined 73 DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i) 74 75 ! 76 !-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt' 77 rifs = rif_wall(k,j,i,wall_index) 78 79 u_i = a * u(k,j,i) + c1 * 0.25 * & 80 ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) 81 82 v_i = b * v(k,j,i) + c2 * 0.25 * & 83 ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) 84 85 ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( & 86 a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) & 87 + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) & 88 ) 89 pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + & 90 b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) ) 91 92 pts = pt_i - hom(k,1,4,0) 93 wspts = ws * pts 94 95 ! 96 !-- (2) Compute wall-parallel absolute velocity vel_total 97 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 ) 98 99 ! 100 !-- (3) Compute wall friction velocity us_wall 101 IF ( rifs >= 0.0 ) THEN 102 103 ! 104 !-- Stable stratification (and neutral) 105 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 106 5.0 * rifs * ( zp - z0(j,i) ) / zp & 107 ) 108 ELSE 109 110 ! 111 !-- Unstable stratification 112 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 113 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) )) 114 115 ! 116 !-- If a borderline case occurs, the formula for stable 117 !-- stratification must be used anyway, or else a zero 118 !-- division would occur in the argument of the logarithm. 119 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 120 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 121 5.0 * rifs * ( zp - z0(j,i) ) / zp & 122 ) 123 ELSE 124 us_wall = kappa * vel_total / ( & 125 LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + & 126 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 127 ) 128 ENDIF 129 130 ENDIF 131 132 ! 133 !-- (4) Compute zp/L (corresponds to neutral Richardson flux 134 !-- number rifs) 135 rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * & 136 ( us_wall**3 + 1E-30 ) ) 137 138 ! 139 !-- Limit the value range of the Richardson numbers. 140 !-- This is necessary for very small velocities (u,w --> 0), 141 !-- because the absolute value of rif can then become very 142 !-- large, which in consequence would result in very large 143 !-- shear stresses and very small momentum fluxes (both are 144 !-- generally unrealistic). 145 IF ( rifs < rif_min ) rifs = rif_min 146 IF ( rifs > rif_max ) rifs = rif_max 147 148 ! 149 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 150 IF ( rifs >= 0.0 ) THEN 151 152 ! 153 !-- Stable stratification (and neutral) 154 wall_flux(k,j,i) = kappa * & 155 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 156 ( LOG( zp / z0(j,i) ) + & 157 5.0 * rifs * ( zp - z0(j,i) ) / zp & 158 ) 159 ELSE 160 161 ! 162 !-- Unstable stratification 163 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 164 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) )) 165 166 ! 167 !-- If a borderline case occurs, the formula for stable 168 !-- stratification must be used anyway, or else a zero 169 !-- division would occur in the argument of the logarithm. 170 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 171 wall_flux(k,j,i) = kappa * & 172 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 173 ( LOG( zp / z0(j,i) ) + & 174 5.0 * rifs * ( zp - z0(j,i) ) / zp & 175 ) 176 ELSE 177 wall_flux(k,j,i) = kappa * & 178 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 179 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) )& 180 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 181 ) 182 ENDIF 183 184 ENDIF 185 wall_flux(k,j,i) = -wall_flux(k,j,i) * ABS(wall_flux(k,j,i)) 186 187 ! 188 !-- store rifs for next time step 189 rif_wall(k,j,i,wall_index) = rifs 190 191 ENDDO 192 193 ENDIF 194 195 ENDDO 196 ENDDO 197 198 END SUBROUTINE wall_fluxes 199 200 201 202 !------------------------------------------------------------------------------! 203 ! Call for all grid point i,j 204 !------------------------------------------------------------------------------! 205 SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 ) 206 207 USE arrays_3d 208 USE control_parameters 209 USE grid_variables 210 USE indices 211 USE statistics 212 USE user 213 214 IMPLICIT NONE 215 216 INTEGER :: i, j, k, nzb_w, nzt_w, wall_index 217 REAL :: a, b, c1, c2, h1, h2, zp 218 219 REAL :: pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts 220 221 REAL, DIMENSION(nzb:nzt+1) :: wall_flux 222 223 224 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 225 wall_flux = 0.0 226 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 227 228 ! 229 !-- All subsequent variables are computed for the respective location where 230 !-- the relevant variable is defined 231 DO k = nzb_w, nzt_w 232 233 ! 234 !-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt' 235 rifs = rif_wall(k,j,i,wall_index) 236 237 u_i = a * u(k,j,i) + c1 * 0.25 * & 238 ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) 239 240 v_i = b * v(k,j,i) + c2 * 0.25 * & 241 ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) 242 243 ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( & 244 a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) & 245 + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) & 246 ) 247 pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + b * pt(k,j-1,i) & 248 + ( c1 + c2 ) * pt(k+1,j,i) ) 249 250 pts = pt_i - hom(k,1,4,0) 251 wspts = ws * pts 252 253 ! 254 !-- (2) Compute wall-parallel absolute velocity vel_total 255 vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 ) 256 257 ! 258 !-- (3) Compute wall friction velocity us_wall 259 IF ( rifs >= 0.0 ) THEN 260 261 ! 262 !-- Stable stratification (and neutral) 90 263 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 91 264 5.0 * rifs * ( zp - z0(j,i) ) / zp & 92 265 ) 93 266 ELSE 94 us_wall = kappa * vel_total / ( & 267 268 ! 269 !-- Unstable stratification 270 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 271 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 272 273 ! 274 !-- If a borderline case occurs, the formula for stable stratification 275 !-- must be used anyway, or else a zero division would occur in the 276 !-- argument of the logarithm. 277 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 278 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 279 5.0 * rifs * ( zp - z0(j,i) ) / zp & 280 ) 281 ELSE 282 us_wall = kappa * vel_total / ( & 95 283 LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + & 96 284 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 97 ) 285 ) 286 ENDIF 287 98 288 ENDIF 99 289 100 ENDIF 101 102 ! 103 !-- (4) Compute zp/L (corresponds to neutral Richardson flux number 104 !-- rifs) 105 rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * ( us_wall**3 + 1E-30 ) ) 106 107 ! 108 !-- Limit the value range of the Richardson numbers. 109 !-- This is necessary for very small velocities (u,w --> 0), because 110 !-- the absolute value of rif can then become very large, which in 111 !-- consequence would result in very large shear stresses and very 112 !-- small momentum fluxes (both are generally unrealistic). 113 IF ( rifs < rif_min ) rifs = rif_min 114 IF ( rifs > rif_max ) rifs = rif_max 115 116 ! 117 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 118 IF ( rifs >= 0.0 ) THEN 119 120 ! 121 !-- Stable stratification (and neutral) 122 wall_flux(k) = kappa * & 123 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 124 ( LOG( zp / z0(j,i) ) + & 125 5.0 * rifs * ( zp - z0(j,i) ) / zp & 126 ) 127 ELSE 128 129 ! 130 !-- Unstable stratification 131 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 132 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 133 134 ! 135 !-- If a borderline case occurs, the formula for stable stratification 136 !-- must be used anyway, or else a zero division would occur in the 137 !-- argument of the logarithm. 138 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 290 ! 291 !-- (4) Compute zp/L (corresponds to neutral Richardson flux number 292 !-- rifs) 293 rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) ) 294 295 ! 296 !-- Limit the value range of the Richardson numbers. 297 !-- This is necessary for very small velocities (u,w --> 0), because 298 !-- the absolute value of rif can then become very large, which in 299 !-- consequence would result in very large shear stresses and very 300 !-- small momentum fluxes (both are generally unrealistic). 301 IF ( rifs < rif_min ) rifs = rif_min 302 IF ( rifs > rif_max ) rifs = rif_max 303 304 ! 305 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 306 IF ( rifs >= 0.0 ) THEN 307 308 ! 309 !-- Stable stratification (and neutral) 139 310 wall_flux(k) = kappa * & 140 311 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 141 ( LOG( zp / z0(j,i) ) +&142 5.0 * rifs * ( zp - z0(j,i) ) / zp&312 ( LOG( zp / z0(j,i) ) + & 313 5.0 * rifs * ( zp - z0(j,i) ) / zp & 143 314 ) 144 315 ELSE 145 wall_flux(k) = kappa * & 146 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 147 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) & 148 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 149 ) 316 317 ! 318 !-- Unstable stratification 319 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 320 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 321 322 ! 323 !-- If a borderline case occurs, the formula for stable stratification 324 !-- must be used anyway, or else a zero division would occur in the 325 !-- argument of the logarithm. 326 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 327 wall_flux(k) = kappa * & 328 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 329 ( LOG( zp / z0(j,i) ) + & 330 5.0 * rifs * ( zp - z0(j,i) ) / zp & 331 ) 332 ELSE 333 wall_flux(k) = kappa * & 334 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 335 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) )& 336 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 337 ) 338 ENDIF 339 150 340 ENDIF 151 152 ENDIF 153 wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) ) 154 155 ! 156 !-- store rifs for next time step 157 rif_wall(k,j,i,wall_index) = rifs 158 159 ENDDO 160 161 END SUBROUTINE wall_fluxes 162 163 164 165 SUBROUTINE wall_fluxes_e( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 ) 166 !------------------------------------------------------------------------------! 167 ! Actual revisions: 168 ! ----------------- 169 ! 170 ! 171 ! Former revisions: 172 ! ----------------- 173 ! Initial version (2007/03/07) 174 ! 341 wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) ) 342 343 ! 344 !-- store rifs for next time step 345 rif_wall(k,j,i,wall_index) = rifs 346 347 ENDDO 348 349 END SUBROUTINE wall_fluxes_ij 350 351 352 353 !------------------------------------------------------------------------------! 354 ! Call for all grid points 355 !------------------------------------------------------------------------------! 356 SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall ) 357 358 !------------------------------------------------------------------------------! 175 359 ! Description: 176 360 ! ------------ … … 180 364 !------------------------------------------------------------------------------! 181 365 182 USE arrays_3d 183 USE control_parameters 184 USE grid_variables 185 USE indices 186 USE statistics 187 USE user 188 189 IMPLICIT NONE 190 191 INTEGER :: i, j, k, kk, nzb_w, nzt_w, wall_index 192 REAL :: a, b, c1, c2, h1, h2, vel_zp, zp 193 194 REAL :: rifs 195 196 REAL, DIMENSION(nzb:nzt+1) :: wall_flux 197 198 199 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 200 wall_flux = 0.0 201 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 202 203 ! 204 !-- All subsequent variables are computed for the respective location where 205 !-- the relevant variable is defined 206 DO k = nzb_w, nzt_w 207 208 ! 209 !-- (1) Compute rifs 210 IF ( k == nzb_w ) THEN 211 kk = nzb_w 212 ELSE 213 kk = k-1 214 ENDIF 215 rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + & 216 a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + & 217 c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) & 218 ) 219 220 ! 221 !-- Skip (2) to (4) of wall_fluxes, because here rifs is already available 222 !-- from (1) 223 224 ! 225 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 226 vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + & 227 b * ( v(k,j,i) + v(k,j+1,i) ) + & 228 (c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) & 229 ) 230 231 IF ( rifs >= 0.0 ) THEN 232 233 ! 234 !-- Stable stratification (and neutral) 235 wall_flux(k) = kappa * vel_zp / & 236 ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp ) 237 ELSE 238 239 ! 240 !-- Unstable stratification 241 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 242 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 243 244 ! 245 !-- If a borderline case occurs, the formula for stable stratification 246 !-- must be used anyway, or else a zero division would occur in the 247 !-- argument of the logarithm. 248 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 249 wall_flux(k) = kappa * vel_zp / & 250 ( LOG( zp / z0(j,i) ) + & 251 5.0 * rifs * ( zp - z0(j,i) ) / zp & 252 ) 253 ELSE 254 wall_flux(k) = kappa * vel_zp / & 366 USE arrays_3d 367 USE control_parameters 368 USE grid_variables 369 USE indices 370 USE statistics 371 USE user 372 373 IMPLICIT NONE 374 375 INTEGER :: i, j, k, kk, wall_index 376 REAL :: a, b, c1, c2, h1, h2, vel_zp, zp 377 378 REAL :: rifs 379 380 REAL, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1) :: wall 381 REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wall_flux 382 383 384 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 385 wall_flux = 0.0 386 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 387 388 DO i = nxl, nxr 389 DO j = nys, nyn 390 391 IF ( wall(j,i) /= 0.0 ) THEN 392 ! 393 !-- All subsequent variables are computed for the respective 394 !-- location where the relevant variable is defined 395 DO k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2 396 397 ! 398 !-- (1) Compute rifs 399 IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN 400 kk = nzb_diff_s_inner(j,i)-1 401 ELSE 402 kk = k-1 403 ENDIF 404 rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + & 405 a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + & 406 c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) & 407 ) 408 409 ! 410 !-- Skip (2) to (4) of wall_fluxes, because here rifs is 411 !-- already available from (1) 412 413 ! 414 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 415 vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + & 416 b * ( v(k,j,i) + v(k,j+1,i) ) + & 417 (c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) & 418 ) 419 420 IF ( rifs >= 0.0 ) THEN 421 422 ! 423 !-- Stable stratification (and neutral) 424 wall_flux(k,j,i) = kappa * vel_zp / & 425 ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp ) 426 ELSE 427 428 ! 429 !-- Unstable stratification 430 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 431 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) )) 432 433 ! 434 !-- If a borderline case occurs, the formula for stable 435 !-- stratification must be used anyway, or else a zero 436 !-- division would occur in the argument of the logarithm. 437 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 438 wall_flux(k,j,i) = kappa * vel_zp / & 439 ( LOG( zp / z0(j,i) ) + & 440 5.0 * rifs * ( zp - z0(j,i) ) / zp & 441 ) 442 ELSE 443 wall_flux(k,j,i) = kappa * vel_zp / & 255 444 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) & 256 445 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 257 446 ) 447 ENDIF 448 449 ENDIF 450 wall_flux(k,j,i) = wall_flux(k,j,i) * ABS( wall_flux(k,j,i) ) 451 452 ! 453 !-- Store rifs for next time step 454 rif_wall(k,j,i,wall_index) = rifs 455 456 ENDDO 457 458 ENDIF 459 460 ENDDO 461 ENDDO 462 463 END SUBROUTINE wall_fluxes_e 464 465 466 467 !------------------------------------------------------------------------------! 468 ! Call for grid point i,j 469 !------------------------------------------------------------------------------! 470 SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 ) 471 472 USE arrays_3d 473 USE control_parameters 474 USE grid_variables 475 USE indices 476 USE statistics 477 USE user 478 479 IMPLICIT NONE 480 481 INTEGER :: i, j, k, kk, nzb_w, nzt_w, wall_index 482 REAL :: a, b, c1, c2, h1, h2, vel_zp, zp 483 484 REAL :: rifs 485 486 REAL, DIMENSION(nzb:nzt+1) :: wall_flux 487 488 489 zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) 490 wall_flux = 0.0 491 wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) 492 493 ! 494 !-- All subsequent variables are computed for the respective location where 495 !-- the relevant variable is defined 496 DO k = nzb_w, nzt_w 497 498 ! 499 !-- (1) Compute rifs 500 IF ( k == nzb_w ) THEN 501 kk = nzb_w 502 ELSE 503 kk = k-1 258 504 ENDIF 259 260 ENDIF 261 wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) ) 262 263 ! 264 !-- store rifs for next time step 265 rif_wall(k,j,i,wall_index) = rifs 266 267 ENDDO 268 269 END SUBROUTINE wall_fluxes_e 505 rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + & 506 a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + & 507 c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) & 508 ) 509 510 ! 511 !-- Skip (2) to (4) of wall_fluxes, because here rifs is already available 512 !-- from (1) 513 514 ! 515 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 516 vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + & 517 b * ( v(k,j,i) + v(k,j+1,i) ) + & 518 (c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) & 519 ) 520 521 IF ( rifs >= 0.0 ) THEN 522 523 ! 524 !-- Stable stratification (and neutral) 525 wall_flux(k) = kappa * vel_zp / & 526 ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp ) 527 ELSE 528 529 ! 530 !-- Unstable stratification 531 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 532 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 533 534 ! 535 !-- If a borderline case occurs, the formula for stable stratification 536 !-- must be used anyway, or else a zero division would occur in the 537 !-- argument of the logarithm. 538 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 539 wall_flux(k) = kappa * vel_zp / & 540 ( LOG( zp / z0(j,i) ) + & 541 5.0 * rifs * ( zp - z0(j,i) ) / zp & 542 ) 543 ELSE 544 wall_flux(k) = kappa * vel_zp / & 545 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) & 546 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 547 ) 548 ENDIF 549 550 ENDIF 551 wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) ) 552 553 ! 554 !-- Store rifs for next time step 555 rif_wall(k,j,i,wall_index) = rifs 556 557 ENDDO 558 559 END SUBROUTINE wall_fluxes_e_ij 560 561 END MODULE wall_fluxes_mod
Note: See TracChangeset
for help on using the changeset viewer.
