Ignore:
Timestamp:
Jun 14, 2016 12:18:18 PM (5 years ago)
Author:
suehring
Message:

filter topography by filling holes resolved by only one grid point ;initialization of wall_flags_0 and wall_flags_00 moved to advec_ws

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/init_grid.f90

    r1932 r1942  
    1919! Current revisions:
    2020! -----------------
    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. 
    2223!
    2324! Former revisions:
     
    171172 SUBROUTINE init_grid
    172173 
     174    USE advec_ws,                                                              &
     175        ONLY:  ws_init_flags
    173176
    174177    USE arrays_3d,                                                             &
     
    185188               io_group, inflow_l, inflow_n, inflow_r, inflow_s,               &
    186189               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,       &
    189191               outflow_r, outflow_s, psolver, scalar_advec, topography,        &
    190192               topography_grid_convention, use_surface_fluxes, use_top_fluxes, &
     
    197199       
    198200    USE indices,                                                               &
    199         ONLY:  flags, nbgp, nx, nxl, nxlg, nxlu, nxl_mg, nxr, nxrg, nxr_mg,    &
    200                ny, nyn, nyng, nyn_mg, nys, nysv, 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,              &
    201203               nzb_diff, nzb_diff_s_inner, nzb_diff_s_outer, nzb_diff_u,       &
    202204               nzb_diff_v, nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner,     &
     
    214216    IMPLICIT NONE
    215217
    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
    244248
    245249    INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::                               &
     
    259263    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_tmp    !< dummy to calculate topography indices on u- and v-grid
    260264
    261     LOGICAL  :: flag_set = .FALSE.  !< steering variable for advection flags
     265    LOGICAL  ::  hole = .FALSE.  !< flag to check if any holes resolved by only 1 grid point were filled
    262266
    263267    REAL(wp) ::  dx_l          !< grid spacing along x on different multigrid level
     
    744748
    745749          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
    747804!
    748805!--       Add cyclic or Neumann boundary conditions (additional layers are for
     
    12371294          ENDIF
    12381295
    1239 !
    1240 !--       Test output of flag arrays
    1241 !          i = nxl_l
    1242 !          WRITE (9,*)  ' '
    1243 !          WRITE (9,*)  '*** mg level ', l, ' ***', mg_switch_to_pe0_level
    1244 !          WRITE (9,*)  '    inc=', inc, '  i =', nxl_l
    1245 !          WRITE (9,*)  '    nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l
    1246 !          DO  k = nzt_l+1, nzb, -1
    1247 !             WRITE (9,'(194(1X,I2))')  ( flags(k,j,i), j = nys_l-1, nyn_l+1 )
    1248 !          ENDDO
    1249 
    12501296          inc = inc * 2
    12511297
     
    12541300    ENDIF
    12551301!
    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.
    12571305    ALLOCATE( wall_flags_0(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                     &
    12581306              wall_flags_00(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    12591307    wall_flags_0  = 0
    12601308    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
    16751314    ENDIF
    16761315
Note: See TracChangeset for help on using the changeset viewer.