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