Changeset 2232 for palm/trunk/SOURCE/diffusion_s.f90
 Timestamp:
 May 30, 2017 5:47:52 PM (4 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/diffusion_s.f90
r2119 r2232 20 20 ! Current revisions: 21 21 !  22 ! 22 ! Adjustments to new topography and surface concept 23 23 ! 24 24 ! Former revisions: … … 110 110 !> Call for all grid points 111 111 !! 112 SUBROUTINE diffusion_s( s, s_flux_b, s_flux_t, wall_s_flux ) 112 SUBROUTINE diffusion_s( s, s_flux_def_h_up, s_flux_def_h_down, & 113 s_flux_t, & 114 s_flux_lsm_h_up, s_flux_usm_h_up, & 115 s_flux_def_v_north, s_flux_def_v_south, & 116 s_flux_def_v_east, s_flux_def_v_west, & 117 s_flux_lsm_v_north, s_flux_lsm_v_south, & 118 s_flux_lsm_v_east, s_flux_lsm_v_west, & 119 s_flux_usm_v_north, s_flux_usm_v_south, & 120 s_flux_usm_v_east, s_flux_usm_v_west ) 113 121 114 122 USE arrays_3d, & … … 119 127 120 128 USE grid_variables, & 121 ONLY: ddx2, ddy2 , fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y129 ONLY: ddx2, ddy2 122 130 123 131 USE indices, & 124 132 ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, & 125 nz b_diff_s_inner, nzb_s_inner, nzb_s_outer, nzt, nzt_diff133 nzt, wall_flags_0 126 134 127 135 USE kinds 128 136 137 USE surface_mod, & 138 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 139 surf_usm_v 140 129 141 IMPLICIT NONE 130 142 131 INTEGER(iwp) :: i !< 132 INTEGER(iwp) :: j !< 133 INTEGER(iwp) :: k !< 134 REAL(wp) :: wall_s_flux(0:4) !< 135 REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_b, s_flux_t !< 143 INTEGER(iwp) :: i !< running index x direction 144 INTEGER(iwp) :: j !< running index y direction 145 INTEGER(iwp) :: k !< running index z direction 146 INTEGER(iwp) :: m !< running index surface elements 147 INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)gridpoint 148 INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)gridpoint 149 150 REAL(wp) :: flag !< flag to mask topography grid points 151 REAL(wp) :: mask_bottom !< flag to mask vertical upwardfacing surface 152 REAL(wp) :: mask_east !< flag to mask vertical surface east of the grid point 153 REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point 154 REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point 155 REAL(wp) :: mask_west !< flag to mask vertical surface west of the grid point 156 REAL(wp) :: mask_top !< flag to mask vertical downwardfacing surface 157 158 REAL(wp), DIMENSION(1:surf_def_v(0)%ns) :: s_flux_def_v_north !< flux at northfacing vertical defaulttype surfaces 159 REAL(wp), DIMENSION(1:surf_def_v(1)%ns) :: s_flux_def_v_south !< flux at southfacing vertical defaulttype surfaces 160 REAL(wp), DIMENSION(1:surf_def_v(2)%ns) :: s_flux_def_v_east !< flux at eastfacing vertical defaulttype surfaces 161 REAL(wp), DIMENSION(1:surf_def_v(3)%ns) :: s_flux_def_v_west !< flux at westfacing vertical defaulttype surfaces 162 REAL(wp), DIMENSION(1:surf_def_h(0)%ns) :: s_flux_def_h_up !< flux at horizontal upwardfacing defaulttype surfaces 163 REAL(wp), DIMENSION(1:surf_def_h(1)%ns) :: s_flux_def_h_down !< flux at horizontal donwwardfacing defaulttype surfaces 164 REAL(wp), DIMENSION(1:surf_lsm_h%ns) :: s_flux_lsm_h_up !< flux at horizontal upwardfacing naturaltype surfaces 165 REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) :: s_flux_lsm_v_north !< flux at northfacing vertical naturaltype surfaces 166 REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) :: s_flux_lsm_v_south !< flux at southfacing vertical naturaltype surfaces 167 REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) :: s_flux_lsm_v_east !< flux at eastfacing vertical naturaltype surfaces 168 REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) :: s_flux_lsm_v_west !< flux at westfacing vertical naturaltype surfaces 169 REAL(wp), DIMENSION(1:surf_usm_h%ns) :: s_flux_usm_h_up !< flux at horizontal upwardfacing urbantype surfaces 170 REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) :: s_flux_usm_v_north !< flux at northfacing vertical urbantype surfaces 171 REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) :: s_flux_usm_v_south !< flux at southfacing vertical urbantype surfaces 172 REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) :: s_flux_usm_v_east !< flux at eastfacing vertical urbantype surfaces 173 REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) :: s_flux_usm_v_west !< flux at westfacing vertical urbantype surfaces 174 REAL(wp), DIMENSION(1:surf_def_h(2)%ns) :: s_flux_t !< flux at model top 136 175 #if defined( __nopointer ) 137 176 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s !< … … 144 183 ! 145 184 ! Compute horizontal diffusion 146 DO k = nzb_s_outer(j,i)+1, nzt 185 DO k = nzb+1, nzt 186 ! 187 ! Predetermine flag to mask topography and wallbounded grid points 188 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 189 ! 190 ! Predetermine flag to mask wallbounded grid points, equivalent to 191 ! former s_outer array 192 mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i1), 0 ) ) 193 mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 0 ) ) 194 mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j1,i), 0 ) ) 195 mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 0 ) ) 147 196 148 197 tend(k,j,i) = tend(k,j,i) & 149 198 + 0.5_wp * ( & 150 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)s(k,j,i) ) & 151  ( kh(k,j,i) + kh(k,j,i1) ) * ( s(k,j,i)s(k,j,i1) ) & 152 ) * ddx2 & 199 mask_east * ( kh(k,j,i) + kh(k,j,i+1) ) & 200 * ( s(k,j,i+1)  s(k,j,i) ) & 201  mask_west * ( kh(k,j,i) + kh(k,j,i1) ) & 202 * ( s(k,j,i)  s(k,j,i1) ) & 203 ) * ddx2 * flag & 153 204 + 0.5_wp * ( & 154 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)s(k,j,i) ) & 155  ( kh(k,j,i) + kh(k,j1,i) ) * ( s(k,j,i)s(k,j1,i) ) & 156 ) * ddy2 157 ENDDO 158 159 ! 160 ! Apply prescribed horizontal wall heatflux where necessary 161 IF ( ( wall_w_x(j,i) /= 0.0_wp ) .OR. ( wall_w_y(j,i) /= 0.0_wp ) & 162 ) THEN 163 DO k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i) 164 205 mask_north * ( kh(k,j,i) + kh(k,j+1,i) ) & 206 * ( s(k,j+1,i)  s(k,j,i) ) & 207  mask_south * ( kh(k,j,i) + kh(k,j1,i) ) & 208 * ( s(k,j,i)  s(k,j1,i) ) & 209 ) * ddy2 * flag 210 ENDDO 211 212 ! 213 ! Apply prescribed horizontal wall heatflux where necessary. First, 214 ! determine start and end index for respective (j,i)index. Please 215 ! note, in the flat case following loop will not be entered, as 216 ! surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces 217 ! so far. 218 ! First, for defaulttype surfaces 219 ! Northfacing vertical defaulttype surfaces 220 surf_s = surf_def_v(0)%start_index(j,i) 221 surf_e = surf_def_v(0)%end_index(j,i) 222 DO m = surf_s, surf_e 223 k = surf_def_v(0)%k(m) 224 tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy2 225 ENDDO 226 ! 227 ! Southfacing vertical defaulttype surfaces 228 surf_s = surf_def_v(1)%start_index(j,i) 229 surf_e = surf_def_v(1)%end_index(j,i) 230 DO m = surf_s, surf_e 231 k = surf_def_v(1)%k(m) 232 tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy2 233 ENDDO 234 ! 235 ! Eastfacing vertical defaulttype surfaces 236 surf_s = surf_def_v(2)%start_index(j,i) 237 surf_e = surf_def_v(2)%end_index(j,i) 238 DO m = surf_s, surf_e 239 k = surf_def_v(2)%k(m) 240 tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx2 241 ENDDO 242 ! 243 ! Westfacing vertical defaulttype surfaces 244 surf_s = surf_def_v(3)%start_index(j,i) 245 surf_e = surf_def_v(3)%end_index(j,i) 246 DO m = surf_s, surf_e 247 k = surf_def_v(3)%k(m) 248 tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx2 249 ENDDO 250 ! 251 ! Now, for naturaltype surfaces. 252 ! Northfacing 253 surf_s = surf_lsm_v(0)%start_index(j,i) 254 surf_e = surf_lsm_v(0)%end_index(j,i) 255 DO m = surf_s, surf_e 256 k = surf_lsm_v(0)%k(m) 257 tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy2 258 ENDDO 259 ! 260 ! Southfacing 261 surf_s = surf_lsm_v(1)%start_index(j,i) 262 surf_e = surf_lsm_v(1)%end_index(j,i) 263 DO m = surf_s, surf_e 264 k = surf_lsm_v(1)%k(m) 265 tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy2 266 ENDDO 267 ! 268 ! Eastfacing 269 surf_s = surf_lsm_v(2)%start_index(j,i) 270 surf_e = surf_lsm_v(2)%end_index(j,i) 271 DO m = surf_s, surf_e 272 k = surf_lsm_v(2)%k(m) 273 tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx2 274 ENDDO 275 ! 276 ! Westfacing 277 surf_s = surf_lsm_v(3)%start_index(j,i) 278 surf_e = surf_lsm_v(3)%end_index(j,i) 279 DO m = surf_s, surf_e 280 k = surf_lsm_v(3)%k(m) 281 tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx2 282 ENDDO 283 ! 284 ! Now, for urbantype surfaces. 285 ! Northfacing 286 surf_s = surf_usm_v(0)%start_index(j,i) 287 surf_e = surf_usm_v(0)%end_index(j,i) 288 DO m = surf_s, surf_e 289 k = surf_usm_v(0)%k(m) 290 tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy2 291 ENDDO 292 ! 293 ! Southfacing 294 surf_s = surf_usm_v(1)%start_index(j,i) 295 surf_e = surf_usm_v(1)%end_index(j,i) 296 DO m = surf_s, surf_e 297 k = surf_usm_v(1)%k(m) 298 tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy2 299 ENDDO 300 ! 301 ! Eastfacing 302 surf_s = surf_usm_v(2)%start_index(j,i) 303 surf_e = surf_usm_v(2)%end_index(j,i) 304 DO m = surf_s, surf_e 305 k = surf_usm_v(2)%k(m) 306 tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx2 307 ENDDO 308 ! 309 ! Westfacing 310 surf_s = surf_usm_v(3)%start_index(j,i) 311 surf_e = surf_usm_v(3)%end_index(j,i) 312 DO m = surf_s, surf_e 313 k = surf_usm_v(3)%k(m) 314 tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx2 315 ENDDO 316 317 ! 318 ! Compute vertical diffusion. In case that surface fluxes have been 319 ! prescribed or computed at bottom and/or top, index k starts/ends at 320 ! nzb+2 or nzt1, respectively. Model top is also mask if top flux 321 ! is given. 322 DO k = nzb+1, nzt 323 ! 324 ! Determine flags to mask topography below and above. Flag 0 is 325 ! used to mask topography in general, and flag 8 implies 326 ! information about use_surface_fluxes. Flag 9 is used to control 327 ! flux at model top. 328 mask_bottom = MERGE( 1.0_wp, 0.0_wp, & 329 BTEST( wall_flags_0(k1,j,i), 8 ) ) 330 mask_top = MERGE( 1.0_wp, 0.0_wp, & 331 BTEST( wall_flags_0(k+1,j,i), 8 ) ) * & 332 MERGE( 1.0_wp, 0.0_wp, & 333 BTEST( wall_flags_0(k+1,j,i), 9 ) ) 334 flag = MERGE( 1.0_wp, 0.0_wp, & 335 BTEST( wall_flags_0(k,j,i), 0 ) ) 336 337 tend(k,j,i) = tend(k,j,i) & 338 + 0.5_wp * ( & 339 ( kh(k,j,i) + kh(k+1,j,i) ) * & 340 ( s(k+1,j,i)s(k,j,i) ) * ddzu(k+1) & 341 * rho_air_zw(k) & 342 * mask_top & 343  ( kh(k,j,i) + kh(k1,j,i) ) * & 344 ( s(k,j,i)s(k1,j,i) ) * ddzu(k) & 345 * rho_air_zw(k1) & 346 * mask_bottom & 347 ) * ddzw(k) * drho_air(k) & 348 * flag 349 ENDDO 350 351 ! 352 ! Vertical diffusion at horizontal walls. 353 IF ( use_surface_fluxes ) THEN 354 ! 355 ! Defaulttype surfaces, upwardfacing 356 surf_s = surf_def_h(0)%start_index(j,i) 357 surf_e = surf_def_h(0)%end_index(j,i) 358 DO m = surf_s, surf_e 359 360 k = surf_def_h(0)%k(m) 361 tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) & 362 * ddzw(k) * drho_air(k) 363 364 ENDDO 365 ! 366 ! Defaulttype surfaces, downwardfacing 367 surf_s = surf_def_h(1)%start_index(j,i) 368 surf_e = surf_def_h(1)%end_index(j,i) 369 DO m = surf_s, surf_e 370 371 k = surf_def_h(1)%k(m) 372 tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) & 373 * ddzw(k) * drho_air(k) 374 375 ENDDO 376 ! 377 ! Naturaltype surfaces, upwardfacing 378 surf_s = surf_lsm_h%start_index(j,i) 379 surf_e = surf_lsm_h%end_index(j,i) 380 DO m = surf_s, surf_e 381 382 k = surf_lsm_h%k(m) 383 tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) & 384 * ddzw(k) * drho_air(k) 385 386 ENDDO 387 ! 388 ! Urbantype surfaces, upwardfacing 389 surf_s = surf_usm_h%start_index(j,i) 390 surf_e = surf_usm_h%end_index(j,i) 391 DO m = surf_s, surf_e 392 393 k = surf_usm_h%k(m) 394 tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) & 395 * ddzw(k) * drho_air(k) 396 397 ENDDO 398 399 ENDIF 400 ! 401 ! Vertical diffusion at the last computational gridpoint along zdirection 402 IF ( use_top_fluxes ) THEN 403 surf_s = surf_def_h(2)%start_index(j,i) 404 surf_e = surf_def_h(2)%end_index(j,i) 405 DO m = surf_s, surf_e 406 407 k = surf_def_h(2)%k(m) 165 408 tend(k,j,i) = tend(k,j,i) & 166 + ( fwxp(j,i) * 0.5_wp * & 167 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)s(k,j,i) ) & 168 + ( 1.0_wp  fwxp(j,i) ) * wall_s_flux(1) & 169 fwxm(j,i) * 0.5_wp * & 170 ( kh(k,j,i) + kh(k,j,i1) ) * ( s(k,j,i)s(k,j,i1) ) & 171 + ( 1.0_wp  fwxm(j,i) ) * wall_s_flux(2) & 172 ) * ddx2 & 173 + ( fwyp(j,i) * 0.5_wp * & 174 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)s(k,j,i) ) & 175 + ( 1.0_wp  fwyp(j,i) ) * wall_s_flux(3) & 176 fwym(j,i) * 0.5_wp * & 177 ( kh(k,j,i) + kh(k,j1,i) ) * ( s(k,j,i)s(k,j1,i) ) & 178 + ( 1.0_wp  fwym(j,i) ) * wall_s_flux(4) & 179 ) * ddy2 409 + (  s_flux_t(m) ) * ddzw(k) * drho_air(k) 180 410 ENDDO 181 411 ENDIF 182 412 183 !184 ! Compute vertical diffusion. In case that surface fluxes have been185 ! prescribed or computed at bottom and/or top, index k starts/ends at186 ! nzb+2 or nzt1, respectively.187 DO k = nzb_diff_s_inner(j,i), nzt_diff188 189 tend(k,j,i) = tend(k,j,i) &190 + 0.5_wp * ( &191 ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)s(k,j,i) ) * ddzu(k+1) &192 * rho_air_zw(k) &193  ( kh(k,j,i) + kh(k1,j,i) ) * ( s(k,j,i)s(k1,j,i) ) * ddzu(k) &194 * rho_air_zw(k1) &195 ) * ddzw(k) * drho_air(k)196 ENDDO197 198 !199 ! Vertical diffusion at the first computational gridpoint along200 ! zdirection201 IF ( use_surface_fluxes ) THEN202 203 k = nzb_s_inner(j,i)+1204 205 tend(k,j,i) = tend(k,j,i) &206 + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) ) &207 * ( s(k+1,j,i)s(k,j,i) ) &208 * ddzu(k+1) &209 * rho_air_zw(k) &210 + s_flux_b(j,i) &211 ) * ddzw(k) * drho_air(k)212 213 ENDIF214 215 !216 ! Vertical diffusion at the last computational gridpoint along217 ! zdirection218 IF ( use_top_fluxes ) THEN219 220 k = nzt221 222 tend(k,j,i) = tend(k,j,i) &223 + (  s_flux_t(j,i) &224  0.5_wp * ( kh(k1,j,i)+kh(k,j,i) )&225 * ( s(k,j,i)s(k1,j,i) ) &226 * ddzu(k) &227 * rho_air_zw(k1) &228 ) * ddzw(k) * drho_air(k)229 230 ENDIF231 232 413 ENDDO 233 414 ENDDO 234 415 235 416 END SUBROUTINE diffusion_s 236 237 417 238 418 !! … … 241 421 !> Call for grid point i,j 242 422 !! 243 SUBROUTINE diffusion_s_ij( i, j, s, s_flux_b, s_flux_t, wall_s_flux ) 423 SUBROUTINE diffusion_s_ij( i, j, s, & 424 s_flux_def_h_up, s_flux_def_h_down, & 425 s_flux_t, & 426 s_flux_lsm_h_up, s_flux_usm_h_up, & 427 s_flux_def_v_north, s_flux_def_v_south, & 428 s_flux_def_v_east, s_flux_def_v_west, & 429 s_flux_lsm_v_north, s_flux_lsm_v_south, & 430 s_flux_lsm_v_east, s_flux_lsm_v_west, & 431 s_flux_usm_v_north, s_flux_usm_v_south, & 432 s_flux_usm_v_east, s_flux_usm_v_west ) 244 433 245 434 USE arrays_3d, & … … 250 439 251 440 USE grid_variables, & 252 ONLY: ddx2, ddy2 , fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y441 ONLY: ddx2, ddy2 253 442 254 443 USE indices, & 255 ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_diff_s_inner, nzb_s_inner, & 256 nzb_s_outer, nzt, nzt_diff 444 ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0 257 445 258 446 USE kinds 259 447 448 USE surface_mod, & 449 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 450 surf_usm_v 451 260 452 IMPLICIT NONE 261 453 262 INTEGER(iwp) :: i !< 263 INTEGER(iwp) :: j !< 264 INTEGER(iwp) :: k !< 265 REAL(wp) :: wall_s_flux(0:4) !< 266 REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_b !< 267 REAL(wp), DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_t !< 454 INTEGER(iwp) :: i !< running index x direction 455 INTEGER(iwp) :: j !< running index y direction 456 INTEGER(iwp) :: k !< running index z direction 457 INTEGER(iwp) :: m !< running index surface elements 458 INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)gridpoint 459 INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)gridpoint 460 461 REAL(wp) :: flag !< flag to mask topography grid points 462 REAL(wp) :: mask_bottom !< flag to mask vertical upwardfacing surface 463 REAL(wp) :: mask_east !< flag to mask vertical surface east of the grid point 464 REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point 465 REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point 466 REAL(wp) :: mask_west !< flag to mask vertical surface west of the grid point 467 REAL(wp) :: mask_top !< flag to mask vertical downwardfacing surface 468 469 REAL(wp), DIMENSION(1:surf_def_v(0)%ns) :: s_flux_def_v_north !< flux at northfacing vertical defaulttype surfaces 470 REAL(wp), DIMENSION(1:surf_def_v(1)%ns) :: s_flux_def_v_south !< flux at southfacing vertical defaulttype surfaces 471 REAL(wp), DIMENSION(1:surf_def_v(2)%ns) :: s_flux_def_v_east !< flux at eastfacing vertical defaulttype surfaces 472 REAL(wp), DIMENSION(1:surf_def_v(3)%ns) :: s_flux_def_v_west !< flux at westfacing vertical defaulttype surfaces 473 REAL(wp), DIMENSION(1:surf_def_h(0)%ns) :: s_flux_def_h_up !< flux at horizontal upwardfacing defaulttype surfaces 474 REAL(wp), DIMENSION(1:surf_def_h(1)%ns) :: s_flux_def_h_down !< flux at horizontal donwwardfacing defaulttype surfaces 475 REAL(wp), DIMENSION(1:surf_lsm_h%ns) :: s_flux_lsm_h_up !< flux at horizontal upwardfacing naturaltype surfaces 476 REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) :: s_flux_lsm_v_north !< flux at northfacing vertical urbantype surfaces 477 REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) :: s_flux_lsm_v_south !< flux at southfacing vertical urbantype surfaces 478 REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) :: s_flux_lsm_v_east !< flux at eastfacing vertical urbantype surfaces 479 REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) :: s_flux_lsm_v_west !< flux at westfacing vertical urbantype surfaces 480 REAL(wp), DIMENSION(1:surf_usm_h%ns) :: s_flux_usm_h_up !< flux at horizontal upwardfacing urbantype surfaces 481 REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) :: s_flux_usm_v_north !< flux at northfacing vertical urbantype surfaces 482 REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) :: s_flux_usm_v_south !< flux at southfacing vertical urbantype surfaces 483 REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) :: s_flux_usm_v_east !< flux at eastfacing vertical urbantype surfaces 484 REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) :: s_flux_usm_v_west !< flux at westfacing vertical urbantype surfaces 485 REAL(wp), DIMENSION(1:surf_def_h(2)%ns) :: s_flux_t !< flux at model top 268 486 #if defined( __nopointer ) 269 487 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s !< … … 274 492 ! 275 493 ! Compute horizontal diffusion 276 DO k = nzb_s_outer(j,i)+1, nzt 494 DO k = nzb+1, nzt 495 ! 496 ! Predetermine flag to mask topography and wallbounded grid points 497 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 498 ! 499 ! Predetermine flag to mask wallbounded grid points, equivalent to 500 ! former s_outer array 501 mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i1), 0 ) ) 502 mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 0 ) ) 503 mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j1,i), 0 ) ) 504 mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 0 ) ) 505 ! 506 ! Finally, determine flag to mask both topography itself as well 507 ! as wallbounded grid points, which will be treated further below 277 508 278 509 tend(k,j,i) = tend(k,j,i) & 279 510 + 0.5_wp * ( & 280 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)s(k,j,i) ) & 281  ( kh(k,j,i) + kh(k,j,i1) ) * ( s(k,j,i)s(k,j,i1) ) & 282 ) * ddx2 & 511 mask_east * ( kh(k,j,i) + kh(k,j,i+1) ) & 512 * ( s(k,j,i+1)  s(k,j,i) ) & 513  mask_west * ( kh(k,j,i) + kh(k,j,i1) ) & 514 * ( s(k,j,i)  s(k,j,i1) ) & 515 ) * ddx2 * flag & 283 516 + 0.5_wp * ( & 284 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)s(k,j,i) ) & 285  ( kh(k,j,i) + kh(k,j1,i) ) * ( s(k,j,i)s(k,j1,i) ) & 286 ) * ddy2 287 ENDDO 288 289 ! 290 ! Apply prescribed horizontal wall heatflux where necessary 291 IF ( ( wall_w_x(j,i) /= 0.0_wp ) .OR. ( wall_w_y(j,i) /= 0.0_wp ) ) & 292 THEN 293 DO k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i) 294 295 tend(k,j,i) = tend(k,j,i) & 296 + ( fwxp(j,i) * 0.5_wp * & 297 ( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)s(k,j,i) ) & 298 + ( 1.0_wp  fwxp(j,i) ) * wall_s_flux(1) & 299 fwxm(j,i) * 0.5_wp * & 300 ( kh(k,j,i) + kh(k,j,i1) ) * ( s(k,j,i)s(k,j,i1) ) & 301 + ( 1.0_wp  fwxm(j,i) ) * wall_s_flux(2) & 302 ) * ddx2 & 303 + ( fwyp(j,i) * 0.5_wp * & 304 ( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)s(k,j,i) ) & 305 + ( 1.0_wp  fwyp(j,i) ) * wall_s_flux(3) & 306 fwym(j,i) * 0.5_wp * & 307 ( kh(k,j,i) + kh(k,j1,i) ) * ( s(k,j,i)s(k,j1,i) ) & 308 + ( 1.0_wp  fwym(j,i) ) * wall_s_flux(4) & 309 ) * ddy2 517 mask_north * ( kh(k,j,i) + kh(k,j+1,i) ) & 518 * ( s(k,j+1,i)  s(k,j,i) ) & 519  mask_south * ( kh(k,j,i) + kh(k,j1,i) ) & 520 * ( s(k,j,i)  s(k,j1,i) ) & 521 ) * ddy2 * flag 522 ENDDO 523 524 ! 525 ! Apply prescribed horizontal wall heatflux where necessary. First, 526 ! determine start and end index for respective (j,i)index. Please 527 ! note, in the flat case following loops will not be entered, as 528 ! surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces 529 ! so far. 530 ! First, for defaulttype surfaces 531 ! Northfacing vertical defaulttype surfaces 532 surf_s = surf_def_v(0)%start_index(j,i) 533 surf_e = surf_def_v(0)%end_index(j,i) 534 DO m = surf_s, surf_e 535 k = surf_def_v(0)%k(m) 536 tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy2 537 ENDDO 538 ! 539 ! Southfacing vertical defaulttype surfaces 540 surf_s = surf_def_v(1)%start_index(j,i) 541 surf_e = surf_def_v(1)%end_index(j,i) 542 DO m = surf_s, surf_e 543 k = surf_def_v(1)%k(m) 544 tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy2 545 ENDDO 546 ! 547 ! Eastfacing vertical defaulttype surfaces 548 surf_s = surf_def_v(2)%start_index(j,i) 549 surf_e = surf_def_v(2)%end_index(j,i) 550 DO m = surf_s, surf_e 551 k = surf_def_v(2)%k(m) 552 tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx2 553 ENDDO 554 ! 555 ! Westfacing vertical defaulttype surfaces 556 surf_s = surf_def_v(3)%start_index(j,i) 557 surf_e = surf_def_v(3)%end_index(j,i) 558 DO m = surf_s, surf_e 559 k = surf_def_v(3)%k(m) 560 tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx2 561 ENDDO 562 ! 563 ! Now, for naturaltype surfaces 564 ! Northfacing 565 surf_s = surf_lsm_v(0)%start_index(j,i) 566 surf_e = surf_lsm_v(0)%end_index(j,i) 567 DO m = surf_s, surf_e 568 k = surf_lsm_v(0)%k(m) 569 tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy2 570 ENDDO 571 ! 572 ! Southfacing 573 surf_s = surf_lsm_v(1)%start_index(j,i) 574 surf_e = surf_lsm_v(1)%end_index(j,i) 575 DO m = surf_s, surf_e 576 k = surf_lsm_v(1)%k(m) 577 tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy2 578 ENDDO 579 ! 580 ! Eastfacing 581 surf_s = surf_lsm_v(2)%start_index(j,i) 582 surf_e = surf_lsm_v(2)%end_index(j,i) 583 DO m = surf_s, surf_e 584 k = surf_lsm_v(2)%k(m) 585 tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx2 586 ENDDO 587 ! 588 ! Westfacing 589 surf_s = surf_lsm_v(3)%start_index(j,i) 590 surf_e = surf_lsm_v(3)%end_index(j,i) 591 DO m = surf_s, surf_e 592 k = surf_lsm_v(3)%k(m) 593 tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx2 594 ENDDO 595 ! 596 ! Now, for urbantype surfaces 597 ! Northfacing 598 surf_s = surf_usm_v(0)%start_index(j,i) 599 surf_e = surf_usm_v(0)%end_index(j,i) 600 DO m = surf_s, surf_e 601 k = surf_usm_v(0)%k(m) 602 tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy2 603 ENDDO 604 ! 605 ! Southfacing 606 surf_s = surf_usm_v(1)%start_index(j,i) 607 surf_e = surf_usm_v(1)%end_index(j,i) 608 DO m = surf_s, surf_e 609 k = surf_usm_v(1)%k(m) 610 tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy2 611 ENDDO 612 ! 613 ! Eastfacing 614 surf_s = surf_usm_v(2)%start_index(j,i) 615 surf_e = surf_usm_v(2)%end_index(j,i) 616 DO m = surf_s, surf_e 617 k = surf_usm_v(2)%k(m) 618 tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx2 619 ENDDO 620 ! 621 ! Westfacing 622 surf_s = surf_usm_v(3)%start_index(j,i) 623 surf_e = surf_usm_v(3)%end_index(j,i) 624 DO m = surf_s, surf_e 625 k = surf_usm_v(3)%k(m) 626 tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx2 627 ENDDO 628 629 630 ! 631 ! Compute vertical diffusion. In case that surface fluxes have been 632 ! prescribed or computed at bottom and/or top, index k starts/ends at 633 ! nzb+2 or nzt1, respectively. Model top is also mask if top flux 634 ! is given. 635 DO k = nzb+1, nzt 636 ! 637 ! Determine flags to mask topography below and above. Flag 0 is 638 ! used to mask topography in general, and flag 8 implies 639 ! information about use_surface_fluxes. Flag 9 is used to control 640 ! flux at model top. 641 mask_bottom = MERGE( 1.0_wp, 0.0_wp, & 642 BTEST( wall_flags_0(k1,j,i), 8 ) ) 643 mask_top = MERGE( 1.0_wp, 0.0_wp, & 644 BTEST( wall_flags_0(k+1,j,i), 8 ) ) * & 645 MERGE( 1.0_wp, 0.0_wp, & 646 BTEST( wall_flags_0(k+1,j,i), 9 ) ) 647 flag = MERGE( 1.0_wp, 0.0_wp, & 648 BTEST( wall_flags_0(k,j,i), 0 ) ) 649 650 tend(k,j,i) = tend(k,j,i) & 651 + 0.5_wp * ( & 652 ( kh(k,j,i) + kh(k+1,j,i) ) * & 653 ( s(k+1,j,i)s(k,j,i) ) * ddzu(k+1) & 654 * rho_air_zw(k) & 655 * mask_top & 656  ( kh(k,j,i) + kh(k1,j,i) ) * & 657 ( s(k,j,i)s(k1,j,i) ) * ddzu(k) & 658 * rho_air_zw(k1) & 659 * mask_bottom & 660 ) * ddzw(k) * drho_air(k) & 661 * flag 662 ENDDO 663 664 ! 665 ! Vertical diffusion at horizontal walls. 666 ! TO DO: Adjust for downward facing walls and mask already in main loop 667 IF ( use_surface_fluxes ) THEN 668 ! 669 ! Defaulttype surfaces, upwardfacing 670 surf_s = surf_def_h(0)%start_index(j,i) 671 surf_e = surf_def_h(0)%end_index(j,i) 672 DO m = surf_s, surf_e 673 674 k = surf_def_h(0)%k(m) 675 676 tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) & 677 * ddzw(k) * drho_air(k) 678 ENDDO 679 ! 680 ! Defaulttype surfaces, downwardfacing 681 surf_s = surf_def_h(1)%start_index(j,i) 682 surf_e = surf_def_h(1)%end_index(j,i) 683 DO m = surf_s, surf_e 684 685 k = surf_def_h(1)%k(m) 686 687 tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) & 688 * ddzw(k) * drho_air(k) 689 ENDDO 690 ! 691 ! Naturaltype surfaces, upwardfacing 692 surf_s = surf_lsm_h%start_index(j,i) 693 surf_e = surf_lsm_h%end_index(j,i) 694 DO m = surf_s, surf_e 695 k = surf_lsm_h%k(m) 696 697 tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) & 698 * ddzw(k) * drho_air(k) 699 ENDDO 700 ! 701 ! Urbantype surfaces, upwardfacing 702 surf_s = surf_usm_h%start_index(j,i) 703 surf_e = surf_usm_h%end_index(j,i) 704 DO m = surf_s, surf_e 705 k = surf_usm_h%k(m) 706 707 tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) & 708 * ddzw(k) * drho_air(k) 310 709 ENDDO 311 710 ENDIF 312 313 !314 ! Compute vertical diffusion. In case that surface fluxes have been315 ! prescribed or computed at bottom and/or top, index k starts/ends at316 ! nzb+2 or nzt1, respectively.317 DO k = nzb_diff_s_inner(j,i), nzt_diff318 319 tend(k,j,i) = tend(k,j,i) &320 + 0.5_wp * ( &321 ( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)s(k,j,i) ) * ddzu(k+1) &322 * rho_air_zw(k) &323  ( kh(k,j,i) + kh(k1,j,i) ) * ( s(k,j,i)s(k1,j,i) ) * ddzu(k) &324 * rho_air_zw(k1) &325 ) * ddzw(k) * drho_air(k)326 ENDDO327 328 !329 ! Vertical diffusion at the first computational gridpoint along zdirection330 IF ( use_surface_fluxes ) THEN331 332 k = nzb_s_inner(j,i)+1333 334 tend(k,j,i) = tend(k,j,i) + ( 0.5_wp * ( kh(k,j,i)+kh(k+1,j,i) ) &335 * ( s(k+1,j,i)s(k,j,i) ) &336 * ddzu(k+1) &337 * rho_air_zw(k) &338 + s_flux_b(j,i) &339 ) * ddzw(k) * drho_air(k)340 341 ENDIF342 343 711 ! 344 712 ! Vertical diffusion at the last computational gridpoint along zdirection 345 713 IF ( use_top_fluxes ) THEN 346 347 k = nzt 348 349 tend(k,j,i) = tend(k,j,i) + (  s_flux_t(j,i) & 350  0.5_wp * ( kh(k1,j,i)+kh(k,j,i) ) & 351 * ( s(k,j,i)s(k1,j,i) ) & 352 * ddzu(k) & 353 * rho_air_zw(k1) & 354 ) * ddzw(k) * drho_air(k) 355 714 surf_s = surf_def_h(2)%start_index(j,i) 715 surf_e = surf_def_h(2)%end_index(j,i) 716 DO m = surf_s, surf_e 717 718 k = surf_def_h(2)%k(m) 719 tend(k,j,i) = tend(k,j,i) & 720 + (  s_flux_t(m) ) * ddzw(k) * drho_air(k) 721 ENDDO 356 722 ENDIF 357 723
Note: See TracChangeset
for help on using the changeset viewer.