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

Timestamp:

Mar 28, 2012 6:44:41 AM (12 years ago)

Author:

raasch

Message:

bugfix for timestep calculation in case of Galilei transformation

File:

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

palm/trunk/SOURCE/timestep.f90

@@ -3 +3 @@
 !------------------------------------------------------------------------------!
 ! Current revisions:
-! -----------------
-! 
+! ------------------
+! bugfix for timestep calculation in case of Galilei transformation,
+! special treatment in case of mirror velocity boundary condition removed
 !
 ! Former revisions:
@@ -56 +57 @@
     IMPLICIT NONE
 
-    INTEGER ::  i, j, k
+    INTEGER ::  i, j, k, u_max_cfl_ijk(3), v_max_cfl_ijk(3)
 
     REAL ::  div, dt_diff, dt_diff_l, dt_plant_canopy,                 &
@@ -62 +63 @@
              dt_plant_canopy_u, dt_plant_canopy_v, dt_plant_canopy_w,  & 
              dt_u, dt_v, dt_w, lad_max, percent_change,                &
-             u_gtrans_l, vabs_max, value, v_gtrans_l
+             u_gtrans_l, u_max_cfl, vabs_max, value, v_gtrans_l, v_max_cfl
 
     REAL, DIMENSION(2)         ::  uv_gtrans, uv_gtrans_l
@@ -70 +71 @@
 
     CALL cpu_log( log_point(12), 'calculate_timestep', 'start' )
-
-!
-!-- Determine the maxima of the velocity components.
-    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', &
-                         u_max, u_max_ijk )
-    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'abs', &
-                         v_max, v_max_ijk )
-    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', &
-                         w_max, w_max_ijk )
-
-!
-!-- In case maxima of the horizontal velocity components have been found at the
-!-- bottom boundary (k=nzb), the corresponding maximum at level k=1 is chosen
-!-- if the Dirichlet-boundary condition ('mirror') has been set. This is 
-!-- necessary, because otherwise in case of Galilei-transform a far too large
-!-- velocity (having the respective opposite sign) would be used for the time
-!-- step determination (almost double the mean flow velocity).
-    IF ( ibc_uv_b == 0 )  THEN
-       IF ( u_max_ijk(1) == nzb )  THEN
-          u_max        = -u_max
-          u_max_ijk(1) = nzb + 1
-       ENDIF
-       IF ( v_max_ijk(1) == nzb )  THEN
-          v_max        = -v_max
-          v_max_ijk(1) = nzb + 1
-       ENDIF
-    ENDIF
 
 !
@@ -142 +116 @@
     ENDIF
 
+!
+!-- Determine the maxima of the velocity components.
+    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', 0.0, &
+                         u_max, u_max_ijk )
+    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'abs', 0.0, &
+                         v_max, v_max_ijk )
+    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', 0.0, &
+                         w_max, w_max_ijk )
+
+!
+!-- In case of Galilei transformation, the horizontal velocity maxima have
+!-- to be calculated from the transformed horizontal velocities
+    IF ( galilei_transformation )  THEN
+       CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'absoff', &
+                            u_gtrans, u_max_cfl, u_max_cfl_ijk )
+       CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'absoff', &
+                            v_gtrans, v_max_cfl, v_max_cfl_ijk )
+    ELSE
+       u_max_cfl = u_max
+       v_max_cfl = v_max
+       u_max_cfl_ijk = u_max_ijk
+       v_max_cfl_ijk = v_max_ijk
+    ENDIF
+
+
     IF ( .NOT. dt_fixed )  THEN
 !
@@ -147 +146 @@
 !
 !--    For each component, compute the maximum time step according to the 
-!--    FCL-criterion.
-       dt_u = dx / ( ABS( u_max - u_gtrans ) + 1.0E-10 )
-       dt_v = dy / ( ABS( v_max - v_gtrans ) + 1.0E-10 )
+!--    CFL-criterion.
+       dt_u = dx / ( ABS( u_max_cfl ) + 1.0E-10 )
+       dt_v = dy / ( ABS( v_max_cfl ) + 1.0E-10 )
        dt_w = dzu(MAX( 1, w_max_ijk(1) )) / ( ABS( w_max ) + 1.0E-10 )
 
@@ -267 +266 @@
 #if defined( __parallel )
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( dt_plant_canopy_l, dt_plant_canopy, 1, MPI_REAL,  &
+          CALL MPI_ALLREDUCE( dt_plant_canopy_l, dt_plant_canopy, 1, MPI_REAL, &
                               MPI_MIN, comm2d, ierr )
 #else
@@ -318 +317 @@
                '&dt_diff         = ', dt_diff, ' s',                          &
                '&dt_plant_canopy = ', dt_plant_canopy, ' s',                  &
-               '&u_max   = ', u_max, ' m/s   k=', u_max_ijk(1),               &
+               '&u_max_cfl   = ', u_max_cfl, ' m/s   k=', u_max_cfl_ijk(1),   &
                '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                    &
-               '&v_max   = ', v_max, ' m/s   k=', v_max_ijk(1),               &
+               '&v_max_cfl   = ', v_max_cfl, ' m/s   k=', v_max_cfl_ijk(1),   &
                '  j=', v_max_ijk(2), '  i=', v_max_ijk(3),                    &
-               '&w_max   = ', w_max, ' m/s   k=', w_max_ijk(1),               &
+               '&w_max       = ', w_max, ' m/s   k=', w_max_ijk(1),           &
                '  j=', w_max_ijk(2), '  i=', w_max_ijk(3)
           CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )