Changeset 187 for palm/trunk/SOURCE
- Timestamp:
- Aug 6, 2008 4:25:09 PM (16 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 4 edited
-
CURRENT_MODIFICATIONS (modified) (3 diffs)
-
prandtl_fluxes.f90 (modified) (7 diffs)
-
production_e.f90 (modified) (9 diffs)
-
wall_fluxes.f90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/CURRENT_MODIFICATIONS
r184 r187 42 42 Changed: 43 43 ------- 44 Modification of the integrated version of the profile function for momentum 45 for unstable stratification; more consistent flux definitions. (wall_fluxes, 46 production_e) 47 44 48 Strict grid matching along z is not needed for mg-solver. (check_parameters) 45 49 … … 79 83 (init_1d_model) 80 84 81 advec_s_ups, advec_u_ups, advec_v_ups, advec_w_ups, calc_spectra, check_open, check_parameters, cpu_statistics, init_1d_model, init_3d_model, modules, palm, parin, poisfft, read_var_list, read_3d_binary, transpose, write_var_list, write_3d_binary85 advec_s_ups, advec_u_ups, advec_v_ups, advec_w_ups, calc_spectra, check_open, check_parameters, cpu_statistics, init_1d_model, init_3d_model, modules, palm, parin, poisfft, production_e, read_var_list, read_3d_binary, transpose, wall_fluxes, write_var_list, write_3d_binary 82 86 83 87 84 88 Errors: 85 89 ------ 90 Bugfix: change definition of us_wall from 1D to 2D. Tests showed that this decreases 91 u* by some 10% and increases TKE and momentum fluxes by some 10% because friction 92 was underestimated and momentum fluxes were wrongly calculated due to the bug. 93 (prandtl_fluxes, wall_fluxes) 94 86 95 Bugfix: calculation of horizontal fluxes at vertical walls (diffusion_s) 87 96 … … 106 115 107 116 Introduce prefix_chr to ensure unique dvrp_file path. 117 108 118 small bugfixes for user_interface sample code (comments): 109 119 - initialize ustvst with 0.0 as it is now computed only until nxr and nyn 110 120 - two ALLOCATE statements moved from user_read_restart_data back to user_init 111 121 - remove 'READ (13) u2_av' statement in user_read_restart_data 122 112 123 Bugfix: remove IF statement in plant_canopy_model_ij (plant_canopy_model) 124 113 125 Bugfix: divide sums(k,8) (e) and sums(k,34) (e*) by ngp_2dh_s_inner(k,sr) 114 126 (like other scalars) (flow_statistics) 127 115 128 Bugfix: dopr_time_count was written on the binary file, which caused that 116 129 NetCDF files newly created by restart files (no append of existing files!) 117 130 contained uneccessary time levels. (read_3d_binary, write_3d_binary) 131 118 132 Bugfix: extra '*' removed in user_statistics sample code (user_interface) 133 119 134 Bugfix: a stop command was missing in some cases of the parallel branch (local_stop) 135 120 136 Bugfix in volume flow control for non-cyclic boundary conditions (pres) 137 121 138 Bugfix: misplaced #endif directives (combine_plot_fields) 122 139 140 check_parameters, diffusion_s, flow_statistics, init_dvrp, init_3d_model, local_stop, plant_canopy_model, poismg, prandtl_fluxes, pres, read_3d_binary, user_interface, wall_fluxes, write_3d_binary 123 141 124 check_parameters, diffusion_s, flow_statistics, init_dvrp, init_3d_model, local_stop, plant_canopy_model, poismg, pres, read_3d_binary, user_interface, write_3d_binary125 -
palm/trunk/SOURCE/prandtl_fluxes.f90
r110 r187 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Bugfix: change definition of us_wall from 1D to 2D: 7 ! Modification of the evaluation of the vertical turbulent momentum 8 ! fluxes u'w' and v'w'; the first usws that is computed corresponds 9 ! to -u'w'/u* and not as priorily assumed to (-u'w')**0.5, the first 10 ! vsws that is computed corresponds to -v'w'/u* and not as priorily 11 ! assumed to (-v'w')**0.5. Therefore, the intermediate result for 12 ! usws has to be multiplied by -u* instead by itself in order to 13 ! get u'w'. Accordingly, the intermediate result for vsws has to be 14 ! multiplied by -u* instead by itself in order to get v'w'. As u* 15 ! is calculated for the position of a scalar an additional 16 ! interpolation of u* to the position of u and v, respectively, 17 ! is necessary. As u* is not determined for the ghost points on 18 ! each PE, an additional exchange of information from neighbouring 19 ! PEs is necessary.! 20 ! Change: Modification of the integrated version of the profile function for 21 ! momentum for unstable stratification 7 22 ! 8 23 ! Former revisions: … … 81 96 !-- Unstable stratification 82 97 a = SQRT( 1.0 - 16.0 * rif(j,i) ) 83 b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) 84 ! 85 !-- If a borderline case occurs, the formula for stable 86 !-- stratification must be used anyway, or else a zero division 87 !-- would occur in the argument of the logarithm 88 IF ( a == 1.0 .OR. b == 1.0 ) THEN 89 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( & 90 LOG( z_p / z0(j,i) ) + & 91 5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p & 92 ) 93 ELSE 94 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( & 95 LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) & 96 ) 97 ENDIF 98 b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p ) 99 100 ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / ( & 101 LOG( z_p / z0(j,i) ) - & 102 2.0 * LOG( ( 1.0 + a ) / ( 1.0 + b ) ) ) 98 103 ENDIF 99 104 … … 168 173 ! 169 174 !-- Unstable stratification 170 a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) ) 171 b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) ) 172 ! 173 !-- If a borderline case occurs, the formula for stable stratification 174 !-- must be used anyway, or else a zero division would occur in the 175 !-- argument of the logarithm. 176 IF ( a == 1.0 .OR. b == 1.0 ) THEN 177 us(j,i) = kappa * uv_total / ( & 178 LOG( z_p / z0(j,i) ) + & 179 5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p & 180 ) 181 ELSE 182 us(j,i) = kappa * uv_total / ( & 183 LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + & 184 2.0 * ( ATAN( b ) - ATAN( a ) ) & 185 ) 186 ENDIF 175 a = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) ) 176 b = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) ) 177 178 us(j,i) = kappa * uv_total / ( & 179 LOG( z_p / z0(j,i) ) - & 180 LOG( ( 1.0 + a )**2 * ( 1.0 + a**2 ) / ( & 181 ( 1.0 + b )**2 * ( 1.0 + b**2 ) ) ) + & 182 2.0 * ( ATAN( a ) - ATAN( b ) ) & 183 ) 187 184 ENDIF 188 185 ENDDO 189 186 ENDDO 190 187 188 ! 189 !-- Values of us at ghost point locations are needed for the evaluation of usws 190 !-- and vsws. 191 CALL exchange_horiz_2d( us ) 191 192 ! 192 193 !-- Compute u'w' for the total model domain. … … 212 213 ! 213 214 !-- Unstable stratification 214 a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) ) 215 b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) ) 216 ! 217 !-- If a borderline case occurs, the formula for stable stratification 218 !-- must be used anyway, or else a zero division would occur in the 219 !-- argument of the logarithm. 220 IF ( a == 1.0 .OR. B == 1.0 ) THEN 221 usws(j,i) = kappa * u(k+1,j,i) / ( & 222 LOG( z_p / z0(j,i) ) + & 223 5.0 * rifm * ( z_p - z0(j,i) ) / z_p & 215 a = SQRT( SQRT( 1.0 - 16.0 * rifm ) ) 216 b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) ) 217 218 usws(j,i) = kappa * u(k+1,j,i) / ( & 219 LOG( z_p / z0(j,i) ) - & 220 LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / ( & 221 (1.0 + b )**2 * ( 1.0 + b**2 ) ) ) + & 222 2.0 * ( ATAN( a ) - ATAN( b ) ) & 224 223 ) 225 ELSE226 usws(j,i) = kappa * u(k+1,j,i) / ( &227 LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &228 2.0 * ( ATAN( b ) - ATAN( a ) ) &229 )230 ENDIF231 224 ENDIF 232 usws(j,i) = -usws(j,i) * ABS( usws(j,i) )225 usws(j,i) = -usws(j,i) * 0.5 * ( us(j,i-1) + us(j,i) ) 233 226 ENDDO 234 227 ENDDO … … 257 250 ! 258 251 !-- Unstable stratification 259 a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) ) 260 b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) ) 261 ! 262 !-- If a borderline case occurs, the formula for stable stratification 263 !-- must be used anyway, or else a zero division would occur in the 264 !-- argument of the logarithm. 265 IF ( a == 1.0 .OR. b == 1.0 ) THEN 266 vsws(j,i) = kappa * v(k+1,j,i) / ( & 267 LOG( z_p / z0(j,i) ) + & 268 5.0 * rifm * ( z_p - z0(j,i) ) / z_p & 252 a = SQRT( SQRT( 1.0 - 16.0 * rifm ) ) 253 b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) ) 254 255 vsws(j,i) = kappa * v(k+1,j,i) / ( & 256 LOG( z_p / z0(j,i) ) - & 257 LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / ( & 258 (1.0 + b )**2 * ( 1.0 + b**2 ) ) ) + & 259 2.0 * ( ATAN( a ) - ATAN( b ) ) & 269 260 ) 270 ELSE271 vsws(j,i) = kappa * v(k+1,j,i) / ( &272 LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &273 2.0 * ( ATAN( b ) - ATAN( a ) ) &274 )275 ENDIF276 261 ENDIF 277 vsws(j,i) = -vsws(j,i) * ABS( vsws(j,i) )262 vsws(j,i) = -vsws(j,i) * 0.5 * ( us(j-1,i) + us(j,i) ) 278 263 ENDDO 279 264 ENDDO … … 317 302 ! 318 303 !-- Unstable stratification 319 a = SQRT( 1.0 - 16.0 * rif(j,i) ) 320 b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) 321 ! 322 !-- If a borderline case occurs, the formula for stable 323 !-- stratification must be used anyway, or else a zero division 324 !-- would occur in the argument of the logarithm. 325 IF ( a == 1.0 .OR. b == 1.0 ) THEN 326 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( & 327 LOG( z_p / z0(j,i) ) + & 328 5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p & 329 ) 330 ELSE 331 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( & 332 LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) & 333 ) 334 ENDIF 304 a = SQRT( 1.0 - 16.0 * rif(j,i) ) 305 b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p ) 306 307 qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / ( & 308 LOG( z_p / z0(j,i) ) - & 309 2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) ) 335 310 ENDIF 336 311 … … 341 316 342 317 ! 343 !-- Exchange the boundaries for u* and the momentum fluxes (fluxes only for 344 !-- completeness's sake). 345 CALL exchange_horiz_2d( us ) 318 !-- Exchange the boundaries for the momentum fluxes (only for sake of 319 !-- completeness) 346 320 CALL exchange_horiz_2d( usws ) 347 321 CALL exchange_horiz_2d( vsws ) -
palm/trunk/SOURCE/production_e.f90
r139 r187 4 4 ! Actual revisions: 5 5 ! ----------------- 6 ! 6 ! Change: add 'minus' sign to fluxes obtained from subroutine wall_fluxes_e for 7 ! consistency with subroutine wall_fluxes 7 8 ! 8 9 ! Former revisions: … … 164 165 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 165 166 usvs, 1.0, 0.0, 0.0, 0.0 ) 166 dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)167 ! dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)167 dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i) 168 ! dudy = - wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i) 168 169 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 169 170 wsvs, 0.0, 0.0, 1.0, 0.0 ) 170 dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)171 ! dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)171 dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i) 172 ! dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i) 172 173 ELSE 173 174 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & … … 180 181 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 181 182 vsus, 0.0, 1.0, 0.0, 0.0 ) 182 dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)183 ! dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)183 dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i) 184 ! dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i) 184 185 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 185 186 wsus, 0.0, 0.0, 0.0, 1.0 ) 186 dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)187 ! dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)187 dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i) 188 ! dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i) 188 189 ELSE 189 190 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & … … 220 221 221 222 IF ( wall_e_y(j,i) /= 0.0 ) THEN 222 dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)223 ! dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)224 dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)225 ! dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)223 dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i) 224 ! dudy = - wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i) 225 dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i) 226 ! dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i) 226 227 ELSE 227 228 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & … … 232 233 233 234 IF ( wall_e_x(j,i) /= 0.0 ) THEN 234 dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)235 ! dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)236 dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)237 ! dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)235 dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i) 236 ! dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i) 237 dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i) 238 ! dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i) 238 239 ELSE 239 240 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & … … 629 630 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 630 631 usvs, 1.0, 0.0, 0.0, 0.0 ) 631 dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)632 dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i) 632 633 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 633 634 wsvs, 0.0, 0.0, 1.0, 0.0 ) 634 dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)635 dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i) 635 636 ELSE 636 637 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & … … 643 644 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 644 645 vsus, 0.0, 1.0, 0.0, 0.0 ) 645 dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)646 dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i) 646 647 CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, & 647 648 wsus, 0.0, 0.0, 0.0, 1.0 ) 648 dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)649 dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i) 649 650 ELSE 650 651 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & … … 679 680 680 681 IF ( wall_e_y(j,i) /= 0.0 ) THEN 681 dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)682 dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)682 dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i) 683 dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i) 683 684 ELSE 684 685 dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - & … … 689 690 690 691 IF ( wall_e_x(j,i) /= 0.0 ) THEN 691 dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)692 dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)692 dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i) 693 dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i) 693 694 ELSE 694 695 dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - & -
palm/trunk/SOURCE/wall_fluxes.f90
r75 r187 3 3 ! Actual revisions: 4 4 ! ----------------- 5 ! 5 ! Bugfix: change definition of us_wall from 1D to 2D: 6 ! Modification of the evaluation of the vertical turbulent momentum 7 ! fluxes u'w' and v'w'; the first usws that is computed corresponds 8 ! to -u'w'/u* and not as priorily assumed to (-u'w')**0.5, the first 9 ! vsws that is computed corresponds to -v'w'/u* and not as priorily 10 ! assumed to (-v'w')**0.5. Therefore, the intermediate result for 11 ! usws has to be multiplied by -u* instead by itself in order to 12 ! get u'w'. Accordingly, the intermediate result for vsws has to be 13 ! multiplied by -u* instead by itself in order to get v'w'. 14 ! This requires the calculation of us_wall (and vel_total, u_i, v_i, ws) 15 ! also in wall_fluxes_e. 16 ! Bugfix: storage of rifs to rifs_wall in wall_fluxes_e removed 17 ! Change: add 'minus' sign to fluxes produced by subroutine wall_fluxes_e for 18 ! consistency with subroutine wall_fluxes 19 ! Change: Modification of the integrated version of the profile function for 20 ! momentum for unstable stratification 6 21 ! 7 22 ! Former revisions: … … 69 84 ! 70 85 !-- All subsequent variables are computed for the respective 71 !-- location where the re levant variable is defined86 !-- location where the respective flux is defined. 72 87 DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i) 73 88 … … 109 124 ! 110 125 !-- Unstable stratification 111 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 112 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) )) 113 114 ! 115 !-- If a borderline case occurs, the formula for stable 116 !-- stratification must be used anyway, or else a zero 117 !-- division would occur in the argument of the logarithm. 118 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 119 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 120 5.0 * rifs * ( zp - z0(j,i) ) / zp & 121 ) 122 ELSE 123 us_wall = kappa * vel_total / ( & 124 LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + & 125 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 126 ) 127 ENDIF 128 126 h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 127 h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) ) 128 129 us_wall = kappa * vel_total / ( & 130 LOG( zp / z0(j,i) ) - & 131 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 132 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 133 2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) & 134 ) 129 135 ENDIF 130 136 … … 160 166 ! 161 167 !-- Unstable stratification 162 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 163 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) )) 164 165 ! 166 !-- If a borderline case occurs, the formula for stable 167 !-- stratification must be used anyway, or else a zero 168 !-- division would occur in the argument of the logarithm. 169 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 170 wall_flux(k,j,i) = kappa * & 171 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 172 ( LOG( zp / z0(j,i) ) + & 173 5.0 * rifs * ( zp - z0(j,i) ) / zp & 174 ) 175 ELSE 176 wall_flux(k,j,i) = kappa * & 177 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 178 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) )& 179 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 180 ) 181 ENDIF 182 168 h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 169 h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) ) 170 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 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 175 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 176 2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) & 177 ) 183 178 ENDIF 184 wall_flux(k,j,i) = -wall_flux(k,j,i) * ABS(wall_flux(k,j,i))179 wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall 185 180 186 181 ! … … 226 221 ! 227 222 !-- All subsequent variables are computed for the respective location where 228 !-- the re levant variable is defined223 !-- the respective flux is defined. 229 224 DO k = nzb_w, nzt_w 230 225 … … 266 261 ! 267 262 !-- Unstable stratification 268 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 269 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 270 271 ! 272 !-- If a borderline case occurs, the formula for stable stratification 273 !-- must be used anyway, or else a zero division would occur in the 274 !-- argument of the logarithm. 275 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 276 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 277 5.0 * rifs * ( zp - z0(j,i) ) / zp & 278 ) 279 ELSE 280 us_wall = kappa * vel_total / ( & 281 LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + & 282 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 283 ) 284 ENDIF 285 263 h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 264 h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) ) 265 266 us_wall = kappa * vel_total / ( & 267 LOG( zp / z0(j,i) ) - & 268 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 269 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 270 2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) & 271 ) 286 272 ENDIF 287 273 … … 315 301 ! 316 302 !-- Unstable stratification 317 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 318 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 319 320 ! 321 !-- If a borderline case occurs, the formula for stable stratification 322 !-- must be used anyway, or else a zero division would occur in the 323 !-- argument of the logarithm. 324 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 325 wall_flux(k) = kappa * & 326 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 327 ( LOG( zp / z0(j,i) ) + & 328 5.0 * rifs * ( zp - z0(j,i) ) / zp & 329 ) 330 ELSE 331 wall_flux(k) = kappa * & 332 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & 333 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) )& 334 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 335 ) 336 ENDIF 337 303 h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 304 h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) ) 305 306 wall_flux(k) = kappa * & 307 ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / ( & 308 LOG( zp / z0(j,i) ) - & 309 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 310 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 311 2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) & 312 ) 338 313 ENDIF 339 wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )314 wall_flux(k) = -wall_flux(k) * us_wall 340 315 341 316 ! … … 371 346 372 347 INTEGER :: i, j, k, kk, wall_index 373 REAL :: a, b, c1, c2, h1, h2, vel_zp, zp 348 REAL :: a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, & 349 ws, zp 374 350 375 351 REAL :: rifs … … 388 364 IF ( wall(j,i) /= 0.0 ) THEN 389 365 ! 390 !-- All subsequent variables are computed for the respective 391 !-- location where the relevant variable is defined 366 !-- All subsequent variables are computed for scalar locations. 392 367 DO k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2 393 394 ! 395 !-- (1) Compute rifs 368 ! 369 !-- (1) Compute rifs, u_i, v_i, and ws 396 370 IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN 397 371 kk = nzb_diff_s_inner(j,i)-1 … … 404 378 ) 405 379 406 ! 407 !-- Skip (2) to (4) of wall_fluxes, because here rifs is 408 !-- already available from (1) 409 380 u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) ) 381 v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) ) 382 ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) ) 383 ! 384 !-- (2) Compute wall-parallel absolute velocity vel_total and 385 !-- interpolate appropriate velocity component vel_zp. 386 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 ) 387 vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws ) 388 ! 389 !-- (3) Compute wall friction velocity us_wall 390 IF ( rifs >= 0.0 ) THEN 391 392 ! 393 !-- Stable stratification (and neutral) 394 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 395 5.0 * rifs * ( zp - z0(j,i) ) / zp & 396 ) 397 ELSE 398 399 ! 400 !-- Unstable stratification 401 h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 402 h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) ) 403 404 us_wall = kappa * vel_total / ( & 405 LOG( zp / z0(j,i) ) - & 406 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 407 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 408 2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) & 409 ) 410 ENDIF 411 412 ! 413 !-- Skip step (4) of wall_fluxes, because here rifs is already 414 !-- available from (1) 410 415 ! 411 416 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 412 vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + &413 b * ( v(k,j,i) + v(k,j+1,i) ) + &414 (c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) &415 )416 417 417 418 IF ( rifs >= 0.0 ) THEN … … 425 426 ! 426 427 !-- Unstable stratification 427 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 428 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) )) 429 430 ! 431 !-- If a borderline case occurs, the formula for stable 432 !-- stratification must be used anyway, or else a zero 433 !-- division would occur in the argument of the logarithm. 434 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 435 wall_flux(k,j,i) = kappa * vel_zp / & 436 ( LOG( zp / z0(j,i) ) + & 437 5.0 * rifs * ( zp - z0(j,i) ) / zp & 438 ) 439 ELSE 440 wall_flux(k,j,i) = kappa * vel_zp / & 441 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) & 442 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 443 ) 444 ENDIF 445 428 h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 429 h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) ) 430 431 wall_flux(k,j,i) = kappa * vel_zp / ( & 432 LOG( zp / z0(j,i) ) - & 433 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 434 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 435 2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) & 436 ) 446 437 ENDIF 447 wall_flux(k,j,i) = wall_flux(k,j,i) * ABS( wall_flux(k,j,i) ) 448 449 ! 450 !-- Store rifs for next time step 451 rif_wall(k,j,i,wall_index) = rifs 438 wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall 452 439 453 440 ENDDO … … 476 463 477 464 INTEGER :: i, j, k, kk, nzb_w, nzt_w, wall_index 478 REAL :: a, b, c1, c2, h1, h2, vel_zp, zp 465 REAL :: a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, & 466 ws, zp 479 467 480 468 REAL :: rifs … … 488 476 489 477 ! 490 !-- All subsequent variables are computed for the respective location where 491 !-- the relevant variable is defined 478 !-- All subsequent variables are computed for scalar locations. 492 479 DO k = nzb_w, nzt_w 493 480 494 481 ! 495 !-- (1) Compute rifs 482 !-- (1) Compute rifs, u_i, v_i, and ws 496 483 IF ( k == nzb_w ) THEN 497 484 kk = nzb_w … … 504 491 ) 505 492 506 ! 507 !-- Skip (2) to (4) of wall_fluxes, because here rifs is already available 508 !-- from (1) 509 493 u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) ) 494 v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) ) 495 ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) ) 496 ! 497 !-- (2) Compute wall-parallel absolute velocity vel_total and 498 !-- interpolate appropriate velocity component vel_zp. 499 vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 ) 500 vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws ) 501 ! 502 !-- (3) Compute wall friction velocity us_wall 503 IF ( rifs >= 0.0 ) THEN 504 505 ! 506 !-- Stable stratification (and neutral) 507 us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & 508 5.0 * rifs * ( zp - z0(j,i) ) / zp & 509 ) 510 ELSE 511 512 ! 513 !-- Unstable stratification 514 h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 515 h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) ) 516 517 us_wall = kappa * vel_total / ( & 518 LOG( zp / z0(j,i) ) - & 519 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 520 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 521 2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) & 522 ) 523 ENDIF 524 525 ! 526 !-- Skip step (4) of wall_fluxes, because here rifs is already 527 !-- available from (1) 510 528 ! 511 529 !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') 530 !-- First interpolate the velocity (this is different from 531 !-- subroutine wall_fluxes because fluxes in subroutine 532 !-- wall_fluxes_e are defined at scalar locations). 512 533 vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + & 513 534 b * ( v(k,j,i) + v(k,j+1,i) ) + & … … 525 546 ! 526 547 !-- Unstable stratification 527 h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 528 h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) 529 530 ! 531 !-- If a borderline case occurs, the formula for stable stratification 532 !-- must be used anyway, or else a zero division would occur in the 533 !-- argument of the logarithm. 534 IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN 535 wall_flux(k) = kappa * vel_zp / & 536 ( LOG( zp / z0(j,i) ) + & 537 5.0 * rifs * ( zp - z0(j,i) ) / zp & 538 ) 539 ELSE 540 wall_flux(k) = kappa * vel_zp / & 541 ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) & 542 + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & 543 ) 544 ENDIF 545 548 h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) ) 549 h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) ) 550 551 wall_flux(k) = kappa * vel_zp / ( & 552 LOG( zp / z0(j,i) ) - & 553 LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / ( & 554 ( 1.0 + h2 )**2 * ( 1.0 + h2**2 ) ) ) + & 555 2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) & 556 ) 546 557 ENDIF 547 wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) ) 548 549 ! 550 !-- Store rifs for next time step 551 rif_wall(k,j,i,wall_index) = rifs 558 wall_flux(k) = - wall_flux(k) * us_wall 552 559 553 560 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
