Changeset 4109 for palm/trunk/SOURCE/init_grid.f90

Timestamp:

Jul 22, 2019 5:00:34 PM (5 years ago)

Author:

suehring

Message:

Control discretization of advection term: separate initialization of WS advection flags for momentum and scalars. In this context, resort the bits and do some minor formatting; Make initialization of scalar-advection flags more flexible, i.e. introduce an arguemnt list to indicate non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols); Introduce extended 'degradation zones', where horizontal advection of passive scalars is discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3 grid points). Even though no building is within the numerical stencil, first-order scheme is used. At fourth and fifth grid point the order of the horizontal advection scheme is successively upgraded. These extended degradation zones are used to avoid stationary numerical oscillations, which are responsible for high concentration maxima that may appear under shear-free stable conditions. Therefore, an additional 3D interger array used to store control flags is introduced; Change interface for scalar advection routine; Bugfix, avoid uninitialized value sk_num in vector version of WS scalar advection; Chemistry: Decycling boundary conditions are only set at the ghost points not on the prognostic grid points; Land-surface model: Relax checks for non-consistent initialization in case static or dynamic input is provided. For example, soil_temperature or deep_soil_temperature is not mandatory any more if dynamic input is available. Also, improper settings of x_type in namelist are only checked if no static file is available.

File:

• palm/trunk/SOURCE/init_grid.f90 (modified) (8 diffs)

palm/trunk/SOURCE/init_grid.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! - Separate initialization of advection flags for momentum and scalars.
+! - Change subroutine interface for ws_init_flags_scalar to pass boundary flags
 ! 
 ! Former revisions:
@@ -373 +374 @@
  
     USE advec_ws,                                                              &
-        ONLY:  ws_init_flags
+        ONLY:  ws_init_flags_momentum, ws_init_flags_scalar
 
     USE arrays_3d,                                                             &
@@ -380 +381 @@
     USE control_parameters,                                                    &
         ONLY:  bc_lr_cyc, bc_ns_cyc,                                           &
+               bc_dirichlet_l,                                                 &
+               bc_dirichlet_n,                                                 &
+               bc_dirichlet_r,                                                 &
+               bc_dirichlet_s,                                                 &
+               bc_radiation_l,                                                 &
+               bc_radiation_n,                                                 &
+               bc_radiation_r,                                                 &
+               bc_radiation_s,                                                 &
                constant_flux_layer, dz, dz_max, dz_stretch_factor,             &
                dz_stretch_factor_array, dz_stretch_level, dz_stretch_level_end,&
@@ -392 +401 @@
         
     USE indices,                                                               &
-        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz,   &
+        ONLY:  advc_flags_m,                                                   &
+               advc_flags_s,                                                   &
+               nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz,   &
                nzb, nzb_diff, nzb_diff_s_inner, nzb_diff_s_outer,              &
                nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner,                 &
@@ -845 +856 @@
 !-- Calculate wall flag arrays for the multigrid method.
 !-- Please note, wall flags are only applied in the non-optimized version.
-    IF ( psolver == 'multigrid_noopt' )  CALL poismg_noopt_init 
+    IF ( psolver == 'multigrid_noopt' )  CALL poismg_noopt_init
 
 !
 !-- Init flags for ws-scheme to degrade order of the numerics near walls, i.e. 
