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

Last change on this file since 4808 was 4768, checked in by suehring, 4 years ago

Enable 3D data output also with 64-bit precision

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