Changeset 1942 for palm/trunk/SOURCE/init_grid.f90
- Timestamp:
- Jun 14, 2016 12:18:18 PM (9 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/init_grid.f90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/init_grid.f90
r1932 r1942 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Topography filter implemented to fill holes resolved by only one grid point. 22 ! Initialization of flags for ws-scheme moved to advec_ws. 22 23 ! 23 24 ! Former revisions: … … 171 172 SUBROUTINE init_grid 172 173 174 USE advec_ws, & 175 ONLY: ws_init_flags 173 176 174 177 USE arrays_3d, & … … 185 188 io_group, inflow_l, inflow_n, inflow_r, inflow_s, & 186 189 masking_method, maximum_grid_level, message_string, & 187 momentum_advec, nest_domain, nest_bound_l, nest_bound_n, & 188 nest_bound_r, nest_bound_s, ocean, outflow_l, outflow_n, & 190 momentum_advec, nest_domain, ocean, outflow_l, outflow_n, & 189 191 outflow_r, outflow_s, psolver, scalar_advec, topography, & 190 192 topography_grid_convention, use_surface_fluxes, use_top_fluxes, & … … 197 199 198 200 USE indices, & 199 ONLY: flags, nbgp, nx, nxl, nxlg, nxl u, nxl_mg, nxr, nxrg, nxr_mg,&200 ny, nyn, nyng, nyn_mg, nys, nys v, nys_mg, nysg, nz, nzb,&201 ONLY: flags, nbgp, nx, nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, & 202 ny, nyn, nyng, nyn_mg, nys, nys_mg, nysg, nz, nzb, & 201 203 nzb_diff, nzb_diff_s_inner, nzb_diff_s_outer, nzb_diff_u, & 202 204 nzb_diff_v, nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner, & … … 214 216 IMPLICIT NONE 215 217 216 INTEGER(iwp) :: bh !< temporary vertical index of building height 217 INTEGER(iwp) :: blx !< grid point number of building size along x 218 INTEGER(iwp) :: bly !< grid point number of building size along y 219 INTEGER(iwp) :: bxl !< index for left building wall 220 INTEGER(iwp) :: bxr !< index for right building wall 221 INTEGER(iwp) :: byn !< index for north building wall 222 INTEGER(iwp) :: bys !< index for south building wall 223 INTEGER(iwp) :: ch !< temporary vertical index for canyon height 224 INTEGER(iwp) :: cwx !< grid point number of canyon size along x 225 INTEGER(iwp) :: cwy !< grid point number of canyon size along y 226 INTEGER(iwp) :: cxl !< index for left canyon wall 227 INTEGER(iwp) :: cxr !< index for right canyon wall 228 INTEGER(iwp) :: cyn !< index for north canyon wall 229 INTEGER(iwp) :: cys !< index for south canyon wall 230 INTEGER(iwp) :: gls !< number of lateral ghost points at total model domain boundaries required for multigrid solver 231 INTEGER(iwp) :: i !< index variable along x 232 INTEGER(iwp) :: ii !< loop variable for reading topography file 233 INTEGER(iwp) :: inc !< incremental parameter for coarsening grid level 234 INTEGER(iwp) :: j !< index variable along y 235 INTEGER(iwp) :: k !< index variable along z 236 INTEGER(iwp) :: l !< loop variable 237 INTEGER(iwp) :: nxl_l !< index of left PE boundary for multigrid level 238 INTEGER(iwp) :: nxr_l !< index of right PE boundary for multigrid level 239 INTEGER(iwp) :: nyn_l !< index of north PE boundary for multigrid level 240 INTEGER(iwp) :: nys_l !< index of south PE boundary for multigrid level 241 INTEGER(iwp) :: nzb_si !< dummy index for local nzb_s_inner 242 INTEGER(iwp) :: nzt_l !< index of top PE boundary for multigrid level 243 INTEGER(iwp) :: vi !< dummy for vertical influence 218 INTEGER(iwp) :: bh !< temporary vertical index of building height 219 INTEGER(iwp) :: blx !< grid point number of building size along x 220 INTEGER(iwp) :: bly !< grid point number of building size along y 221 INTEGER(iwp) :: bxl !< index for left building wall 222 INTEGER(iwp) :: bxr !< index for right building wall 223 INTEGER(iwp) :: byn !< index for north building wall 224 INTEGER(iwp) :: bys !< index for south building wall 225 INTEGER(iwp) :: ch !< temporary vertical index for canyon height 226 INTEGER(iwp) :: cwx !< grid point number of canyon size along x 227 INTEGER(iwp) :: cwy !< grid point number of canyon size along y 228 INTEGER(iwp) :: cxl !< index for left canyon wall 229 INTEGER(iwp) :: cxr !< index for right canyon wall 230 INTEGER(iwp) :: cyn !< index for north canyon wall 231 INTEGER(iwp) :: cys !< index for south canyon wall 232 INTEGER(iwp) :: gls !< number of lateral ghost points at total model domain boundaries required for multigrid solver 233 INTEGER(iwp) :: i !< index variable along x 234 INTEGER(iwp) :: ii !< loop variable for reading topography file 235 INTEGER(iwp) :: inc !< incremental parameter for coarsening grid level 236 INTEGER(iwp) :: j !< index variable along y 237 INTEGER(iwp) :: k !< index variable along z 238 INTEGER(iwp) :: l !< loop variable 239 INTEGER(iwp) :: nxl_l !< index of left PE boundary for multigrid level 240 INTEGER(iwp) :: nxr_l !< index of right PE boundary for multigrid level 241 INTEGER(iwp) :: nyn_l !< index of north PE boundary for multigrid level 242 INTEGER(iwp) :: nys_l !< index of south PE boundary for multigrid level 243 INTEGER(iwp) :: nzb_si !< dummy index for local nzb_s_inner 244 INTEGER(iwp) :: nzt_l !< index of top PE boundary for multigrid level 245 INTEGER(iwp) :: num_hole !< number of holes (in topography) resolved by only one grid point 246 INTEGER(iwp) :: num_wall !< number of surrounding vertical walls for a single grid point 247 INTEGER(iwp) :: vi !< dummy for vertical influence 244 248 245 249 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: & … … 259 263 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_tmp !< dummy to calculate topography indices on u- and v-grid 260 264 261 LOGICAL :: flag_set = .FALSE. !< steering variable for advection flags265 LOGICAL :: hole = .FALSE. !< flag to check if any holes resolved by only 1 grid point were filled 262 266 263 267 REAL(wp) :: dx_l !< grid spacing along x on different multigrid level … … 744 748 745 749 DEALLOCATE ( topo_height ) 746 750 ! 751 !-- Filter topography, i.e. fill holes resolved by only one grid point. 752 !-- Such holes are suspected to lead to velocity blow-ups as continuity 753 !-- equation on discrete grid cannot be fulfilled in such case. 754 !-- For now, check only for holes and fill them to the lowest height level 755 !-- of the directly adjoining grid points along x- and y- direction. 756 !-- Before checking for holes, set lateral boundary conditions for 757 !-- topography. After hole-filling, boundary conditions must be set again! 758 IF ( bc_ns_cyc ) THEN 759 nzb_local(-1,:) = nzb_local(ny,:) 760 nzb_local(ny+1,:) = nzb_local(0,:) 761 ELSE 762 nzb_local(-1,:) = nzb_local(0,:) 763 nzb_local(ny+1,:) = nzb_local(ny,:) 764 ENDIF 765 766 IF ( bc_lr_cyc ) THEN 767 nzb_local(:,-1) = nzb_local(:,nx) 768 nzb_local(:,nx+1) = nzb_local(:,0) 769 ELSE 770 nzb_local(:,-1) = nzb_local(:,0) 771 nzb_local(:,nx+1) = nzb_local(:,nx) 772 ENDIF 773 774 num_hole = 0 775 DO i = 0, nx 776 DO j = 0, ny 777 778 num_wall = 0 779 780 IF ( nzb_local(j-1,i) > nzb_local(j,i) ) & 781 num_wall = num_wall + 1 782 IF ( nzb_local(j+1,i) > nzb_local(j,i) ) & 783 num_wall = num_wall + 1 784 IF ( nzb_local(j,i-1) > nzb_local(j,i) ) & 785 num_wall = num_wall + 1 786 IF ( nzb_local(j,i+1) > nzb_local(j,i) ) & 787 num_wall = num_wall + 1 788 789 IF ( num_wall == 4 ) THEN 790 hole = .TRUE. 791 nzb_local(j,i) = MIN( nzb_local(j-1,i), nzb_local(j+1,i), & 792 nzb_local(j,i-1), nzb_local(j,i+1) ) 793 num_hole = num_hole + 1 794 ENDIF 795 ENDDO 796 ENDDO 797 ! 798 !-- Create an informative message if any hole was removed. 799 IF ( hole ) THEN 800 WRITE( message_string, * ) num_hole, 'hole(s) resolved by only '//& 801 'one grid point were filled' 802 CALL message( 'init_grid', 'PA0430', 0, 0, 0, 6, 0 ) 803 ENDIF 747 804 ! 748 805 !-- Add cyclic or Neumann boundary conditions (additional layers are for … … 1237 1294 ENDIF 1238 1295 1239 !1240 !-- Test output of flag arrays1241 ! i = nxl_l1242 ! WRITE (9,*) ' '1243 ! WRITE (9,*) '*** mg level ', l, ' ***', mg_switch_to_pe0_level1244 ! WRITE (9,*) ' inc=', inc, ' i =', nxl_l1245 ! WRITE (9,*) ' nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l1246 ! DO k = nzt_l+1, nzb, -11247 ! WRITE (9,'(194(1X,I2))') ( flags(k,j,i), j = nys_l-1, nyn_l+1 )1248 ! ENDDO1249 1250 1296 inc = inc * 2 1251 1297 … … 1254 1300 ENDIF 1255 1301 ! 1256 !-- Allocate flags needed for masking walls. 1302 !-- Allocate flags needed for masking walls. Even though these flags are only 1303 !-- required in the ws-scheme, the arrays need to be allocated as they are 1304 !-- used in OpenACC directives. 1257 1305 ALLOCATE( wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 1258 1306 wall_flags_00(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 1259 1307 wall_flags_0 = 0 1260 1308 wall_flags_00 = 0 1261 1262 IF ( scalar_advec == 'ws-scheme' .OR. & 1263 scalar_advec == 'ws-scheme-mono' ) THEN 1264 ! 1265 !-- Set flags to steer the degradation of the advection scheme in advec_ws 1266 !-- near topography, inflow- and outflow boundaries as well as bottom and 1267 !-- top of model domain. wall_flags_0 remains zero for all non-prognostic 1268 !-- grid points. 1269 DO i = nxl, nxr 1270 DO j = nys, nyn 1271 DO k = nzb_s_inner(j,i)+1, nzt 1272 ! 1273 !-- scalar - x-direction 1274 !-- WS1 (0), WS3 (1), WS5 (2) 1275 IF ( ( k <= nzb_s_inner(j,i+1) .OR. k <= nzb_s_inner(j,i+2) & 1276 .OR. k <= nzb_s_inner(j,i-1) ) & 1277 .OR. ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1278 .AND. i == nxl ) .OR. & 1279 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1280 .AND. i == nxr ) ) & 1281 THEN 1282 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 ) 1283 ELSEIF ( ( k <= nzb_s_inner(j,i+3) .AND. k > nzb_s_inner(j,i+1)& 1284 .AND. k > nzb_s_inner(j,i+2)& 1285 .AND. k > nzb_s_inner(j,i-1)& 1286 ) .OR. & 1287 ( k <= nzb_s_inner(j,i-2) .AND. k > nzb_s_inner(j,i+1)& 1288 .AND. k > nzb_s_inner(j,i+2)& 1289 .AND. k > nzb_s_inner(j,i-1)& 1290 ) & 1291 .OR. & 1292 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1293 .AND. i == nxr-1 ) .OR. & 1294 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1295 .AND. i == nxlu ) ) & 1296 THEN 1297 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 ) 1298 ELSEIF ( k > nzb_s_inner(j,i+1) .AND. k > nzb_s_inner(j,i+2) & 1299 .AND. k > nzb_s_inner(j,i+3) .AND. k > nzb_s_inner(j,i-1) & 1300 .AND. k > nzb_s_inner(j,i-2) ) & 1301 THEN 1302 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 ) 1303 ENDIF 1304 ! 1305 !-- scalar - y-direction 1306 !-- WS1 (3), WS3 (4), WS5 (5) 1307 IF ( ( k <= nzb_s_inner(j+1,i) .OR. k <= nzb_s_inner(j+2,i) & 1308 .OR. k <= nzb_s_inner(j-1,i) ) & 1309 .OR. & 1310 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1311 .AND. j == nys ) .OR. & 1312 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1313 .AND. j == nyn ) ) & 1314 THEN 1315 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 ) 1316 ! 1317 !-- WS3 1318 ELSEIF ( ( k <= nzb_s_inner(j+3,i) .AND. k > nzb_s_inner(j+1,i)& 1319 .AND. k > nzb_s_inner(j+2,i)& 1320 .AND. k > nzb_s_inner(j-1,i)& 1321 ) .OR. & 1322 ( k <= nzb_s_inner(j-2,i) .AND. k > nzb_s_inner(j+1,i)& 1323 .AND. k > nzb_s_inner(j+2,i)& 1324 .AND. k > nzb_s_inner(j-1,i)& 1325 ) & 1326 .OR. & 1327 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1328 .AND. j == nysv ) .OR. & 1329 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1330 .AND. j == nyn-1 ) ) & 1331 THEN 1332 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 4 ) 1333 ! 1334 !-- WS5 1335 ELSEIF ( k > nzb_s_inner(j+1,i) .AND. k > nzb_s_inner(j+2,i) & 1336 .AND. k > nzb_s_inner(j+3,i) .AND. k > nzb_s_inner(j-1,i) & 1337 .AND. k > nzb_s_inner(j-2,i) ) & 1338 THEN 1339 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 ) 1340 ENDIF 1341 ! 1342 !-- scalar - z-direction 1343 !-- WS1 (6), WS3 (7), WS5 (8) 1344 flag_set = .FALSE. 1345 IF ( k == nzb_s_inner(j,i) + 1 .OR. k == nzt ) THEN 1346 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 6 ) 1347 flag_set = .TRUE. 1348 ELSEIF ( k == nzb_s_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN 1349 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 7 ) 1350 flag_set = .TRUE. 1351 ELSEIF ( k > nzb_s_inner(j,i) .AND. .NOT. flag_set ) THEN 1352 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 8 ) 1353 ENDIF 1354 1355 ENDDO 1356 ENDDO 1357 ENDDO 1358 ENDIF 1359 1360 IF ( momentum_advec == 'ws-scheme' ) THEN 1361 ! 1362 !-- Set wall_flags_0 to steer the degradation of the advection scheme in advec_ws 1363 !-- near topography, inflow- and outflow boundaries as well as bottom and 1364 !-- top of model domain. wall_flags_0 remains zero for all non-prognostic 1365 !-- grid points. 1366 DO i = nxl, nxr 1367 DO j = nys, nyn 1368 DO k = nzb+1, nzt 1369 ! 1370 !-- At first, set flags to WS1. 1371 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 1372 !-- in order to handle the left/south flux. 1373 !-- near vertical walls. 1374 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 ) 1375 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 ) 1376 ! 1377 !-- u component - x-direction 1378 !-- WS1 (9), WS3 (10), WS5 (11) 1379 IF ( k <= nzb_u_inner(j,i+1) .OR. & 1380 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1381 .AND. i <= nxlu ) .OR. & 1382 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1383 .AND. i == nxr ) ) & 1384 THEN 1385 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 ) 1386 ELSEIF ( ( k <= nzb_u_inner(j,i+2) .AND. & 1387 k > nzb_u_inner(j,i+1) ) .OR. & 1388 k <= nzb_u_inner(j,i-1) & 1389 .OR. & 1390 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1391 .AND. i == nxr-1 ) .OR. & 1392 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1393 .AND. i == nxlu+1) ) & 1394 THEN 1395 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 10 ) 1396 ! 1397 !-- Clear flag for WS1 1398 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 9 ) 1399 ELSEIF ( k > nzb_u_inner(j,i+1) .AND. k > nzb_u_inner(j,i+2) & 1400 .AND. k > nzb_u_inner(j,i-1) ) & 1401 THEN 1402 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 11 ) 1403 ! 1404 !-- Clear flag for WS1 1405 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 9 ) 1406 ENDIF 1407 ! 1408 !-- u component - y-direction 1409 !-- WS1 (12), WS3 (13), WS5 (14) 1410 IF ( k <= nzb_u_inner(j+1,i) .OR. & 1411 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1412 .AND. j == nys ) .OR. & 1413 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1414 .AND. j == nyn ) ) & 1415 THEN 1416 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 ) 1417 ELSEIF ( ( k <= nzb_u_inner(j+2,i) .AND. & 1418 k > nzb_u_inner(j+1,i) ) .OR. & 1419 k <= nzb_u_inner(j-1,i) & 1420 .OR. & 1421 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1422 .AND. j == nysv ) .OR. & 1423 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1424 .AND. j == nyn-1 ) ) & 1425 THEN 1426 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 ) 1427 ! 1428 !-- Clear flag for WS1 1429 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 12 ) 1430 ELSEIF ( k > nzb_u_inner(j+1,i) .AND. k > nzb_u_inner(j+2,i) & 1431 .AND. k > nzb_u_inner(j-1,i) ) & 1432 THEN 1433 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 14 ) 1434 ! 1435 !-- Clear flag for WS1 1436 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 12 ) 1437 ENDIF 1438 ! 1439 !-- u component - z-direction 1440 !-- WS1 (15), WS3 (16), WS5 (17) 1441 flag_set = .FALSE. 1442 IF ( k == nzb_u_inner(j,i) + 1 .OR. k == nzt ) THEN 1443 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 15 ) 1444 flag_set = .TRUE. 1445 ELSEIF ( k == nzb_u_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN 1446 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 16 ) 1447 flag_set = .TRUE. 1448 ELSEIF ( k > nzb_u_inner(j,i) .AND. .NOT. flag_set ) THEN 1449 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 17 ) 1450 ENDIF 1451 1452 ENDDO 1453 ENDDO 1454 ENDDO 1455 1456 DO i = nxl, nxr 1457 DO j = nys, nyn 1458 DO k = nzb+1, nzt 1459 ! 1460 !-- At first, set flags to WS1. 1461 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 1462 !-- in order to handle the left/south flux. 1463 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 ) 1464 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 ) 1465 ! 1466 !-- v component - x-direction 1467 !-- WS1 (18), WS3 (19), WS5 (20) 1468 IF ( k <= nzb_v_inner(j,i+1) .OR. & 1469 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1470 .AND. i == nxl ) .OR. & 1471 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1472 .AND. i == nxr ) ) & 1473 THEN 1474 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 ) 1475 ! 1476 !-- WS3 1477 ELSEIF ( ( k <= nzb_v_inner(j,i+2) .AND. & 1478 k > nzb_v_inner(j,i+1) ) .OR. & 1479 k <= nzb_v_inner(j,i-1) & 1480 .OR. & 1481 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1482 .AND. i == nxr-1 ) .OR. & 1483 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1484 .AND. i == nxlu ) ) & 1485 THEN 1486 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 ) 1487 ! 1488 !-- Clear flag for WS1 1489 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 18 ) 1490 ELSEIF ( k > nzb_v_inner(j,i+1) .AND. k > nzb_v_inner(j,i+2) & 1491 .AND. k > nzb_v_inner(j,i-1) ) & 1492 THEN 1493 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 20 ) 1494 ! 1495 !-- Clear flag for WS1 1496 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 18 ) 1497 ENDIF 1498 ! 1499 !-- v component - y-direction 1500 !-- WS1 (21), WS3 (22), WS5 (23) 1501 IF ( k <= nzb_v_inner(j+1,i) .OR. & 1502 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1503 .AND. j <= nysv ) .OR. & 1504 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1505 .AND. j == nyn ) ) & 1506 THEN 1507 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 ) 1508 ELSEIF ( ( k <= nzb_v_inner(j+2,i) .AND. & 1509 k > nzb_v_inner(j+1,i) ) .OR. & 1510 k <= nzb_v_inner(j-1,i) & 1511 .OR. & 1512 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1513 .AND. j == nysv+1) .OR. & 1514 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1515 .AND. j == nyn-1 ) ) & 1516 THEN 1517 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 ) 1518 ! 1519 !-- Clear flag for WS1 1520 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 21 ) 1521 ELSEIF ( k > nzb_v_inner(j+1,i) .AND. k > nzb_v_inner(j+2,i) & 1522 .AND. k > nzb_v_inner(j-1,i) ) & 1523 THEN 1524 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 ) 1525 ! 1526 !-- Clear flag for WS1 1527 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 21 ) 1528 ENDIF 1529 ! 1530 !-- v component - z-direction 1531 !-- WS1 (24), WS3 (25), WS5 (26) 1532 flag_set = .FALSE. 1533 IF ( k == nzb_v_inner(j,i) + 1 .OR. k == nzt ) THEN 1534 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 24 ) 1535 flag_set = .TRUE. 1536 ELSEIF ( k == nzb_v_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN 1537 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 25 ) 1538 flag_set = .TRUE. 1539 ELSEIF ( k > nzb_v_inner(j,i) .AND. .NOT. flag_set ) THEN 1540 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 26 ) 1541 ENDIF 1542 1543 ENDDO 1544 ENDDO 1545 ENDDO 1546 DO i = nxl, nxr 1547 DO j = nys, nyn 1548 DO k = nzb+1, nzt 1549 ! 1550 !-- At first, set flags to WS1. 1551 !-- Since fluxes are swapped in advec_ws.f90, this is necessary to 1552 !-- in order to handle the left/south flux. 1553 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 ) 1554 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 ) 1555 ! 1556 !-- w component - x-direction 1557 !-- WS1 (27), WS3 (28), WS5 (29) 1558 IF ( k <= nzb_w_inner(j,i+1) .OR. & 1559 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1560 .AND. i == nxl ) .OR. & 1561 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1562 .AND. i == nxr ) ) & 1563 THEN 1564 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 ) 1565 ELSEIF ( ( k <= nzb_w_inner(j,i+2) .AND. & 1566 k > nzb_w_inner(j,i+1) ) .OR. & 1567 k <= nzb_w_inner(j,i-1) & 1568 .OR. & 1569 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1570 .AND. i == nxr-1 ) .OR. & 1571 ( ( inflow_l .OR. outflow_l .OR. nest_bound_l ) & 1572 .AND. i == nxlu ) ) & 1573 THEN 1574 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 ) 1575 ! 1576 !-- Clear flag for WS1 1577 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 27 ) 1578 ELSEIF ( k > nzb_w_inner(j,i+1) .AND. k > nzb_w_inner(j,i+2) & 1579 .AND. k > nzb_w_inner(j,i-1) ) & 1580 THEN 1581 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i),29 ) 1582 ! 1583 !-- Clear flag for WS1 1584 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 27 ) 1585 ENDIF 1586 ! 1587 !-- w component - y-direction 1588 !-- WS1 (30), WS3 (31), WS5 (32) 1589 IF ( k <= nzb_w_inner(j+1,i) .OR. & 1590 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1591 .AND. j == nys ) .OR. & 1592 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1593 .AND. j == nyn ) ) & 1594 THEN 1595 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 ) 1596 ELSEIF ( ( k <= nzb_w_inner(j+2,i) .AND. & 1597 k > nzb_w_inner(j+1,i) ) .OR. & 1598 k <= nzb_w_inner(j-1,i) & 1599 .OR. & 1600 ( ( inflow_s .OR. outflow_s .OR. nest_bound_s ) & 1601 .AND. j == nysv ) .OR. & 1602 ( ( inflow_n .OR. outflow_n .OR. nest_bound_n ) & 1603 .AND. j == nyn-1 ) ) & 1604 THEN 1605 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 ) 1606 ! 1607 !-- Clear flag for WS1 1608 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 30 ) 1609 ELSEIF ( k > nzb_w_inner(j+1,i) .AND. k > nzb_w_inner(j+2,i) & 1610 .AND. k > nzb_w_inner(j-1,i) ) & 1611 THEN 1612 wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 0 ) 1613 ! 1614 !-- Clear flag for WS1 1615 wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 30 ) 1616 ENDIF 1617 ! 1618 !-- w component - z-direction 1619 !-- WS1 (33), WS3 (34), WS5 (35) 1620 flag_set = .FALSE. 1621 IF ( k == nzb_w_inner(j,i) .OR. k == nzb_w_inner(j,i) + 1 & 1622 .OR. k == nzt ) THEN 1623 ! 1624 !-- Please note, at k == nzb_w_inner(j,i) a flag is explictely 1625 !-- set, although this is not a prognostic level. However, 1626 !-- contrary to the advection of u,v and s this is necessary 1627 !-- because flux_t(nzb_w_inner(j,i)) is used for the tendency 1628 !-- at k == nzb_w_inner(j,i)+1. 1629 wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 1 ) 1630 flag_set = .TRUE. 1631 ELSEIF ( k == nzb_w_inner(j,i) + 2 .OR. k == nzt - 1 ) THEN 1632 wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 2 ) 1633 flag_set = .TRUE. 1634 ELSEIF ( k > nzb_w_inner(j,i) .AND. .NOT. flag_set ) THEN 1635 wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 3 ) 1636 ENDIF 1637 1638 ENDDO 1639 ENDDO 1640 ENDDO 1641 1642 ENDIF 1643 1644 ! 1645 !-- Exchange 3D integer wall_flags. 1646 IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' & 1647 .OR. scalar_advec == 'ws-scheme-mono' ) THEN 1648 ! 1649 !-- Exchange ghost points for advection flags 1650 CALL exchange_horiz_int( wall_flags_0, nbgp ) 1651 CALL exchange_horiz_int( wall_flags_00, nbgp ) 1652 ! 1653 !-- Set boundary flags at inflow and outflow boundary in case of 1654 !-- non-cyclic boundary conditions. 1655 IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) THEN 1656 wall_flags_0(:,:,nxl-1) = wall_flags_0(:,:,nxl) 1657 wall_flags_00(:,:,nxl-1) = wall_flags_00(:,:,nxl) 1658 ENDIF 1659 1660 IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) THEN 1661 wall_flags_0(:,:,nxr+1) = wall_flags_0(:,:,nxr) 1662 wall_flags_00(:,:,nxr+1) = wall_flags_00(:,:,nxr) 1663 ENDIF 1664 1665 IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) THEN 1666 wall_flags_0(:,nyn+1,:) = wall_flags_0(:,nyn,:) 1667 wall_flags_00(:,nyn+1,:) = wall_flags_00(:,nyn,:) 1668 ENDIF 1669 1670 IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) THEN 1671 wall_flags_0(:,nys-1,:) = wall_flags_0(:,nys,:) 1672 wall_flags_00(:,nys-1,:) = wall_flags_00(:,nys,:) 1673 ENDIF 1674 1309 ! 1310 !-- Init flags for ws-scheme to degrade order near walls 1311 IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' .OR.& 1312 scalar_advec == 'ws-scheme-mono' ) THEN 1313 CALL ws_init_flags 1675 1314 ENDIF 1676 1315
Note: See TracChangeset
for help on using the changeset viewer.
