wrfv2_fire/dyn_em/adapt_timestep_em.F

   1 RECURSIVE 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, ONLY : wrf_dm_maxval, wrf_dm_minval, wrf_dm_mintile_double, wrf_dm_tile_val_int, wrf_dm_maxtile_real
  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   TYPE(WRFU_TimeInterval)                    :: parent_dtInterval
  37   INTEGER                                    :: num_small_steps
  38   integer                                    :: tile
  39   LOGICAL                                    :: stepping_to_bc
  40   INTEGER                                    :: bc_time, output_time
  41   double precision                           :: dt = 0
  42   INTEGER, PARAMETER                         :: precision = 100
  43   INTEGER                                    :: dt_num, dt_den, dt_whole
  44   REAL                                       :: factor
  45   INTEGER                                    :: num, den
  46   TYPE(WRFU_TimeInterval)                    :: last_dtInterval
  47   REAL                                       :: real_time
  48
  49   !
  50   ! If use_last2 is true, this routine will use the time step
  51   !   from 2 steps ago to compute the next time step.  This
  52   !   is used along with step_to_output and step_to_bc
  53
  54   use_last2 = .FALSE.
  55
  56   !
  57   ! If this step has already been adapted, no need to do it again.
  58   ! time step can already be adapted when adaptation_domain is
  59   ! enabled.
  60   !
  61
  62   if (grid%last_step_updated == grid%itimestep) then
  63      return
  64   else
  65      grid%last_step_updated = grid%itimestep
  66   endif
  67
  68   !
  69   ! For nests, set adapt_this_step_using_child to parent's value
  70   !
  71   if (grid%id .ne. 1) then
  72      grid%adapt_this_step_using_child = grid%parents(1)%ptr%adapt_this_step_using_child;
  73   endif
  74
  75   !
  76   ! For nests, if we're not adapting using child nest, we only want to change
  77   !    nests' time steps when the time is conincident with the parent's time.
  78   !    So, if dtbc is not zero, simply return and leave the last time step in
  79   !    place.
  80   !
  81
  82 !  if ((grid%id .ne. 1) .and. (.not. grid%adapt_this_step_using_child)) then
  83 !     if (abs(grid%dtbc) > 0.0001) then
  84 !        return
  85 !     endif
  86 !  endif
  87
  88   last_dtInterval = grid%last_dtInterval
  89
  90   !
  91   ! Get current time
  92   !
  93
  94   tmpTimeInterval = domain_get_current_time ( grid ) - &
  95        domain_get_start_time ( grid )
  96
  97   !
  98   ! Calculate current time in seconds since beginning of model run.
  99   !   Unfortunately, ESMF does not seem to have a way to return
 100   !   floating point seconds based on a TimeInterval.  So, we will
 101   !   calculate it here--but, this is not clean!!
 102   !
 103   curr_secs = real_time(tmpTimeInterval)
 104
 105   !
 106   ! Calculate the maximum allowable increase in the time step given
 107   !   the user input max_step_increase_pct value and the nest ratio.
 108   !
 109   max_increase_factor = 1. + grid%max_step_increase_pct / 100.
 110
 111   !
 112   ! If this is the first time step of the model run (indicated by current_time
 113   !    eq start_time), then set the time step to the input starting_time_step.
 114   !
 115   ! Else, calculate the time step based on cfl.
 116   !
 117   if (domain_get_current_time ( grid ) .eq. domain_get_start_time ( grid )) then
 118      CALL WRFU_TimeIntervalSet(dtInterval, Sn=grid%starting_time_step, Sd=1)
 119      curr_secs = 0
 120      CALL WRFU_TimeIntervalSet(last_dtInterval, Sn=0, Sd=1)
 121
 122   else
 123      if (grid%max_cfl_val < 0.001) then
 124         !
 125         ! If the max_cfl_val is small, then we increase dtInterval the maximum
 126         !    amount allowable.
 127         !
 128         num = INT(max_increase_factor * precision + 0.5)
 129         den = precision
 130         dtInterval = last_dtInterval * num / den
 131
 132      else
 133         !
 134         ! If the max_cfl_val is greater than the user input target cfl, we
 135         !    reduce the time step,
 136         ! else, we increase it.
 137         !
 138         if (grid%max_cfl_val .gt. grid%target_cfl) then
 139            !
 140            ! If we are reducing the time step, we go below target cfl by half the
 141            !   difference between max and target
 142            !   This tends to keep the model more stable.
 143            !
 144
 145            factor = ( (grid%target_cfl - 0.5 * &
 146                 (grid%max_cfl_val - grid%target_cfl) ) / grid%max_cfl_val )
 147            num = INT(factor * precision + 0.5)
 148            den = precision
 149
 150            dtInterval = last_dtInterval * num / den
 151
 152         else
 153            !
 154            ! Linearly increase dtInterval (we'll limit below)
 155            !
 156
 157            factor = grid%target_cfl / grid%max_cfl_val
 158            num = INT(factor * precision + 0.5)
 159            den = precision
 160            dtInterval = last_dtInterval * num / den
 161         endif
 162
 163      endif
 164
 165   endif
 166
 167   ! Limit the increase of dtInterval to the specified input limit
 168
 169   num = NINT( max_increase_factor * precision )
 170   den = precision
 171   tmpTimeInterval = last_dtInterval * num / den
 172   if ( (domain_get_current_time ( grid ) .ne. domain_get_start_time ( grid )) &
 173        .and. (dtInterval .gt. tmpTimeInterval ) ) then
 174      dtInterval = tmpTimeInterval
 175   endif
 176
 177   !
 178   ! Here, we round off dtInterval to nearest 1/100.  This prevents
 179   !    the denominator from getting too large and causing overflow.
 180   !
 181   dt = real_time(dtInterval)
 182   num = NINT(dt * precision)
 183   den = precision
 184   CALL WRFU_TimeIntervalSet(dtInterval, Sn=num, Sd=den)
 185
 186   ! Limit the maximum dtInterval based on user input
 187
 188   CALL WRFU_TimeIntervalSet(tmpTimeInterval, Sn=grid%max_time_step, Sd=1)
 189   if (dtInterval .gt. tmpTimeInterval ) then
 190      dtInterval = tmpTimeInterval
 191   endif
 192
 193   ! Limit the minimum dtInterval based on user input.
 194
 195   CALL WRFU_TimeIntervalSet(tmpTimeInterval, Sn=grid%min_time_step, Sd=1)
 196   if (dtInterval .lt. tmpTimeInterval ) then
 197      dtInterval = tmpTimeInterval
 198   endif
 199
 200   !
 201   ! Now, if this is a nest, and we are adapting based upon parent,
 202   !   we round down the time step to the nearest
 203   !   value that divides evenly into the parent time step.
 204   ! If this is a nest, and we are adapting based upon the child (i.e., the
 205   !   nest), we update the parent timestep to the next smallest multiple
 206   !   timestep.
 207   !
 208   if (grid%nested) then
 209
 210      dt = real_time(dtInterval)
 211
 212      if (.not. grid%adapt_this_step_using_child) then
 213
 214         ! We'll calculate real numbers to get the number of small steps:
 215
 216         num_small_steps = CEILING( grid%parents(1)%ptr%dt / dt )
 217
 218 #ifdef DM_PARALLEL
 219         call wrf_dm_maxval(num_small_steps, idex, jdex)
 220 #endif
 221         dtInterval = domain_get_time_step(grid%parents(1)%ptr) / &
 222              num_small_steps
 223      else
 224
 225         num_small_steps = FLOOR( grid%parents(1)%ptr%dt / dt )
 226
 227 #ifdef DM_PARALLEL
 228         call wrf_dm_minval(num_small_steps, idex, jdex)
 229 #endif
 230         if (num_small_steps < 1) then
 231            num_small_steps = 1
 232         endif
 233
 234      endif
 235   endif
 236
 237
 238   !
 239   ! Setup the values for several variables from the tile with the
 240   !   minimum dt.
 241   !
 242   dt = real_time(dtInterval)
 243
 244 #ifdef DM_PARALLEL
 245   call wrf_dm_mintile_double(dt, tile)
 246   CALL WRFU_TimeIntervalGet(dtInterval,Sn=dt_num,Sd=dt_den,S=dt_whole)
 247   call wrf_dm_tile_val_int(dt_num, tile)
 248   call wrf_dm_tile_val_int(dt_den, tile)
 249   call wrf_dm_tile_val_int(dt_whole, tile)
 250   CALL WRFU_TimeIntervalSet(dtInterval, Sn = dt_whole*dt_den + dt_num, Sd = dt_den)
 251
 252   call wrf_dm_maxtile_real(grid%max_cfl_val, tile)
 253   call wrf_dm_maxtile_real(grid%max_vert_cfl, tile)
 254   call wrf_dm_maxtile_real(grid%max_horiz_cfl, tile)
 255 #endif
 256
 257   if ((grid%nested) .and. (grid%adapt_this_step_using_child)) then
 258
 259      grid%dt = real_time(dtInterval)
 260
 261      ! Set parent step here.
 262      grid%parents(1)%ptr%dt = grid%dt * num_small_steps
 263      parent_dtInterval = dtInterval * num_small_steps
 264
 265      !
 266      ! Update the parent clock based on the new time step
 267      !
 268      CALL WRFU_ClockSet ( grid%parents(1)%ptr%domain_clock,        &
 269           timeStep=parent_dtInterval, &
 270           rc=rc )
 271
 272   endif
 273
 274
 275   !
 276   ! Assure that we fall on a BC time.  Due to a bug in WRF, the time
 277   !   step must fall on the boundary times.  Only modify the dtInterval
 278   !   when this is not the first time step on this domain.
 279   !
 280
 281   grid%stepping_to_time = .FALSE.
 282   time_to_bc = grid%interval_seconds - grid%dtbc
 283   num = INT(time_to_bc * precision + 0.5)
 284   den = precision
 285   CALL WRFU_TimeIntervalSet(tmpTimeInterval, Sn=num, Sd=den)
 286
 287   if ( ( tmpTimeInterval .LT. dtInterval * 2 ) .and. &
 288        ( tmpTimeInterval .GT. dtInterval ) ) then
 289      dtInterval = tmpTimeInterval / 2
 290
 291      use_last2 = .TRUE.
 292      stepping_to_bc = .true.
 293      grid%stepping_to_time = .TRUE.
 294
 295   elseif (tmpTimeInterval .LE. dtInterval) then
 296
 297      bc_time = NINT ( (curr_secs + time_to_bc) / ( grid%interval_seconds ) ) &
 298           * ( grid%interval_seconds )
 299      CALL WRFU_TimeIntervalSet(tmpTimeInterval, S=bc_time)
 300      dtInterval = tmpTimeInterval - &
 301           (domain_get_current_time(grid) - domain_get_start_time(grid))
 302
 303      use_last2 = .TRUE.
 304      stepping_to_bc = .true.
 305      grid%stepping_to_time = .TRUE.
 306   else
 307      stepping_to_bc = .false.
 308   endif
 309
 310   !
 311   ! If the user has requested that we step to output, then
 312   !   assure that we fall on an output time.  We look out two time steps to
 313   !   avoid having a very short time step.  Very short time steps can cause model
 314   !   instability.
 315   !
 316
 317   if ((grid%step_to_output_time) .and. (.not. stepping_to_bc) .and. &
 318        (.not. grid%nested)) then
 319
 320      time_to_output = grid%history_interval*60 - &
 321           mod(curr_secs, grid%history_interval*60.0)
 322      num = INT(time_to_output * precision + 0.5)
 323      den = precision
 324      call WRFU_TimeIntervalSet(tmpTimeInterval, Sn=num, Sd=den)
 325
 326      if ( ( tmpTimeInterval .LT. dtInterval * 2 ) .and. &
 327           ( tmpTimeInterval .GT. dtInterval ) ) then
 328         dtInterval = tmpTimeInterval / 2
 329         use_last2 = .TRUE.
 330         grid%stepping_to_time = .TRUE.
 331
 332      elseif (tmpTimeInterval .LE. dtInterval) then
 333         !
 334         ! We will do some tricks here to assure that we fall exactly on an
 335         !   output time.  Without the tricks, round-off error causes problems!
 336         !
 337
 338         !
 339         ! Calculate output time.  We round to nearest history time to assure
 340         !    we don't have any rounding error.
 341         !
 342         output_time = NINT ( (curr_secs + time_to_output) /  &
 343              ( grid%history_interval * 60 ) ) * (grid%history_interval * 60)
 344         CALL WRFU_TimeIntervalSet(tmpTimeInterval, S=output_time)
 345         dtInterval = tmpTimeInterval - &
 346              (domain_get_current_time(grid) - domain_get_start_time(grid))
 347
 348         use_last2 = .TRUE.
 349         grid%stepping_to_time = .TRUE.
 350      endif
 351   endif
 352
 353   !
 354   ! Now, set adapt_this_step_using_child only if we are not stepping to an
 355   !   output time, or, it's not the start of the model run.
 356   ! Note:  adapt_this_step_using_child is updated just before recursive call to
 357   !    adapt_timestep--see end of this function.
 358   !
 359
 360   if (grid%id == 1) then
 361      if ((grid%adaptation_domain > 1) .and. &
 362           (grid%max_dom == 2) .and. &
 363           (.not. grid%stepping_to_time) .and. &
 364           (domain_get_current_time(grid) .ne. &
 365           domain_get_start_time(grid)) &
 366           ) then
 367
 368         grid%adapt_this_step_using_child = .TRUE.
 369      else
 370         grid%adapt_this_step_using_child = .FALSE.
 371      endif
 372   endif
 373
 374
 375   if (use_last2) then
 376      grid%last_dtInterval = last_dtInterval
 377   else
 378      grid%last_dtInterval = dtInterval
 379   endif
 380
 381   grid%dt = real_time(dtInterval)
 382
 383   grid%last_max_vert_cfl = grid%max_vert_cfl
 384
 385   !
 386   ! Update the clock based on the new time step
 387   !
 388   CALL WRFU_ClockSet ( grid%domain_clock,        &
 389        timeStep=dtInterval, &
 390        rc=rc )
 391
 392   !
 393   ! If we're are adapting based on the child time step,
 394   !    we call the child from here.  This assures that
 395   !    child and parent are updated in sync.
 396   ! Note: This is not necessary when we are adapting based
 397   !    upon parent.
 398   !
 399   if ((grid%id == 1) .and. (grid%adapt_this_step_using_child)) then
 400      !
 401      ! Finally, check if we can adapt using child.  If we are
 402      !   stepping to an output time, we cannot adapt based upon
 403      !   child.  So, we reset the variable before calling the child.
 404      ! This covers the case that, within this parent time-step that
 405      !   we just calculated, we are stepping to an output time.
 406      !
 407      if (grid%stepping_to_time) then
 408         grid%adapt_this_step_using_child = .FALSE.
 409      endif
 410      call adapt_timestep(grid%nests(1)%ptr, config_flags)
 411   endif
 412
 413   !
 414   ! Lateral boundary weight recomputation based on time step.
 415   !
 416   if (grid%id == 1) then
 417      CALL lbc_fcx_gcx ( grid%fcx , grid%gcx , grid%spec_bdy_width , &
 418           grid%spec_zone , grid%relax_zone , grid%dt , config_flags%spec_exp , &
 419           config_flags%specified , config_flags%nested )
 420   endif
 421
 422 END SUBROUTINE adapt_timestep
 423
 424 FUNCTION real_time( timeinterval ) RESULT ( out_time )
 425
 426   USE module_domain
 427
 428   IMPLICIT NONE
 429
 430 ! This function returns a floating point time from an input time interval
 431 !
 432 ! Unfortunately, the ESMF did not provide this functionality.
 433 !
 434 ! Be careful with the output because, due to rounding, the time is only
 435 !   approximate.
 436 !
 437 ! Todd Hutchinson, WSI
 438 ! 4/17/2007
 439
 440 ! !RETURN VALUE:
 441       REAL :: out_time
 442       INTEGER :: dt_num, dt_den, dt_whole
 443
 444 ! !ARGUMENTS:
 445       TYPE(WRFU_TimeInterval), intent(INOUT) :: timeinterval
 446
 447       CALL WRFU_TimeIntervalGet(timeinterval,Sn=dt_num,Sd=dt_den,S=dt_whole)
 448       if (ABS(dt_den) < 1) then
 449          out_time = dt_whole
 450       else
 451          out_time = dt_whole + dt_num / REAL(dt_den)
 452       endif
 453 END FUNCTION
 454
 455
 456 FUNCTION real_time_r8( timeinterval ) RESULT ( out_time )
 457
 458   USE module_domain
 459
 460   IMPLICIT NONE
 461
 462 ! This function returns a double precision floating point time from an input time interval
 463 !
 464 ! Unfortunately, the ESMF did not provide this functionality.
 465 !
 466 ! Be careful with the output because, due to rounding, the time is only
 467 !   approximate.
 468 !
 469 ! Todd Hutchinson, WSI 4/17/2007
 470 ! Converted to r8, William.Gustafson@pnl.gov; 8-May-2008
 471
 472 ! !RETURN VALUE:
 473       REAL(KIND=8) :: out_time
 474       INTEGER(KIND=8) :: dt_whole
 475       INTEGER :: dt_num, dt_den
 476
 477 ! !ARGUMENTS:
 478       TYPE(WRFU_TimeInterval), intent(INOUT) :: timeinterval
 479
 480       CALL WRFU_TimeIntervalGet(timeinterval,Sn=dt_num,Sd=dt_den,S_i8=dt_whole)
 481       if (ABS(dt_den) < 1) then
 482          out_time = REAL(dt_whole)
 483       else
 484          out_time = REAL(dt_whole) + REAL(dt_num,8)/REAL(dt_den,8)
 485       endif
 486 END FUNCTION real_time_r8