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