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

Last change on this file since 4365 was 4360, checked in by suehring, 4 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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