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

Last change on this file since 4598 was 4517, checked in by raasch, 12 months ago

added restart with MPI-IO for reading local arrays

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