wrfv2_fire/dyn_em/adapt_timestep_em.F

   1 SUBROUTINE adapt_timestep(grid, config_flags)
   2
   3 !--------------------------------------------------------------------------
   4 !<DESCRIPTION>
   5 !<pre>
   6 !
   7 ! This routine sets the time step based on the cfl condition.  It's used to
   8 !   dynamically adapt the timestep as the model runs.
   9 !
  10 ! T. Hutchinson, WSI
  11 ! March 2007
  12 !
  13 !</pre>
  14 !</DESCRIPTION>
  15 !--------------------------------------------------------------------------
  16
  17 ! Driver layer modules
  18   USE module_domain
  19   USE module_configure
  20   USE module_dm
  21   USE module_bc_em
  22
  23   IMPLICIT NONE
  24
  25   TYPE(domain) , TARGET , INTENT(INOUT)      :: grid
  26   TYPE (grid_config_rec_type) , INTENT(IN)   :: config_flags
  27
  28   LOGICAL                                    :: use_last2
  29   REAL                                       :: curr_secs
  30   REAL                                       :: max_increase_factor
  31   REAL                                       :: time_to_output, &
  32                                                 time_to_bc
  33   INTEGER                                    :: idex=0, jdex=0
  34   INTEGER                                    :: rc
  35   TYPE(WRFU_TimeInterval)                    :: tmpTimeInterval, dtInterval
  36   INTEGER                                    :: num_small_steps
  37   integer                                    :: tile
  38   LOGICAL                                    :: stepping_to_bc
  39   INTEGER                                    :: bc_time, output_time
  40   double precision                           :: dt
  41   INTEGER, PARAMETER                         :: precision = 100
  42   INTEGER                                    :: dt_num, dt_den, dt_whole
  43   REAL                                       :: factor
  44   INTEGER                                    :: num, den
  45   TYPE(WRFU_TimeInterval)                    :: last_dtInterval
  46   REAL                                       :: real_time
  47
  48   !
  49   ! If use_last2 is true, this routine will use the time step
  50   !   from 2 steps ago to compute the next time step.  This
  51   !   is used along with step_to_output and step_to_bc
  52
  53   use_last2 = .FALSE.
  54
  55   !
  56   ! For nests, we only want to change nests' time steps when the time is
  57   !    conincident with the parent's time.  So, if dtbc is not
  58   !    zero, simply return and leave the last time step in place.
  59   !
  60   if (config_flags%nested) then
  61      if (abs(grid%dtbc) > 0.0001) then
  62         return
  63      endif
  64   endif
  65
  66   last_dtInterval = grid%last_dtInterval
  67
  68   !
  69   ! Get current time
  70   !
  71
  72   tmpTimeInterval = domain_get_current_time ( grid ) - &
  73        domain_get_start_time ( grid )
  74
  75   !
  76   ! Calculate current time in seconds since beginning of model run.
  77   !   Unfortunately, ESMF does not seem to have a way to return
  78   !   floating point seconds based on a TimeInterval.  So, we will
  79   !   calculate it here--but, this is not clean!!
  80   !
  81   curr_secs = real_time(tmpTimeInterval)
  82
  83   !
  84   ! Calculate the maximum allowable increase in the time step given
  85   !   the user input max_step_increase_pct value and the nest ratio.
  86   !
  87   if (config_flags%nested) then
  88      max_increase_factor = 1. + &
  89           grid%parent_time_step_ratio * grid%max_step_increase_pct / 100.
  90   else
  91      max_increase_factor = 1. + grid%max_step_increase_pct / 100.
  92   endif
  93
  94   !
  95   ! If this is the first time step of the model run (indicated by current_time
  96   !    eq start_time), then set the time step to the input starting_time_step.
  97   !
  98   ! Else, calculate the time step based on cfl.
  99   !
 100   if (domain_get_current_time ( grid ) .eq. domain_get_start_time ( grid )) then
 101      CALL WRFU_TimeIntervalSet(dtInterval, Sn=grid%starting_time_step, Sd=1)
 102      curr_secs = 0
 103      CALL WRFU_TimeIntervalSet(last_dtInterval, Sn=0, Sd=1)
 104
 105   else
 106      if (grid%max_cfl_val < 0.001) then
 107         !
 108         ! If the max_cfl_val is small, then we increase dtInterval the maximum
 109         !    amount allowable.
 110         !
 111         num = INT(max_increase_factor * precision + 0.5)
 112         den = precision
 113         dtInterval = last_dtInterval * num / den
 114
 115      else
 116         !
 117         ! If the max_cfl_val is greater than the user input target cfl, we
 118         !    reduce the time step,
 119         ! else, we increase it.
 120         !
 121         if (grid%max_cfl_val .gt. grid%target_cfl) then
 122            !
 123            ! If we are reducing the time step, we go below target cfl by half the
 124            !   difference between max and target (or 75% of the last time step,
 125            !   which ever is greater).
 126            !   This tends to keep the model more stable.
 127            !
 128
 129            factor = MAX ( 0.75 , ( (grid%target_cfl - 0.5 * &
 130                 (grid%max_cfl_val - grid%target_cfl) ) / grid%max_cfl_val ) )
 131            num = INT(factor * precision + 0.5)
 132            den = precision
 133
 134            dtInterval = last_dtInterval * num / den
 135
 136         else
 137            !
 138            ! Linearly increase dtInterval (we'll limit below)
 139            !
 140
 141            factor = grid%target_cfl / grid%max_cfl_val
 142            num = INT(factor * precision + 0.5)
 143            den = precision
 144            dtInterval = last_dtInterval * num / den
 145         endif
 146
 147      endif
 148
 149   endif
 150
 151   ! Limit the increase of dtInterval to the specified input limit
 152
 153   num = NINT( max_increase_factor * precision )
 154   den = precision
 155   tmpTimeInterval = last_dtInterval * num / den
 156   if ( (domain_get_current_time ( grid ) .ne. domain_get_start_time ( grid )) &
 157        .and. (dtInterval .gt. tmpTimeInterval ) ) then
 158      dtInterval = tmpTimeInterval
 159   endif
 160
 161   !
 162   ! Here, we round off dtInterval to nearest 1/100.  This prevents
 163   !    the denominator from getting too large and causing overflow.
 164   !
 165   dt = real_time(dtInterval)
 166   num = NINT(dt * precision)
 167   den = precision
 168   CALL WRFU_TimeIntervalSet(dtInterval, Sn=num, Sd=den)
 169
 170   ! Limit the maximum dtInterval based on user input
 171
 172   CALL WRFU_TimeIntervalSet(tmpTimeInterval, Sn=grid%max_time_step, Sd=1)
 173   if (dtInterval .gt. tmpTimeInterval ) then
 174      dtInterval = tmpTimeInterval
 175   endif
 176
 177   ! Limit the minimum dtInterval based on user input.
 178
 179   CALL WRFU_TimeIntervalSet(tmpTimeInterval, Sn=grid%min_time_step, Sd=1)
 180   if (dtInterval .lt. tmpTimeInterval ) then
 181      dtInterval = tmpTimeInterval
 182   endif
 183
 184   !
 185   ! Now, if this is a nest, we round down the time step to the nearest
 186   !   value that divides evenly into the parent time step.
 187   !
 188   if (config_flags%nested) then
 189      ! We'll calculate real numbers to get the number of small steps:
 190
 191      dt = real_time(dtInterval)
 192
 193      num_small_steps = CEILING( grid%parents(1)%ptr%dt / dt )
 194
 195 #ifdef DM_PARALLEL
 196      call wrf_dm_maxval(num_small_steps, idex, jdex)
 197 #endif
 198      dtInterval = domain_get_time_step(grid%parents(1)%ptr) / &
 199           num_small_steps
 200   endif
 201
 202   !
 203   ! Setup the values for several variables from the tile with the
 204   !   minimum dt.
 205   !
 206   dt = real_time(dtInterval)
 207
 208 #ifdef DM_PARALLEL
 209   call wrf_dm_mintile_double(dt, tile)
 210   CALL WRFU_TimeIntervalGet(dtInterval,Sn=dt_num,Sd=dt_den,S=dt_whole)
 211   call wrf_dm_tile_val_int(dt_num, tile)
 212   call wrf_dm_tile_val_int(dt_den, tile)
 213   call wrf_dm_tile_val_int(dt_whole, tile)
 214   CALL WRFU_TimeIntervalSet(dtInterval, Sn = dt_whole*dt_den + dt_num, Sd = dt_den)
 215
 216   call wrf_dm_maxtile_real(grid%max_cfl_val, tile)
 217   call wrf_dm_maxtile_real(grid%max_vert_cfl, tile)
 218   call wrf_dm_maxtile_real(grid%max_horiz_cfl, tile)
 219 #endif
 220
 221   !
 222   ! Assure that we fall on a BC time.  Due to a bug in WRF, the time
 223   !   step must fall on the boundary times.
 224   !
 225
 226   time_to_bc = grid%interval_seconds - grid%dtbc
 227   num = INT(time_to_bc * precision + 0.5)
 228   den = precision
 229   CALL WRFU_TimeIntervalSet(tmpTimeInterval, Sn=num, Sd=den)
 230
 231   if ( ( tmpTimeInterval .LT. dtInterval * 2 ) .and. &
 232        ( tmpTimeInterval .GT. dtInterval ) ) then
 233      dtInterval = tmpTimeInterval / 2
 234
 235      use_last2 = .TRUE.
 236      stepping_to_bc = .true.
 237
 238   elseif (tmpTimeInterval .LE. dtInterval) then
 239
 240      bc_time = NINT ( (curr_secs + time_to_bc) / ( grid%interval_seconds ) ) &
 241           * ( grid%interval_seconds )
 242      CALL WRFU_TimeIntervalSet(tmpTimeInterval, S=bc_time)
 243      dtInterval = tmpTimeInterval - &
 244           (domain_get_current_time(grid) - domain_get_start_time(grid))
 245
 246      use_last2 = .TRUE.
 247      stepping_to_bc = .true.
 248   else
 249      stepping_to_bc = .false.
 250   endif
 251
 252   !
 253   ! If the user has requested that we step to output, then
 254   !   assure that we fall on an output time.  We look out two time steps to
 255   !   avoid having a very short time step.  Very short time steps can cause model
 256   !   instability.
 257   !
 258
 259   if ((grid%step_to_output_time) .and. (.not. stepping_to_bc) .and. &
 260        (.not. config_flags%nested)) then
 261
 262      time_to_output = grid%history_interval*60 - &
 263           mod(curr_secs, grid%history_interval*60.0)
 264      num = INT(time_to_output * precision + 0.5)
 265      den = precision
 266      call WRFU_TimeIntervalSet(tmpTimeInterval, Sn=num, Sd=den)
 267
 268      if ( ( tmpTimeInterval .LT. dtInterval * 2 ) .and. &
 269           ( tmpTimeInterval .GT. dtInterval ) ) then
 270         dtInterval = tmpTimeInterval / 2
 271         use_last2 = .TRUE.
 272
 273      elseif (tmpTimeInterval .LE. dtInterval) then
 274         !
 275         ! We will do some tricks here to assure that we fall exactly on an
 276         !   output time.  Without the tricks, round-off error causes problems!
 277         !
 278
 279         !
 280         ! Calculate output time.  We round to nearest history time to assure
 281         !    we don't have any rounding error.
 282         !
 283         output_time = NINT ( (curr_secs + time_to_output) /  &
 284              ( grid%history_interval * 60 ) ) * (grid%history_interval * 60)
 285         CALL WRFU_TimeIntervalSet(tmpTimeInterval, S=output_time)
 286         dtInterval = tmpTimeInterval - &
 287              (domain_get_current_time(grid) - domain_get_start_time(grid))
 288
 289         use_last2 = .TRUE.
 290      endif
 291   endif
 292
 293   if (use_last2) then
 294      grid%last_dtInterval = last_dtInterval
 295   else
 296      grid%last_dtInterval = dtInterval
 297   endif
 298
 299   grid%dt = real_time(dtInterval)
 300
 301   grid%last_max_vert_cfl = grid%max_vert_cfl
 302
 303   !
 304   ! Update the clock based on the new time step
 305   !
 306   CALL WRFU_ClockSet ( grid%domain_clock,        &
 307        timeStep=dtInterval, &
 308        rc=rc )
 309
 310   !
 311   ! Lateral boundary weight recomputation based on time step.
 312   !
 313   CALL lbc_fcx_gcx ( grid%fcx , grid%gcx , grid%spec_bdy_width , &
 314                      grid%spec_zone , grid%relax_zone , grid%dt , config_flags%spec_exp , &
 315                      config_flags%specified , config_flags%nested )
 316
 317 END SUBROUTINE adapt_timestep
 318
 319 FUNCTION real_time( timeinterval ) RESULT ( out_time )
 320
 321   USE module_domain
 322
 323   IMPLICIT NONE
 324
 325 ! This function returns a floating point time from an input time interval
 326 !
 327 ! Unfortunately, the ESMF did not provide this functionality.
 328 !
 329 ! Be careful with the output because, due to rounding, the time is only
 330 !   approximate.
 331 !
 332 ! Todd Hutchinson, WSI
 333 ! 4/17/2007
 334
 335 ! !RETURN VALUE:
 336       REAL :: out_time
 337       INTEGER :: dt_num, dt_den, dt_whole
 338
 339 ! !ARGUMENTS:
 340       TYPE(WRFU_TimeInterval), intent(INOUT) :: timeinterval
 341
 342       CALL WRFU_TimeIntervalGet(timeinterval,Sn=dt_num,Sd=dt_den,S=dt_whole)
 343       if (ABS(dt_den) < 1) then
 344          out_time = dt_whole
 345       else
 346          out_time = dt_whole + dt_num / REAL(dt_den)
 347       endif
 348 END FUNCTION