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

Last change on this file since 4505 was 4505, checked in by schwenkel, 15 months ago

Add flag for saturation check

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