Changeset 1942 for palm/trunk/SOURCE/init_grid.f90
 Timestamp:
 Jun 14, 2016 12:18:18 PM (6 years ago)
 File:

 1 edited
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 wsscheme 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 vgrid 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 blowups 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 holefilling, 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(j1,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,i1) > 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(j1,i), nzb_local(j+1,i), & 792 nzb_local(j,i1), 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_l1, 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 wsscheme, 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 == 'wsscheme' .OR. & 1263 scalar_advec == 'wsschememono' ) 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 nonprognostic 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  xdirection 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,i1) ) & 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,i1)& 1286 ) .OR. & 1287 ( k <= nzb_s_inner(j,i2) .AND. k > nzb_s_inner(j,i+1)& 1288 .AND. k > nzb_s_inner(j,i+2)& 1289 .AND. k > nzb_s_inner(j,i1)& 1290 ) & 1291 .OR. & 1292 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1293 .AND. i == nxr1 ) .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,i1) & 1300 .AND. k > nzb_s_inner(j,i2) ) & 1301 THEN 1302 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 ) 1303 ENDIF 1304 ! 1305 ! scalar  ydirection 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(j1,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(j1,i)& 1321 ) .OR. & 1322 ( k <= nzb_s_inner(j2,i) .AND. k > nzb_s_inner(j+1,i)& 1323 .AND. k > nzb_s_inner(j+2,i)& 1324 .AND. k > nzb_s_inner(j1,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 == nyn1 ) ) & 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(j1,i) & 1337 .AND. k > nzb_s_inner(j2,i) ) & 1338 THEN 1339 wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 ) 1340 ENDIF 1341 ! 1342 ! scalar  zdirection 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 == 'wsscheme' ) 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 nonprognostic 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  xdirection 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,i1) & 1389 .OR. & 1390 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1391 .AND. i == nxr1 ) .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,i1) ) & 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  ydirection 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(j1,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 == nyn1 ) ) & 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(j1,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  zdirection 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  xdirection 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,i1) & 1480 .OR. & 1481 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1482 .AND. i == nxr1 ) .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,i1) ) & 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  ydirection 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(j1,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 == nyn1 ) ) & 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(j1,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  zdirection 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  xdirection 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,i1) & 1568 .OR. & 1569 ( ( inflow_r .OR. outflow_r .OR. nest_bound_r ) & 1570 .AND. i == nxr1 ) .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,i1) ) & 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  ydirection 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(j1,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 == nyn1 ) ) & 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(j1,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  zdirection 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 == 'wsscheme' .OR. scalar_advec == 'wsscheme' & 1647 .OR. scalar_advec == 'wsschememono' ) 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 ! noncyclic boundary conditions. 1655 IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) THEN 1656 wall_flags_0(:,:,nxl1) = wall_flags_0(:,:,nxl) 1657 wall_flags_00(:,:,nxl1) = 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(:,nys1,:) = wall_flags_0(:,nys,:) 1672 wall_flags_00(:,nys1,:) = wall_flags_00(:,nys,:) 1673 ENDIF 1674 1309 ! 1310 ! Init flags for wsscheme to degrade order near walls 1311 IF ( momentum_advec == 'wsscheme' .OR. scalar_advec == 'wsscheme' .OR.& 1312 scalar_advec == 'wsschememono' ) THEN 1313 CALL ws_init_flags 1675 1314 ENDIF 1676 1315
Note: See TracChangeset
for help on using the changeset viewer.