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

Last change on this file since 4359 was 4359, checked in by suehring, 21 months ago

Fix for last commit

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