source: palm/trunk/SOURCE/timestep.f90 @ 1484

Last change on this file since 1484 was 1484, checked in by kanani, 9 years ago

New:
---
Subroutine init_plant_canopy added to module plant_canopy_model_mod. (plant_canopy_model)
Alternative method for lad-profile construction added, also, new parameters added.
(header, package_parin, plant_canopy_model, read_var_list, write_var_list)
plant_canopy_model-dependency added to several subroutines. (Makefile)
New package/namelist canopy_par for canopy-related parameters added. (package_parin)

Changed:
---
Code structure of the plant canopy model changed, all canopy-model related code
combined to module plant_canopy_model_mod. (check_parameters, init_3d_model,
modules, timestep)
Module plant_canopy_model_mod added in USE-lists of some subroutines. (check_parameters,
header, init_3d_model, package_parin, read_var_list, user_init_plant_canopy, write_var_list)
Canopy initialization moved to new subroutine init_plant_canopy. (check_parameters,
init_3d_model, plant_canopy_model)
Calculation of canopy timestep-criterion removed, instead, the canopy
drag is now directly limited in the calculation of the canopy tendency terms.
(plant_canopy_model, timestep)
Some parameters renamed. (check_parameters, header, init_plant_canopy,
plant_canopy_model, read_var_list, write_var_list)
Unnecessary 3d-arrays removed. (init_plant_canopy, plant_canopy_model, user_init_plant_canopy)
Parameter checks regarding canopy initialization added. (check_parameters)
All canopy steering parameters moved from namelist inipar to canopy_par. (package_parin, parin)
Some redundant MPI communication removed. (init_plant_canopy)

Bugfix:
---
Missing KIND-attribute for REAL constant added. (check_parameters)
DO-WHILE-loop for lad-profile output restricted. (header)
Removed double-listing of use_upstream_for_tke in ONLY-list of module
control_parameters. (prognostic_equations)

  • Property svn:keywords set to Id
