Changeset 4335 for palm/trunk

Timestamp:

Dec 12, 2019 4:39:05 PM (4 years ago)

Author:

suehring

Message:

Plant canopy: implement fix for LAD on building edges also for the vectorl-optimized branch

Location:

palm/trunk/SOURCE

Files:

• plant_canopy_model_mod.f90 (modified) (6 diffs)
• synthetic_turbulence_generator_mod.f90 (modified) (2 diffs)

palm/trunk/SOURCE/plant_canopy_model_mod.f90

@@ -27 +27 @@
 ! -----------------
 ! $Id$
+! Fix for LAD at building edges also implemented in vector branch.
+! 
+! 4331 2019-12-10 18:25:02Z suehring
 ! Typo corrected
 ! 
@@ -1566 +1569 @@
        INTEGER(iwp) ::  kk        !< running index for flat lad arrays
 
+       LOGICAL ::  building_edge_e !< control flag indicating an eastward-facing building edge
+       LOGICAL ::  building_edge_n !< control flag indicating a north-facing building edge
+       LOGICAL ::  building_edge_s !< control flag indicating a south-facing building edge
+       LOGICAL ::  building_edge_w !< control flag indicating a westward-facing building edge
+
        REAL(wp) ::  ddt_3d    !< inverse of the LES timestep (dt_3d)
        REAL(wp) ::  lad_local !< local lad value
@@ -1586 +1594 @@
                 DO  j = nys, nyn
 !
+!--                Set control flags indicating east- and westward-orientated 
+!--                building edges. Note, building_egde_w is set from the perspective
+!--                of the potential rooftop grid point, while building_edge_e is
+!--                is set from the perspective of the non-building grid point.
+                   building_edge_w = ANY( BTEST( wall_flags_static_0(:,j,i),   6 ) )&
+                        .AND.  .NOT. ANY( BTEST( wall_flags_static_0(:,j,i-1), 6 ) )
+                   building_edge_e = ANY( BTEST( wall_flags_static_0(:,j,i-1), 6 ) )&
+                        .AND.  .NOT. ANY( BTEST( wall_flags_static_0(:,j,i),   6 ) )
+!
 !--                Determine topography-top index on u-grid
                    k_wall = topo_top_ind(j,i,1)
@@ -1593 +1610 @@
 !
 !--                   In order to create sharp boundaries of the plant canopy, 
-!--                   the lad on the u-grid at index (k,j,i) is equal to 
-!--                   lad_s(k,j,i), rather than being interpolated from the 
-!--                   surrounding lad_s, because this would yield smaller lad 
-!--                   at the canopy boundaries than inside of the canopy.
+!--                   the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
+!--                   rather than being interpolated from the surrounding lad_s, 
+!--                   because this would yield smaller lad at the canopy boundaries 
+!--                   than inside of the canopy.
 !--                   For the same reason, the lad at the rightmost(i+1)canopy 
-!--                   boundary on the u-grid equals lad_s(k,j,i).
+!--                   boundary on the u-grid equals lad_s(k,j,i), which is considered
+!--                   in the next if-statement. Note, at left-sided building edges 
+!--                   this is not applied, here the LAD is equals the LAD at grid
+!--                   point (k,j,i), in order to avoid that LAD is mistakenly mapped 
+!--                   on top of a roof where (usually) is no LAD is defined. 
                       lad_local = lad_s(kk,j,i)
