Changeset 1942


Ignore:
Timestamp:
Jun 14, 2016 12:18:18 PM (8 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

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

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

    r1877 r1942  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Initialization of flags for ws-scheme moved from init_grid.
    2222!
    2323! Former revisions:
     
    167167    PUBLIC   advec_s_ws, advec_s_ws_acc, advec_u_ws, advec_u_ws_acc,          &
    168168             advec_v_ws, advec_v_ws_acc, advec_w_ws, advec_w_ws_acc,          &
    169              ws_init, ws_statistics
     169             ws_init, ws_init_flags, ws_statistics
    170170
    171171    INTERFACE ws_init
    172172       MODULE PROCEDURE ws_init
    173173    END INTERFACE ws_init
     174
     175    INTERFACE ws_init_flags
     176       MODULE PROCEDURE ws_init_flags
     177    END INTERFACE ws_init_flags
    174178
    175179    INTERFACE ws_statistics
     
    364368
    365369    END SUBROUTINE ws_init
     370
     371!------------------------------------------------------------------------------!
     372! Description:
     373! ------------
     374!> Initialization of flags for WS-scheme used to degrade the order of the scheme
     375!> near walls.
     376!------------------------------------------------------------------------------!
     377    SUBROUTINE ws_init_flags
     378
     379       USE control_parameters,                                                 &
     380           ONLY:  inflow_l, inflow_n, inflow_r, inflow_s, momentum_advec,      &
     381                  nest_bound_l, nest_bound_n, nest_bound_r, nest_bound_s,      &
     382                  outflow_l, outflow_n, outflow_r, outflow_s, scalar_advec
     383
     384       USE indices,                                                            &
     385           ONLY:  nbgp, nxl, nxlu, nxr, nyn, nys, nysv, nzb, nzb_s_inner,      &
     386                  nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt, wall_flags_0,    &
     387                  wall_flags_00
     388
     389       USE kinds
     390
     391       IMPLICIT NONE
     392
     393       INTEGER(iwp) ::  i  !< index variable along x
     394       INTEGER(iwp) ::  j  !< index variable along y
     395       INTEGER(iwp) ::  k  !< index variable along z
     396
     397       LOGICAL      ::  flag_set !< steering variable for advection flags
     398   
     399
     400       IF ( scalar_advec == 'ws-scheme' .OR.                                   &
     401            scalar_advec == 'ws-scheme-mono' )  THEN
     402!
     403!--       Set flags to steer the degradation of the advection scheme in advec_ws
     404!--       near topography, inflow- and outflow boundaries as well as bottom and
     405!--       top of model domain. wall_flags_0 remains zero for all non-prognostic
     406!--       grid points.
     407          DO  i = nxl, nxr
     408             DO  j = nys, nyn
     409                DO  k = nzb_s_inner(j,i)+1, nzt
     410!
     411!--                scalar - x-direction
     412!--                WS1 (0), WS3 (1), WS5 (2)
     413                   IF ( ( k <= nzb_s_inner(j,i+1) .OR. k <= nzb_s_inner(j,i+2) &   
     414                     .OR. k <= nzb_s_inner(j,i-1) )                            &
     415                       .OR. ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )    &
     416                              .AND. i == nxl   )    .OR.                       &
     417                            ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )    &
     418                              .AND. i == nxr   ) )                             &
     419                   THEN
     420                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 0 )
     421                   ELSEIF ( ( k <= nzb_s_inner(j,i+3) .AND. k > nzb_s_inner(j,i+1)&
     422                                                      .AND. k > nzb_s_inner(j,i+2)&
     423                                                      .AND. k > nzb_s_inner(j,i-1)&
     424                            )                       .OR.                          &
     425                            ( k <= nzb_s_inner(j,i-2) .AND. k > nzb_s_inner(j,i+1)&
     426                                                      .AND. k > nzb_s_inner(j,i+2)&
     427                                                      .AND. k > nzb_s_inner(j,i-1)&
     428                            )                                                     &
     429                                                    .OR.                          &
     430                            ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )       &
     431                              .AND. i == nxr-1 )    .OR.                          &
     432                            ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )       &
     433                              .AND. i == nxlu  ) )                                &
     434                   THEN
     435                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 1 )
     436                   ELSEIF ( k > nzb_s_inner(j,i+1) .AND. k > nzb_s_inner(j,i+2)&
     437                      .AND. k > nzb_s_inner(j,i+3) .AND. k > nzb_s_inner(j,i-1)&
     438                      .AND. k > nzb_s_inner(j,i-2) )                           &
     439                   THEN
     440                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 2 )
     441                   ENDIF
     442!
     443!--                scalar - y-direction
     444!--                WS1 (3), WS3 (4), WS5 (5)
     445                   IF ( ( k <= nzb_s_inner(j+1,i)  .OR. k <= nzb_s_inner(j+2,i)   &   
     446                                                   .OR. k <= nzb_s_inner(j-1,i) ) &
     447                                                    .OR.                          &
     448                            ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )       &
     449                              .AND. j == nys   )    .OR.                          &
     450                            ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )       &
     451                              .AND. j == nyn   ) )                                &
     452                   THEN
     453                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 3 )
     454!
     455!--                WS3
     456                   ELSEIF ( ( k <= nzb_s_inner(j+3,i) .AND. k > nzb_s_inner(j+1,i)&
     457                                                      .AND. k > nzb_s_inner(j+2,i)&
     458                                                      .AND. k > nzb_s_inner(j-1,i)&
     459                            )                           .OR.                      &
     460                            ( k <= nzb_s_inner(j-2,i) .AND. k > nzb_s_inner(j+1,i)&
     461                                                      .AND. k > nzb_s_inner(j+2,i)&
     462                                                      .AND. k > nzb_s_inner(j-1,i)&
     463                            )                                                     &
     464                                                        .OR.                      &
     465                            ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )       &
     466                              .AND. j == nysv  )    .OR.                          &
     467                            ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )       &
     468                              .AND. j == nyn-1 ) )                                &
     469                   THEN
     470                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 4 )
     471!
     472!--                WS5
     473                   ELSEIF ( k > nzb_s_inner(j+1,i) .AND. k > nzb_s_inner(j+2,i)&
     474                      .AND. k > nzb_s_inner(j+3,i) .AND. k > nzb_s_inner(j-1,i)&
     475                      .AND. k > nzb_s_inner(j-2,i) )                           &
     476                   THEN
     477                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 )
     478                   ENDIF
     479!
     480!--                scalar - z-direction
     481!--                WS1 (6), WS3 (7), WS5 (8)
     482                   flag_set = .FALSE.
     483                   IF ( k == nzb_s_inner(j,i) + 1 .OR. k == nzt )  THEN
     484                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 6 )
     485                      flag_set = .TRUE.
     486                   ELSEIF ( k == nzb_s_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
     487                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 7 )
     488                      flag_set = .TRUE.
     489                   ELSEIF ( k > nzb_s_inner(j,i) .AND. .NOT. flag_set )  THEN
     490                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 8 )
     491                   ENDIF
     492
     493                ENDDO
     494             ENDDO
     495          ENDDO
     496       ENDIF
     497
     498       IF ( momentum_advec == 'ws-scheme' )  THEN
     499!
     500!--       Set wall_flags_0 to steer the degradation of the advection scheme in advec_ws
     501!--       near topography, inflow- and outflow boundaries as well as bottom and
     502!--       top of model domain. wall_flags_0 remains zero for all non-prognostic
     503!--       grid points.
     504          DO  i = nxl, nxr
     505             DO  j = nys, nyn
     506                DO  k = nzb+1, nzt
     507
     508!--                At first, set flags to WS1.
     509!--                Since fluxes are swapped in advec_ws.f90, this is necessary to
     510!--                in order to handle the left/south flux.
     511!--                near vertical walls.
     512                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 )
     513                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 )
     514!
     515!--                u component - x-direction
     516!--                WS1 (9), WS3 (10), WS5 (11)
     517                   IF ( k <= nzb_u_inner(j,i+1)     .OR.                       &
     518                            ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )    &
     519                              .AND. i <= nxlu  )    .OR.                       &
     520                            ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )    &
     521                              .AND. i == nxr   ) )                             &
     522                   THEN
     523                       wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 9 )
     524                   ELSEIF ( ( k <= nzb_u_inner(j,i+2) .AND.                    &
     525                              k >  nzb_u_inner(j,i+1) ) .OR.                   &
     526                              k <= nzb_u_inner(j,i-1)                          &
     527                                                        .OR.                   &
     528                            ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )    &
     529                              .AND. i == nxr-1 )    .OR.                       &
     530                            ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )    &
     531                              .AND. i == nxlu+1) )                             &
     532                   THEN
     533                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 10 )
     534!
     535!--                   Clear flag for WS1
     536                      wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 9 )
     537                   ELSEIF ( k > nzb_u_inner(j,i+1) .AND. k > nzb_u_inner(j,i+2)   &
     538                                                   .AND. k > nzb_u_inner(j,i-1) ) &
     539                   THEN   
     540                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 11 )
     541!
     542!--                   Clear flag for WS1
     543                      wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 9 )
     544                   ENDIF
     545!
     546!--                u component - y-direction
     547!--                WS1 (12), WS3 (13), WS5 (14)
     548                   IF ( k <= nzb_u_inner(j+1,i) .OR.                           &
     549                            ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )    &
     550                              .AND. j == nys   )    .OR.                       &
     551                            ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )    &
     552                              .AND. j == nyn   ) )                             &
     553                   THEN
     554                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 12 )
     555                   ELSEIF ( ( k <= nzb_u_inner(j+2,i) .AND.                    &
     556                              k >  nzb_u_inner(j+1,i) ) .OR.                   &
     557                              k <= nzb_u_inner(j-1,i)                          &
     558                                                        .OR.                   &
     559                            ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )    &
     560                              .AND. j == nysv  )    .OR.                       &
     561                            ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )    &
     562                              .AND. j == nyn-1 ) )                             &
     563                   THEN
     564                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 13 )
     565!
     566!--                   Clear flag for WS1
     567                      wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 12 )
     568                   ELSEIF ( k > nzb_u_inner(j+1,i) .AND. k > nzb_u_inner(j+2,i)   &
     569                                                   .AND. k > nzb_u_inner(j-1,i) ) &
     570                   THEN
     571                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 14 )
     572!
     573!--                   Clear flag for WS1
     574                      wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 12 )
     575                   ENDIF
     576!
     577!--                u component - z-direction
     578!--                WS1 (15), WS3 (16), WS5 (17)
     579                   flag_set = .FALSE.
     580                   IF ( k == nzb_u_inner(j,i) + 1 .OR. k == nzt )  THEN
     581                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 15 )
     582                      flag_set = .TRUE.
     583                   ELSEIF ( k == nzb_u_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
     584                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 16 )
     585                      flag_set = .TRUE.
     586                   ELSEIF ( k > nzb_u_inner(j,i) .AND. .NOT. flag_set )  THEN
     587                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 17 )
     588                   ENDIF
     589
     590                ENDDO
     591             ENDDO
     592          ENDDO
     593
     594          DO  i = nxl, nxr
     595             DO  j = nys, nyn
     596                DO  k = nzb+1, nzt
     597!
     598!--                At first, set flags to WS1.
     599!--                Since fluxes are swapped in advec_ws.f90, this is necessary to
     600!--                in order to handle the left/south flux.
     601                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 )
     602                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 )
     603!
     604!--                v component - x-direction
     605!--                WS1 (18), WS3 (19), WS5 (20)
     606                   IF ( k <= nzb_v_inner(j,i+1) .OR.                           &
     607                            ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )    &
     608                              .AND. i == nxl   )    .OR.                       &
     609                            ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )    &
     610                              .AND. i == nxr   ) )                             &
     611                  THEN
     612                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 18 )
     613!
     614!--                WS3
     615                   ELSEIF ( ( k <= nzb_v_inner(j,i+2) .AND.                    &
     616                              k >  nzb_v_inner(j,i+1) ) .OR.                   &
     617                              k <= nzb_v_inner(j,i-1)                          &
     618                                                    .OR.                       &
     619                            ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )    &
     620                              .AND. i == nxr-1 )    .OR.                       &
     621                            ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )    &
     622                              .AND. i == nxlu  ) )                             &
     623                   THEN
     624                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 19 )
     625!
     626!--                   Clear flag for WS1
     627                      wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 18 )
     628                   ELSEIF ( k > nzb_v_inner(j,i+1) .AND. k > nzb_v_inner(j,i+2)   &
     629                                                   .AND. k > nzb_v_inner(j,i-1) ) &
     630                   THEN
     631                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 20 )
     632!
     633!--                   Clear flag for WS1
     634                      wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 18 )
     635                   ENDIF
     636!
     637!--                v component - y-direction
     638!--                WS1 (21), WS3 (22), WS5 (23)
     639                   IF ( k <= nzb_v_inner(j+1,i) .OR.                           &
     640                            ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )    &
     641                              .AND. j <= nysv  )    .OR.                       &
     642                            ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )    &
     643                              .AND. j == nyn   ) )                             &
     644                   THEN
     645                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 21 )
     646                   ELSEIF ( ( k <= nzb_v_inner(j+2,i) .AND.                    &
     647                              k >  nzb_v_inner(j+1,i) ) .OR.                   &
     648                              k <= nzb_v_inner(j-1,i)                          &
     649                                                        .OR.                   &
     650                            ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )    &
     651                              .AND. j == nysv+1)    .OR.                       &
     652                            ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )    &
     653                              .AND. j == nyn-1 ) )                             &
     654                   THEN
     655                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 22 )
     656!
     657!--                   Clear flag for WS1
     658                      wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 21 )
     659                   ELSEIF ( k > nzb_v_inner(j+1,i) .AND. k > nzb_v_inner(j+2,i)   &
     660                                                   .AND. k > nzb_v_inner(j-1,i) ) &
     661                   THEN
     662                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 23 )
     663!
     664!--                   Clear flag for WS1
     665                      wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 21 )
     666                   ENDIF
     667!
     668!--                v component - z-direction
     669!--                WS1 (24), WS3 (25), WS5 (26)
     670                   flag_set = .FALSE.
     671                   IF ( k == nzb_v_inner(j,i) + 1 .OR. k == nzt )  THEN
     672                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 24 )
     673                      flag_set = .TRUE.
     674                   ELSEIF ( k == nzb_v_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
     675                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 25 )
     676                      flag_set = .TRUE.
     677                   ELSEIF ( k > nzb_v_inner(j,i) .AND. .NOT. flag_set )  THEN
     678                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 26 )
     679                   ENDIF
     680
     681                ENDDO
     682             ENDDO
     683          ENDDO
     684          DO  i = nxl, nxr
     685             DO  j = nys, nyn
     686                DO  k = nzb+1, nzt
     687!
     688!--                At first, set flags to WS1.
     689!--                Since fluxes are swapped in advec_ws.f90, this is necessary to
     690!--                in order to handle the left/south flux.
     691                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 )
     692                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 )
     693!
     694!--                w component - x-direction
     695!--                WS1 (27), WS3 (28), WS5 (29)
     696                   IF ( k <= nzb_w_inner(j,i+1) .OR.                           &
     697                            ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )    &
     698                              .AND. i == nxl   )    .OR.                       &
     699                            ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )    &
     700                              .AND. i == nxr   ) )                             &
     701                   THEN
     702                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 27 )
     703                   ELSEIF ( ( k <= nzb_w_inner(j,i+2) .AND.                    &
     704                              k >  nzb_w_inner(j,i+1) ) .OR.                   &
     705                              k <= nzb_w_inner(j,i-1)                          &
     706                                                        .OR.                   &
     707                            ( ( inflow_r .OR. outflow_r .OR. nest_bound_r )    &
     708                              .AND. i == nxr-1 )    .OR.                       &
     709                            ( ( inflow_l .OR. outflow_l .OR. nest_bound_l )    &
     710                              .AND. i == nxlu  ) )                             &
     711                   THEN
     712                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 28 )
     713!   
     714!--                   Clear flag for WS1
     715                      wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 27 )
     716                   ELSEIF ( k > nzb_w_inner(j,i+1) .AND. k > nzb_w_inner(j,i+2)   &
     717                                                   .AND. k > nzb_w_inner(j,i-1) ) &
     718                   THEN
     719                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i),29 )
     720!   
     721!--                   Clear flag for WS1
     722                      wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 27 )
     723                   ENDIF
     724!
     725!--                w component - y-direction
     726!--                WS1 (30), WS3 (31), WS5 (32)
     727                   IF ( k <= nzb_w_inner(j+1,i) .OR.                           &
     728                            ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )    &
     729                              .AND. j == nys   )    .OR.                       &
     730                            ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )    &
     731                              .AND. j == nyn   ) )                             &
     732                   THEN
     733                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 30 )
     734                   ELSEIF ( ( k <= nzb_w_inner(j+2,i) .AND.                    &
     735                              k >  nzb_w_inner(j+1,i) ) .OR.                   &
     736                               k <= nzb_w_inner(j-1,i)                         &
     737                                                        .OR.                   &
     738                            ( ( inflow_s .OR. outflow_s .OR. nest_bound_s )    &
     739                              .AND. j == nysv  )    .OR.                       &
     740                            ( ( inflow_n .OR. outflow_n .OR. nest_bound_n )    &
     741                              .AND. j == nyn-1 ) )                             &
     742                   THEN
     743                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 )
     744!
     745!--                   Clear flag for WS1
     746                      wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 30 )
     747                   ELSEIF ( k > nzb_w_inner(j+1,i) .AND. k > nzb_w_inner(j+2,i)   &
     748                                                   .AND. k > nzb_w_inner(j-1,i) ) &
     749                   THEN
     750                      wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 0 )
     751!
     752!--                   Clear flag for WS1
     753                      wall_flags_0(k,j,i) = IBCLR( wall_flags_0(k,j,i), 30 )
     754                   ENDIF
     755!
     756!--                w component - z-direction
     757!--                WS1 (33), WS3 (34), WS5 (35)
     758                   flag_set = .FALSE.
     759                   IF ( k == nzb_w_inner(j,i) .OR. k == nzb_w_inner(j,i) + 1   &
     760                                              .OR. k == nzt )  THEN
     761!
     762!--                   Please note, at k == nzb_w_inner(j,i) a flag is explictely
     763!--                   set, although this is not a prognostic level. However,
     764!--                   contrary to the advection of u,v and s this is necessary
     765!--                   because flux_t(nzb_w_inner(j,i)) is used for the tendency
     766!--                   at k == nzb_w_inner(j,i)+1.
     767                      wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 1 )
     768                      flag_set = .TRUE.
     769                   ELSEIF ( k == nzb_w_inner(j,i) + 2 .OR. k == nzt - 1 )  THEN
     770                      wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 2 )
     771                      flag_set = .TRUE.
     772                   ELSEIF ( k > nzb_w_inner(j,i) .AND. .NOT. flag_set )  THEN
     773                      wall_flags_00(k,j,i) = IBSET( wall_flags_00(k,j,i), 3 )
     774                   ENDIF
     775
     776                ENDDO
     777             ENDDO
     778          ENDDO
     779
     780       ENDIF
     781
     782
     783!
     784!--    Exchange 3D integer wall_flags.
     785       IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme'  &
     786       .OR. scalar_advec == 'ws-scheme-mono' )  THEN 
     787!
     788!--       Exchange ghost points for advection flags
     789          CALL exchange_horiz_int( wall_flags_0,  nbgp )
     790          CALL exchange_horiz_int( wall_flags_00, nbgp )
     791!
     792!--       Set boundary flags at inflow and outflow boundary in case of
     793!--       non-cyclic boundary conditions.
     794         IF ( inflow_l .OR. outflow_l .OR. nest_bound_l )  THEN
     795             wall_flags_0(:,:,nxl-1)  = wall_flags_0(:,:,nxl)
     796             wall_flags_00(:,:,nxl-1) = wall_flags_00(:,:,nxl)
     797         ENDIF
     798
     799         IF ( inflow_r .OR. outflow_r .OR. nest_bound_r )  THEN
     800            wall_flags_0(:,:,nxr+1)  = wall_flags_0(:,:,nxr)
     801            wall_flags_00(:,:,nxr+1) = wall_flags_00(:,:,nxr)
     802          ENDIF
     803
     804          IF ( inflow_n .OR. outflow_n .OR. nest_bound_n )  THEN
     805             wall_flags_0(:,nyn+1,:)  = wall_flags_0(:,nyn,:)
     806             wall_flags_00(:,nyn+1,:) = wall_flags_00(:,nyn,:)
     807          ENDIF
     808
     809          IF ( inflow_s .OR. outflow_s  .OR. nest_bound_s )  THEN
     810             wall_flags_0(:,nys-1,:)  = wall_flags_0(:,nys,:)
     811             wall_flags_00(:,nys-1,:) = wall_flags_00(:,nys,:)
     812          ENDIF
     813 
     814       ENDIF
     815
     816
     817    END SUBROUTINE ws_init_flags
    366818
    367819
  • 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.