Changeset 866


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

bugfix for timestep calculation in case of Galilei transformation

Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

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

    r863 r866  
    44! Current revisions:
    55! -----------------
     6! use only 60% of the geostrophic wind as translation speed in case of Galilean
     7! transformation and use_ug_for_galilei_tr = .T. in order to mimimize the
     8! timestep
    69!
    710! Former revisions:
     
    13281331            ug_vertical_gradient_level(1) == 0.0 .AND. &
    13291332            vg_vertical_gradient_level(1) == 0.0 )  THEN
    1330           u_gtrans = ug_surface
    1331           v_gtrans = vg_surface
     1333          u_gtrans = ug_surface * 0.6
     1334          v_gtrans = vg_surface * 0.6
    13321335       ELSEIF ( use_ug_for_galilei_tr .AND.                &
    13331336                ug_vertical_gradient_level(1) /= 0.0 )  THEN
  • palm/trunk/SOURCE/global_min_max.f90

    r668 r866  
    1  SUBROUTINE global_min_max( i1, i2, j1, j2, k1, k2, ar, mode, value, &
     1 SUBROUTINE global_min_max( i1, i2, j1, j2, k1, k2, ar, mode, offset, value, &
    22                            value_ijk, value1, value1_ijk )
    33
    44!------------------------------------------------------------------------------!
    55! Current revisions:
    6 ! -----------------
     6! ------------------
     7! new mode "absoff" accounts for an offset in the respective array
    78!
    89! Former revisions:
     
    4344                          fmin_ijk_l(3), value_ijk(3)
    4445    INTEGER, OPTIONAL ::  value1_ijk(3)
    45     REAL              ::  value, &
     46    REAL              ::  offset, value, &
    4647                          ar(i1:i2,j1:j2,k1:k2)
    4748#if defined( __ibm )
     
    6263!--    Determine the local minimum
    6364       fmin_ijk_l = MINLOC( ar )
    64        fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - nbgp    ! MINLOC assumes lowerbound = 1
     65       fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - nbgp ! MINLOC assumes lowerbound = 1
    6566       fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - nbgp
    6667       fmin_ijk_l(3) = k1 + fmin_ijk_l(3) - 1
     
    103104!--    Determine the local maximum
    104105       fmax_ijk_l = MAXLOC( ar )
    105        fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - nbgp    ! MAXLOC assumes lowerbound = 1
     106       fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - nbgp ! MAXLOC assumes lowerbound = 1
    106107       fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - nbgp
    107108       fmax_ijk_l(3) = k1 + fmax_ijk_l(3) - 1
     
    198199
    199200!
     201!-- Determine absolute maximum of ( array - offset )
     202    IF ( mode == 'absoff' )  THEN
     203
     204!
     205!--    Determine the local absolut maximum
     206       fmax_l(1)     = 0.0
     207       fmax_ijk_l(1) =  i1
     208       fmax_ijk_l(2) =  j1
     209       fmax_ijk_l(3) =  k1
     210       DO  k = k1, k2
     211          DO  j = j1, j2
     212!
     213!--          Attention: the lowest gridpoint is excluded here, because there
     214!--          ---------  is no advection at nzb=0 and mode 'absoff' is only
     215!--                     used for calculating u,v extrema for CFL-criteria
     216             DO  i = i1+1, i2
     217                IF ( ABS( ar(i,j,k) - offset ) > fmax_l(1) )  THEN
     218                   fmax_l(1) = ABS( ar(i,j,k) - offset )
     219                   fmax_ijk_l(1) = i
     220                   fmax_ijk_l(2) = j
     221                   fmax_ijk_l(3) = k
     222                ENDIF
     223             ENDDO
     224          ENDDO
     225       ENDDO
     226
     227!
     228!--    Set a flag in case that the determined value is negative.
     229!--    A constant offset has to be subtracted in order to handle the special
     230!--    case i=0 correctly
     231       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0 )  THEN
     232          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
     233       ENDIF
     234
     235#if defined( __parallel )
     236       fmax_l(2)  = myid
     237       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     238       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
     239                           ierr )
     240
     241!
     242!--    Determine the global absolut maximum. Result stored on PE0.
     243       id_fmax = fmax(2)
     244       IF ( id_fmax /= 0 )  THEN
     245          IF ( myid == 0 )  THEN
     246             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
     247                            status, ierr )
     248          ELSEIF ( myid == id_fmax )  THEN
     249             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
     250          ENDIF
     251       ELSE
     252          fmax_ijk = fmax_ijk_l
     253       ENDIF
     254!
     255!--    Send the indices of the just determined absolut maximum to other PEs
     256       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
     257#else
     258       fmax(1)  = fmax_l(1)
     259       fmax_ijk = fmax_ijk_l
     260#endif
     261
     262    ENDIF
     263
     264!
    200265!-- Determine output parameters
    201266    SELECT CASE( mode )
     
    218283          value1_ijk = fmax_ijk
    219284
    220        CASE( 'abs' )
     285       CASE( 'abs', 'absoff' )
    221286
    222287          value     = fmax(1)
  • palm/trunk/SOURCE/modules.f90

    r863 r866  
    44! Current revisions:
    55! -----------------
     6! interface for global_min_max changed
    67!
    78! Former revisions:
     
    978979    INTERFACE
    979980
    980        SUBROUTINE global_min_max ( i1, i2, j1, j2, k1, k2, feld, mode, wert, &
    981                                    wert_ijk, wert1, wert1_ijk )
     981       SUBROUTINE global_min_max ( i1, i2, j1, j2, k1, k2, feld, mode, offset, &
     982                                   wert, wert_ijk, wert1, wert1_ijk )
    982983          CHARACTER (LEN=*), INTENT(IN) ::  mode
    983984          INTEGER, INTENT(IN)           ::  i1, i2, j1, j2, k1, k2
    984985          INTEGER                       ::  wert_ijk(3)
    985986          INTEGER, OPTIONAL             ::  wert1_ijk(3)
    986           REAL                          ::  wert
     987          REAL                          ::  offset, wert
    987988          REAL, OPTIONAL                ::  wert1
    988989          REAL, INTENT(IN)              ::  feld(i1:i2,j1:j2,k1:k2)
  • 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.