source: palm/trunk/SOURCE/dynamics_mod.f90 @ 4691

Last change on this file since 4691 was 4627, checked in by raasch, 4 years ago

bugfix for r4646

  • Property svn:keywords set to Id
File size: 68.0 KB
Line 
1!> @file dynamics_mod.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2020 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: dynamics_mod.f90 4627 2020-07-26 10:14:44Z suehring $
26! bugfix for r4626
27!
28! 4626 2020-07-26 09:49:48Z raasch
29! file re-formatted to follow the PALM coding standard
30!
31! 4517 2020-05-03 14:29:30Z raasch
32! added restart with MPI-IO for reading local arrays
33!
34! 4505 2020-04-20 15:37:15Z schwenkel
35! Add flag for saturation check
36!
37! 4495 2020-04-13 20:11:20Z resler
38! restart data handling with MPI-IO added
39!
40! 4360 2020-01-07 11:25:50Z suehring
41! Bugfix for last commit.
42!
43! 4359 2019-12-30 13:36:50Z suehring
44! Refine post-initialization check for realistically inital values of mixing ratio. Give an error
45! message for faulty initial values, but only a warning in a restart run.
46!
47! 4347 2019-12-18 13:18:33Z suehring
48! Implement post-initialization check for realistically inital values of mixing ratio
49!
50! 4281 2019-10-29 15:15:39Z schwenkel
51! Moved boundary conditions in dynamics module
52!
53! 4097 2019-07-15 11:59:11Z suehring
54! Avoid overlong lines - limit is 132 characters per line
55!
56! 4047 2019-06-21 18:58:09Z knoop
57! Initial introduction of the dynamics module with only dynamics_swap_timelevel implemented
58!
59!
60! Description:
61! ------------
62!> This module contains the dynamics of PALM.
63!--------------------------------------------------------------------------------------------------!
64 MODULE dynamics_mod
65
66
67    USE arrays_3d,                                                                                 &
68        ONLY:  c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l,                      &
69               dzu,                                                                                &
70               exner,                                                                              &
71               hyp,                                                                                &
72               pt, pt_1, pt_2, pt_init, pt_p,                                                      &
73               q, q_1, q_2, q_p,                                                                   &
74               s, s_1, s_2, s_p,                                                                   &
75               u, u_1, u_2, u_init, u_p, u_m_l, u_m_n, u_m_r, u_m_s,                               &
76               v, v_1, v_2, v_p, v_init, v_m_l, v_m_n, v_m_r, v_m_s,                               &
77               w, w_1, w_2, w_p, w_m_l, w_m_n, w_m_r, w_m_s
78
79    USE basic_constants_and_equations_mod,                                                         &
80        ONLY:  magnus,                                                                             &
81               rd_d_rv
82
83    USE control_parameters,                                                                        &
84        ONLY:  bc_dirichlet_l,                                                                     &
85               bc_dirichlet_s,                                                                     &
86               bc_radiation_l,                                                                     &
87               bc_radiation_n,                                                                     &
88               bc_radiation_r,                                                                     &
89               bc_radiation_s,                                                                     &
90               bc_pt_t_val,                                                                        &
91               bc_q_t_val,                                                                         &
92               bc_s_t_val,                                                                         &
93               check_realistic_q,                                                                  &
94               child_domain,                                                                       &
95               coupling_mode,                                                                      &
96               dt_3d,                                                                              &
97               humidity,                                                                           &
98               ibc_pt_b,                                                                           &
99               ibc_pt_t,                                                                           &
100               ibc_q_b,                                                                            &
101               ibc_q_t,                                                                            &
102               ibc_s_b,                                                                            &
103               ibc_s_t,                                                                            &
104               ibc_uv_b,                                                                           &
105               ibc_uv_t,                                                                           &
106               initializing_actions,                                                               &
107               intermediate_timestep_count,                                                        &
108               length,                                                                             &
109               message_string,                                                                     &
110               nesting_offline,                                                                    &
111               neutral,                                                                            &
112               nudging,                                                                            &
113               passive_scalar,                                                                     &
114               restart_string,                                                                     &
115               tsc,                                                                                &
116               use_cmax
117
118    USE grid_variables,                                                                            &
119        ONLY:  ddx,                                                                                &
120               ddy,                                                                                &
121               dx,                                                                                 &
122               dy
123
124    USE indices,                                                                                   &
125        ONLY:  nbgp,                                                                               &
126               nx,                                                                                 &
127               nxl,                                                                                &
128               nxlg,                                                                               &
129               nxr,                                                                                &
130               nxrg,                                                                               &
131               ny,                                                                                 &
132               nys,                                                                                &
133               nysg,                                                                               &
134               nyn,                                                                                &
135               nyng,                                                                               &
136               nzb,                                                                                &
137               nzt
138
139    USE kinds
140
141    USE pegrid
142
143    USE pmc_interface,                                                                             &
144        ONLY : nesting_mode
145
146!    USE restart_data_mpi_io_mod,                                                                   &
147!        ONLY:
148
149    USE surface_mod,                                                                               &
150        ONLY :  bc_h
151
152
153    IMPLICIT NONE
154
155    LOGICAL ::  dynamics_module_enabled = .FALSE.   !<
156
157    SAVE
158
159    PRIVATE
160
161!
162!-- Public functions
163    PUBLIC                                                                                         &
164       dynamics_parin,                                                                             &
165       dynamics_check_parameters,                                                                  &
166       dynamics_check_data_output_ts,                                                              &
167       dynamics_check_data_output_pr,                                                              &
168       dynamics_check_data_output,                                                                 &
169       dynamics_init_masks,                                                                        &
170       dynamics_define_netcdf_grid,                                                                &
171       dynamics_init_arrays,                                                                       &
172       dynamics_init,                                                                              &
173       dynamics_init_checks,                                                                       &
174       dynamics_header,                                                                            &
175       dynamics_actions,                                                                           &
176       dynamics_non_advective_processes,                                                           &
177       dynamics_exchange_horiz,                                                                    &
178       dynamics_prognostic_equations,                                                              &
179       dynamics_boundary_conditions,                                                               &
180       dynamics_swap_timelevel,                                                                    &
181       dynamics_3d_data_averaging,                                                                 &
182       dynamics_data_output_2d,                                                                    &
183       dynamics_data_output_3d,                                                                    &
184       dynamics_statistics,                                                                        &
185       dynamics_rrd_global,                                                                        &
186       dynamics_rrd_local,                                                                         &
187       dynamics_wrd_global,                                                                        &
188       dynamics_wrd_local,                                                                         &
189       dynamics_last_actions
190
191!
192!-- Public parameters, constants and initial values
193    PUBLIC                                                                                         &
194       dynamics_module_enabled
195
196    INTERFACE dynamics_parin
197       MODULE PROCEDURE dynamics_parin
198    END INTERFACE dynamics_parin
199
200    INTERFACE dynamics_check_parameters
201       MODULE PROCEDURE dynamics_check_parameters
202    END INTERFACE dynamics_check_parameters
203
204    INTERFACE dynamics_check_data_output_ts
205       MODULE PROCEDURE dynamics_check_data_output_ts
206    END INTERFACE dynamics_check_data_output_ts
207
208    INTERFACE dynamics_check_data_output_pr
209       MODULE PROCEDURE dynamics_check_data_output_pr
210    END INTERFACE dynamics_check_data_output_pr
211
212    INTERFACE dynamics_check_data_output
213       MODULE PROCEDURE dynamics_check_data_output
214    END INTERFACE dynamics_check_data_output
215
216    INTERFACE dynamics_init_masks
217       MODULE PROCEDURE dynamics_init_masks
218    END INTERFACE dynamics_init_masks
219
220    INTERFACE dynamics_define_netcdf_grid
221       MODULE PROCEDURE dynamics_define_netcdf_grid
222    END INTERFACE dynamics_define_netcdf_grid
223
224    INTERFACE dynamics_init_arrays
225       MODULE PROCEDURE dynamics_init_arrays
226    END INTERFACE dynamics_init_arrays
227
228    INTERFACE dynamics_init
229       MODULE PROCEDURE dynamics_init
230    END INTERFACE dynamics_init
231
232    INTERFACE dynamics_init_checks
233       MODULE PROCEDURE dynamics_init_checks
234    END INTERFACE dynamics_init_checks
235
236    INTERFACE dynamics_header
237       MODULE PROCEDURE dynamics_header
238    END INTERFACE dynamics_header
239
240    INTERFACE dynamics_actions
241       MODULE PROCEDURE dynamics_actions
242       MODULE PROCEDURE dynamics_actions_ij
243    END INTERFACE dynamics_actions
244
245    INTERFACE dynamics_non_advective_processes
246       MODULE PROCEDURE dynamics_non_advective_processes
247       MODULE PROCEDURE dynamics_non_advective_processes_ij
248    END INTERFACE dynamics_non_advective_processes
249
250    INTERFACE dynamics_exchange_horiz
251       MODULE PROCEDURE dynamics_exchange_horiz
252    END INTERFACE dynamics_exchange_horiz
253
254    INTERFACE dynamics_prognostic_equations
255       MODULE PROCEDURE dynamics_prognostic_equations
256       MODULE PROCEDURE dynamics_prognostic_equations_ij
257    END INTERFACE dynamics_prognostic_equations
258
259    INTERFACE dynamics_boundary_conditions
260       MODULE PROCEDURE dynamics_boundary_conditions
261    END INTERFACE dynamics_boundary_conditions
262
263    INTERFACE dynamics_swap_timelevel
264       MODULE PROCEDURE dynamics_swap_timelevel
265    END INTERFACE dynamics_swap_timelevel
266
267    INTERFACE dynamics_3d_data_averaging
268       MODULE PROCEDURE dynamics_3d_data_averaging
269    END INTERFACE dynamics_3d_data_averaging
270
271    INTERFACE dynamics_data_output_2d
272       MODULE PROCEDURE dynamics_data_output_2d
273    END INTERFACE dynamics_data_output_2d
274
275    INTERFACE dynamics_data_output_3d
276       MODULE PROCEDURE dynamics_data_output_3d
277    END INTERFACE dynamics_data_output_3d
278
279    INTERFACE dynamics_statistics
280       MODULE PROCEDURE dynamics_statistics
281    END INTERFACE dynamics_statistics
282
283    INTERFACE dynamics_rrd_global
284       MODULE PROCEDURE dynamics_rrd_global_ftn
285       MODULE PROCEDURE dynamics_rrd_global_mpi
286    END INTERFACE dynamics_rrd_global
287
288    INTERFACE dynamics_rrd_local
289       MODULE PROCEDURE dynamics_rrd_local_ftn
290       MODULE PROCEDURE dynamics_rrd_local_mpi
291    END INTERFACE dynamics_rrd_local
292
293    INTERFACE dynamics_wrd_global
294       MODULE PROCEDURE dynamics_wrd_global
295    END INTERFACE dynamics_wrd_global
296
297    INTERFACE dynamics_wrd_local
298       MODULE PROCEDURE dynamics_wrd_local
299    END INTERFACE dynamics_wrd_local
300
301    INTERFACE dynamics_last_actions
302       MODULE PROCEDURE dynamics_last_actions
303    END INTERFACE dynamics_last_actions
304
305
306 CONTAINS
307
308
309!--------------------------------------------------------------------------------------------------!
310! Description:
311! ------------
312!> Read module-specific namelist
313!--------------------------------------------------------------------------------------------------!
314 SUBROUTINE dynamics_parin
315
316
317    CHARACTER (LEN=80)  ::  line  !< dummy string that contains the current line of the parameter file
318
319    NAMELIST /dynamics_parameters/                                                                 &
320       dynamics_module_enabled
321
322
323    line = ' '
324!
325!-- Try to find module-specific namelist
326    REWIND ( 11 )
327    line = ' '
328    DO   WHILE ( INDEX( line, '&dynamics_parameters' ) == 0 )
329       READ ( 11, '(A)', END=12 )  line
330    ENDDO
331    BACKSPACE ( 11 )
332
333!-- Set default module switch to true
334    dynamics_module_enabled = .TRUE.
335
336!-- Read user-defined namelist
337    READ ( 11, dynamics_parameters, ERR = 10 )
338
339    GOTO 12
340
34110  BACKSPACE( 11 )
342    READ( 11 , '(A)') line
343    CALL parin_fail_message( 'dynamics_parameters', line )
344
34512  CONTINUE
346
347 END SUBROUTINE dynamics_parin
348
349
350!--------------------------------------------------------------------------------------------------!
351! Description:
352! ------------
353!> Check control parameters and deduce further quantities.
354!--------------------------------------------------------------------------------------------------!
355 SUBROUTINE dynamics_check_parameters
356
357
358 END SUBROUTINE dynamics_check_parameters
359
360
361!--------------------------------------------------------------------------------------------------!
362! Description:
363! ------------
364!> Set module-specific timeseries units and labels
365!--------------------------------------------------------------------------------------------------!
366 SUBROUTINE dynamics_check_data_output_ts( dots_max, dots_num, dots_label, dots_unit )
367
368    INTEGER(iwp),      INTENT(IN)     ::  dots_max
369
370    CHARACTER (LEN=*), DIMENSION(dots_max), INTENT(INOUT)  :: dots_label
371    CHARACTER (LEN=*), DIMENSION(dots_max), INTENT(INOUT)  :: dots_unit
372
373    INTEGER(iwp),      INTENT(INOUT)  ::  dots_num
374
375!
376!-- Next line is to avoid compiler warning about unused variables. Please remove.
377    IF ( dots_num == 0  .OR.  dots_label(1)(1:1) == ' '  .OR.  dots_unit(1)(1:1) == ' ' )  CONTINUE
378
379
380 END SUBROUTINE dynamics_check_data_output_ts
381
382
383!--------------------------------------------------------------------------------------------------!
384! Description:
385! ------------
386!> Set the unit of module-specific profile output quantities. For those variables not recognized,
387!> the parameter unit is set to "illegal", which tells the calling routine that the output variable
388!> is not defined and leads to a program abort.
389!--------------------------------------------------------------------------------------------------!
390 SUBROUTINE dynamics_check_data_output_pr( variable, var_count, unit, dopr_unit )
391
392
393    CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
394    CHARACTER (LEN=*) ::  unit     !<
395    CHARACTER (LEN=*) ::  variable !<
396
397    INTEGER(iwp) ::  var_count     !<
398
399!
400!-- Next line is to avoid compiler warning about unused variables. Please remove.
401    IF ( unit(1:1) == ' '  .OR.  dopr_unit(1:1) == ' '  .OR.  var_count == 0 )  CONTINUE
402
403    SELECT CASE ( TRIM( variable ) )
404
405!       CASE ( 'var_name' )
406
407       CASE DEFAULT
408          unit = 'illegal'
409
410    END SELECT
411
412
413 END SUBROUTINE dynamics_check_data_output_pr
414
415
416!--------------------------------------------------------------------------------------------------!
417! Description:
418! ------------
419!> Set the unit of module-specific output quantities. For those variables not recognized,
420!> the parameter unit is set to "illegal", which tells the calling routine that the output variable
421!< is not defined and leads to a program abort.
422!--------------------------------------------------------------------------------------------------!
423 SUBROUTINE dynamics_check_data_output( variable, unit )
424
425
426    CHARACTER (LEN=*) ::  unit     !<
427    CHARACTER (LEN=*) ::  variable !<
428
429    SELECT CASE ( TRIM( variable ) )
430
431!       CASE ( 'u2' )
432
433       CASE DEFAULT
434          unit = 'illegal'
435
436    END SELECT
437
438
439 END SUBROUTINE dynamics_check_data_output
440
441
442!--------------------------------------------------------------------------------------------------!
443!
444! Description:
445! ------------
446!> Initialize module-specific masked output
447!--------------------------------------------------------------------------------------------------!
448 SUBROUTINE dynamics_init_masks( variable, unit )
449
450
451    CHARACTER (LEN=*) ::  unit     !<
452    CHARACTER (LEN=*) ::  variable !<
453
454
455    SELECT CASE ( TRIM( variable ) )
456
457!       CASE ( 'u2' )
458
459       CASE DEFAULT
460          unit = 'illegal'
461
462    END SELECT
463
464
465 END SUBROUTINE dynamics_init_masks
466
467
468!--------------------------------------------------------------------------------------------------!
469! Description:
470! ------------
471!> Initialize module-specific arrays
472!--------------------------------------------------------------------------------------------------!
473 SUBROUTINE dynamics_init_arrays
474
475
476 END SUBROUTINE dynamics_init_arrays
477
478
479!--------------------------------------------------------------------------------------------------!
480! Description:
481! ------------
482!> Execution of module-specific initializing actions
483!--------------------------------------------------------------------------------------------------!
484 SUBROUTINE dynamics_init
485
486
487 END SUBROUTINE dynamics_init
488
489
490!--------------------------------------------------------------------------------------------------!
491! Description:
492! ------------
493!> Perform module-specific post-initialization checks
494!--------------------------------------------------------------------------------------------------!
495 SUBROUTINE dynamics_init_checks
496
497    INTEGER(iwp) ::  i !< loop index in x-direction
498    INTEGER(iwp) ::  j !< loop index in y-direction
499    INTEGER(iwp) ::  k !< loop index in z-direction
500
501    LOGICAL      ::  realistic_q = .TRUE. !< flag indicating realistic mixing ratios
502
503    REAL(wp)     ::  e_s !< saturation water vapor pressure
504    REAL(wp)     ::  q_s !< saturation mixing ratio
505    REAL(wp)     ::  t_l !< actual temperature
506
507!
508!-- Check for realistic initial mixing ratio. This must be in a realistic phyiscial range and must
509!-- not exceed the saturation mixing ratio by more than 2 percent. Please note, the check is
510!-- performed for each grid point (not just for a vertical profile), in order to cover also
511!-- three-dimensional initialization. Note, this check gives an error only for the initial run not
512!-- for a restart run. In case there are no cloud physics considered, the mixing ratio can exceed
513!-- the saturation moisture. This case a warning is given.
514    IF ( humidity  .AND.  .NOT. neutral  .AND.  check_realistic_q )  THEN
515       DO  i = nxl, nxr
516          DO  j = nys, nyn
517             DO  k = nzb+1, nzt
518!
519!--             Calculate actual temperature, water vapor saturation pressure, and based on this
520!--             the saturation mixing ratio.
521                t_l = exner(k) * pt(k,j,i)
522                e_s = magnus( t_l )
523                q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
524
525                IF ( q(k,j,i) > 1.02_wp * q_s )  realistic_q = .FALSE.
526             ENDDO
527          ENDDO
528       ENDDO
529!
530!--    Since the check is performed locally, merge the logical flag from all mpi ranks,
531!--    in order to do not print the error message multiple times.
532#if defined( __parallel )
533       CALL MPI_ALLREDUCE( MPI_IN_PLACE, realistic_q, 1, MPI_LOGICAL, MPI_LAND, comm2d, ierr)
534#endif
535
536       IF ( .NOT. realistic_q  .AND.                                                               &
537            TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
538          message_string = 'The initial mixing ratio exceeds the saturation mixing ratio.'
539          CALL message( 'dynamic_init_checks', 'PA0697', 2, 2, 0, 6, 0 )
540       ELSEIF ( .NOT. realistic_q  .AND.                                                           &
541                TRIM( initializing_actions ) == 'read_restart_data' )  THEN
542          message_string = 'The mixing ratio exceeds the saturation mixing ratio.'
543          CALL message( 'dynamic_init_checks', 'PA0697', 0, 1, 0, 6, 0 )
544       ENDIF
545    ENDIF
546
547 END SUBROUTINE dynamics_init_checks
548
549
550!--------------------------------------------------------------------------------------------------!
551! Description:
552! ------------
553!> Set the grids on which module-specific output quantities are defined. Allowed values for
554!> grid_x are "x" and "xu", for grid_y "y" and "yv", and for grid_z "zu" and "zw".
555!--------------------------------------------------------------------------------------------------!
556 SUBROUTINE dynamics_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
557
558
559    CHARACTER (LEN=*) ::  grid_x     !<
560    CHARACTER (LEN=*) ::  grid_y     !<
561    CHARACTER (LEN=*) ::  grid_z     !<
562    CHARACTER (LEN=*) ::  variable   !<
563
564    LOGICAL ::  found   !<
565
566
567    SELECT CASE ( TRIM( variable ) )
568
569!       CASE ( 'u2' )
570
571       CASE DEFAULT
572          found  = .FALSE.
573          grid_x = 'none'
574          grid_y = 'none'
575          grid_z = 'none'
576
577    END SELECT
578
579
580 END SUBROUTINE dynamics_define_netcdf_grid
581
582
583!--------------------------------------------------------------------------------------------------!
584! Description:
585! ------------
586!> Print a header with module-specific information.
587!--------------------------------------------------------------------------------------------------!
588 SUBROUTINE dynamics_header( io )
589
590
591    INTEGER(iwp) ::  io   !<
592
593!
594!-- If no module-specific variables are read from the namelist-file, no information will be printed.
595    IF ( .NOT. dynamics_module_enabled )  THEN
596       WRITE ( io, 100 )
597       RETURN
598    ENDIF
599
600!
601!-- Printing the information.
602    WRITE ( io, 110 )
603
604!
605!-- Format-descriptors
606100 FORMAT (//' *** dynamic module disabled'/)
607110 FORMAT (//1X,78('#')                                                                           &
608            //' User-defined variables and actions:'/                                              &
609              ' -----------------------------------'//)
610
611 END SUBROUTINE dynamics_header
612
613
614!--------------------------------------------------------------------------------------------------!
615! Description:
616! ------------
617!> Execute module-specific actions for all grid points
618!--------------------------------------------------------------------------------------------------!
619 SUBROUTINE dynamics_actions( location )
620
621
622    CHARACTER (LEN=*) ::  location !<
623
624!    INTEGER(iwp) ::  i !<
625!    INTEGER(iwp) ::  j !<
626!    INTEGER(iwp) ::  k !<
627
628!
629!-- Here the user-defined actions follow
630!-- No calls for single grid points are allowed at locations before and after the timestep, since
631!-- these calls are not within an i,j-loop
632    SELECT CASE ( location )
633
634       CASE ( 'before_timestep' )
635
636
637       CASE ( 'before_prognostic_equations' )
638
639
640       CASE ( 'after_integration' )
641
642
643       CASE ( 'after_timestep' )
644
645
646       CASE ( 'u-tendency' )
647
648
649       CASE ( 'v-tendency' )
650
651
652       CASE ( 'w-tendency' )
653
654
655       CASE ( 'pt-tendency' )
656
657
658       CASE ( 'sa-tendency' )
659
660
661       CASE ( 'e-tendency' )
662
663
664       CASE ( 'q-tendency' )
665
666
667       CASE ( 's-tendency' )
668
669
670       CASE DEFAULT
671          CONTINUE
672
673    END SELECT
674
675 END SUBROUTINE dynamics_actions
676
677
678!--------------------------------------------------------------------------------------------------!
679! Description:
680! ------------
681!> Execute module-specific actions for grid point i,j
682!--------------------------------------------------------------------------------------------------!
683 SUBROUTINE dynamics_actions_ij( i, j, location )
684
685
686    CHARACTER (LEN=*) ::  location
687
688    INTEGER(iwp) ::  i
689    INTEGER(iwp) ::  j
690
691!
692!-- Here the user-defined actions follow
693    SELECT CASE ( location )
694
695       CASE ( 'u-tendency' )
696
697!
698!--       Next line is to avoid compiler warning about unused variables. Please remove.
699          IF ( i +  j < 0 )  CONTINUE
700
701       CASE ( 'v-tendency' )
702
703
704       CASE ( 'w-tendency' )
705
706
707       CASE ( 'pt-tendency' )
708
709
710       CASE ( 'sa-tendency' )
711
712
713       CASE ( 'e-tendency' )
714
715
716       CASE ( 'q-tendency' )
717
718
719       CASE ( 's-tendency' )
720
721
722       CASE DEFAULT
723          CONTINUE
724
725    END SELECT
726
727 END SUBROUTINE dynamics_actions_ij
728
729
730!--------------------------------------------------------------------------------------------------!
731! Description:
732! ------------
733!> Compute module-specific non-advective processes for all grid points
734!--------------------------------------------------------------------------------------------------!
735 SUBROUTINE dynamics_non_advective_processes
736
737
738
739 END SUBROUTINE dynamics_non_advective_processes
740
741
742!--------------------------------------------------------------------------------------------------!
743! Description:
744! ------------
745!> Compute module-specific non-advective processes for grid points i,j
746!--------------------------------------------------------------------------------------------------!
747 SUBROUTINE dynamics_non_advective_processes_ij( i, j )
748
749
750    INTEGER(iwp) ::  i                 !<
751    INTEGER(iwp) ::  j                 !<
752
753!
754!--    Next line is just to avoid compiler warnings about unused variables. You may remove it.
755       IF ( i + j < 0 )  CONTINUE
756
757
758 END SUBROUTINE dynamics_non_advective_processes_ij
759
760
761!--------------------------------------------------------------------------------------------------!
762! Description:
763! ------------
764!> Perform module-specific horizontal boundary exchange
765!--------------------------------------------------------------------------------------------------!
766 SUBROUTINE dynamics_exchange_horiz
767
768
769
770 END SUBROUTINE dynamics_exchange_horiz
771
772
773!--------------------------------------------------------------------------------------------------!
774! Description:
775! ------------
776!> Compute module-specific prognostic equations for all grid points
777!--------------------------------------------------------------------------------------------------!
778 SUBROUTINE dynamics_prognostic_equations
779
780
781
782 END SUBROUTINE dynamics_prognostic_equations
783
784
785!--------------------------------------------------------------------------------------------------!
786! Description:
787! ------------
788!> Compute module-specific prognostic equations for grid point i,j
789!--------------------------------------------------------------------------------------------------!
790 SUBROUTINE dynamics_prognostic_equations_ij( i, j, i_omp_start, tn )
791
792
793    INTEGER(iwp), INTENT(IN) ::  i            !< grid index in x-direction
794    INTEGER(iwp), INTENT(IN) ::  i_omp_start  !< first loop index of i-loop in prognostic_equations
795    INTEGER(iwp), INTENT(IN) ::  j            !< grid index in y-direction
796    INTEGER(iwp), INTENT(IN) ::  tn           !< task number of openmp task
797
798!
799!-- Next line is just to avoid compiler warnings about unused variables. You may remove it.
800    IF ( i + j + i_omp_start + tn < 0 )  CONTINUE
801
802 END SUBROUTINE dynamics_prognostic_equations_ij
803
804
805!--------------------------------------------------------------------------------------------------!
806! Description:
807! ------------
808!> Compute boundary conditions of dynamics model
809!--------------------------------------------------------------------------------------------------!
810 SUBROUTINE dynamics_boundary_conditions
811
812    IMPLICIT NONE
813
814    INTEGER(iwp) ::  i  !< grid index x direction
815    INTEGER(iwp) ::  j  !< grid index y direction
816    INTEGER(iwp) ::  k  !< grid index z direction
817    INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
818    INTEGER(iwp) ::  m  !< running index surface elements
819
820    REAL(wp)    ::  c_max !< maximum phase velocity allowed by CFL criterion, used for outflow boundary condition
821    REAL(wp)    ::  denom !< horizontal gradient of velocity component normal to the outflow boundary
822
823!
824!-- Bottom boundary
825    IF ( ibc_uv_b == 1 )  THEN
826       u_p(nzb,:,:) = u_p(nzb+1,:,:)
827       v_p(nzb,:,:) = v_p(nzb+1,:,:)
828    ENDIF
829!
830!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
831!-- of downward-facing surfaces.
832    DO  l = 0, 1
833       !$OMP PARALLEL DO PRIVATE( i, j, k )
834       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
835       !$ACC PRESENT(bc_h, w_p)
836       DO  m = 1, bc_h(l)%ns
837          i = bc_h(l)%i(m)
838          j = bc_h(l)%j(m)
839          k = bc_h(l)%k(m)
840          w_p(k+bc_h(l)%koff,j,i) = 0.0_wp
841       ENDDO
842    ENDDO
843
844!
845!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
846    IF ( ibc_uv_t == 0 )  THEN
847        !$ACC KERNELS PRESENT(u_p, v_p, u_init, v_init)
848        u_p(nzt+1,:,:) = u_init(nzt+1)
849        v_p(nzt+1,:,:) = v_init(nzt+1)
850        !$ACC END KERNELS
851    ELSEIF ( ibc_uv_t == 1 )  THEN
852        u_p(nzt+1,:,:) = u_p(nzt,:,:)
853        v_p(nzt+1,:,:) = v_p(nzt,:,:)
854    ENDIF
855
856!
857!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
858    IF ( .NOT.  child_domain  .AND.  .NOT.  nesting_offline  .AND.                                 &
859         TRIM(coupling_mode) /= 'vnested_fine' )  THEN
860       !$ACC KERNELS PRESENT(w_p)
861       w_p(nzt:nzt+1,:,:) = 0.0_wp  !< nzt is not a prognostic level (but cf. pres)
862       !$ACC END KERNELS
863    ENDIF
864
865!
866!-- Temperature at bottom and top boundary.
867!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by the sea surface temperature
868!-- of the coupled ocean model.
869!-- Dirichlet
870    IF ( .NOT. neutral )  THEN
871       IF ( ibc_pt_b == 0 )  THEN
872          DO  l = 0, 1
873             !$OMP PARALLEL DO PRIVATE( i, j, k )
874             DO  m = 1, bc_h(l)%ns
875                i = bc_h(l)%i(m)
876                j = bc_h(l)%j(m)
877                k = bc_h(l)%k(m)
878                pt_p(k+bc_h(l)%koff,j,i) = pt(k+bc_h(l)%koff,j,i)
879             ENDDO
880          ENDDO
881!
882!--    Neumann, zero-gradient
883       ELSEIF ( ibc_pt_b == 1 )  THEN
884          DO  l = 0, 1
885             !$OMP PARALLEL DO PRIVATE( i, j, k )
886             !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
887             !$ACC PRESENT(bc_h, pt_p)
888             DO  m = 1, bc_h(l)%ns
889                i = bc_h(l)%i(m)
890                j = bc_h(l)%j(m)
891                k = bc_h(l)%k(m)
892                pt_p(k+bc_h(l)%koff,j,i) = pt_p(k,j,i)
893             ENDDO
894          ENDDO
895       ENDIF
896
897!
898!--    Temperature at top boundary
899       IF ( ibc_pt_t == 0 )  THEN
900           pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
901!
902!--        In case of nudging adjust top boundary to pt which is
903!--        read in from NUDGING-DATA
904           IF ( nudging )  THEN
905              pt_p(nzt+1,:,:) = pt_init(nzt+1)
906           ENDIF
907       ELSEIF ( ibc_pt_t == 1 )  THEN
908           pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
909       ELSEIF ( ibc_pt_t == 2 )  THEN
910           !$ACC KERNELS PRESENT(pt_p, dzu)
911           pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
912           !$ACC END KERNELS
913       ENDIF
914    ENDIF
915!
916!-- Boundary conditions for total water content, bottom and top boundary (see also temperature)
917    IF ( humidity )  THEN
918!
919!--    Surface conditions for constant_humidity_flux
920!--    Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate
921!--    belongs to the atmospheric grid point, therefore, set q_p at k-1
922       IF ( ibc_q_b == 0 ) THEN
923
924          DO  l = 0, 1
925             !$OMP PARALLEL DO PRIVATE( i, j, k )
926             DO  m = 1, bc_h(l)%ns
927                i = bc_h(l)%i(m)
928                j = bc_h(l)%j(m)
929                k = bc_h(l)%k(m)
930                q_p(k+bc_h(l)%koff,j,i) = q(k+bc_h(l)%koff,j,i)
931             ENDDO
932          ENDDO
933
934       ELSE
935
936          DO  l = 0, 1
937             !$OMP PARALLEL DO PRIVATE( i, j, k )
938             DO  m = 1, bc_h(l)%ns
939                i = bc_h(l)%i(m)
940                j = bc_h(l)%j(m)
941                k = bc_h(l)%k(m)
942                q_p(k+bc_h(l)%koff,j,i) = q_p(k,j,i)
943             ENDDO
944          ENDDO
945       ENDIF
946!
947!--    Top boundary
948       IF ( ibc_q_t == 0 ) THEN
949          q_p(nzt+1,:,:) = q(nzt+1,:,:)
950       ELSEIF ( ibc_q_t == 1 ) THEN
951          q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
952       ENDIF
953    ENDIF
954!
955!-- Boundary conditions for scalar, bottom and top boundary (see also temperature)
956    IF ( passive_scalar )  THEN
957!
958!--    Surface conditions for constant_humidity_flux
959!--    Run loop over all non-natural and natural walls. Note, in wall-datatype the k coordinate
960!--    belongs to the atmospheric grid point, therefore, set s_p at k-1
961       IF ( ibc_s_b == 0 ) THEN
962
963          DO  l = 0, 1
964             !$OMP PARALLEL DO PRIVATE( i, j, k )
965             DO  m = 1, bc_h(l)%ns
966                i = bc_h(l)%i(m)
967                j = bc_h(l)%j(m)
968                k = bc_h(l)%k(m)
969                s_p(k+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i)
970             ENDDO
971          ENDDO
972
973       ELSE
974
975          DO  l = 0, 1
976             !$OMP PARALLEL DO PRIVATE( i, j, k )
977             DO  m = 1, bc_h(l)%ns
978                i = bc_h(l)%i(m)
979                j = bc_h(l)%j(m)
980                k = bc_h(l)%k(m)
981                s_p(k+bc_h(l)%koff,j,i) = s_p(k,j,i)
982             ENDDO
983          ENDDO
984       ENDIF
985!
986!--    Top boundary condition
987       IF ( ibc_s_t == 0 )  THEN
988          s_p(nzt+1,:,:) = s(nzt+1,:,:)
989       ELSEIF ( ibc_s_t == 1 )  THEN
990          s_p(nzt+1,:,:) = s_p(nzt,:,:)
991       ELSEIF ( ibc_s_t == 2 )  THEN
992          s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
993       ENDIF
994
995    ENDIF
996!
997!-- In case of inflow or nest boundary at the south boundary the boundary for v is at nys and in
998!-- case of inflow or nest boundary at the left boundary the boundary for u is at nxl. Since in
999!-- prognostic_equations (cache optimized version) these levels are handled as a prognostic level,
1000!-- boundary values have to be restored here.
1001    IF ( bc_dirichlet_s )  THEN
1002       v_p(:,nys,:) = v_p(:,nys-1,:)
1003    ELSEIF ( bc_dirichlet_l ) THEN
1004       u_p(:,:,nxl) = u_p(:,:,nxl-1)
1005    ENDIF
1006
1007!
1008!-- The same restoration for u at i=nxl and v at j=nys as above must be made in case of nest
1009!-- boundaries. This must not be done in case of vertical nesting mode as in that case the lateral
1010!-- boundaries are actually cyclic.
1011!-- Lateral oundary conditions for TKE and dissipation are set in tcm_boundary_conds.
1012    IF ( nesting_mode /= 'vertical'  .OR.  nesting_offline )  THEN
1013       IF ( bc_dirichlet_s )  THEN
1014          v_p(:,nys,:) = v_p(:,nys-1,:)
1015       ENDIF
1016       IF ( bc_dirichlet_l )  THEN
1017          u_p(:,:,nxl) = u_p(:,:,nxl-1)
1018       ENDIF
1019    ENDIF
1020
1021!
1022!-- Lateral boundary conditions for scalar quantities at the outflow.
1023!-- Lateral oundary conditions for TKE and dissipation are set in tcm_boundary_conds.
1024    IF ( bc_radiation_s )  THEN
1025       pt_p(:,nys-1,:)     = pt_p(:,nys,:)
1026       IF ( humidity )  THEN
1027          q_p(:,nys-1,:) = q_p(:,nys,:)
1028       ENDIF
1029       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
1030    ELSEIF ( bc_radiation_n )  THEN
1031       pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
1032       IF ( humidity )  THEN
1033          q_p(:,nyn+1,:) = q_p(:,nyn,:)
1034       ENDIF
1035       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
1036    ELSEIF ( bc_radiation_l )  THEN
1037       pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
1038       IF ( humidity )  THEN
1039          q_p(:,:,nxl-1) = q_p(:,:,nxl)
1040       ENDIF
1041       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
1042    ELSEIF ( bc_radiation_r )  THEN
1043       pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
1044       IF ( humidity )  THEN
1045          q_p(:,:,nxr+1) = q_p(:,:,nxr)
1046       ENDIF
1047       IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
1048    ENDIF
1049
1050!
1051!-- Radiation boundary conditions for the velocities at the respective outflow.
1052!-- The phase velocity is either assumed to the maximum phase velocity that ensures numerical
1053!-- stability (CFL-condition) or calculated after Orlanski(1976) and averaged along the outflow
1054!-- boundary.
1055    IF ( bc_radiation_s )  THEN
1056
1057       IF ( use_cmax )  THEN
1058          u_p(:,-1,:) = u(:,0,:)
1059          v_p(:,0,:)  = v(:,1,:)
1060          w_p(:,-1,:) = w(:,0,:)
1061       ELSEIF ( .NOT. use_cmax )  THEN
1062
1063          c_max = dy / dt_3d
1064
1065          c_u_m_l = 0.0_wp
1066          c_v_m_l = 0.0_wp
1067          c_w_m_l = 0.0_wp
1068
1069          c_u_m = 0.0_wp
1070          c_v_m = 0.0_wp
1071          c_w_m = 0.0_wp
1072
1073!
1074!--       Calculate the phase speeds for u, v, and w, first local and then average along the outflow
1075!--       boundary.
1076          DO  k = nzb+1, nzt+1
1077             DO  i = nxl, nxr
1078
1079                denom = u_m_s(k,0,i) - u_m_s(k,1,i)
1080
1081                IF ( denom /= 0.0_wp )  THEN
1082                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
1083                   IF ( c_u(k,i) < 0.0_wp )  THEN
1084                      c_u(k,i) = 0.0_wp
1085                   ELSEIF ( c_u(k,i) > c_max )  THEN
1086                      c_u(k,i) = c_max
1087                   ENDIF
1088                ELSE
1089                   c_u(k,i) = c_max
1090                ENDIF
1091
1092                denom = v_m_s(k,1,i) - v_m_s(k,2,i)
1093
1094                IF ( denom /= 0.0_wp )  THEN
1095                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
1096                   IF ( c_v(k,i) < 0.0_wp )  THEN
1097                      c_v(k,i) = 0.0_wp
1098                   ELSEIF ( c_v(k,i) > c_max )  THEN
1099                      c_v(k,i) = c_max
1100                   ENDIF
1101                ELSE
1102                   c_v(k,i) = c_max
1103                ENDIF
1104
1105                denom = w_m_s(k,0,i) - w_m_s(k,1,i)
1106
1107                IF ( denom /= 0.0_wp )  THEN
1108                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
1109                   IF ( c_w(k,i) < 0.0_wp )  THEN
1110                      c_w(k,i) = 0.0_wp
1111                   ELSEIF ( c_w(k,i) > c_max )  THEN
1112                      c_w(k,i) = c_max
1113                   ENDIF
1114                ELSE
1115                   c_w(k,i) = c_max
1116                ENDIF
1117
1118                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
1119                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
1120                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
1121
1122             ENDDO
1123          ENDDO
1124
1125#if defined( __parallel )
1126          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
1127          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx,   &
1128                              ierr )
1129          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
1130          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx,   &
1131                              ierr )
1132          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
1133          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx,   &
1134                              ierr )
1135#else
1136          c_u_m = c_u_m_l
1137          c_v_m = c_v_m_l
1138          c_w_m = c_w_m_l
1139#endif
1140
1141          c_u_m = c_u_m / (nx+1)
1142          c_v_m = c_v_m / (nx+1)
1143          c_w_m = c_w_m / (nx+1)
1144
1145!
1146!--       Save old timelevels for the next timestep
1147          IF ( intermediate_timestep_count == 1 )  THEN
1148             u_m_s(:,:,:) = u(:,0:1,:)
1149             v_m_s(:,:,:) = v(:,1:2,:)
1150             w_m_s(:,:,:) = w(:,0:1,:)
1151          ENDIF
1152
1153!
1154!--       Calculate the new velocities
1155          DO  k = nzb+1, nzt+1
1156             DO  i = nxlg, nxrg
1157                u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) *                              &
1158                                          ( u(k,-1,i) - u(k,0,i) ) * ddy
1159
1160                v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v_m(k) *                              &
1161                                          ( v(k,0,i) - v(k,1,i) ) * ddy
1162
1163                w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) *                              &
1164                                          ( w(k,-1,i) - w(k,0,i) ) * ddy
1165             ENDDO
1166          ENDDO
1167
1168!
1169!--       Bottom boundary at the outflow
1170          IF ( ibc_uv_b == 0 )  THEN
1171             u_p(nzb,-1,:) = 0.0_wp
1172             v_p(nzb,0,:)  = 0.0_wp
1173          ELSE
1174             u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
1175             v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
1176          ENDIF
1177          w_p(nzb,-1,:) = 0.0_wp
1178
1179!
1180!--       Top boundary at the outflow
1181          IF ( ibc_uv_t == 0 )  THEN
1182             u_p(nzt+1,-1,:) = u_init(nzt+1)
1183             v_p(nzt+1,0,:)  = v_init(nzt+1)
1184          ELSE
1185             u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
1186             v_p(nzt+1,0,:)  = v_p(nzt,0,:)
1187          ENDIF
1188          w_p(nzt:nzt+1,-1,:) = 0.0_wp
1189
1190       ENDIF
1191
1192    ENDIF
1193
1194    IF ( bc_radiation_n )  THEN
1195
1196       IF ( use_cmax )  THEN
1197          u_p(:,ny+1,:) = u(:,ny,:)
1198          v_p(:,ny+1,:) = v(:,ny,:)
1199          w_p(:,ny+1,:) = w(:,ny,:)
1200       ELSEIF ( .NOT. use_cmax )  THEN
1201
1202          c_max = dy / dt_3d
1203
1204          c_u_m_l = 0.0_wp
1205          c_v_m_l = 0.0_wp
1206          c_w_m_l = 0.0_wp
1207
1208          c_u_m = 0.0_wp
1209          c_v_m = 0.0_wp
1210          c_w_m = 0.0_wp
1211
1212!
1213!--       Calculate the phase speeds for u, v, and w, first local and then average along the outflow
1214!--       boundary.
1215          DO  k = nzb+1, nzt+1
1216             DO  i = nxl, nxr
1217
1218                denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
1219
1220                IF ( denom /= 0.0_wp )  THEN
1221                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
1222                   IF ( c_u(k,i) < 0.0_wp )  THEN
1223                      c_u(k,i) = 0.0_wp
1224                   ELSEIF ( c_u(k,i) > c_max )  THEN
1225                      c_u(k,i) = c_max
1226                   ENDIF
1227                ELSE
1228                   c_u(k,i) = c_max
1229                ENDIF
1230
1231                denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
1232
1233                IF ( denom /= 0.0_wp )  THEN
1234                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
1235                   IF ( c_v(k,i) < 0.0_wp )  THEN
1236                      c_v(k,i) = 0.0_wp
1237                   ELSEIF ( c_v(k,i) > c_max )  THEN
1238                      c_v(k,i) = c_max
1239                   ENDIF
1240                ELSE
1241                   c_v(k,i) = c_max
1242                ENDIF
1243
1244                denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
1245
1246                IF ( denom /= 0.0_wp )  THEN
1247                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
1248                   IF ( c_w(k,i) < 0.0_wp )  THEN
1249                      c_w(k,i) = 0.0_wp
1250                   ELSEIF ( c_w(k,i) > c_max )  THEN
1251                      c_w(k,i) = c_max
1252                   ENDIF
1253                ELSE
1254                   c_w(k,i) = c_max
1255                ENDIF
1256
1257                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
1258                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
1259                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
1260
1261             ENDDO
1262          ENDDO
1263
1264#if defined( __parallel )
1265          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
1266          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx,   &
1267                              ierr )
1268          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
1269          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx,   &
1270                              ierr )
1271          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
1272          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dx,   &
1273                              ierr )
1274#else
1275          c_u_m = c_u_m_l
1276          c_v_m = c_v_m_l
1277          c_w_m = c_w_m_l
1278#endif
1279
1280          c_u_m = c_u_m / (nx+1)
1281          c_v_m = c_v_m / (nx+1)
1282          c_w_m = c_w_m / (nx+1)
1283
1284!
1285!--       Save old timelevels for the next timestep
1286          IF ( intermediate_timestep_count == 1 )  THEN
1287                u_m_n(:,:,:) = u(:,ny-1:ny,:)
1288                v_m_n(:,:,:) = v(:,ny-1:ny,:)
1289                w_m_n(:,:,:) = w(:,ny-1:ny,:)
1290          ENDIF
1291
1292!
1293!--       Calculate the new velocities
1294          DO  k = nzb+1, nzt+1
1295             DO  i = nxlg, nxrg
1296                u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) *                          &
1297                                              ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
1298
1299                v_p(k,ny+1,i) = v(k,ny+1,i)  - dt_3d * tsc(2) * c_v_m(k) *                         &
1300                                               ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
1301
1302                w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) *                          &
1303                                              ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
1304             ENDDO
1305          ENDDO
1306
1307!
1308!--       Bottom boundary at the outflow
1309          IF ( ibc_uv_b == 0 )  THEN
1310             u_p(nzb,ny+1,:) = 0.0_wp
1311             v_p(nzb,ny+1,:) = 0.0_wp
1312          ELSE
1313             u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
1314             v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
1315          ENDIF
1316          w_p(nzb,ny+1,:) = 0.0_wp
1317
1318!
1319!--       Top boundary at the outflow
1320          IF ( ibc_uv_t == 0 )  THEN
1321             u_p(nzt+1,ny+1,:) = u_init(nzt+1)
1322             v_p(nzt+1,ny+1,:) = v_init(nzt+1)
1323          ELSE
1324             u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
1325             v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
1326          ENDIF
1327          w_p(nzt:nzt+1,ny+1,:) = 0.0_wp
1328
1329       ENDIF
1330
1331    ENDIF
1332
1333    IF ( bc_radiation_l )  THEN
1334
1335       IF ( use_cmax )  THEN
1336          u_p(:,:,0)  = u(:,:,1)
1337          v_p(:,:,-1) = v(:,:,0)
1338          w_p(:,:,-1) = w(:,:,0)
1339       ELSEIF ( .NOT. use_cmax )  THEN
1340
1341          c_max = dx / dt_3d
1342
1343          c_u_m_l = 0.0_wp
1344          c_v_m_l = 0.0_wp
1345          c_w_m_l = 0.0_wp
1346
1347          c_u_m = 0.0_wp
1348          c_v_m = 0.0_wp
1349          c_w_m = 0.0_wp
1350
1351!
1352!--       Calculate the phase speeds for u, v, and w, first local and then average along the outflow
1353!--       boundary.
1354          DO  k = nzb+1, nzt+1
1355             DO  j = nys, nyn
1356
1357                denom = u_m_l(k,j,1) - u_m_l(k,j,2)
1358
1359                IF ( denom /= 0.0_wp )  THEN
1360                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
1361                   IF ( c_u(k,j) < 0.0_wp )  THEN
1362                      c_u(k,j) = 0.0_wp
1363                   ELSEIF ( c_u(k,j) > c_max )  THEN
1364                      c_u(k,j) = c_max
1365                   ENDIF
1366                ELSE
1367                   c_u(k,j) = c_max
1368                ENDIF
1369
1370                denom = v_m_l(k,j,0) - v_m_l(k,j,1)
1371
1372                IF ( denom /= 0.0_wp )  THEN
1373                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
1374                   IF ( c_v(k,j) < 0.0_wp )  THEN
1375                      c_v(k,j) = 0.0_wp
1376                   ELSEIF ( c_v(k,j) > c_max )  THEN
1377                      c_v(k,j) = c_max
1378                   ENDIF
1379                ELSE
1380                   c_v(k,j) = c_max
1381                ENDIF
1382
1383                denom = w_m_l(k,j,0) - w_m_l(k,j,1)
1384
1385                IF ( denom /= 0.0_wp )  THEN
1386                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
1387                   IF ( c_w(k,j) < 0.0_wp )  THEN
1388                      c_w(k,j) = 0.0_wp
1389                   ELSEIF ( c_w(k,j) > c_max )  THEN
1390                      c_w(k,j) = c_max
1391                   ENDIF
1392                ELSE
1393                   c_w(k,j) = c_max
1394                ENDIF
1395
1396                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
1397                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
1398                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
1399
1400             ENDDO
1401          ENDDO
1402
1403#if defined( __parallel )
1404          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
1405          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy,   &
1406                              ierr )
1407          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
1408          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy,   &
1409                              ierr )
1410          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
1411          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy,   &
1412                              ierr )
1413#else
1414          c_u_m = c_u_m_l
1415          c_v_m = c_v_m_l
1416          c_w_m = c_w_m_l
1417#endif
1418
1419          c_u_m = c_u_m / (ny+1)
1420          c_v_m = c_v_m / (ny+1)
1421          c_w_m = c_w_m / (ny+1)
1422
1423!
1424!--       Save old timelevels for the next timestep
1425          IF ( intermediate_timestep_count == 1 )  THEN
1426                u_m_l(:,:,:) = u(:,:,1:2)
1427                v_m_l(:,:,:) = v(:,:,0:1)
1428                w_m_l(:,:,:) = w(:,:,0:1)
1429          ENDIF
1430
1431!
1432!--       Calculate the new velocities
1433          DO  k = nzb+1, nzt+1
1434             DO  j = nysg, nyng
1435                u_p(k,j,0)  = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) *                               &
1436                                         ( u(k,j,0) - u(k,j,1) ) * ddx
1437
1438                v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) *                              &
1439                                          ( v(k,j,-1) - v(k,j,0) ) * ddx
1440
1441                w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) *                              &
1442                                          ( w(k,j,-1) - w(k,j,0) ) * ddx
1443             ENDDO
1444          ENDDO
1445
1446!
1447!--       Bottom boundary at the outflow
1448          IF ( ibc_uv_b == 0 )  THEN
1449             u_p(nzb,:,0)  = 0.0_wp
1450             v_p(nzb,:,-1) = 0.0_wp
1451          ELSE
1452             u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
1453             v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
1454          ENDIF
1455          w_p(nzb,:,-1) = 0.0_wp
1456
1457!
1458!--       Top boundary at the outflow
1459          IF ( ibc_uv_t == 0 )  THEN
1460             u_p(nzt+1,:,0)  = u_init(nzt+1)
1461             v_p(nzt+1,:,-1) = v_init(nzt+1)
1462          ELSE
1463             u_p(nzt+1,:,0)  = u_p(nzt,:,0)
1464             v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
1465          ENDIF
1466          w_p(nzt:nzt+1,:,-1) = 0.0_wp
1467
1468       ENDIF
1469
1470    ENDIF
1471
1472    IF ( bc_radiation_r )  THEN
1473
1474       IF ( use_cmax )  THEN
1475          u_p(:,:,nx+1) = u(:,:,nx)
1476          v_p(:,:,nx+1) = v(:,:,nx)
1477          w_p(:,:,nx+1) = w(:,:,nx)
1478       ELSEIF ( .NOT. use_cmax )  THEN
1479
1480          c_max = dx / dt_3d
1481
1482          c_u_m_l = 0.0_wp
1483          c_v_m_l = 0.0_wp
1484          c_w_m_l = 0.0_wp
1485
1486          c_u_m = 0.0_wp
1487          c_v_m = 0.0_wp
1488          c_w_m = 0.0_wp
1489
1490!
1491!--       Calculate the phase speeds for u, v, and w, first local and then average along the outflow
1492!--       boundary.
1493          DO  k = nzb+1, nzt+1
1494             DO  j = nys, nyn
1495
1496                denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
1497
1498                IF ( denom /= 0.0_wp )  THEN
1499                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
1500                   IF ( c_u(k,j) < 0.0_wp )  THEN
1501                      c_u(k,j) = 0.0_wp
1502                   ELSEIF ( c_u(k,j) > c_max )  THEN
1503                      c_u(k,j) = c_max
1504                   ENDIF
1505                ELSE
1506                   c_u(k,j) = c_max
1507                ENDIF
1508
1509                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
1510
1511                IF ( denom /= 0.0_wp )  THEN
1512                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
1513                   IF ( c_v(k,j) < 0.0_wp )  THEN
1514                      c_v(k,j) = 0.0_wp
1515                   ELSEIF ( c_v(k,j) > c_max )  THEN
1516                      c_v(k,j) = c_max
1517                   ENDIF
1518                ELSE
1519                   c_v(k,j) = c_max
1520                ENDIF
1521
1522                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
1523
1524                IF ( denom /= 0.0_wp )  THEN
1525                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
1526                   IF ( c_w(k,j) < 0.0_wp )  THEN
1527                      c_w(k,j) = 0.0_wp
1528                   ELSEIF ( c_w(k,j) > c_max )  THEN
1529                      c_w(k,j) = c_max
1530                   ENDIF
1531                ELSE
1532                   c_w(k,j) = c_max
1533                ENDIF
1534
1535                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
1536                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
1537                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
1538
1539             ENDDO
1540          ENDDO
1541
1542#if defined( __parallel )
1543          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
1544          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy,   &
1545                              ierr )
1546          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
1547          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy,   &
1548                              ierr )
1549          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
1550          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, MPI_SUM, comm1dy,   &
1551                              ierr )
1552#else
1553          c_u_m = c_u_m_l
1554          c_v_m = c_v_m_l
1555          c_w_m = c_w_m_l
1556#endif
1557
1558          c_u_m = c_u_m / (ny+1)
1559          c_v_m = c_v_m / (ny+1)
1560          c_w_m = c_w_m / (ny+1)
1561
1562!
1563!--       Save old timelevels for the next timestep
1564          IF ( intermediate_timestep_count == 1 )  THEN
1565                u_m_r(:,:,:) = u(:,:,nx-1:nx)
1566                v_m_r(:,:,:) = v(:,:,nx-1:nx)
1567                w_m_r(:,:,:) = w(:,:,nx-1:nx)
1568          ENDIF
1569
1570!
1571!--       Calculate the new velocities
1572          DO  k = nzb+1, nzt+1
1573             DO  j = nysg, nyng
1574                u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) *                          &
1575                                              ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
1576
1577                v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) *                          &
1578                                              ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
1579
1580                w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) *                          &
1581                                              ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
1582             ENDDO
1583          ENDDO
1584
1585!
1586!--       Bottom boundary at the outflow
1587          IF ( ibc_uv_b == 0 )  THEN
1588             u_p(nzb,:,nx+1) = 0.0_wp
1589             v_p(nzb,:,nx+1) = 0.0_wp
1590          ELSE
1591             u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
1592             v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
1593          ENDIF
1594          w_p(nzb,:,nx+1) = 0.0_wp
1595
1596!
1597!--       Top boundary at the outflow
1598          IF ( ibc_uv_t == 0 )  THEN
1599             u_p(nzt+1,:,nx+1) = u_init(nzt+1)
1600             v_p(nzt+1,:,nx+1) = v_init(nzt+1)
1601          ELSE
1602             u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
1603             v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
1604          ENDIF
1605          w_p(nzt:nzt+1,:,nx+1) = 0.0_wp
1606
1607       ENDIF
1608
1609    ENDIF
1610
1611 END SUBROUTINE dynamics_boundary_conditions
1612!--------------------------------------------------------------------------------------------------!
1613! Description:
1614! ------------
1615!> Swap timelevels of module-specific array pointers
1616!--------------------------------------------------------------------------------------------------!
1617 SUBROUTINE dynamics_swap_timelevel ( mod_count )
1618
1619
1620    INTEGER, INTENT(IN) :: mod_count
1621
1622
1623    SELECT CASE ( mod_count )
1624
1625       CASE ( 0 )
1626
1627          u  => u_1;   u_p  => u_2
1628          v  => v_1;   v_p  => v_2
1629          w  => w_1;   w_p  => w_2
1630          IF ( .NOT. neutral )  THEN
1631             pt => pt_1;  pt_p => pt_2
1632          ENDIF
1633          IF ( humidity )  THEN
1634             q => q_1;    q_p => q_2
1635          ENDIF
1636          IF ( passive_scalar )  THEN
1637             s => s_1;    s_p => s_2
1638          ENDIF
1639
1640       CASE ( 1 )
1641
1642          u  => u_2;   u_p  => u_1
1643          v  => v_2;   v_p  => v_1
1644          w  => w_2;   w_p  => w_1
1645          IF ( .NOT. neutral )  THEN
1646             pt => pt_2;  pt_p => pt_1
1647          ENDIF
1648          IF ( humidity )  THEN
1649             q => q_2;    q_p => q_1
1650          ENDIF
1651          IF ( passive_scalar )  THEN
1652             s => s_2;    s_p => s_1
1653          ENDIF
1654
1655    END SELECT
1656
1657 END SUBROUTINE dynamics_swap_timelevel
1658
1659
1660!--------------------------------------------------------------------------------------------------!
1661! Description:
1662! ------------
1663!> Sum up and time-average module-specific output quantities as well as allocate the array necessary
1664!> for storing the average.
1665!--------------------------------------------------------------------------------------------------!
1666 SUBROUTINE dynamics_3d_data_averaging( mode, variable )
1667
1668
1669    CHARACTER (LEN=*) ::  mode    !<
1670    CHARACTER (LEN=*) :: variable !<
1671
1672
1673    IF ( mode == 'allocate' )  THEN
1674
1675       SELECT CASE ( TRIM( variable ) )
1676
1677!          CASE ( 'u2' )
1678
1679          CASE DEFAULT
1680             CONTINUE
1681
1682       END SELECT
1683
1684    ELSEIF ( mode == 'sum' )  THEN
1685
1686       SELECT CASE ( TRIM( variable ) )
1687
1688!          CASE ( 'u2' )
1689
1690          CASE DEFAULT
1691             CONTINUE
1692
1693       END SELECT
1694
1695    ELSEIF ( mode == 'average' )  THEN
1696
1697       SELECT CASE ( TRIM( variable ) )
1698
1699!          CASE ( 'u2' )
1700
1701       END SELECT
1702
1703    ENDIF
1704
1705
1706 END SUBROUTINE dynamics_3d_data_averaging
1707
1708
1709!--------------------------------------------------------------------------------------------------!
1710! Description:
1711! ------------
1712!> Resorts the module-specific output quantity with indices (k,j,i) to a temporary array with
1713!> indices (i,j,k) and sets the grid on which it is defined.
1714!> Allowed values for grid are "zu" and "zw".
1715!--------------------------------------------------------------------------------------------------!
1716 SUBROUTINE dynamics_data_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do,     &
1717                                     nzt_do, fill_value )
1718
1719
1720    CHARACTER (LEN=*)             ::  grid       !<
1721    CHARACTER (LEN=*), INTENT(IN) ::  mode       !< either 'xy', 'xz' or 'yz'
1722    CHARACTER (LEN=*)             ::  variable   !<
1723
1724    INTEGER(iwp) ::  av     !< flag to control data output of instantaneous or time-averaged data
1725!    INTEGER(iwp) ::  i      !< grid index along x-direction
1726!    INTEGER(iwp) ::  j      !< grid index along y-direction
1727!    INTEGER(iwp) ::  k      !< grid index along z-direction
1728!    INTEGER(iwp) ::  m      !< running index surface elements
1729    INTEGER(iwp) ::  nzb_do !< lower limit of the domain (usually nzb)
1730    INTEGER(iwp) ::  nzt_do !< upper limit of the domain (usually nzt+1)
1731
1732    LOGICAL      ::  found !<
1733    LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
1734
1735    REAL(wp), INTENT(IN) ::  fill_value
1736
1737    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
1738
1739!
1740!-- Next line is just to avoid compiler warnings about unused variables. You may remove it.
1741    IF ( two_d .AND. av + LEN( mode ) + local_pf(nxl,nys,nzb_do) + fill_value < 0.0 )  CONTINUE
1742
1743    found = .TRUE.
1744
1745    SELECT CASE ( TRIM( variable ) )
1746
1747!       CASE ( 'u2_xy', 'u2_xz', 'u2_yz' )
1748
1749       CASE DEFAULT
1750          found = .FALSE.
1751          grid  = 'none'
1752
1753    END SELECT
1754
1755
1756 END SUBROUTINE dynamics_data_output_2d
1757
1758
1759!--------------------------------------------------------------------------------------------------!
1760! Description:
1761! ------------
1762!> Resorts the module-specific output quantity with indices (k,j,i) to a temporary array with
1763!> indices (i,j,k).
1764!--------------------------------------------------------------------------------------------------!
1765 SUBROUTINE dynamics_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
1766
1767
1768    CHARACTER (LEN=*) ::  variable !<
1769
1770    INTEGER(iwp) ::  av    !<
1771!    INTEGER(iwp) ::  i     !<
1772!    INTEGER(iwp) ::  j     !<
1773!    INTEGER(iwp) ::  k     !<
1774    INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
1775    INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
1776
1777    LOGICAL      ::  found !<
1778
1779    REAL(wp), INTENT(IN) ::  fill_value    !< value for the _FillValue attribute
1780
1781    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
1782
1783!
1784!-- Next line is to avoid compiler warning about unused variables. Please remove.
1785    IF ( av + local_pf(nxl,nys,nzb_do) + fill_value < 0.0 )  CONTINUE
1786
1787
1788    found = .TRUE.
1789
1790    SELECT CASE ( TRIM( variable ) )
1791
1792!       CASE ( 'u2' )
1793
1794       CASE DEFAULT
1795          found = .FALSE.
1796
1797    END SELECT
1798
1799
1800 END SUBROUTINE dynamics_data_output_3d
1801
1802
1803!--------------------------------------------------------------------------------------------------!
1804! Description:
1805! ------------
1806!> Calculation of module-specific statistics, i.e. horizontally averaged profiles and time series.
1807!> This is called for every statistic region sr, but at least for the region "total domain" (sr=0).
1808!--------------------------------------------------------------------------------------------------!
1809 SUBROUTINE dynamics_statistics( mode, sr, tn )
1810
1811
1812    CHARACTER (LEN=*) ::  mode   !<
1813!    INTEGER(iwp) ::  i    !<
1814!    INTEGER(iwp) ::  j    !<
1815!    INTEGER(iwp) ::  k    !<
1816    INTEGER(iwp) ::  sr   !<
1817    INTEGER(iwp) ::  tn   !<
1818
1819!
1820!-- Next line is to avoid compiler warning about unused variables. Please remove.
1821    IF ( sr == 0  .OR.  tn == 0 )  CONTINUE
1822
1823    IF ( mode == 'profiles' )  THEN
1824
1825    ELSEIF ( mode == 'time_series' )  THEN
1826
1827    ENDIF
1828
1829 END SUBROUTINE dynamics_statistics
1830
1831
1832!--------------------------------------------------------------------------------------------------!
1833! Description:
1834! ------------
1835!> Read module-specific global restart data (Fortran binary format).
1836!--------------------------------------------------------------------------------------------------!
1837 SUBROUTINE dynamics_rrd_global_ftn( found )
1838
1839    LOGICAL, INTENT(OUT)  ::  found
1840
1841
1842    found = .TRUE.
1843
1844
1845    SELECT CASE ( restart_string(1:length) )
1846
1847       CASE ( 'global_paramter' )
1848!          READ ( 13 )  global_parameter
1849
1850       CASE DEFAULT
1851
1852          found = .FALSE.
1853
1854    END SELECT
1855
1856
1857 END SUBROUTINE dynamics_rrd_global_ftn
1858
1859
1860!--------------------------------------------------------------------------------------------------!
1861! Description:
1862! ------------
1863!> Read module-specific global restart data (MPI-IO).
1864!--------------------------------------------------------------------------------------------------!
1865 SUBROUTINE dynamics_rrd_global_mpi
1866
1867!    CALL rrd_mpi_io( 'global_parameter', global_parameter )
1868    CONTINUE
1869
1870 END SUBROUTINE dynamics_rrd_global_mpi
1871
1872
1873!--------------------------------------------------------------------------------------------------!
1874! Description:
1875! ------------
1876!> Read module-specific local restart data arrays (Fortran binary format).
1877!> Subdomain index limits on file are given by nxl_on_file, etc.
1878!> Indices nxlc, etc. indicate the range of gridpoints to be mapped from the subdomain on file (f)
1879!> to the subdomain of the current PE (c). They have been calculated in routine rrd_local.
1880!--------------------------------------------------------------------------------------------------!
1881 SUBROUTINE dynamics_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf,     &
1882                                    nync, nyn_on_file, nysf, nysc, nys_on_file, tmp_2d, tmp_3d,    &
1883                                    found )
1884
1885
1886    INTEGER(iwp) ::  k               !<
1887    INTEGER(iwp) ::  nxlc            !<
1888    INTEGER(iwp) ::  nxlf            !<
1889    INTEGER(iwp) ::  nxl_on_file     !<
1890    INTEGER(iwp) ::  nxrc            !<
1891    INTEGER(iwp) ::  nxrf            !<
1892    INTEGER(iwp) ::  nxr_on_file     !<
1893    INTEGER(iwp) ::  nync            !<
1894    INTEGER(iwp) ::  nynf            !<
1895    INTEGER(iwp) ::  nyn_on_file     !<
1896    INTEGER(iwp) ::  nysc            !<
1897    INTEGER(iwp) ::  nysf            !<
1898    INTEGER(iwp) ::  nys_on_file     !<
1899
1900    LOGICAL, INTENT(OUT)  ::  found
1901
1902    REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d   !<
1903    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
1904
1905!
1906!-- Next line is to avoid compiler warning about unused variables. Please remove.
1907    IF ( k + nxlc + nxlf + nxrc + nxrf + nync + nynf + nysc + nysf +                               &
1908         tmp_2d(nys_on_file,nxl_on_file) +                                                         &
1909         tmp_3d(nzb,nys_on_file,nxl_on_file) < 0.0 )  CONTINUE
1910!
1911!-- Here the reading of user-defined restart data follows:
1912!-- Sample for user-defined output
1913
1914    found = .TRUE.
1915
1916    SELECT CASE ( restart_string(1:length) )
1917
1918!       CASE ( 'u2_av' )
1919
1920       CASE DEFAULT
1921
1922          found = .FALSE.
1923
1924    END SELECT
1925
1926 END SUBROUTINE dynamics_rrd_local_ftn
1927
1928
1929!--------------------------------------------------------------------------------------------------!
1930! Description:
1931! ------------
1932!> Read module-specific local restart data arrays (MPI-IO).
1933!--------------------------------------------------------------------------------------------------!
1934 SUBROUTINE dynamics_rrd_local_mpi
1935
1936    IMPLICIT NONE
1937
1938!    LOGICAL ::  array_found  !<
1939
1940
1941!    CALL rd_mpi_io_check_array( 'u2_av' , found = array_found )
1942!    IF ( array_found )  THEN
1943!       IF ( .NOT. ALLOCATED( u2_av ) )  ALLOCATE( u2_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1944!       CALL rrd_mpi_io( 'u2_av', u2_av )
1945!    ENDIF
1946
1947    CONTINUE
1948
1949 END SUBROUTINE dynamics_rrd_local_mpi
1950
1951
1952
1953!--------------------------------------------------------------------------------------------------!
1954! Description:
1955! ------------
1956!> Writes global module-specific restart data into binary file(s) for restart runs.
1957!--------------------------------------------------------------------------------------------------!
1958 SUBROUTINE dynamics_wrd_global
1959
1960
1961 END SUBROUTINE dynamics_wrd_global
1962
1963
1964!--------------------------------------------------------------------------------------------------!
1965! Description:
1966! ------------
1967!> Writes processor specific and module-specific restart data into binary file(s) for restart runs.
1968!--------------------------------------------------------------------------------------------------!
1969 SUBROUTINE dynamics_wrd_local
1970
1971
1972 END SUBROUTINE dynamics_wrd_local
1973
1974
1975!--------------------------------------------------------------------------------------------------!
1976! Description:
1977! ------------
1978!> Execute module-specific actions at the very end of the program.
1979!--------------------------------------------------------------------------------------------------!
1980 SUBROUTINE dynamics_last_actions
1981
1982
1983 END SUBROUTINE dynamics_last_actions
1984
1985 END MODULE dynamics_mod
Note: See TracBrowser for help on using the repository browser.