Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

adapt_timestep_em.F @ 409

Last change on this file since 409 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
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_h3600 + grid%history_interval_d86400
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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format