-                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp )&
-                      THEN
-                         lad_local = lad_s(kk,j,i-1)
-                      ENDIF
+                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j,i-1) > 0.0_wp  &
+                           .AND.  .NOT. building_edge_w )  lad_local = lad_s(kk,j,i-1)
+!
+!--                   In order to avoid that LAD is mistakenly considered at right-
+!--                   sided building edges (here the topography-top index for the
+!--                   u-component at index j,i is still on the building while the
+!--                   topography top for the scalar isn't), LAD is taken from grid
+!--                   point (j,i-1).
+                      IF ( lad_local > 0.0_wp .AND. lad_s(kk,j,i-1) == 0.0_wp  &
+                           .AND.  building_edge_e )  lad_local = lad_s(kk,j,i-1)
 
                       pre_tend = 0.0_wp
@@ -1650 +1677 @@
                 DO  j = nysv, nyn
 !
+!--                Set control flags indicating north- and southward-orientated 
+!--                building edges. Note, building_egde_s is set from the perspective
+!--                of the potential rooftop grid point, while building_edge_n is
+!--                is set from the perspective of the non-building grid point.
+                   building_edge_s = ANY( BTEST( wall_flags_static_0(:,j,i),   6 ) )&
+                        .AND.  .NOT. ANY( BTEST( wall_flags_static_0(:,j-1,i), 6 ) )
+                   building_edge_n = ANY( BTEST( wall_flags_static_0(:,j-1,i), 6 ) )&
+                        .AND.  .NOT. ANY( BTEST( wall_flags_static_0(:,j,i),   6 ) )
+!
 !--                Determine topography-top index on v-grid
                    k_wall = topo_top_ind(j,i,2)
@@ -1658 +1694 @@
 !
 !--                   In order to create sharp boundaries of the plant canopy, 
-!--                   the lad on the v-grid at index (k,j,i) is equal to 
-!--                   lad_s(k,j,i), rather than being interpolated from the 
-!--                   surrounding lad_s, because this would yield smaller lad 
-!--                   at the canopy boundaries than inside of the canopy.
-!--                   For the same reason, the lad at the northmost(j+1) canopy 
-!--                   boundary on the v-grid equals lad_s(k,j,i).
+!--                   the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
+!--                   rather than being interpolated from the surrounding lad_s, 
+!--                   because this would yield smaller lad at the canopy boundaries 
+!--                   than inside of the canopy.
+!--                   For the same reason, the lad at the northmost(j+1)canopy 
+!--                   boundary on the v-grid equals lad_s(k,j,i), which is considered
+!--                   in the next if-statement. Note, at left-sided building edges 
+!--                   this is not applied, here the LAD is equals the LAD at grid
+!--                   point (k,j,i), in order to avoid that LAD is mistakenly mapped 
+!--                   on top of a roof where (usually) is no LAD is defined. 
                       lad_local = lad_s(kk,j,i)
-                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp )&
-                      THEN
-                         lad_local = lad_s(kk,j-1,i)
-                      ENDIF
+                      IF ( lad_local == 0.0_wp .AND. lad_s(kk,j-1,i) > 0.0_wp &
+                       .AND.  .NOT. building_edge_s )  lad_local = lad_s(kk,j-1,i)
+!
+!--                   In order to avoid that LAD is mistakenly considered at right-
+!--                   sided building edges (here the topography-top index for the
+!--                   u-component at index j,i is still on the building while the
+!--                   topography top for the scalar isn't), LAD is taken from grid
+!--                   point (j,i-1). 
+                      IF ( lad_local > 0.0_wp .AND. lad_s(kk,j-1,i) == 0.0_wp  &
+                       .AND.  building_edge_n )  lad_local = lad_s(kk,j-1,i)
 
                       pre_tend = 0.0_wp

palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Commentation of last commit
+! 
+! 4332 2019-12-10 19:44:12Z suehring
 ! Limit initial velocity seeds in restart runs, if not the seed calculation 
 ! may become unstable. Further, minor bugfix in initial velocity seed 
@@ -939 +942 @@
           DO  i = nxlg, nxrg
              DO  k = nzb, nzt+1
-
+!
+!--             In case the correlation coefficients are very small, the 
+!--             velocity seeds may become very large finally creating 
+!--             numerical instabilities in the synthetic turbulence generator.
+!--             Empirically, a value of 10E-8 seems to be sufficient.
                 IF  ( a11(k) > 10E-8_wp )  THEN
                    fu_xz(k,i) = ( u(k,j,i) - u_init(k) ) / a11(k)