Changeset 2232 for palm/trunk/SOURCE/boundary_conds.f90
- Timestamp:
- May 30, 2017 5:47:52 PM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/boundary_conds.f90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/boundary_conds.f90
r2119 r2232 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Set boundary conditions on topography top using flag method. 23 23 ! 24 24 ! Former revisions: … … 168 168 USE indices, & 169 169 ONLY: nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, & 170 nzb, nz b_s_inner, nzb_w_inner, nzt170 nzb, nzt, wall_flags_0 171 171 172 172 USE kinds … … 177 177 ONLY : nesting_mode 178 178 179 USE surface_mod, & 180 ONLY : bc_h 179 181 180 182 IMPLICIT NONE 181 183 182 INTEGER(iwp) :: i !< 183 INTEGER(iwp) :: j !< 184 INTEGER(iwp) :: k !< 184 INTEGER(iwp) :: i !< grid index x direction 185 INTEGER(iwp) :: j !< grid index y direction 186 INTEGER(iwp) :: k !< grid index z direction 187 INTEGER(iwp) :: kb !< variable to set respective boundary value, depends on facing. 188 INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls 189 INTEGER(iwp) :: m !< running index surface elements 185 190 186 191 REAL(wp) :: c_max !< … … 194 199 v_p(nzb,:,:) = v_p(nzb+1,:,:) 195 200 ENDIF 196 197 DO i = nxlg, nxrg 198 DO j = nysg, nyng 199 w_p(nzb_w_inner(j,i),j,i) = 0.0_wp 201 ! 202 !-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case 203 !-- of downward-facing surfaces. 204 DO l = 0, 1 205 ! 206 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set, 207 !-- for downward-facing surfaces at topography bottom (k+1). 208 kb = MERGE( -1, 1, l == 0 ) 209 !$OMP PARALLEL DO PRIVATE( i, j, k ) 210 DO m = 1, bc_h(l)%ns 211 i = bc_h(l)%i(m) 212 j = bc_h(l)%j(m) 213 k = bc_h(l)%k(m) 214 w_p(k+kb,j,i) = 0.0_wp 200 215 ENDDO 201 216 ENDDO … … 216 231 217 232 ! 218 !-- Temperature at bottom boundary.233 !-- Temperature at bottom and top boundary. 219 234 !-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by 220 235 !-- the sea surface temperature of the coupled ocean model. 236 !-- Dirichlet 221 237 IF ( ibc_pt_b == 0 ) THEN 222 DO i = nxlg, nxrg 223 DO j = nysg, nyng 224 pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i) 238 DO l = 0, 1 239 ! 240 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set, 241 !-- for downward-facing surfaces at topography bottom (k+1). 242 kb = MERGE( -1, 1, l == 0 ) 243 !$OMP PARALLEL DO PRIVATE( i, j, k ) 244 DO m = 1, bc_h(l)%ns 245 i = bc_h(l)%i(m) 246 j = bc_h(l)%j(m) 247 k = bc_h(l)%k(m) 248 pt_p(k+kb,j,i) = pt(k+kb,j,i) 225 249 ENDDO 226 250 ENDDO 251 ! 252 !-- Neumann, zero-gradient 227 253 ELSEIF ( ibc_pt_b == 1 ) THEN 228 DO i = nxlg, nxrg 229 DO j = nysg, nyng 230 pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i) 254 DO l = 0, 1 255 ! 256 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set, 257 !-- for downward-facing surfaces at topography bottom (k+1). 258 kb = MERGE( -1, 1, l == 0 ) 259 !$OMP PARALLEL DO PRIVATE( i, j, k ) 260 DO m = 1, bc_h(l)%ns 261 i = bc_h(l)%i(m) 262 j = bc_h(l)%j(m) 263 k = bc_h(l)%k(m) 264 pt_p(k+kb,j,i) = pt_p(k,j,i) 231 265 ENDDO 232 266 ENDDO … … 253 287 !-- Generally Neumann conditions with de/dz=0 are assumed 254 288 IF ( .NOT. constant_diffusion ) THEN 255 DO i = nxlg, nxrg 256 DO j = nysg, nyng 257 e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i) 289 290 DO l = 0, 1 291 ! 292 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set, 293 !-- for downward-facing surfaces at topography bottom (k+1). 294 kb = MERGE( -1, 1, l == 0 ) 295 !$OMP PARALLEL DO PRIVATE( i, j, k ) 296 DO m = 1, bc_h(l)%ns 297 i = bc_h(l)%i(m) 298 j = bc_h(l)%j(m) 299 k = bc_h(l)%k(m) 300 e_p(k+kb,j,i) = e_p(k,j,i) 258 301 ENDDO 259 302 ENDDO 303 260 304 IF ( .NOT. nest_domain ) THEN 261 305 e_p(nzt+1,:,:) = e_p(nzt,:,:) … … 268 312 ! 269 313 !-- Bottom boundary: Neumann condition because salinity flux is always 270 !-- given 271 DO i = nxlg, nxrg 272 DO j = nysg, nyng 273 sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i) 314 !-- given. 315 DO l = 0, 1 316 ! 317 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set, 318 !-- for downward-facing surfaces at topography bottom (k+1). 319 kb = MERGE( -1, 1, l == 0 ) 320 !$OMP PARALLEL DO PRIVATE( i, j, k ) 321 DO m = 1, bc_h(l)%ns 322 i = bc_h(l)%i(m) 323 j = bc_h(l)%j(m) 324 k = bc_h(l)%k(m) 325 sa_p(k+kb,j,i) = sa_p(k,j,i) 274 326 ENDDO 275 327 ENDDO 276 277 328 ! 278 329 !-- Top boundary: Dirichlet or Neumann … … 291 342 ! 292 343 !-- Surface conditions for constant_humidity_flux 344 !-- Run loop over all non-natural and natural walls. Note, in wall-datatype 345 !-- the k coordinate belongs to the atmospheric grid point, therefore, set 346 !-- q_p at k-1 293 347 IF ( ibc_q_b == 0 ) THEN 294 DO i = nxlg, nxrg 295 DO j = nysg, nyng 296 q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i) 348 349 DO l = 0, 1 350 ! 351 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set, 352 !-- for downward-facing surfaces at topography bottom (k+1). 353 kb = MERGE( -1, 1, l == 0 ) 354 !$OMP PARALLEL DO PRIVATE( i, j, k ) 355 DO m = 1, bc_h(l)%ns 356 i = bc_h(l)%i(m) 357 j = bc_h(l)%j(m) 358 k = bc_h(l)%k(m) 359 q_p(k+kb,j,i) = q(k+kb,j,i) 297 360 ENDDO 298 361 ENDDO 362 299 363 ELSE 300 DO i = nxlg, nxrg 301 DO j = nysg, nyng 302 q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i) 364 !$OMP PARALLEL DO PRIVATE( i, j, k ) 365 DO m = 1, bc_h(0)%ns 366 i = bc_h(0)%i(m) 367 j = bc_h(0)%j(m) 368 k = bc_h(0)%k(m) 369 q_p(k-1,j,i) = q_p(k,j,i) 370 ENDDO 371 372 DO l = 0, 1 373 ! 374 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set, 375 !-- for downward-facing surfaces at topography bottom (k+1). 376 kb = MERGE( -1, 1, l == 0 ) 377 !$OMP PARALLEL DO PRIVATE( i, j, k ) 378 DO m = 1, bc_h(l)%ns 379 i = bc_h(l)%i(m) 380 j = bc_h(l)%j(m) 381 k = bc_h(l)%k(m) 382 q_p(k+kb,j,i) = q_p(k,j,i) 303 383 ENDDO 304 384 ENDDO … … 315 395 ! 316 396 !-- Surface conditions rain water (Dirichlet) 317 DO i = nxlg, nxrg 318 DO j = nysg, nyng 319 qr_p(nzb_s_inner(j,i),j,i) = 0.0_wp 320 nr_p(nzb_s_inner(j,i),j,i) = 0.0_wp 321 ENDDO 397 !-- Run loop over all non-natural and natural walls. Note, in wall-datatype 398 !-- the k coordinate belongs to the atmospheric grid point, therefore, set 399 !-- qr_p and nr_p at k-1 400 !$OMP PARALLEL DO PRIVATE( i, j, k ) 401 DO m = 1, bc_h(0)%ns 402 i = bc_h(0)%i(m) 403 j = bc_h(0)%j(m) 404 k = bc_h(0)%k(m) 405 qr_p(k-1,j,i) = 0.0_wp 406 nr_p(k-1,j,i) = 0.0_wp 322 407 ENDDO 323 408 ! … … 334 419 ! 335 420 !-- Surface conditions for constant_humidity_flux 421 !-- Run loop over all non-natural and natural walls. Note, in wall-datatype 422 !-- the k coordinate belongs to the atmospheric grid point, therefore, set 423 !-- s_p at k-1 336 424 IF ( ibc_s_b == 0 ) THEN 337 DO i = nxlg, nxrg 338 DO j = nysg, nyng 339 s_p(nzb_s_inner(j,i),j,i) = s(nzb_s_inner(j,i),j,i) 425 426 DO l = 0, 1 427 ! 428 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set, 429 !-- for downward-facing surfaces at topography bottom (k+1). 430 kb = MERGE( -1, 1, l == 0 ) 431 !$OMP PARALLEL DO PRIVATE( i, j, k ) 432 DO m = 1, bc_h(l)%ns 433 i = bc_h(l)%i(m) 434 j = bc_h(l)%j(m) 435 k = bc_h(l)%k(m) 436 s_p(k+kb,j,i) = s(k+kb,j,i) 340 437 ENDDO 341 438 ENDDO 439 342 440 ELSE 343 DO i = nxlg, nxrg 344 DO j = nysg, nyng 345 s_p(nzb_s_inner(j,i),j,i) = s_p(nzb_s_inner(j,i)+1,j,i) 441 !$OMP PARALLEL DO PRIVATE( i, j, k ) 442 DO m = 1, bc_h(0)%ns 443 i = bc_h(0)%i(m) 444 j = bc_h(0)%j(m) 445 k = bc_h(0)%k(m) 446 s_p(k-1,j,i) = s_p(k,j,i) 447 ENDDO 448 449 DO l = 0, 1 450 ! 451 !-- Set kb, for upward-facing surfaces value at topography top (k-1) is set, 452 !-- for downward-facing surfaces at topography bottom (k+1). 453 kb = MERGE( -1, 1, l == 0 ) 454 !$OMP PARALLEL DO PRIVATE( i, j, k ) 455 DO m = 1, bc_h(l)%ns 456 i = bc_h(l)%i(m) 457 j = bc_h(l)%j(m) 458 k = bc_h(l)%k(m) 459 s_p(k+kb,j,i) = s_p(k,j,i) 346 460 ENDDO 347 461 ENDDO … … 394 508 IF ( outflow_s ) THEN 395 509 pt_p(:,nys-1,:) = pt_p(:,nys,:) 396 IF ( .NOT. constant_diffusion 510 IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:) 397 511 IF ( humidity ) THEN 398 512 q_p(:,nys-1,:) = q_p(:,nys,:)
Note: See TracChangeset
for help on using the changeset viewer.
