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

Last change on this file since 4358 was 4358, checked in by suehring, 16 months ago

Bugfix, time coordinate is relative to origin_time rather than to 00:00:00 UTC; Refine post-initialization check for realistically inital values of mixing ratio. Give an error message for faulty initial values, but only a warning in a restart run.

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