Changeset 4281 for palm/trunk/SOURCE


Ignore:
Timestamp:
Oct 29, 2019 3:15:39 PM (5 years ago)
Author:
schwenkel
Message:

Moved boundary_conds to dynamics module

Location:
palm/trunk/SOURCE
Files:
1 deleted
4 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r4270 r4281  
    2525# -----------------
    2626# $Id$
     27# delete boundary_conds, added missing dependencies
     28#
     29# 4270 2019-10-23 10:46:20Z monakurppa
    2730# - Implement offline nesting for salsa and add dependency of nesting_offl_mod
    2831#   on salsa_mod
     
    127130        basic_constants_and_equations_mod.f90 \
    128131        biometeorology_mod.f90 \
    129         boundary_conds.f90 \
    130132        buoyancy.f90 \
    131133        calc_mean_profile.f90 \
     
    345347        palm_date_time_mod.o \
    346348        radiation_model_mod.o
    347 boundary_conds.o: \
    348         bulk_cloud_model_mod.o \
    349         chemistry_model_mod.o \
    350         mod_kinds.o \
    351         modules.o \
    352         salsa_mod.o \
    353         surface_mod.o \
    354         turbulence_closure_mod.o
    355349bulk_cloud_model_mod.o: \
    356350        basic_constants_and_equations_mod.o \
     
    553547dynamics_mod.o: \
    554548        mod_kinds.o \
     549        surface_mod.o \
     550        pmc_interface_mod.o \
    555551        modules.o
    556552exchange_horiz.o: \
  • palm/trunk/SOURCE/dynamics_mod.f90

    r4097 r4281  
    2525! -----------------
    2626! $Id$
     27! Moved boundary conditions in dynamics module
     28!
     29! 4097 2019-07-15 11:59:11Z suehring
    2730! Avoid overlong lines - limit is 132 characters per line
    2831!
     
    3942
    4043    USE arrays_3d, &
    41         ONLY:  pt, pt_1, pt_2, pt_p, &
     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, &
    4247               q, q_1, q_2, q_p, &
    4348               s, s_1, s_2, s_p, &
    44                u, u_1, u_2, u_p, &
    45                v, v_1, v_2, v_p, &
    46                w, w_1, w_2, w_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
    4752
    4853    USE control_parameters, &
    49         ONLY:  length, &
     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, &
    5078               restart_string, &
    5179               humidity, &
    5280               neutral, &
    53                passive_scalar
     81               passive_scalar, &
     82               tsc, &
     83               use_cmax
     84
     85    USE grid_variables, &
     86        ONLY:  ddx, &
     87               ddy, &
     88               dx, &
     89               dy
    5490
    5591    USE indices, &
    5692        ONLY:  nbgp, &
     93               nx, &
    5794               nxl, &
     95               nxlg, &
    5896               nxr, &
     97               nxrg, &
     98               ny, &
    5999               nys, &
     100               nysg, &
    60101               nyn, &
     102               nyng, &
    61103               nzb, &
    62104               nzt
    63105
    64106    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
    65116
    66117    IMPLICIT NONE
     
    90141       dynamics_exchange_horiz, &
    91142       dynamics_prognostic_equations, &
     143       dynamics_boundary_conditions, &
    92144       dynamics_swap_timelevel, &
    93145       dynamics_3d_data_averaging, &
     
    169221    END INTERFACE dynamics_prognostic_equations
    170222
     223    INTERFACE dynamics_boundary_conditions
     224       MODULE PROCEDURE dynamics_boundary_conditions
     225    END INTERFACE dynamics_boundary_conditions
     226
    171227    INTERFACE dynamics_swap_timelevel
    172228       MODULE PROCEDURE dynamics_swap_timelevel
     
    657713
    658714
     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
    6591529!------------------------------------------------------------------------------!
    6601530! Description:
  • palm/trunk/SOURCE/module_interface.f90

    r4272 r4281  
    2525! -----------------
    2626! $Id$
     27! Added dynamics boundary conditions
     28!
     29! 4272 2019-10-23 15:18:57Z schwenkel
    2730! Further modularization of boundary conditions: moved boundary conditions to
    2831! respective modules
     
    179182               dynamics_exchange_horiz, &
    180183               dynamics_prognostic_equations, &
     184               dynamics_boundary_conditions, &
    181185               dynamics_swap_timelevel, &
    182186               dynamics_3d_data_averaging, &
     
    197201               tcm_actions, &
    198202               tcm_prognostic_equations, &
     203               tcm_boundary_conds, &
    199204               tcm_swap_timelevel, &
    200205               tcm_3d_data_averaging, &
     
    12991304    IF ( debug_output_timestep )  CALL debug_message( 'module-specific boundary_conditions', 'start' )
    13001305
     1306    CALL dynamics_boundary_conditions
     1307    CALL tcm_boundary_conds
     1308
    13011309    IF ( bulk_cloud_model    )  CALL bcm_boundary_conditions
    13021310    IF ( air_chemistry       )  CALL chem_boundary_conditions
  • palm/trunk/SOURCE/time_integration.f90

    r4276 r4281  
    2525! -----------------
    2626! $Id$
     27! Moved boundary conditions to module interface
     28!
     29! 4276 2019-10-28 16:03:29Z schwenkel
    2730! Further modularization of lpm code components
    2831!
     
    730733!
    731734!--       Boundary conditions for the prognostic quantities (except of the
    732 !--       velocities at the outflow in case of a non-cyclic lateral wall)
    733           CALL boundary_conds
    734 !
    735 !--       Boundary conditions for module-specific variables
     735!--       velocities at the outflow in case of a non-cyclic lateral wall) and
     736!--       boundary conditions for module-specific variables
    736737          CALL module_interface_boundary_conditions
    737738!
Note: See TracChangeset for help on using the changeset viewer.