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

Last change on this file since 4343 was 4281, checked in by schwenkel, 5 years ago

Moved boundary_conds to dynamics module

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