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

Last change on this file since 4347 was 4347, checked in by suehring, 20 months ago

Check for realistically initial values of mixing ratio implemented. Mixing ratio must not exceed its saturation value, else surface fluxes in the land-surface scheme may become unrealistic. Test cases updated

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