File size: 17.3 KB
Line 
1 SUBROUTINE timestep
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! Changes due to new module structure of the plant canopy model:
23!   calculations and parameters related to the plant canopy model removed
24!   (the limitation of the canopy drag, i.e. that the canopy drag itself should
25!   not change the sign of the velocity components, is now assured for in the
26!   calculation of the canopy tendency terms in subroutine plant_canopy_model)
27!
28! Former revisions:
29! -----------------
30! $Id: timestep.f90 1484 2014-10-21 10:53:05Z kanani $
31!
32! 1342 2014-03-26 17:04:47Z kanani
33! REAL constants defined as wp-kind
34!
35! 1322 2014-03-20 16:38:49Z raasch
36! REAL functions provided with KIND-attribute
37!
38! 1320 2014-03-20 08:40:49Z raasch
39! ONLY-attribute added to USE-statements,
40! kind-parameters added to all INTEGER and REAL declaration statements,
41! kinds are defined in new module kinds,
42! old module precision_kind is removed,
43! revision history before 2012 removed,
44! comment fields (!:) to be used for variable explanations added to
45! all variable declaration statements
46!
47! 1257 2013-11-08 15:18:40Z raasch
48! openacc porting
49! bugfix for calculation of advective timestep in case of vertically stretched
50! grids
51!
52! 1092 2013-02-02 11:24:22Z raasch
53! unused variables removed
54!
55! 1053 2012-11-13 17:11:03Z hoffmann
56! timestep is reduced in two-moment cloud scheme according to the maximum
57! terminal velocity of rain drops
58!
59! 1036 2012-10-22 13:43:42Z raasch
60! code put under GPL (PALM 3.9)
61!
62! 1001 2012-09-13 14:08:46Z raasch
63! all actions concerning leapfrog scheme removed
64!
65! 978 2012-08-09 08:28:32Z fricke
66! restriction of the outflow damping layer in the diffusion criterion removed
67!
68! 866 2012-03-28 06:44:41Z raasch
69! bugfix for timestep calculation in case of Galilei transformation,
70! special treatment in case of mirror velocity boundary condition removed
71!
72! Revision 1.1  1997/08/11 06:26:19  raasch
73! Initial revision
74!
75!
76! Description:
77! ------------
78! Compute the time step under consideration of the FCL and diffusion criterion.
79!------------------------------------------------------------------------------!
80
81    USE arrays_3d,                                                             &
82        ONLY:  dzu, dzw, kh, km, u, v, w
83
84    USE cloud_parameters,                                                      &
85        ONLY:  dt_precipitation
86
87    USE control_parameters,                                                    &
88        ONLY:  cfl_factor, coupling_mode, dt_3d, dt_fixed, dt_max,             &
89               galilei_transformation, old_dt, message_string,                 &
90               stop_dt, terminate_coupled, terminate_coupled_remote,           &
91               timestep_reason, u_gtrans, use_ug_for_galilei_tr, v_gtrans
92
93    USE cpulog,                                                                &
94        ONLY:  cpu_log, log_point
95
96    USE grid_variables,                                                        &
97        ONLY:  dx, dx2, dy, dy2
98
99    USE indices,                                                               &
100        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
101
102    USE interfaces
103
104    USE kinds
105
106    USE pegrid
107
108    USE statistics,                                                            &
109        ONLY:  flow_statistics_called, hom, u_max, u_max_ijk, v_max, v_max_ijk,&
110               w_max, w_max_ijk
111
112    IMPLICIT NONE
113
114    INTEGER(iwp) ::  i !:
115    INTEGER(iwp) ::  j !:
116    INTEGER(iwp) ::  k !:
117
118    REAL(wp) ::  div               !:
119    REAL(wp) ::  dt_diff           !:
120    REAL(wp) ::  dt_diff_l         !:
121    REAL(wp) ::  dt_u              !:
122    REAL(wp) ::  dt_u_l            !:
123    REAL(wp) ::  dt_v              !:
124    REAL(wp) ::  dt_v_l            !:
125    REAL(wp) ::  dt_w              !:
126    REAL(wp) ::  dt_w_l            !:
127    REAL(wp) ::  u_gtrans_l        !:
128    REAL(wp) ::  u_max_l           !:
129    REAL(wp) ::  u_min_l           !:
130    REAL(wp) ::  value             !:
131    REAL(wp) ::  v_gtrans_l        !:
132    REAL(wp) ::  v_max_l           !:
133    REAL(wp) ::  v_min_l           !:
134    REAL(wp) ::  w_max_l           !:
135    REAL(wp) ::  w_min_l           !:
136 
137    REAL(wp), DIMENSION(2)         ::  uv_gtrans   !:
138    REAL(wp), DIMENSION(2)         ::  uv_gtrans_l !:
139    REAL(wp), DIMENSION(3)         ::  reduce      !:
140    REAL(wp), DIMENSION(3)         ::  reduce_l    !:
141    REAL(wp), DIMENSION(nzb+1:nzt) ::  dxyz2_min   !: 
142
143
144
145    CALL cpu_log( log_point(12), 'calculate_timestep', 'start' )
146
147!
148!-- In case of Galilei-transform not using the geostrophic wind as translation
149!-- velocity, compute the volume-averaged horizontal velocity components, which
150!-- will then be subtracted from the horizontal wind for the time step and
151!-- horizontal advection routines.
152    IF ( galilei_transformation  .AND. .NOT.  use_ug_for_galilei_tr )  THEN
153       IF ( flow_statistics_called )  THEN
154!
155!--       Horizontal averages already existent, just need to average them
156!--       vertically.
157          u_gtrans = 0.0_wp
158          v_gtrans = 0.0_wp
159          DO  k = nzb+1, nzt
160             u_gtrans = u_gtrans + hom(k,1,1,0)
161             v_gtrans = v_gtrans + hom(k,1,2,0)
162          ENDDO
163          u_gtrans = u_gtrans / REAL( nzt - nzb, KIND=wp )
164          v_gtrans = v_gtrans / REAL( nzt - nzb, KIND=wp )
165       ELSE
166!
167!--       Averaging over the entire model domain.
168          u_gtrans_l = 0.0_wp
169          v_gtrans_l = 0.0_wp
170          !$acc parallel present( u, v )
171          DO  i = nxl, nxr
172             DO  j = nys, nyn
173                DO  k = nzb+1, nzt
174                   u_gtrans_l = u_gtrans_l + u(k,j,i)
175                   v_gtrans_l = v_gtrans_l + v(k,j,i)
176                ENDDO
177             ENDDO
178          ENDDO
179          !$acc end parallel
180          uv_gtrans_l(1) = u_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb), KIND=wp )
181          uv_gtrans_l(2) = v_gtrans_l / REAL( (nxr-nxl+1)*(nyn-nys+1)*(nzt-nzb), KIND=wp )
182#if defined( __parallel )
183          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
184          CALL MPI_ALLREDUCE( uv_gtrans_l, uv_gtrans, 2, MPI_REAL, MPI_SUM, &
185                              comm2d, ierr )
186          u_gtrans = uv_gtrans(1) / REAL( numprocs, KIND=wp )
187          v_gtrans = uv_gtrans(2) / REAL( numprocs, KIND=wp )
188#else
189          u_gtrans = uv_gtrans_l(1)
190          v_gtrans = uv_gtrans_l(2)
191#endif
192       ENDIF
193    ENDIF
194
195!
196!-- Determine the maxima of the velocity components, including their
197!-- grid index positions.
198#if defined( __openacc )
199    IF ( dt_fixed )  THEN  ! otherwise do it further below for better cache usage
200       u_max_l = -999999.9_wp
201       u_min_l =  999999.9_wp
202       v_max_l = -999999.9_wp
203       v_min_l =  999999.9_wp
204       w_max_l = -999999.9_wp
205       w_min_l =  999999.9_wp
206       !$acc parallel present( u, v, w )
207       DO  i = nxl, nxr
208          DO  j = nys, nyn
209             DO  k = nzb+1, nzt
210                u_max_l = MAX( u_max_l, u(k,j,i) )
211                u_min_l = MIN( u_min_l, u(k,j,i) )
212                v_max_l = MAX( v_max_l, v(k,j,i) )
213                v_min_l = MIN( v_min_l, v(k,j,i) )
214                w_max_l = MAX( w_max_l, w(k,j,i) )
215                w_min_l = MIN( w_min_l, w(k,j,i) )
216             ENDDO
217          ENDDO
218       ENDDO
219       !$acc end parallel
220#if defined( __parallel )
221       reduce_l(1) = u_max_l
222       reduce_l(2) = v_max_l
223       reduce_l(3) = w_max_l
224       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
225       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MAX, comm2d, ierr )
226       u_max = reduce(1)
227       v_max = reduce(2)
228       w_max = reduce(3)
229       reduce_l(1) = u_min_l
230       reduce_l(2) = v_min_l
231       reduce_l(3) = w_min_l
232       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
233       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
234       IF ( ABS( reduce(1) ) > u_max )  u_max = reduce(1)
235       IF ( ABS( reduce(2) ) > v_max )  v_max = reduce(2)
236       IF ( ABS( reduce(3) ) > w_max )  w_max = reduce(3)
237#else
238       IF ( ABS( u_min_l ) > u_max_l )  THEN
239          u_max = u_min_l
240       ELSE
241          u_max = u_max_l
242       ENDIF
243       IF ( ABS( v_min_l ) > v_max_l )  THEN
244          v_max = v_min_l
245       ELSE
246          v_max = v_max_l
247       ENDIF
248       IF ( ABS( w_min_l ) > w_max_l )  THEN
249          w_max = w_min_l
250       ELSE
251          w_max = w_max_l
252       ENDIF
253#endif
254    ENDIF
255#else
256    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u, 'abs', 0.0_wp, &
257                         u_max, u_max_ijk )
258    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v, 'abs', 0.0_wp, &
259                         v_max, v_max_ijk )
260    CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w, 'abs', 0.0_wp, &
261                         w_max, w_max_ijk )
262#endif
263
264    IF ( .NOT. dt_fixed )  THEN
265#if defined( __openacc )
266!
267!--    Variable time step:
268!--    Calculate the maximum time step according to the CFL-criterion,
269!--    individually for each velocity component
270       dt_u_l  =  999999.9_wp
271       dt_v_l  =  999999.9_wp
272       dt_w_l  =  999999.9_wp
273       u_max_l = -999999.9_wp
274       u_min_l =  999999.9_wp
275       v_max_l = -999999.9_wp
276       v_min_l =  999999.9_wp
277       w_max_l = -999999.9_wp
278       w_min_l =  999999.9_wp
279       !$acc parallel loop collapse(3) present( u, v, w )
280       DO  i = nxl, nxr
281          DO  j = nys, nyn
282             DO  k = nzb+1, nzt
283                dt_u_l  = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) )
284                dt_v_l  = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) )
285                dt_w_l  = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
286                u_max_l = MAX( u_max_l, u(k,j,i) )
287                u_min_l = MIN( u_min_l, u(k,j,i) )
288                v_max_l = MAX( v_max_l, v(k,j,i) )
289                v_min_l = MIN( v_min_l, v(k,j,i) )
290                w_max_l = MAX( w_max_l, w(k,j,i) )
291                w_min_l = MIN( w_min_l, w(k,j,i) )
292             ENDDO
293          ENDDO
294       ENDDO
295       !$acc end parallel
296
297#if defined( __parallel )
298       reduce_l(1) = dt_u_l
299       reduce_l(2) = dt_v_l
300       reduce_l(3) = dt_w_l
301       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
302       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
303       dt_u = reduce(1)
304       dt_v = reduce(2)
305       dt_w = reduce(3)
306
307       reduce_l(1) = u_max_l
308       reduce_l(2) = v_max_l
309       reduce_l(3) = w_max_l
310       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
311       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MAX, comm2d, ierr )
312       u_max = reduce(1)
313       v_max = reduce(2)
314       w_max = reduce(3)
315       reduce_l(1) = u_min_l
316       reduce_l(2) = v_min_l
317       reduce_l(3) = w_min_l
318       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
319       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
320       IF ( ABS( reduce(1) ) > u_max )  u_max = reduce(1)
321       IF ( ABS( reduce(2) ) > v_max )  v_max = reduce(2)
322       IF ( ABS( reduce(3) ) > w_max )  w_max = reduce(3)
323#else
324       dt_u = dt_u_l
325       dt_v = dt_v_l
326       dt_w = dt_w_l
327
328       IF ( ABS( u_min_l ) > u_max_l )  THEN
329          u_max = u_min_l
330       ELSE
331          u_max = u_max_l
332       ENDIF
333       IF ( ABS( v_min_l ) > v_max_l )  THEN
334          v_max = v_min_l
335       ELSE
336          v_max = v_max_l
337       ENDIF
338       IF ( ABS( w_min_l ) > w_max_l )  THEN
339          w_max = w_min_l
340       ELSE
341          w_max = w_max_l
342       ENDIF
343#endif
344
345#else
346!
347!--    Variable time step:
348!--    Calculate the maximum time step according to the CFL-criterion,
349!--    individually for each velocity component
350       dt_u_l = 999999.9_wp
351       dt_v_l = 999999.9_wp
352       dt_w_l = 999999.9_wp
353       DO  i = nxl, nxr
354          DO  j = nys, nyn
355             DO  k = nzb+1, nzt
356                dt_u_l = MIN( dt_u_l, ( dx     / ( ABS( u(k,j,i) - u_gtrans ) + 1.0E-10_wp ) ) )
357                dt_v_l = MIN( dt_v_l, ( dy     / ( ABS( v(k,j,i) - v_gtrans ) + 1.0E-10_wp ) ) )
358                dt_w_l = MIN( dt_w_l, ( dzu(k) / ( ABS( w(k,j,i) )            + 1.0E-10_wp ) ) )
359             ENDDO
360          ENDDO
361       ENDDO
362
363#if defined( __parallel )
364       reduce_l(1) = dt_u_l
365       reduce_l(2) = dt_v_l
366       reduce_l(3) = dt_w_l
367       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
368       CALL MPI_ALLREDUCE( reduce_l, reduce, 3, MPI_REAL, MPI_MIN, comm2d, ierr )
369       dt_u = reduce(1)
370       dt_v = reduce(2)
371       dt_w = reduce(3)
372#else
373       dt_u = dt_u_l
374       dt_v = dt_v_l
375       dt_w = dt_w_l
376#endif
377
378#endif
379
380!
381!--    Compute time step according to the diffusion criterion.
382!--    First calculate minimum grid spacing which only depends on index k
383!--    Note: also at k=nzb+1 a full grid length is being assumed, although
384!--          in the Prandtl-layer friction term only dz/2 is used.
385!--          Experience from the old model seems to justify this.
386       dt_diff_l = 999999.0_wp
387
388       DO  k = nzb+1, nzt
389           dxyz2_min(k) = MIN( dx2, dy2, dzw(k)*dzw(k) ) * 0.125_wp
390       ENDDO
391
392!$OMP PARALLEL private(i,j,k,value) reduction(MIN: dt_diff_l)
393!$OMP DO
394       !$acc parallel loop collapse(3) present( kh, km )
395       DO  i = nxl, nxr
396          DO  j = nys, nyn
397             DO  k = nzb+1, nzt
398                dt_diff_l = MIN( dt_diff_l, dxyz2_min(k) / &
399                                       ( MAX( kh(k,j,i), km(k,j,i) ) + 1E-20_wp ) )
400             ENDDO
401          ENDDO
402       ENDDO
403       !$acc end parallel
404!$OMP END PARALLEL
405#if defined( __parallel )
406       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
407       CALL MPI_ALLREDUCE( dt_diff_l, dt_diff, 1, MPI_REAL, MPI_MIN, comm2d, &
408                           ierr )
409#else
410       dt_diff = dt_diff_l
411#endif
412
413!
414!--    The time step is the minimum of the 3-4 components and the diffusion time
415!--    step minus a reduction (cfl_factor) to be on the safe side.
416!--    The time step must not exceed the maximum allowed value.
417       dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w,   &
418                                 dt_precipitation )
419       dt_3d = MIN( dt_3d, dt_max )
420
421!
422!--    Remember the restricting time step criterion for later output.
423       IF ( MIN( dt_u, dt_v, dt_w ) < dt_diff )  THEN
424          timestep_reason = 'A'
425       ELSE
426          timestep_reason = 'D'
427       ENDIF
428
429!
430!--    Set flag if the time step becomes too small.
431       IF ( dt_3d < ( 0.00001_wp * dt_max ) )  THEN
432          stop_dt = .TRUE.
433
434          WRITE( message_string, * ) 'Time step has reached minimum limit.',  &
435               '&dt              = ', dt_3d, ' s  Simulation is terminated.', &
436               '&old_dt          = ', old_dt, ' s',                           &
437               '&dt_u            = ', dt_u, ' s',                             &
438               '&dt_v            = ', dt_v, ' s',                             &
439               '&dt_w            = ', dt_w, ' s',                             &
440               '&dt_diff         = ', dt_diff, ' s',                          &
441               '&u_max           = ', u_max, ' m/s   k=', u_max_ijk(1),       &
442               '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                    &
443               '&v_max           = ', v_max, ' m/s   k=', v_max_ijk(1),       &
444               '  j=', v_max_ijk(2), '  i=', v_max_ijk(3),                    &
445               '&w_max           = ', w_max, ' m/s   k=', w_max_ijk(1),       &
446               '  j=', w_max_ijk(2), '  i=', w_max_ijk(3)
447          CALL message( 'timestep', 'PA0312', 0, 1, 0, 6, 0 )
448!
449!--       In case of coupled runs inform the remote model of the termination
450!--       and its reason, provided the remote model has not already been
451!--       informed of another termination reason (terminate_coupled > 0) before.
452#if defined( __parallel )
453          IF ( coupling_mode /= 'uncoupled' .AND. terminate_coupled == 0 )  THEN
454             terminate_coupled = 2
455             IF ( myid == 0 ) THEN
456                CALL MPI_SENDRECV( &
457                     terminate_coupled,        1, MPI_INTEGER, target_id,  0, &
458                     terminate_coupled_remote, 1, MPI_INTEGER, target_id,  0, &
459                     comm_inter, status, ierr )
460             ENDIF
461             CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr)
462          ENDIF
463#endif
464       ENDIF
465
466!
467!--    Ensure a smooth value (two significant digits) of the timestep.
468       div = 1000.0_wp
469       DO  WHILE ( dt_3d < div )
470          div = div / 10.0_wp
471       ENDDO
472       dt_3d = NINT( dt_3d * 100.0_wp / div ) * div / 100.0_wp
473
474!
475!--    Adjust the time step
476       old_dt = dt_3d
477
478    ENDIF
479
480    CALL cpu_log( log_point(12), 'calculate_timestep', 'stop' )
481
482 END SUBROUTINE timestep
Note: See TracBrowser for help on using the repository browser.