-!-- to decrease the numerical stencil appropriately.
-    IF ( momentum_advec == 'ws-scheme'  .OR.  scalar_advec == 'ws-scheme' )    &
-       CALL ws_init_flags
+!-- to decrease the numerical stencil appropriately. The order of the scheme
+!-- is degraded near solid walls as well as near non-cyclic inflow and outflow
+!-- boundaries. Do this separately for momentum and scalars.
+    IF ( momentum_advec == 'ws-scheme' )  THEN
+       ALLOCATE( advc_flags_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       CALL ws_init_flags_momentum
+    ENDIF
+    IF ( scalar_advec == 'ws-scheme'   )  THEN
+       ALLOCATE( advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       advc_flags_s = 0
+       
+       CALL ws_init_flags_scalar( bc_dirichlet_l  .OR.  bc_radiation_l,        &
+                                  bc_dirichlet_n  .OR.  bc_radiation_n,        &
+                                  bc_dirichlet_r  .OR.  bc_radiation_r,        &
+                                  bc_dirichlet_s  .OR.  bc_radiation_s,        &
+                                  advc_flags_s )
+    ENDIF
 
 !
@@ -861 +886 @@
        DO  j = nys, nyn
           DO  k = nzb, nzt + 1
-             k_top = MAX( k_top, MERGE( k, 0,                                  &
-                                        .NOT. BTEST( topo(k,j,i), 0 ) ) )
+             k_top = MAX( k_top, MERGE( k, 0, .NOT. BTEST( topo(k,j,i), 0 ) ) )
           ENDDO
        ENDDO
@@ -2496 +2520 @@
     USE control_parameters,                                                    &
         ONLY:  bc_lr_cyc, bc_ns_cyc, constant_flux_layer, land_surface,        &
-               use_surface_fluxes, use_top_fluxes, urban_surface
+               scalar_advec, use_surface_fluxes, use_top_fluxes, urban_surface
 
     USE indices,                                                               &
@@ -2724 +2748 @@
           wall_flags_0(nzt+1,j,i) = IBSET( wall_flags_0(nzt+1,j,i), 22 )
           wall_flags_0(nzt+1,j,i) = IBSET( wall_flags_0(nzt+1,j,i), 23 )
+!
+!--       Set flags indicating that topography is close by in horizontal 
+!--       direction, i.e. flags that infold the topography. These will be used
+!--       to set advection flags for passive scalars, where due to large 
+!--       gradients near buildings stationary numerical oscillations can produce 
+!--       unrealistically high concentrations. This is only necessary if 
+!--       WS-scheme is applied for scalar advection. Note, these flags will be 
+!--       only used for passive scalars such as chemical species or aerosols.
+          IF ( scalar_advec == 'ws-scheme' )  THEN
+             DO k = nzb, nzt
+                IF ( BTEST( wall_flags_0(k,j,i), 0 )  .AND. (                   &
+                     ANY( .NOT. BTEST( wall_flags_0(k,j-3:j+3,i-1), 0 ) )  .OR.&
+                     ANY( .NOT. BTEST( wall_flags_0(k,j-3:j+3,i-2), 0 ) )  .OR.&
+                     ANY( .NOT. BTEST( wall_flags_0(k,j-3:j+3,i-3), 0 ) )  .OR.&
+                     ANY( .NOT. BTEST( wall_flags_0(k,j-3:j+3,i+1), 0 ) )  .OR.&
+                     ANY( .NOT. BTEST( wall_flags_0(k,j-3:j+3,i+2), 0 ) )  .OR.&
+                     ANY( .NOT. BTEST( wall_flags_0(k,j-3:j+3,i+3), 0 ) )  .OR.&
+                     ANY( .NOT. BTEST( wall_flags_0(k,j-1,i-3:i+3), 0 ) )  .OR.&
+                     ANY( .NOT. BTEST( wall_flags_0(k,j-2,i-3:i+3), 0 ) )  .OR.&
+                     ANY( .NOT. BTEST( wall_flags_0(k,j-3,i-3:i+3), 0 ) )  .OR.&
+                     ANY( .NOT. BTEST( wall_flags_0(k,j+1,i-3:i+3), 0 ) )  .OR.&
+                     ANY( .NOT. BTEST( wall_flags_0(k,j+2,i-3:i+3), 0 ) )  .OR.&
+                     ANY( .NOT. BTEST( wall_flags_0(k,j+3,i-3:i+3), 0 ) )      &
+                                                            ) )                &
+                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 31 )
+                     
+             ENDDO
+          ENDIF
        ENDDO
     ENDDO