Changeset 316 for palm

Timestamp:

May 13, 2009 3:31:23 PM (15 years ago)

Author:

letzel

Message:

• Additional timestep criterion in case of simulations with plant canopy (modules, timestep)

Location:

palm/trunk/SOURCE

Files:

• CURRENT_MODIFICATIONS (modified) (2 diffs)
• modules.f90 (modified) (2 diffs)
• timestep.f90 (modified) (5 diffs)

palm/trunk/SOURCE/CURRENT_MODIFICATIONS

@@ -1 +1 @@
 New:
 ---
+Additional timestep criterion in case of simulations with plant canopy (modules, timestep)
+
 Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1 
 (check_parameters)
@@ -35 +37 @@
 write_var_list)
 
-check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, modules, mrun, package_parin, parin, prognostic_equations, read_var_list, surface_coupler, time_integration, user_check_parameters, user_data_output_dvrp, user_header, user_init_grid, write_var_list
+check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, modules, mrun, package_parin, parin, prognostic_equations, read_var_list, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_dvrp, user_header, user_init_grid, write_var_list
 
 New: set_particle_attributes, set_slicer_attributes_dvrp

palm/trunk/SOURCE/modules.f90

@@ -5 +5 @@
 ! Actual revisions:
 ! -----------------
+! +cdc_max, lad_u_max, lad_v_max, lad_w_max, 
+! cdc_max_ijk, lad_u_max_ijk, lad_v_max_ijk, lad_w_max_ijk
 ! typographical error in dots_unit fixed
 ! +clip_dvrp_*, cluster_size, color_interval, dvrpsize_interval, dvrp_overlap,
@@ -1217 +1219 @@
     CHARACTER (LEN=40) ::  region(0:9)
     INTEGER ::  pr_palm = 80, statistic_regions = 0, var_ts = 100
-    INTEGER ::  u_max_ijk(3), v_max_ijk(3), w_max_ijk(3)
+    INTEGER ::  cdc_max_ijk(3), lad_u_max_ijk(3), lad_v_max_ijk(3),  &
+                lad_w_max_ijk(3), u_max_ijk(3), v_max_ijk(3), w_max_ijk(3)
     LOGICAL ::  flow_statistics_called = .FALSE.
-    REAL ::     u_max, v_max, w_max
+    REAL ::     cdc_max, lad_u_max, lad_v_max, lad_w_max, u_max, v_max, w_max
     REAL, DIMENSION(:), ALLOCATABLE       ::  sums_divnew_l, sums_divold_l
     REAL, DIMENSION(:,:), ALLOCATABLE     ::  sums, sums_wsts_bc_l, ts_value

palm/trunk/SOURCE/timestep.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
+! Additional timestep criterion in case of simulations with plant canopy
+!
 ! Output of messages replaced by message handling routine.
 !
@@ -45 +47 @@
     INTEGER ::  i, j, k
 
-    REAL ::  div, dt_diff, dt_diff_l, dt_u, dt_v, dt_w, percent_change, &
-             u_gtrans_l, value, v_gtrans_l
+    REAL ::  div, dt_diff, dt_diff_l, dt_plant_canopy,  & 
+             dt_u, dt_v, dt_w, lad_max, percent_change, &
+             u_gtrans_l, vabs_max, value, v_gtrans_l
 
     REAL, DIMENSION(2)         ::  uv_gtrans, uv_gtrans_l
@@ -176 +179 @@
 
 !
-!--    The time step is the minimum of the 3 components and the diffusion time 
+!--    Additional timestep criterion with plant canopies: 
+!--    it is not allowed to extract more than the available momentum
+       IF ( plant_canopy ) THEN
+!
+!--       Determine the maxima of the leaf drag coefficient and the leaf 
+!--       area density.
+          CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, cdc, &
+                               'max', cdc_max, cdc_max_ijk )
+          CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, lad_u, &
+                               'max', lad_u_max, lad_u_max_ijk )
+          CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, lad_v, &
+                               'max', lad_v_max, lad_v_max_ijk )
+          CALL global_min_max( nzb, nzt+1, nys-1, nyn+1, nxl-1, nxr+1, lad_w, &
+                               'max', lad_w_max, lad_w_max_ijk )
+!
+!--       Determine the potential maximum of the absolute value of the velocity
+          vabs_max = SQRT( u_max**2 + v_max**2 + w_max**2 )
+          lad_max  = MAX( lad_u_max, lad_v_max, lad_w_max )
+
+          dt_plant_canopy = 1.0 / ( cdc_max * lad_max * vabs_max )
+
+       ELSE
+!
+!--       Use dt_diff as dummy value to avoid extra IF branches further below
+          dt_plant_canopy = dt_diff
+
+       ENDIF
+
+!
+!--    The time step is the minimum of the 3-4 components and the diffusion time
 !--    step minus a reduction to be on the safe side. Factor 0.5 is necessary 
 !--    since the leap-frog scheme always progresses by 2 * delta t.
@@ -183 +215 @@
 !--    The time step must not exceed the maximum allowed value.
        IF ( timestep_scheme(1:5) == 'runge' )  THEN
-          dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w )
-       ELSE
-          dt_3d = cfl_factor * 0.5 * MIN( dt_diff, dt_u, dt_v, dt_w )
+          dt_3d = cfl_factor * MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w )
+       ELSE
+          dt_3d = cfl_factor * 0.5 *  &
+                  MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w )
        ENDIF
        dt_3d = MIN( dt_3d, dt_max )
@@ -191 +224 @@
 !
 !--    Remember the restricting time step criterion for later output.
-       IF ( dt_diff > MIN( dt_u, dt_v, dt_w ) )  THEN
+       IF ( MIN( dt_u, dt_v, dt_w ) < MIN( dt_diff, dt_plant_canopy ) )  THEN
           timestep_reason = 'A'
+       ELSEIF ( dt_plant_canopy < dt_diff )  THEN
+          timestep_reason = 'C'
        ELSE
           timestep_reason = 'D'