Changeset 316 for palm/trunk/SOURCE/timestep.f90

Timestamp:

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

Author:

letzel

Message:

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

File:

• palm/trunk/SOURCE/timestep.f90 (modified) (5 diffs)

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'