Ignore:
Timestamp:
Mar 28, 2012 6:44:41 AM (9 years ago)
Author:
raasch
Message:

bugfix for timestep calculation in case of Galilei transformation

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/timestep.f90

    r708 r866  
    33!------------------------------------------------------------------------------!
    44! Current revisions:
    5 ! -----------------
    6 !
     5! ------------------
     6! bugfix for timestep calculation in case of Galilei transformation,
     7! special treatment in case of mirror velocity boundary condition removed
    78!
    89! Former revisions:
     
    5657    IMPLICIT NONE
    5758
    58     INTEGER ::  i, j, k
     59    INTEGER ::  i, j, k, u_max_cfl_ijk(3), v_max_cfl_ijk(3)
    5960
    6061    REAL ::  div, dt_diff, dt_diff_l, dt_plant_canopy,                 &
     
    6263             dt_plant_canopy_u, dt_plant_canopy_v, dt_plant_canopy_w,  &
    6364             dt_u, dt_v, dt_w, lad_max, percent_change,                &
    64              u_gtrans_l, vabs_max, value, v_gtrans_l
     65             u_gtrans_l, u_max_cfl, vabs_max, value, v_gtrans_l, v_max_cfl
    6566
    6667    REAL, DIMENSION(2)         ::  uv_gtrans, uv_gtrans_l
     
    7071
    7172    CALL cpu_log( log_point(12), 'calculate_timestep', 'start' )
    72 
    73 !
    74 !-- Determine the maxima of the velocity components.
    75     CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', &
    76                          u_max, u_max_ijk )
    77     CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'abs', &
    78                          v_max, v_max_ijk )
    79     CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', &
    80                          w_max, w_max_ijk )
    81 
    82 !
    83 !-- In case maxima of the horizontal velocity components have been found at the
    84 !-- bottom boundary (k=nzb), the corresponding maximum at level k=1 is chosen
    85 !-- if the Dirichlet-boundary condition ('mirror') has been set. This is
    86 !-- necessary, because otherwise in case of Galilei-transform a far too large
    87 !-- velocity (having the respective opposite sign) would be used for the time
    88 !-- step determination (almost double the mean flow velocity).
    89     IF ( ibc_uv_b == 0 )  THEN
    90        IF ( u_max_ijk(1) == nzb )  THEN
    91           u_max        = -u_max
    92           u_max_ijk(1) = nzb + 1
    93        ENDIF
    94        IF ( v_max_ijk(1) == nzb )  THEN
    95           v_max        = -v_max
    96           v_max_ijk(1) = nzb + 1
    97        ENDIF
    98     ENDIF
    9973
    10074!
     
    142116    ENDIF
    143117
     118!
     119!-- Determine the maxima of the velocity components.
     120    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', 0.0, &
     121                         u_max, u_max_ijk )
     122    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'abs', 0.0, &
     123                         v_max, v_max_ijk )
     124    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', 0.0, &
     125                         w_max, w_max_ijk )
     126
     127!
     128!-- In case of Galilei transformation, the horizontal velocity maxima have
     129!-- to be calculated from the transformed horizontal velocities
     130    IF ( galilei_transformation )  THEN
     131       CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'absoff', &
     132                            u_gtrans, u_max_cfl, u_max_cfl_ijk )
     133       CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'absoff', &
     134                            v_gtrans, v_max_cfl, v_max_cfl_ijk )
     135    ELSE
     136       u_max_cfl = u_max
     137       v_max_cfl = v_max
     138       u_max_cfl_ijk = u_max_ijk
     139       v_max_cfl_ijk = v_max_ijk
     140    ENDIF
     141
     142
    144143    IF ( .NOT. dt_fixed )  THEN
    145144!
     
    147146!
    148147!--    For each component, compute the maximum time step according to the
    149 !--    FCL-criterion.
    150        dt_u = dx / ( ABS( u_max - u_gtrans ) + 1.0E-10 )
    151        dt_v = dy / ( ABS( v_max - v_gtrans ) + 1.0E-10 )
     148!--    CFL-criterion.
     149       dt_u = dx / ( ABS( u_max_cfl ) + 1.0E-10 )
     150       dt_v = dy / ( ABS( v_max_cfl ) + 1.0E-10 )
    152151       dt_w = dzu(MAX( 1, w_max_ijk(1) )) / ( ABS( w_max ) + 1.0E-10 )
    153152
     
    267266#if defined( __parallel )
    268267          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    269           CALL MPI_ALLREDUCE( dt_plant_canopy_l, dt_plant_canopy, 1, MPI_REAL,  &
     268          CALL MPI_ALLREDUCE( dt_plant_canopy_l, dt_plant_canopy, 1, MPI_REAL, &
    270269                              MPI_MIN, comm2d, ierr )
    271270#else
     
    318317               '&dt_diff         = ', dt_diff, ' s',                          &
    319318               '&dt_plant_canopy = ', dt_plant_canopy, ' s',                  &
    320                '&u_max   = ', u_max, ' m/s   k=', u_max_ijk(1),               &
     319               '&u_max_cfl   = ', u_max_cfl, ' m/s   k=', u_max_cfl_ijk(1),   &
    321320               '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                    &
    322                '&v_max   = ', v_max, ' m/s   k=', v_max_ijk(1),               &
     321               '&v_max_cfl   = ', v_max_cfl, ' m/s   k=', v_max_cfl_ijk(1),   &
    323322               '  j=', v_max_ijk(2), '  i=', v_max_ijk(3),                    &
    324                '&w_max   = ', w_max, ' m/s   k=', w_max_ijk(1),               &
     323               '&w_max       = ', w_max, ' m/s   k=', w_max_ijk(1),           &
    325324               '  j=', w_max_ijk(2), '  i=', w_max_ijk(3)
    326325          CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
Note: See TracChangeset for help on using the changeset viewer.