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

Last change on this file since 4740 was 4731, checked in by schwenkel, 4 years ago

Move exchange_horiz from time_integration to modules

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