Changeset 2232 for palm/trunk/SOURCE/lpm_init.f90

Timestamp:

May 30, 2017 5:47:52 PM (7 years ago)

Author:

suehring

Message:

Adjustments according new topography and surface-modelling concept implemented

File:

• palm/trunk/SOURCE/lpm_init.f90 (modified) (13 diffs)

palm/trunk/SOURCE/lpm_init.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
+! Adjustments according to new topography realization
 ! 
 ! 
@@ -146 +147 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  de_dx, de_dy, de_dz, zu, zw, z0
+        ONLY:  de_dx, de_dy, de_dz, zu, zw
 
     USE control_parameters,                                                    &
@@ -157 +158 @@
     USE indices,                                                               &
         ONLY:  nx, nxl, nxlg, nxrg, nxr, ny, nyn, nys, nyng, nysg, nz, nzb,    &
-               nzb_w_inner, nzt
+               nzb_max, nzt, wall_flags_0
 
     USE kinds
@@ -196 +197 @@
         ONLY:  random_function
 
+    USE surface_mod,                                                           &
+        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
+
     IMPLICIT NONE
 
@@ -287 +291 @@
        number_of_particle_groups = max_number_of_particle_groups
     ENDIF
+!
+!-- Check if downward-facing walls exist. This case, reflection boundary
+!-- conditions (as well as subgrid-scale velocities) may do not work
+!-- propably (not realized so far). 
+    IF ( surf_def_h(1)%ns >= 1 )  THEN
+       WRITE( message_string, * ) 'Overhanging topograpyh do not work '// &   
+                                  'with particles' 
+       CALL message( 'lpm_init', 'PA0212', 0, 1, 0, 6, 0 )
+
+    ENDIF 
 
 !
@@ -355 +369 @@
        
        ALLOCATE ( log_z_z0(0:number_of_sublayers) ) 
-       z_p         = zu(nzb+1) - zw(nzb)
+       z_p = zu(nzb+1) - zw(nzb)
 
 !
@@ -362 +376 @@
 !--    However, sensitivity studies showed that the effect is 
 !--    negligible. 
-       z0_av_local  = SUM( z0(nys:nyn,nxl:nxr) )
+       z0_av_local  = SUM( surf_def_h(0)%z0 ) + SUM( surf_lsm_h%z0 ) +         &
+                      SUM( surf_usm_h%z0 )
        z0_av_global = 0.0_wp
 
@@ -577 +592 @@
     INTEGER(iwp)               ::  j           !< loop variable ( particles per point )
     INTEGER(iwp)               ::  jp          !< index variable along y
+    INTEGER(iwp)               ::  k           !< index variable along z
+    INTEGER(iwp)               ::  k_surf      !< index of surface grid point
     INTEGER(iwp)               ::  kp          !< index variable along z
     INTEGER(iwp)               ::  loop_stride !< loop variable for initialization
@@ -679 +696 @@
 !
 !--                            Determine the grid indices of the particle position
-                              ip = ( tmp_particle%x + 0.5_wp * dx ) * ddx
+                               ip = ( tmp_particle%x + 0.5_wp * dx ) * ddx
                                jp = ( tmp_particle%y + 0.5_wp * dy ) * ddy
                                kp = tmp_particle%z / dz + 1 + offset_ocean_nzt
+!
+!--                            Determine surface level. Therefore, check for
+!--                            upward-facing wall on w-grid. MAXLOC will return 
+!--                            the index of the lowest upward-facing wall.
+                               k_surf = MAXLOC(                                &
+                                             MERGE( 1, 0,                      &
+                                   BTEST( wall_flags_0(nzb:nzb_max,jp,ip), 18 )&
+                                                  ), DIM = 1                   &
+                                              ) - 1
 
                                IF ( seed_follows_topography )  THEN
 !
 !--                               Particle height is given relative to topography
-                                  kp = kp + nzb_w_inner(jp,ip)
-                                  tmp_particle%z = tmp_particle%z +            &
-                                                         zw(nzb_w_inner(jp,ip))
-                                  IF ( kp > nzt )  THEN
+                                  kp = kp + k_surf
+                                  tmp_particle%z = tmp_particle%z + zw(k_surf)
+!--                               Skip particle release if particle position is
+!--                               above model top, or within topography in case
+!--                               of overhanging structures.
+                                  IF ( kp > nzt  .OR.                          &
+                                 .NOT. BTEST( wall_flags_0(kp,jp,ip), 0 ) )  THEN
                                      pos_x = pos_x + pdx(i)
                                      CYCLE xloop
                                   ENDIF
+!
+!--                            Skip particle release if particle position is 
+!--                            below surface, or within topography in case
+!--                            of overhanging structures.
                                ELSEIF ( .NOT. seed_follows_topography .AND.    &
-                                        tmp_particle%z <= zw(nzb_w_inner(jp,ip)) )  THEN
+                                         tmp_particle%z <= zw(k_surf)  .OR.    &
+                                        .NOT. BTEST( wall_flags_0(kp,jp,ip), 0 ) )&
+                               THEN
                                   pos_x = pos_x + pdx(i)
                                   CYCLE xloop                               
@@ -820 +855 @@
                                      pdx(particles(n)%group)
                       particles(n)%x = particles(n)%x +                        &
-                              MERGE( rand_contr, SIGN( dx, rand_contr ), &
+                              MERGE( rand_contr, SIGN( dx, rand_contr ),       &
                                      ABS( rand_contr ) < dx                    &
                                    ) 
@@ -828 +863 @@
                                      pdy(particles(n)%group)
                       particles(n)%y = particles(n)%y +                        &
-                              MERGE( rand_contr, SIGN( dy, rand_contr ), &
+                              MERGE( rand_contr, SIGN( dy, rand_contr ),       &
                                      ABS( rand_contr ) < dy                    &
                                    ) 
@@ -836 +871 @@
                                      pdz(particles(n)%group)
                       particles(n)%z = particles(n)%z +                        &
-                              MERGE( rand_contr, SIGN( dz, rand_contr ), &
+                              MERGE( rand_contr, SIGN( dz, rand_contr ),       &
                                      ABS( rand_contr ) < dz                    &
                                    ) 
@@ -854 +889 @@
                    i = ( particles(n)%x + 0.5_wp * dx ) * ddx
                    j = ( particles(n)%y + 0.5_wp * dy ) * ddy
-                   IF ( particles(n)%z <= zw(nzb_w_inner(j,i)) )  THEN
+                   k =   particles(n)%z / dz + 1 + offset_ocean_nzt
+!
+!--                Check if particle is within topography
+                   IF ( .NOT. BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
                       particles(n)%particle_mask = .FALSE.
                       deleted_particles = deleted_particles + 1
                    ENDIF
+
                 ENDDO
              ENDDO