Changeset 2232 for palm/trunk/SOURCE/diffusion_v.f90
- Timestamp:
- May 30, 2017 5:47:52 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/diffusion_v.f90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/diffusion_v.f90
r2119 r2232 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Adjustments to new topography and surface concept 23 23 ! 24 24 ! Former revisions: … … 90 90 91 91 92 USE wall_fluxes_mod93 94 92 PRIVATE 95 93 PUBLIC diffusion_v … … 111 109 112 110 USE arrays_3d, & 113 ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w, & 114 drho_air, rho_air_zw 111 ONLY: ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw 115 112 116 113 USE control_parameters, & 117 ONLY: constant_top_momentumflux, topography, use_surface_fluxes,&114 ONLY: constant_top_momentumflux, use_surface_fluxes, & 118 115 use_top_fluxes 119 116 120 117 USE grid_variables, & 121 ONLY: ddx, ddy, ddy2 , fxm, fxp, wall_v118 ONLY: ddx, ddy, ddy2 122 119 123 120 USE indices, & 124 ONLY: nxl, nxr, nyn, nys, nysv, nzb, nzb_diff_v, nzb_v_inner, & 125 nzb_v_outer, nzt, nzt_diff 121 ONLY: nxl, nxr, nyn, nys, nysv, nzb, nzt, wall_flags_0 126 122 127 123 USE kinds 128 124 125 USE surface_mod, & 126 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 127 surf_usm_v 128 129 129 IMPLICIT NONE 130 130 131 INTEGER(iwp) :: i !<132 INTEGER(iwp) :: j !<133 INTEGER(iwp) :: k !<134 REAL(wp) :: kmxm !<135 REAL(wp) :: kmxp !<136 REAL(wp) :: kmzm !<137 REAL(wp) :: kmzp !<138 139 REAL(wp) , DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: vsus !<140 141 ! 142 !-- First calculate horizontal momentum flux v'u' at vertical walls, 143 !-- if neccessary 144 IF ( topography /= 'flat' ) THEN145 CALL wall_fluxes( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, nzb_v_inner, &146 nzb_v_outer, wall_v )147 ENDIF131 INTEGER(iwp) :: i !< running index x direction 132 INTEGER(iwp) :: j !< running index y direction 133 INTEGER(iwp) :: k !< running index z direction 134 INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall 135 INTEGER(iwp) :: m !< running index surface elements 136 INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint 137 INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint 138 139 REAL(wp) :: flag !< flag to mask topography grid points 140 REAL(wp) :: kmxm !< 141 REAL(wp) :: kmxp !< 142 REAL(wp) :: kmzm !< 143 REAL(wp) :: kmzp !< 144 REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface 145 REAL(wp) :: mask_east !< flag to mask vertical surface south of the grid point 146 REAL(wp) :: mask_west !< flag to mask vertical surface north of the grid point 147 REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface 148 148 149 149 DO i = nxl, nxr … … 151 151 ! 152 152 !-- Compute horizontal diffusion 153 DO k = nzb_v_outer(j,i)+1, nzt 153 DO k = nzb+1, nzt 154 155 ! 156 !-- Predetermine flag to mask topography and wall-bounded grid points. 157 !-- It is sufficient to masked only east- and west-facing surfaces, which 158 !-- need special treatment for the v-component. 159 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 2 ) ) 160 mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 2 ) ) 161 mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 2 ) ) 154 162 ! 155 163 !-- Interpolate eddy diffusivities on staggered gridpoints 156 kmxp = 0.25_wp * & 157 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 158 kmxm = 0.25_wp * & 159 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 160 161 tend(k,j,i) = tend(k,j,i) & 162 & + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & 163 & + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 164 & - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 165 & - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 166 & ) * ddx & 167 & + 2.0_wp * ( & 168 & km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 169 & - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 170 & ) * ddy2 164 kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 165 kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 166 167 tend(k,j,i) = tend(k,j,i) + ( & 168 mask_east * kmxp * ( & 169 ( v(k,j,i+1) - v(k,j,i) ) * ddx & 170 + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 171 ) & 172 - mask_west * kmxm * ( & 173 ( v(k,j,i) - v(k,j,i-1) ) * ddx & 174 + ( u(k,j,i) - u(k,j-1,i) ) * ddy & 175 ) & 176 ) * ddx * flag & 177 + 2.0_wp * ( & 178 km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 179 - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 180 ) * ddy2 * flag 181 171 182 ENDDO 172 183 173 184 ! 174 !-- Wall functions at the left and right walls, respectively 175 IF ( wall_v(j,i) /= 0.0_wp ) THEN 176 177 DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i) 178 kmxp = 0.25_wp * & 179 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 180 kmxm = 0.25_wp * & 181 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 182 183 tend(k,j,i) = tend(k,j,i) & 184 + 2.0_wp * ( & 185 km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 186 - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 187 ) * ddy2 & 188 + ( fxp(j,i) * ( & 189 kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & 190 + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 191 ) & 192 - fxm(j,i) * ( & 193 kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 194 + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 195 ) & 196 + wall_v(j,i) * vsus(k,j,i) & 197 ) * ddx 198 ENDDO 199 ENDIF 200 201 ! 202 !-- Compute vertical diffusion. In case of simulating a Prandtl 203 !-- layer, index k starts at nzb_v_inner+2. 204 DO k = nzb_diff_v(j,i), nzt_diff 185 !-- Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3) 186 !-- surfaces. Note, in the the flat case, loops won't be entered as 187 !-- start_index > end_index. Furtermore, note, no vertical natural surfaces 188 !-- so far. 189 !-- Default-type surfaces 190 DO l = 2, 3 191 surf_s = surf_def_v(l)%start_index(j,i) 192 surf_e = surf_def_v(l)%end_index(j,i) 193 DO m = surf_s, surf_e 194 k = surf_def_v(l)%k(m) 195 tend(k,j,i) = tend(k,j,i) + & 196 surf_def_v(l)%mom_flux_uv(m) * ddx 197 ENDDO 198 ENDDO 199 ! 200 !-- Natural-type surfaces 201 DO l = 2, 3 202 surf_s = surf_lsm_v(l)%start_index(j,i) 203 surf_e = surf_lsm_v(l)%end_index(j,i) 204 DO m = surf_s, surf_e 205 k = surf_lsm_v(l)%k(m) 206 tend(k,j,i) = tend(k,j,i) + & 207 surf_lsm_v(l)%mom_flux_uv(m) * ddx 208 ENDDO 209 ENDDO 210 ! 211 !-- Urban-type surfaces 212 DO l = 2, 3 213 surf_s = surf_usm_v(l)%start_index(j,i) 214 surf_e = surf_usm_v(l)%end_index(j,i) 215 DO m = surf_s, surf_e 216 k = surf_usm_v(l)%k(m) 217 tend(k,j,i) = tend(k,j,i) + & 218 surf_usm_v(l)%mom_flux_uv(m) * ddx 219 ENDDO 220 ENDDO 221 ! 222 !-- Compute vertical diffusion. In case of simulating a surface layer, 223 !-- respective grid diffusive fluxes are masked (flag 10) within this 224 !-- loop, and added further below, else, simple gradient approach is 225 !-- applied. Model top is also mask if top-momentum flux is given. 226 DO k = nzb+1, nzt 227 ! 228 !-- Determine flags to mask topography below and above. Flag 2 is 229 !-- used to mask topography in general, while flag 8 implies also 230 !-- information about use_surface_fluxes. Flag 9 is used to control 231 !-- momentum flux at model top. 232 mask_bottom = MERGE( 1.0_wp, 0.0_wp, & 233 BTEST( wall_flags_0(k-1,j,i), 8 ) ) 234 mask_top = MERGE( 1.0_wp, 0.0_wp, & 235 BTEST( wall_flags_0(k+1,j,i), 8 ) ) * & 236 MERGE( 1.0_wp, 0.0_wp, & 237 BTEST( wall_flags_0(k+1,j,i), 9 ) ) 238 flag = MERGE( 1.0_wp, 0.0_wp, & 239 BTEST( wall_flags_0(k,j,i), 2 ) ) 205 240 ! 206 241 !-- Interpolate eddy diffusivities on staggered gridpoints … … 213 248 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 214 249 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 215 & ) * rho_air_zw(k) 250 & ) * rho_air_zw(k) * mask_top & 216 251 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 217 252 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 218 & ) * rho_air_zw(k-1) 219 & ) * ddzw(k) * drho_air(k) 253 & ) * rho_air_zw(k-1) * mask_bottom & 254 & ) * ddzw(k) * drho_air(k) * flag 220 255 ENDDO 221 256 … … 228 263 !-- comparison with other (LES) models showed that the values of 229 264 !-- the momentum flux becomes too large in this case. 230 !-- The term containing w(k-1,..) (see above equation) is removed here231 !-- because the vertical velocity is assumed to be zero at the surface.232 265 IF ( use_surface_fluxes ) THEN 233 k = nzb_v_inner(j,i)+1 234 ! 235 !-- Interpolate eddy diffusivities on staggered gridpoints 236 kmzp = 0.25_wp * & 237 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 238 239 tend(k,j,i) = tend(k,j,i) & 240 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 241 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 242 & ) * rho_air_zw(k) & 243 & - ( -vsws(j,i) ) & 244 & ) * ddzw(k) * drho_air(k) 266 ! 267 !-- Default-type surfaces, upward-facing 268 surf_s = surf_def_h(0)%start_index(j,i) 269 surf_e = surf_def_h(0)%end_index(j,i) 270 DO m = surf_s, surf_e 271 k = surf_def_h(0)%k(m) 272 273 tend(k,j,i) = tend(k,j,i) & 274 + ( - ( - surf_def_h(0)%vsws(m) ) & 275 ) * ddzw(k) * drho_air(k) 276 ENDDO 277 ! 278 !-- Default-type surfaces, dowward-facing 279 surf_s = surf_def_h(1)%start_index(j,i) 280 surf_e = surf_def_h(1)%end_index(j,i) 281 DO m = surf_s, surf_e 282 k = surf_def_h(1)%k(m) 283 284 tend(k,j,i) = tend(k,j,i) & 285 + ( - surf_def_h(1)%vsws(m) & 286 ) * ddzw(k) * drho_air(k) 287 ENDDO 288 ! 289 !-- Natural-type surfaces, upward-facing 290 surf_s = surf_lsm_h%start_index(j,i) 291 surf_e = surf_lsm_h%end_index(j,i) 292 DO m = surf_s, surf_e 293 k = surf_lsm_h%k(m) 294 295 tend(k,j,i) = tend(k,j,i) & 296 + ( - ( - surf_lsm_h%vsws(m) ) & 297 ) * ddzw(k) * drho_air(k) 298 299 ENDDO 300 ! 301 !-- Urban-type surfaces, upward-facing 302 surf_s = surf_usm_h%start_index(j,i) 303 surf_e = surf_usm_h%end_index(j,i) 304 DO m = surf_s, surf_e 305 k = surf_usm_h%k(m) 306 307 tend(k,j,i) = tend(k,j,i) & 308 + ( - ( - surf_usm_h%vsws(m) ) & 309 ) * ddzw(k) * drho_air(k) 310 311 ENDDO 245 312 ENDIF 246 247 ! 248 !-- Vertical diffusion at the first gridpoint below the top boundary, 249 !-- if the momentum flux at the top is prescribed by the user 250 IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN 251 k = nzt 252 ! 253 !-- Interpolate eddy diffusivities on staggered gridpoints 254 kmzm = 0.25_wp * & 255 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 256 257 tend(k,j,i) = tend(k,j,i) & 258 & + ( ( -vswst(j,i) ) & 259 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 260 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 261 & ) * rho_air_zw(k-1) & 262 & ) * ddzw(k) * drho_air(k) 313 ! 314 !-- Add momentum flux at model top 315 IF ( use_top_fluxes ) THEN 316 surf_s = surf_def_h(2)%start_index(j,i) 317 surf_e = surf_def_h(2)%end_index(j,i) 318 DO m = surf_s, surf_e 319 320 k = surf_def_h(2)%k(m) 321 322 tend(k,j,i) = tend(k,j,i) & 323 + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k) 324 ENDDO 263 325 ENDIF 264 326 … … 277 339 278 340 USE arrays_3d, & 279 ONLY: ddzu, ddzw, km, tend, u, v, vsws, vswst, w, & 280 drho_air, rho_air_zw 341 ONLY: ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw 281 342 282 343 USE control_parameters, & 283 ONLY: constant_top_momentumflux, use_surface_fluxes, use_top_fluxes 344 ONLY: constant_top_momentumflux, use_surface_fluxes, & 345 use_top_fluxes 284 346 285 347 USE grid_variables, & 286 ONLY: ddx, ddy, ddy2 , fxm, fxp, wall_v348 ONLY: ddx, ddy, ddy2 287 349 288 350 USE indices, & 289 ONLY: nzb, nz b_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff351 ONLY: nzb, nzt, wall_flags_0 290 352 291 353 USE kinds 292 354 355 USE surface_mod, & 356 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 357 surf_usm_v 358 293 359 IMPLICIT NONE 294 360 295 INTEGER(iwp) :: i !< 296 INTEGER(iwp) :: j !< 297 INTEGER(iwp) :: k !< 298 REAL(wp) :: kmxm !< 299 REAL(wp) :: kmxp !< 300 REAL(wp) :: kmzm !< 301 REAL(wp) :: kmzp !< 302 303 REAL(wp), DIMENSION(nzb:nzt+1) :: vsus !< 361 362 INTEGER(iwp) :: i !< running index x direction 363 INTEGER(iwp) :: j !< running index y direction 364 INTEGER(iwp) :: k !< running index z direction 365 INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall 366 INTEGER(iwp) :: m !< running index surface elements 367 INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint 368 INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint 369 370 REAL(wp) :: flag !< flag to mask topography grid points 371 REAL(wp) :: kmxm !< 372 REAL(wp) :: kmxp !< 373 REAL(wp) :: kmzm !< 374 REAL(wp) :: kmzp !< 375 REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface 376 REAL(wp) :: mask_east !< flag to mask vertical surface south of the grid point 377 REAL(wp) :: mask_west !< flag to mask vertical surface north of the grid point 378 REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface 304 379 305 380 ! 306 381 !-- Compute horizontal diffusion 307 DO k = nzb_v_outer(j,i)+1, nzt 382 DO k = nzb+1, nzt 383 ! 384 !-- Predetermine flag to mask topography and wall-bounded grid points. 385 !-- It is sufficient to masked only east- and west-facing surfaces, which 386 !-- need special treatment for the v-component. 387 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 2 ) ) 388 mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 2 ) ) 389 mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 2 ) ) 308 390 ! 309 391 !-- Interpolate eddy diffusivities on staggered gridpoints … … 311 393 kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 312 394 313 tend(k,j,i) = tend(k,j,i) & 314 & + ( kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & 315 & + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 316 & - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 317 & - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 318 & ) * ddx & 319 & + 2.0_wp * ( & 320 & km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 321 & - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 322 & ) * ddy2 395 tend(k,j,i) = tend(k,j,i) + ( & 396 mask_east * kmxp * ( & 397 ( v(k,j,i+1) - v(k,j,i) ) * ddx & 398 + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 399 ) & 400 - mask_west * kmxm * ( & 401 ( v(k,j,i) - v(k,j,i-1) ) * ddx & 402 + ( u(k,j,i) - u(k,j-1,i) ) * ddy & 403 ) & 404 ) * ddx * flag & 405 + 2.0_wp * ( & 406 km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 407 - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 408 ) * ddy2 * flag 323 409 ENDDO 324 410 325 411 ! 326 !-- Wall functions at the left and right walls, respectively 327 IF ( wall_v(j,i) /= 0.0_wp ) THEN 328 329 ! 330 !-- Calculate the horizontal momentum flux v'u' 331 CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i), & 332 vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp ) 333 334 DO k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i) 335 kmxp = 0.25_wp * & 336 ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) ) 337 kmxm = 0.25_wp * & 338 ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) ) 339 340 tend(k,j,i) = tend(k,j,i) & 341 + 2.0_wp * ( & 342 km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) & 343 - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) & 344 ) * ddy2 & 345 + ( fxp(j,i) * ( & 346 kmxp * ( v(k,j,i+1) - v(k,j,i) ) * ddx & 347 + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy & 348 ) & 349 - fxm(j,i) * ( & 350 kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx & 351 + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy & 352 ) & 353 + wall_v(j,i) * vsus(k) & 354 ) * ddx 355 ENDDO 356 ENDIF 357 358 ! 359 !-- Compute vertical diffusion. In case of simulating a Prandtl layer, 360 !-- index k starts at nzb_v_inner+2. 361 DO k = nzb_diff_v(j,i), nzt_diff 412 !-- Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3) 413 !-- surfaces. Note, in the the flat case, loops won't be entered as 414 !-- start_index > end_index. Furtermore, note, no vertical natural surfaces 415 !-- so far. 416 !-- Default-type surfaces 417 DO l = 2, 3 418 surf_s = surf_def_v(l)%start_index(j,i) 419 surf_e = surf_def_v(l)%end_index(j,i) 420 DO m = surf_s, surf_e 421 k = surf_def_v(l)%k(m) 422 tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx 423 ENDDO 424 ENDDO 425 ! 426 !-- Natural-type surfaces 427 DO l = 2, 3 428 surf_s = surf_lsm_v(l)%start_index(j,i) 429 surf_e = surf_lsm_v(l)%end_index(j,i) 430 DO m = surf_s, surf_e 431 k = surf_lsm_v(l)%k(m) 432 tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx 433 ENDDO 434 ENDDO 435 ! 436 !-- Urban-type surfaces 437 DO l = 2, 3 438 surf_s = surf_usm_v(l)%start_index(j,i) 439 surf_e = surf_usm_v(l)%end_index(j,i) 440 DO m = surf_s, surf_e 441 k = surf_usm_v(l)%k(m) 442 tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx 443 ENDDO 444 ENDDO 445 ! 446 !-- Compute vertical diffusion. In case of simulating a surface layer, 447 !-- respective grid diffusive fluxes are masked (flag 8) within this 448 !-- loop, and added further below, else, simple gradient approach is 449 !-- applied. Model top is also mask if top-momentum flux is given. 450 DO k = nzb+1, nzt 451 ! 452 !-- Determine flags to mask topography below and above. Flag 2 is 453 !-- used to mask topography in general, while flag 10 implies also 454 !-- information about use_surface_fluxes. Flag 9 is used to control 455 !-- momentum flux at model top. 456 mask_bottom = MERGE( 1.0_wp, 0.0_wp, & 457 BTEST( wall_flags_0(k-1,j,i), 8 ) ) 458 mask_top = MERGE( 1.0_wp, 0.0_wp, & 459 BTEST( wall_flags_0(k+1,j,i), 8 ) ) * & 460 MERGE( 1.0_wp, 0.0_wp, & 461 BTEST( wall_flags_0(k+1,j,i), 9 ) ) 462 flag = MERGE( 1.0_wp, 0.0_wp, & 463 BTEST( wall_flags_0(k,j,i), 2 ) ) 362 464 ! 363 465 !-- Interpolate eddy diffusivities on staggered gridpoints … … 368 470 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 369 471 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 370 & ) * rho_air_zw(k) 472 & ) * rho_air_zw(k) * mask_top & 371 473 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 372 474 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 373 & ) * rho_air_zw(k-1) 374 & ) * ddzw(k) * drho_air(k) 475 & ) * rho_air_zw(k-1) * mask_bottom & 476 & ) * ddzw(k) * drho_air(k) * flag 375 477 ENDDO 376 478 … … 383 485 !-- other (LES) models showed that the values of the momentum flux becomes 384 486 !-- too large in this case. 385 !-- The term containing w(k-1,..) (see above equation) is removed here386 !-- because the vertical velocity is assumed to be zero at the surface.387 487 IF ( use_surface_fluxes ) THEN 388 k = nzb_v_inner(j,i)+1 389 ! 390 !-- Interpolate eddy diffusivities on staggered gridpoints 391 kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) ) 392 393 tend(k,j,i) = tend(k,j,i) & 394 & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) & 395 & + ( w(k,j,i) - w(k,j-1,i) ) * ddy & 396 & ) * rho_air_zw(k) & 397 & - ( -vsws(j,i) ) & 398 & ) * ddzw(k) * drho_air(k) 488 ! 489 !-- Default-type surfaces, upward-facing 490 surf_s = surf_def_h(0)%start_index(j,i) 491 surf_e = surf_def_h(0)%end_index(j,i) 492 DO m = surf_s, surf_e 493 k = surf_def_h(0)%k(m) 494 495 tend(k,j,i) = tend(k,j,i) & 496 + ( - ( - surf_def_h(0)%vsws(m) ) & 497 ) * ddzw(k) * drho_air(k) 498 ENDDO 499 ! 500 !-- Default-type surfaces, dowward-facing 501 surf_s = surf_def_h(1)%start_index(j,i) 502 surf_e = surf_def_h(1)%end_index(j,i) 503 DO m = surf_s, surf_e 504 k = surf_def_h(1)%k(m) 505 506 tend(k,j,i) = tend(k,j,i) & 507 + ( - surf_def_h(1)%vsws(m) & 508 ) * ddzw(k) * drho_air(k) 509 ENDDO 510 ! 511 !-- Natural-type surfaces, upward-facing 512 surf_s = surf_lsm_h%start_index(j,i) 513 surf_e = surf_lsm_h%end_index(j,i) 514 DO m = surf_s, surf_e 515 k = surf_lsm_h%k(m) 516 517 tend(k,j,i) = tend(k,j,i) & 518 + ( - ( - surf_lsm_h%vsws(m) ) & 519 ) * ddzw(k) * drho_air(k) 520 521 ENDDO 522 ! 523 !-- Urban-type surfaces, upward-facing 524 surf_s = surf_usm_h%start_index(j,i) 525 surf_e = surf_usm_h%end_index(j,i) 526 DO m = surf_s, surf_e 527 k = surf_usm_h%k(m) 528 529 tend(k,j,i) = tend(k,j,i) & 530 + ( - ( - surf_usm_h%vsws(m) ) & 531 ) * ddzw(k) * drho_air(k) 532 533 ENDDO 399 534 ENDIF 400 401 ! 402 !-- Vertical diffusion at the first gridpoint below the top boundary, 403 !-- if the momentum flux at the top is prescribed by the user 404 IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN 405 k = nzt 406 ! 407 !-- Interpolate eddy diffusivities on staggered gridpoints 408 kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) ) 409 410 tend(k,j,i) = tend(k,j,i) & 411 & + ( ( -vswst(j,i) ) & 412 & - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) & 413 & + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy & 414 & ) * rho_air_zw(k-1) & 415 & ) * ddzw(k) * drho_air(k) 535 ! 536 !-- Add momentum flux at model top 537 IF ( use_top_fluxes ) THEN 538 surf_s = surf_def_h(2)%start_index(j,i) 539 surf_e = surf_def_h(2)%end_index(j,i) 540 DO m = surf_s, surf_e 541 542 k = surf_def_h(2)%k(m) 543 544 tend(k,j,i) = tend(k,j,i) & 545 + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k) 546 ENDDO 416 547 ENDIF 417 548
Note: See TracChangeset
for help on using the changeset viewer.
