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

Last change on this file since 4495 was 4495, checked in by raasch, 5 years ago

restart data handling with MPI-IO added, first part

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