source: lmdz_wrf/trunk/WRFV3/dyn_em/adapt_timestep_em.F @ 1531

Last change on this file since 1531 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 17.1 KB
Line 
1RECURSIVE 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
423END SUBROUTINE adapt_timestep
424
425SUBROUTINE 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
480END SUBROUTINE calc_dt
481
482
483FUNCTION 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
512END FUNCTION
513
514FUNCTION 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
544END FUNCTION real_time_r8
Note: See TracBrowser for help on using the repository browser.