Changeset 76 for palm/trunk/SOURCE

Timestamp:

Mar 29, 2007 12:58:32 AM (18 years ago)

Author:

raasch

Message:

changes for Neumann p conditions both at bottom and top

Location:

palm/trunk/SOURCE

Files:

• CURRENT_MODIFICATIONS (modified) (2 diffs)
• header.f90 (modified) (1 diff)
• poisfft.f90 (modified) (5 diffs)
• pres.f90 (modified) (3 diffs)

palm/trunk/SOURCE/CURRENT_MODIFICATIONS

@@ -35 +35 @@
 General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d
 
+The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.
+
 +age_m in particle_type
 
@@ -58 +60 @@
 The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.
 
-advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary
+advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary
 
 

palm/trunk/SOURCE/header.f90

@@ -283 +283 @@
     WRITE ( io, 254 )  nx, ny, nzt+1, MIN( nnx, nx+1 ), MIN( nny, ny+1 ), &
                        MIN( nnz+2, nzt+2 )
-    IF ( nxa == nx  .AND.  nya == ny  .AND.  nza == nz )  THEN
-       WRITE ( io, 255 )
-    ELSE
-       WRITE ( io, 256 )  nnx-(nxa-nx), nny-(nya-ny), nzt+2
+    IF ( numprocs > 1 )  THEN
+       IF ( nxa == nx  .AND.  nya == ny  .AND.  nza == nz )  THEN
+          WRITE ( io, 255 )
+       ELSE
+          WRITE ( io, 256 )  nnx-(nxa-nx), nny-(nya-ny), nzt+2
+       ENDIF
     ENDIF
     IF ( sloping_surface )  WRITE ( io, 260 )  alpha_surface

palm/trunk/SOURCE/poisfft.f90

@@ -4 +4 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Tridiagonal coefficients adjusted for Neumann boundary conditions both at
+! the bottom and the top.
 !
 ! Former revisions:
@@ -393 +394 @@
 !------------------------------------------------------------------------------!
 
+          USE control_parameters
+
           IMPLICIT NONE
 
@@ -443 +446 @@
              ENDDO
           ENDDO
+
+!
+!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
+!--       The respective values of ar should be zero at all k-levels if
+!--       acceleration of horizontally averaged vertical velocity is zero.
+          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
+             IF ( j == 0  .AND.  nxl_z == 0 )  THEN
+#if defined( __parallel )
+                DO  k = 1, nz
+                   ar(nxl_z,j,k) = 0.0
+                ENDDO
+#else
+                DO  k = 1, nz
+                   ar(k,j,nxl_z) = 0.0
+                ENDDO
+#endif
+             ENDIF
+          ENDIF
 
        END SUBROUTINE substi
@@ -1466 +1487 @@
           IMPLICIT NONE
 
-          INTEGER ::  i, j, k
+          INTEGER ::  i, k
 
           REAL, DIMENSION(0:nx,nz)       ::  ar
@@ -1495 +1516 @@
           ENDDO
 
+!
+!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
+!--       The respective values of ar should be zero at all k-levels if
+!--       acceleration of horizontally averaged vertical velocity is zero.
+          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
+             IF ( j == 0 )  THEN
+                DO  k = 1, nz
+                   ar(0,k) = 0.0
+                ENDDO
+             ENDIF
+          ENDIF
+
        END SUBROUTINE substi_1dd
 

palm/trunk/SOURCE/pres.f90

@@ -5 +5 @@
 ! -----------------
 ! Volume flow control for non-cyclic boundary conditions added (currently only
-! for the north boundary!!), 2nd+3rd argument removed from exchange horiz
+! for the north boundary!!), 2nd+3rd argument removed from exchange horiz,
+! mean vertical velocity is removed in case of Neumann boundary conditions
+! both at the bottom and the top
 !
 ! Former revisions:
@@ -45 +47 @@
 
     REAL, DIMENSION(1:2) ::  volume_flow_l, volume_flow_offset
+    REAL, DIMENSION(1:nzt) ::  w_l, w_l_l
 
 
@@ -96 +99 @@
     ENDIF
 
+!
+!-- Remove mean vertical velocity
+    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
+       IF ( simulated_time > 0.0 )  THEN ! otherwise nzb_w_inner is not yet known
+          w_l = 0.0;  w_l_l = 0.0
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb_w_inner(j,i)+1, nzt
+                   w_l_l(k) = w_l_l(k) + w(k,j,i)
+                ENDDO
+             ENDDO
+          ENDDO
+#if defined( __parallel )   
+          CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, &
+                              ierr )
+#else
+          w_l = w_l_l  
+#endif
+          DO  k = 1, nzt
+             w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
+          ENDDO
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  k = nzb_w_inner(j,i)+1, nzt
+                   w(k,j,i) = w(k,j,i) - w_l(k)
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+    ENDIF
 
 !