Changeset 73


Ignore:
Timestamp:
Mar 20, 2007 8:33:14 AM (15 years ago)
Author:
raasch
Message:

preliminary changes for radiation conditions

Location:
palm/trunk/SOURCE
Files:
14 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/CURRENT_MODIFICATIONS

    r72 r73  
    3333Changed:
    3434-------
     35General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d
     36
    3537+age_m in particle_type
    3638
     
    4446Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.
    4547
     48q is not allowed to become negative (prognostic_equations).
     49
    4650__vtk directives removed from main program.
    4751
    4852The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.
    4953
    50 check_parameters, data_output_dvrp, data_output_ptseries, data_output_ts, flow_statistics, header, init_particles, init_3d_model, modules, palm, package_parin, parin, time_integration
     54check_parameters, data_output_dvrp, data_output_ptseries, data_output_ts, exchange_horiz_2d, flow_statistics, header, init_particles, init_pegrid, init_3d_model, modules, palm, package_parin, parin, prognostic_equations, read_3d_binary, time_integration, write_3d_binary
    5155
    5256
     
    5862Bugfix in sample for reading user defined data from restart file (user_init)
    5963
     64Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)
     65
    6066Check for possible negative humidities in the initial humidity profile.
    6167
     
    6470
    6571Makefile
    66 check_parameters, init_1d_model, user_interface
     72check_parameters, init_1d_model, init_3d_model, user_interface
  • palm/trunk/SOURCE/boundary_conds.f90

    r39 r73  
    44! Actual revisions:
    55! -----------------
    6 !
     6! The "main" part sets conditions for time level t+dt insteat of level t,
     7! outflow boundary conditions changed from Neumann to radiation condition
    78!
    89! Former revisions:
     
    5152    INTEGER ::  i, j, k
    5253
     54    REAL    ::  c_max, c_u, c_v, c_w, denom
     55
    5356
    5457    IF ( range == 'main')  THEN
     
    5659!--    Bottom boundary
    5760       IF ( ibc_uv_b == 0 )  THEN
    58           u(nzb,:,:) = -u(nzb+1,:,:)  ! satisfying the Dirichlet condition with
    59           v(nzb,:,:) = -v(nzb+1,:,:)  ! an extra layer below the surface where
    60        ELSE                           ! the u and v component change their sign
    61           u(nzb,:,:) = u(nzb+1,:,:)
    62           v(nzb,:,:) = v(nzb+1,:,:)
     61!
     62!--       Satisfying the Dirichlet condition with an extra layer below the
     63!--       surface where the u and v component change their sign
     64          u_p(nzb,:,:) = -u_p(nzb+1,:,:)
     65          v_p(nzb,:,:) = -v_p(nzb+1,:,:)
     66       ELSE
     67          u_p(nzb,:,:) = u_p(nzb+1,:,:)
     68          v_p(nzb,:,:) = v_p(nzb+1,:,:)
    6369       ENDIF
    6470       DO  i = nxl-1, nxr+1
    6571          DO  j = nys-1, nyn+1
    66              w(nzb_w_inner(j,i),j,i) = 0.0
     72             w_p(nzb_w_inner(j,i),j,i) = 0.0
    6773          ENDDO
    6874       ENDDO
     
    7177!--    Top boundary
    7278       IF ( ibc_uv_t == 0 )  THEN
    73           u(nzt+1,:,:) = ug(nzt+1)
    74           v(nzt+1,:,:) = vg(nzt+1)
     79          u_p(nzt+1,:,:) = ug(nzt+1)
     80          v_p(nzt+1,:,:) = vg(nzt+1)
    7581       ELSE
    76           u(nzt+1,:,:) = u(nzt,:,:)
    77           v(nzt+1,:,:) = v(nzt,:,:)
    78        ENDIF
    79        w(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
     82          u_p(nzt+1,:,:) = u_p(nzt,:,:)
     83          v_p(nzt+1,:,:) = v_p(nzt,:,:)
     84       ENDIF
     85       w_p(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
    8086
    8187!
    8288!--    Temperature at bottom boundary
    8389       IF ( ibc_pt_b == 0 )  THEN
    84           IF ( timestep_scheme(1:5) /= 'runge' )  THEN
    85              DO  i = nxl-1, nxr+1
    86                 DO  j = nys-1, nyn+1
    87                    pt(nzb_s_inner(j,i),j,i) = pt_m(nzb_s_inner(j,i),j,i)
    88                 ENDDO
    89              ENDDO
    90           ELSE
    91              DO  i = nxl-1, nxr+1
    92                 DO  j = nys-1, nyn+1
    93                    pt(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i),j,i) 
    94                                             ! pt_m is not used for Runge-Kutta
    95                 ENDDO
    96              ENDDO
    97           ENDIF
     90          DO  i = nxl-1, nxr+1
     91             DO  j = nys-1, nyn+1
     92                pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
     93             ENDDO
     94          ENDDO
    9895       ELSE
    9996          DO  i = nxl-1, nxr+1
    10097             DO  j = nys-1, nyn+1
    101                 pt(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i)+1,j,i)
     98                pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
    10299             ENDDO
    103100          ENDDO
     
    107104!--    Temperature at top boundary
    108105       IF ( ibc_pt_t == 0 )  THEN
    109           IF ( timestep_scheme(1:5) /= 'runge' )  THEN
    110              pt(nzt+1,:,:) = pt_m(nzt+1,:,:)
    111           ELSE
    112              pt(nzt+1,:,:) = pt_p(nzt+1,:,:)  ! pt_m not used for Runge-Kutta
    113           ENDIF
     106          pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
    114107       ELSEIF ( ibc_pt_t == 1 )  THEN
    115           pt(nzt+1,:,:) = pt(nzt,:,:)
     108          pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
    116109       ELSEIF ( ibc_pt_t == 2 )  THEN
    117           pt(nzt+1,:,:) = pt(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
     110          pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
    118111       ENDIF
    119112
     
    124117          DO  i = nxl-1, nxr+1
    125118             DO  j = nys-1, nyn+1
    126                 e(nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i)
    127              ENDDO
    128           ENDDO
    129           e(nzt+1,:,:) = e(nzt,:,:)
     119                e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
     120             ENDDO
     121          ENDDO
     122          e_p(nzt+1,:,:) = e_p(nzt,:,:)
    130123       ENDIF
    131124
     
    137130!--       Surface conditions for constant_moisture_flux
    138131          IF ( ibc_q_b == 0 ) THEN
    139              IF ( timestep_scheme(1:5) /= 'runge' )  THEN
    140                 DO  i = nxl-1, nxr+1
    141                    DO  j = nys-1, nyn+1
    142                       q(nzb_s_inner(j,i),j,i) = q_m(nzb_s_inner(j,i),j,i)
    143                    ENDDO
     132             DO  i = nxl-1, nxr+1
     133                DO  j = nys-1, nyn+1
     134                   q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
    144135                ENDDO
    145              ELSE
    146                 DO  i = nxl-1, nxr+1
    147                    DO  j = nys-1, nyn+1
    148                       q(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i),j,i)
    149                    ENDDO                      ! q_m is not used for Runge-Kutta
    150                 ENDDO
    151              ENDIF
     136             ENDDO
    152137          ELSE
    153138             DO  i = nxl-1, nxr+1
    154139                DO  j = nys-1, nyn+1
    155                    q(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i)+1,j,i)
     140                   q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
    156141                ENDDO
    157142             ENDDO
     
    159144!
    160145!--       Top boundary
    161           q(nzt+1,:,:) = q(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
     146          q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
    162147       ENDIF
    163148
     
    167152!--    divergence. Dirichlet conditions are used for all other quantities.
    168153       IF ( inflow_s )  THEN
    169           v(:,nys,:) = v(:,nys-1,:)
     154          v_p(:,nys,:) = v_p(:,nys-1,:)
    170155       ELSEIF ( inflow_n ) THEN
    171           v(:,nyn+vynp,:) = v(:,nyn+vynp+1,:)
     156          v_p(:,nyn+vynp,:) = v_p(:,nyn+vynp+1,:)
    172157       ELSEIF ( inflow_l ) THEN
    173           u(:,:,nxl) = u(:,:,nxl-1)
     158          u_p(:,:,nxl) = u_p(:,:,nxl-1)
    174159       ELSEIF ( inflow_r ) THEN
    175           u(:,:,nxr+uxrp) = u(:,:,nxr+uxrp+1)
     160          u_p(:,:,nxr+uxrp) = u_p(:,:,nxr+uxrp+1)
    176161       ENDIF
    177162
     
    179164!--    Lateral boundary conditions for scalar quantities at the outflow
    180165       IF ( outflow_s )  THEN
    181           pt(:,nys-1,:)     = pt(:,nys,:)
    182           IF ( .NOT. constant_diffusion     )  e(:,nys-1,:) = e(:,nys,:)
    183           IF ( moisture .OR. passive_scalar )  q(:,nys-1,:) = q(:,nys,:)
     166          pt_p(:,nys-1,:)     = pt_p(:,nys,:)
     167          IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
     168          IF ( moisture .OR. passive_scalar )  q_p(:,nys-1,:) = q_p(:,nys,:)
    184169       ELSEIF ( outflow_n )  THEN
    185           pt(:,nyn+1,:)     = pt(:,nyn,:)
    186           IF ( .NOT. constant_diffusion     )  e(:,nyn+1,:) = e(:,nyn,:)
    187           IF ( moisture .OR. passive_scalar )  q(:,nyn+1,:) = q(:,nyn,:)
     170          pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
     171          IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
     172          IF ( moisture .OR. passive_scalar )  q_p(:,nyn+1,:) = q_p(:,nyn,:)
    188173       ELSEIF ( outflow_l )  THEN
    189           pt(:,:,nxl-1)     = pt(:,:,nxl)
    190           IF ( .NOT. constant_diffusion     )  e(:,:,nxl-1) = e(:,:,nxl)
    191           IF ( moisture .OR. passive_scalar )  q(:,:,nxl-1) = q(:,:,nxl)     
     174          pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
     175          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
     176          IF ( moisture .OR. passive_scalar )  q_p(:,:,nxl-1) = q_p(:,:,nxl)     
    192177       ELSEIF ( outflow_r )  THEN
    193           pt(:,:,nxr+1)     = pt(:,:,nxr)
    194           IF ( .NOT. constant_diffusion     )  e(:,:,nxr+1) = e(:,:,nxr)
    195           IF ( moisture .OR. passive_scalar )  q(:,:,nxr+1) = q(:,:,nxr)
     178          pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
     179          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
     180          IF ( moisture .OR. passive_scalar )  q_p(:,:,nxr+1) = q_p(:,:,nxr)
    196181       ENDIF
    197182
     
    200185    IF ( range == 'outflow_uvw' ) THEN
    201186!
    202 !--    Horizontal boundary conditions for the velocities at the outflow.
    203 !--    A Neumann condition is used for the wall normal velocity. The vertical
    204 !--    velocity is assumed as zero and a horizontal average along the wall is
    205 !--    used for the wall parallel horizontal velocity. The combination of all
    206 !--    three conditions ensures that the velocity field is free of divergence
    207 !--    at the outflow (uvmean_outflow_l is calculated in pres).
     187!--    Radiation boundary condition for the velocities at the respective outflow
    208188       IF ( outflow_s ) THEN
    209           v(:,nys-1,:) = v(:,nys,:)
    210           w(:,nys-1,:) = 0.0
    211 !
    212 !--       Compute the mean horizontal wind parallel to and within the outflow
    213 !--       wall and use this as boundary condition for u
    214 #if defined( __parallel )
    215           CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
    216                               MPI_REAL, MPI_SUM, comm1dx, ierr )
    217           uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
    218 #else
    219           uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
    220 #endif
    221           DO  k = nzb, nzt+1
    222              u(k,nys-1,:) = uvmean_outflow(k)
    223           ENDDO
    224        ENDIF
    225 
    226        IF ( outflow_n ) THEN
    227           v(:,nyn+vynp+1,:) = v(:,nyn+vynp,:)
    228           w(:,nyn+1,:)      = 0.0
    229 !
    230 !--       Compute the mean horizontal wind parallel to and within the outflow
    231 !--       wall and use this as boundary condition for u
    232 #if defined( __parallel )
    233           CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
    234                               MPI_REAL, MPI_SUM, comm1dx, ierr )
    235           uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
    236 #else
    237           uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
    238 #endif
    239           DO  k = nzb, nzt+1
    240              u(k,nyn+1,:) = uvmean_outflow(k)
    241           ENDDO
     189!          v(:,nys-1,:) = v(:,nys,:)
     190!          w(:,nys-1,:) = 0.0
     191!!
     192!!--       Compute the mean horizontal wind parallel to and within the outflow
     193!!--       wall and use this as boundary condition for u
     194!#if defined( __parallel )
     195!          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
     196!                              MPI_REAL, MPI_SUM, comm1dx, ierr )
     197!          uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
     198!#else
     199!          uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
     200!#endif
     201!          DO  k = nzb, nzt+1
     202!             u(k,nys-1,:) = uvmean_outflow(k)
     203!          ENDDO
     204       ENDIF
     205
     206       IF ( outflow_n  .AND.                                                 &
     207            intermediate_timestep_count == intermediate_timestep_count_max ) &
     208       THEN
     209
     210          c_max = dy / dt_3d
     211
     212          DO i = nxl-1, nxr+1
     213             DO k = nzb+1, nzt+1
     214
     215!
     216!--             First calculate the phase speeds for u,v, and w
     217                denom = u_m_n(k,ny,i,-2) - u_m_n(k,ny-1,i,-2)
     218
     219                IF ( denom /= 0.0 )  THEN
     220                   c_u = -c_max * ( u_m_n(k,ny,i,-1)-u_m_n(k,ny,i,-2) ) / denom
     221                   IF ( c_u < 0.0 )  THEN
     222                      c_u = 0.0
     223                   ELSEIF ( c_u > c_max )  THEN
     224                      c_u = c_max
     225                   ENDIF
     226                ELSE
     227                   c_u = c_max
     228                ENDIF
     229
     230                denom = v_m_n(k,ny,i,-2) - v_m_n(k,ny-1,i,-2)
     231
     232                IF ( denom /= 0.0 )  THEN
     233                   c_v = -c_max * ( v_m_n(k,ny,i,-1)-v_m_n(k,ny,i,-2) ) / denom
     234                   IF ( c_v < 0.0 )  THEN
     235                      c_v = 0.0
     236                   ELSEIF ( c_v > c_max )  THEN
     237                      c_v = c_max
     238                   ENDIF
     239                ELSE
     240                   c_v = c_max
     241                ENDIF
     242
     243                denom = w_m_n(k,ny,i,-2) - w_m_n(k,ny-1,i,-2)
     244
     245                IF ( denom /= 0.0 )  THEN
     246                   c_w = -c_max * ( w_m_n(k,ny,i,-1)-w_m_n(k,ny,i,-2) ) / denom
     247                   IF ( c_w < 0.0 )  THEN
     248                      c_w = 0.0
     249                   ELSEIF ( c_w > c_max )  THEN
     250                      c_w = c_max
     251                   ENDIF
     252                ELSE
     253                   c_w = c_max
     254                ENDIF
     255
     256!
     257!--             Calculate the new velocities
     258                u(k,ny+1,i) = u_m_n(k,ny+1,i,-1) - dt_3d * c_u * &
     259                              ( u_m_n(k,ny+1,i,-1) - u_m_n(k,ny,i,-1) ) * ddy
     260
     261                v(k,ny+1,i) = v_m_n(k,ny+1,i,-1) - dt_3d * c_v * &
     262                              ( v_m_n(k,ny+1,i,-1) - v_m_n(k,ny,i,-1) ) * ddy
     263
     264                w(k,ny+1,i) = w_m_n(k,ny+1,i,-1) - dt_3d * c_w * &
     265                              ( w_m_n(k,ny+1,i,-1) - w_m_n(k,ny,i,-1) ) * ddy
     266
     267!
     268!--             Swap timelevels for the next timestep
     269                u_m_n(k,:,i,-2) = u_m_n(k,:,i,-1)
     270                u_m_n(k,:,i,-1) = u(k,ny-1:ny+1,i)
     271                v_m_n(k,:,i,-2) = v_m_n(k,:,i,-1)
     272                v_m_n(k,:,i,-1) = v(k,ny-1:ny+1,i)
     273                w_m_n(k,:,i,-2) = w_m_n(k,:,i,-1)
     274                w_m_n(k,:,i,-1) = w(k,ny-1:ny+1,i)
     275
     276             ENDDO
     277          ENDDO
     278
     279!
     280!--       Bottom boundary at the outflow
     281          IF ( ibc_uv_b == 0 )  THEN
     282             u(nzb,ny+1,:) = -u(nzb+1,ny+1,:)
     283             v(nzb,ny+1,:) = -v(nzb+1,ny+1,:) 
     284          ELSE                   
     285             u(nzb,ny+1,:) =  u(nzb+1,ny+1,:)
     286             v(nzb,ny+1,:) =  v(nzb+1,ny+1,:)
     287          ENDIF
     288          w(nzb,ny+1,:) = 0.0
     289
     290!
     291!--       Top boundary at the outflow
     292          IF ( ibc_uv_t == 0 )  THEN
     293             u(nzt+1,ny+1,:) = ug(nzt+1)
     294             v(nzt+1,ny+1,:) = vg(nzt+1)
     295          ELSE
     296             u(nzt+1,ny+1,:) = u(nzt,nyn+1,:)
     297             v(nzt+1,ny+1,:) = v(nzt,nyn+1,:)
     298          ENDIF
     299          w(nzt:nzt+1,ny+1,:) = 0.0
     300
    242301       ENDIF
    243302
    244303       IF ( outflow_l ) THEN
    245           u(:,:,nxl-1) = u(:,:,nxl)
    246           w(:,:,nxl-1) = 0.0
     304!          u(:,:,nxl-1) = u(:,:,nxl)
     305!          w(:,:,nxl-1) = 0.0
    247306!
    248307!--       Compute the mean horizontal wind parallel to and within the outflow
    249308!--       wall and use this as boundary condition for v
    250 #if defined( __parallel )
    251           CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
    252                               MPI_REAL, MPI_SUM, comm1dy, ierr )
    253           uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
    254 #else
    255           uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
    256 #endif
    257           DO  k = nzb, nzt+1
    258              v(k,:,nxl-1) = uvmean_outflow(k)
    259           ENDDO   
    260 
    261        ENDIF
    262 
    263        IF ( outflow_r ) THEN
    264           u(:,:,nxr+uxrp+1) = u(:,:,nxr+uxrp)
    265           w(:,:,nxr+1) = 0.0
    266 !
    267 !--       Compute the mean horizontal wind parallel to and within the outflow
    268 !--       wall and use this as boundary condition for v
    269 #if defined( __parallel )
    270           CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
    271                               MPI_REAL, MPI_SUM, comm1dy, ierr )
    272           uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
    273 #else
    274           uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
    275 #endif
    276           DO  k = nzb, nzt+1
    277              v(k,:,nxr+1) = uvmean_outflow(k)
    278           ENDDO   
     309!#if defined( __parallel )
     310!          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
     311!                              MPI_REAL, MPI_SUM, comm1dy, ierr )
     312!          uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
     313!#else
     314!          uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
     315!#endif
     316!          DO  k = nzb, nzt+1
     317!             v(k,:,nxl-1) = uvmean_outflow(k)
     318!          ENDDO   
     319!
     320       ENDIF
     321
     322       IF ( outflow_r  .AND.                                                 &
     323            intermediate_timestep_count == intermediate_timestep_count_max ) &
     324       THEN
     325
     326          c_max = dx / dt_3d
     327
     328          DO j = nys-1, nyn+1
     329             DO k = nzb+1, nzt+1
     330
     331!
     332!--             First calculate the phase speeds for u,v, and w
     333                denom = u_m_r(k,j,nx,-2) - u_m_r(k,j,nx-1,-2)
     334
     335                IF ( denom /= 0.0 )  THEN
     336                   c_u = -c_max * ( u_m_r(k,j,nx,-1)-u_m_r(k,j,nx,-2) ) / denom
     337                   IF ( c_u < 0.0 )  THEN
     338                      c_u = 0.0
     339                   ELSEIF ( c_u > c_max )  THEN
     340                      c_u = c_max
     341                   ENDIF
     342                ELSE
     343                   c_u = c_max
     344                ENDIF
     345
     346                denom = v_m_r(k,j,nx,-2) - v_m_r(k,j,nx-1,-2)
     347
     348                IF ( denom /= 0.0 )  THEN
     349                   c_v = -c_max * ( v_m_r(k,j,nx,-1)-v_m_r(k,j,nx,-2) ) / denom
     350                   IF ( c_v < 0.0 )  THEN
     351                      c_v = 0.0
     352                   ELSEIF ( c_v > c_max )  THEN
     353                      c_v = c_max
     354                   ENDIF
     355                ELSE
     356                   c_v = c_max
     357                ENDIF
     358
     359                denom = w_m_r(k,j,nx,-2) - w_m_r(k,j,nx-1,-2)
     360
     361                IF ( denom /= 0.0 )  THEN
     362                   c_w = -c_max * ( w_m_r(k,j,nx,-1)-w_m_n(k,j,nx,-2) ) / denom
     363                   IF ( c_w < 0.0 )  THEN
     364                      c_w = 0.0
     365                   ELSEIF ( c_w > c_max )  THEN
     366                      c_w = c_max
     367                   ENDIF
     368                ELSE
     369                   c_w = c_max
     370                ENDIF
     371
     372!
     373!--             Calculate the new velocities
     374                u(k,j,nx+1) = u_m_r(k,j,nx+1,-1) - dt_3d * c_u * &
     375                              ( u_m_r(k,j,nx+1,-1) - u_m_r(k,j,nx,-1) ) * ddx
     376
     377                v(k,j,nx+1) = v_m_r(k,j,nx+1,-1) - dt_3d * c_v * &
     378                              ( v_m_r(k,j,nx+1,-1) - v_m_r(k,j,nx,-1) ) * ddx
     379
     380                w(k,j,nx+1) = w_m_r(k,j,nx+1,-1) - dt_3d * c_w * &
     381                              ( w_m_r(k,j,nx+1,-1) - w_m_r(k,j,nx,-1) ) * ddx
     382
     383!
     384!--             Swap timelevels for the next timestep
     385                u_m_r(k,j,:,-2) = u_m_r(k,j,:,-1)
     386                u_m_r(k,j,:,-1) = u(k,j,nx-1:nx+1)
     387                v_m_r(k,j,:,-2) = v_m_r(k,j,:,-1)
     388                v_m_r(k,j,:,-1) = v(k,j,nx-1:nx+1)
     389                w_m_r(k,j,:,-2) = w_m_r(k,j,:,-1)
     390                w_m_r(k,j,:,-1) = w(k,j,nx-1:nx+1)
     391
     392             ENDDO
     393          ENDDO
     394
     395!
     396!--       Bottom boundary at the outflow
     397          IF ( ibc_uv_b == 0 )  THEN
     398             u(nzb,ny+1,:) = -u(nzb+1,ny+1,:)
     399             v(nzb,ny+1,:) = -v(nzb+1,ny+1,:) 
     400          ELSE                   
     401             u(nzb,ny+1,:) =  u(nzb+1,ny+1,:)
     402             v(nzb,ny+1,:) =  v(nzb+1,ny+1,:)
     403          ENDIF
     404          w(nzb,ny+1,:) = 0.0
     405
     406!
     407!--       Top boundary at the outflow
     408          IF ( ibc_uv_t == 0 )  THEN
     409             u(nzt+1,ny+1,:) = ug(nzt+1)
     410             v(nzt+1,ny+1,:) = vg(nzt+1)
     411          ELSE
     412             u(nzt+1,ny+1,:) = u(nzt,nyn+1,:)
     413             v(nzt+1,ny+1,:) = v(nzt,nyn+1,:)
     414          ENDIF
     415          w(nzt:nzt+1,ny+1,:) = 0.0
     416
    279417       ENDIF
    280418
  • palm/trunk/SOURCE/calc_precipitation.f90

    r72 r73  
    7171                tend(k,j,i) = tend(k,j,i) - dqdt_precip
    7272!
    73 !--             Precipitation rate in (kg * 0.001) / m**2 / s (because 1kg
    74 !--             gives 1 mm)
     73!--             Precipitation rate in kg / m**2 / s (= mm/s)
    7574                precipitation_rate(j,i) = precipitation_rate(j,i) + &
    76                                           dqdt_precip * dzw(k) * 0.001
     75                                          dqdt_precip * dzw(k)
    7776
    7877             ENDDO
    7978!
    80 !--          Sum up the precipitation amount (unit kg * 0.001 / m**2)
     79!--          Sum up the precipitation amount, unit kg / m**2 (= mm)
    8180             IF ( intermediate_timestep_count ==         &
    8281                  intermediate_timestep_count_max  .AND. &
  • palm/trunk/SOURCE/check_parameters.f90

    r72 r73  
    88! and pt_reference are checked,
    99! output of precipitation amount/rate and roughnes length + check
    10 ! possible negative humidities are avoided in initial profile
     10! possible negative humidities are avoided in initial profile,
     11! dirichlet/neumann changed to dirichlet/radiation, etc.
    1112!
    1213! Former revisions:
     
    705706!-- Check boundary conditions and set internal variables:
    706707!-- Lateral boundary conditions
    707     IF ( bc_lr /= 'cyclic'  .AND.  bc_lr /= 'dirichlet/neumann'  .AND. &
    708          bc_lr /= 'neumann/dirichlet' )  THEN
     708    IF ( bc_lr /= 'cyclic'  .AND.  bc_lr /= 'dirichlet/radiation'  .AND. &
     709         bc_lr /= 'radiation/dirichlet' )  THEN
    709710       IF ( myid == 0 )  THEN
    710711          PRINT*, '+++ check_parameters:'
     
    713714       CALL local_stop
    714715    ENDIF
    715     IF ( bc_ns /= 'cyclic'  .AND.  bc_ns /= 'dirichlet/neumann'  .AND. &
    716          bc_ns /= 'neumann/dirichlet' )  THEN
     716    IF ( bc_ns /= 'cyclic'  .AND.  bc_ns /= 'dirichlet/radiation'  .AND. &
     717         bc_ns /= 'radiation/dirichlet' )  THEN
    717718       IF ( myid == 0 )  THEN
    718719          PRINT*, '+++ check_parameters:'
     
    22842285    ENDIF
    22852286
    2286     IF ( outflow_l )  THEN
     2287    IF ( bc_lr == 'radiation/dirichlet' )  THEN
    22872288       dist_nxr    = nx - inflow_disturbance_begin
    22882289       dist_nxl(1) = nx - inflow_disturbance_end
    2289     ELSEIF ( outflow_r )  THEN
     2290    ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
    22902291       dist_nxl    = inflow_disturbance_begin
    22912292       dist_nxr(1) = inflow_disturbance_end
    2292     ELSEIF ( outflow_s )  THEN
     2293    ENDIF
     2294    IF ( bc_ns == 'dirichlet/radiation' )  THEN
    22932295       dist_nyn    = ny - inflow_disturbance_begin
    22942296       dist_nys(1) = ny - inflow_disturbance_end
    2295     ELSEIF ( outflow_n )  THEN
     2297    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
    22962298       dist_nys    = inflow_disturbance_begin
    22972299       dist_nyn(1) = inflow_disturbance_end
  • palm/trunk/SOURCE/exchange_horiz_2d.f90

    r4 r73  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary
     7! conditions
    78!
    89! Former revisions:
     
    9899#endif
    99100
     101!
     102!-- Neumann-conditions at inflow/outflow in case of non-cyclic boundary
     103!-- conditions
     104    IF ( inflow_l .OR. outflow_l )  ar(:,nxl-1) = ar(:,nxl)
     105    IF ( inflow_r .OR. outflow_r )  ar(:,nxr+1) = ar(:,nxr)
     106    IF ( inflow_s .OR. outflow_s )  ar(nys-1,:) = ar(nys,:)
     107    IF ( inflow_n .OR. outflow_n )  ar(nyn+1,:) = ar(nyn,:)
     108
    100109    CALL cpu_log( log_point_s(13), 'exchange_horiz_2d', 'stop' )
    101110
  • palm/trunk/SOURCE/init_3d_model.f90

    r72 r73  
    77! Actual revisions:
    88! -----------------
     9! Arrays for radiation boundary conditions are allocated (u_m_l, u_m_r, etc.),
     10! bugfix for cases with the outflow damping layer extending over more than one
     11! subdomain,
    912! New initializing action "by_user" calls user_init_3d_model,
    1013! precipitation_amount/rate, ts_value are allocated, +module netcdf_control,
     
    206209
    207210!
     211!-- Arrays to store velocity data from t-dt needed for radiation boundary
     212!-- conditions
     213    IF ( outflow_l )  THEN
     214       ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1,-2:-1), &
     215                 v_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1,-2:-1), &
     216                 w_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1,-2:-1) )
     217    ENDIF
     218    IF ( outflow_r )  THEN
     219       ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1,-2:-1), &
     220                 v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1,-2:-1), &
     221                 w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1,-2:-1) )
     222    ENDIF
     223    IF ( outflow_s )  THEN
     224       ALLOCATE( u_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1,-2:-1), &
     225                 v_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1,-2:-1), &
     226                 w_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1,-2:-1) )
     227    ENDIF
     228    IF ( outflow_n )  THEN
     229       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1,-2:-1), &
     230                 v_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1,-2:-1), &
     231                 w_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1,-2:-1) )
     232    ENDIF
     233
     234!
    208235!-- Initial assignment of the pointers
    209236    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     
    726753       ENDIF
    727754
     755!
     756!--    Initialize old timelevels needed for radiation boundary conditions
     757       IF ( outflow_l )  THEN
     758          u_m_l(:,:,:,-2) = u(:,:,-1:1)
     759          v_m_l(:,:,:,-2) = v(:,:,-1:1)
     760          w_m_l(:,:,:,-2) = w(:,:,-1:1)
     761          u_m_l(:,:,:,-1) = u(:,:,-1:1)
     762          v_m_l(:,:,:,-1) = v(:,:,-1:1)
     763          w_m_l(:,:,:,-1) = w(:,:,-1:1)
     764       ENDIF
     765       IF ( outflow_r )  THEN
     766          u_m_r(:,:,:,-2) = u(:,:,nx-1:nx+1)
     767          v_m_r(:,:,:,-2) = v(:,:,nx-1:nx+1)
     768          w_m_r(:,:,:,-2) = w(:,:,nx-1:nx+1)
     769          u_m_r(:,:,:,-1) = u(:,:,nx-1:nx+1)
     770          v_m_r(:,:,:,-1) = v(:,:,nx-1:nx+1)
     771          w_m_r(:,:,:,-1) = w(:,:,nx-1:nx+1)
     772       ENDIF
     773       IF ( outflow_s )  THEN
     774          u_m_s(:,:,:,-2) = u(:,-1:1,:)
     775          v_m_s(:,:,:,-2) = v(:,-1:1,:)
     776          w_m_s(:,:,:,-2) = w(:,-1:1,:)
     777          u_m_s(:,:,:,-1) = u(:,-1:1,:)
     778          v_m_s(:,:,:,-1) = v(:,-1:1,:)
     779          w_m_s(:,:,:,-1) = w(:,-1:1,:)
     780       ENDIF
     781       IF ( outflow_n )  THEN
     782          u_m_n(:,:,:,-2) = u(:,ny-1:ny+1,:)
     783          v_m_n(:,:,:,-2) = v(:,ny-1:ny+1,:)
     784          w_m_n(:,:,:,-2) = w(:,ny-1:ny+1,:)
     785          u_m_n(:,:,:,-1) = u(:,ny-1:ny+1,:)
     786          v_m_n(:,:,:,-1) = v(:,ny-1:ny+1,:)
     787          w_m_n(:,:,:,-1) = w(:,ny-1:ny+1,:)
     788       ENDIF
     789
    728790    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' ) &
    729791    THEN
     
    788850!-- non-cyclic lateral boundaries. A linear increase is assumed over the first
    789851!-- half of the width of the damping layer
    790     IF ( bc_lr /= 'cyclic' )  THEN
     852    IF ( bc_lr == 'dirichlet/radiation' )  THEN
    791853
    792854       DO  i = nxl-1, nxr+1
    793 
    794           IF ( outflow_r )  THEN
    795 
    796              IF ( i >= nx - outflow_damping_width )  THEN
    797                 km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
    798                                ( i - ( nx - outflow_damping_width ) ) / &
    799                                REAL( outflow_damping_width/2 )            &
    800                                                 )
    801              ELSE
    802                 km_damp_x(i) = 0.0
    803              ENDIF
    804 
    805           ELSEIF ( outflow_l )  THEN
    806 
    807              IF ( i <= outflow_damping_width )  THEN
    808                 km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
    809                                ( outflow_damping_width - i ) / &
    810                                REAL( outflow_damping_width/2 )            &
    811                                                 )
    812              ELSE
    813                 km_damp_x(i) = 0.0
    814              ENDIF
    815 
    816           ENDIF
    817 
     855          IF ( i >= nx - outflow_damping_width )  THEN
     856             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
     857                            ( i - ( nx - outflow_damping_width ) ) /   &
     858                            REAL( outflow_damping_width/2 )            &
     859                                             )
     860          ELSE
     861             km_damp_x(i) = 0.0
     862          ENDIF
    818863       ENDDO
    819     ENDIF
    820 
    821     IF ( bc_ns /= 'cyclic' )  THEN
     864
     865    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
     866
     867       DO  i = nxl-1, nxr+1
     868          IF ( i <= outflow_damping_width )  THEN
     869             km_damp_x(i) = km_damp_max * MIN( 1.0,                    &
     870                            ( outflow_damping_width - i ) /            &
     871                            REAL( outflow_damping_width/2 )            &
     872                                             )
     873          ELSE
     874             km_damp_x(i) = 0.0
     875          ENDIF
     876       ENDDO
     877
     878    ENDIF
     879
     880    IF ( bc_ns == 'dirichlet/radiation' )  THEN
    822881
    823882       DO  j = nys-1, nyn+1
    824 
    825           IF ( outflow_n )  THEN
    826 
    827              IF ( j >= ny - outflow_damping_width )  THEN
    828                 km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
    829                                ( j - ( ny - outflow_damping_width ) ) / &
    830                                REAL( outflow_damping_width/2 )            &
    831                                                 )
    832              ELSE
    833                 km_damp_y(j) = 0.0
    834              ENDIF
    835 
    836           ELSEIF ( outflow_s )  THEN
    837 
    838              IF ( j <= outflow_damping_width )  THEN
    839                 km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
    840                                ( outflow_damping_width - j ) / &
    841                                REAL( outflow_damping_width/2 )            &
    842                                                 )
    843              ELSE
    844                 km_damp_y(j) = 0.0
    845              ENDIF
    846 
    847           ENDIF
    848 
     883          IF ( j >= ny - outflow_damping_width )  THEN
     884             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
     885                            ( j - ( ny - outflow_damping_width ) ) /   &
     886                            REAL( outflow_damping_width/2 )            &
     887                                             )
     888          ELSE
     889             km_damp_y(j) = 0.0
     890          ENDIF
    849891       ENDDO
     892
     893    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
     894
     895       DO  j = nys-1, nyn+1
     896          IF ( j <= outflow_damping_width )  THEN
     897             km_damp_y(j) = km_damp_max * MIN( 1.0,                    &
     898                            ( outflow_damping_width - j ) /            &
     899                            REAL( outflow_damping_width/2 )            &
     900                                             )
     901          ELSE
     902             km_damp_y(j) = 0.0
     903          ENDIF
     904       ENDDO
     905
    850906    ENDIF
    851907
  • palm/trunk/SOURCE/init_pegrid.f90

    r4 r73  
    44! Actual revisions:
    55! -----------------
    6 !
     6! uxrp, vynp eliminated,
     7! dirichlet/neumann changed to dirichlet/radiation, etc.
    78!
    89! Former revisions:
     
    495496!
    496497!-- For non-cyclic boundaries extend array u (v) by one gridpoint
    497     IF ( bc_lr /= 'cyclic' )  uxrp = 1
    498     IF ( bc_ns /= 'cyclic' )  vynp = 1
     498!    IF ( bc_lr /= 'cyclic' )  uxrp = 1
     499!    IF ( bc_ns /= 'cyclic' )  vynp = 1
    499500
    500501!
     
    743744!-- by one gridpoint on the left/rightmost (northest/southest) processor
    744745    IF ( pleft == MPI_PROC_NULL )  THEN
    745        IF ( bc_lr == 'dirichlet/neumann' )  THEN
     746       IF ( bc_lr == 'dirichlet/radiation' )  THEN
    746747          inflow_l  = .TRUE.
    747        ELSEIF ( bc_lr == 'neumann/dirichlet' )  THEN
     748       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
    748749          outflow_l = .TRUE.
    749750       ENDIF
     
    751752
    752753    IF ( pright == MPI_PROC_NULL )  THEN
    753        IF ( bc_lr == 'dirichlet/neumann' )  THEN
     754       IF ( bc_lr == 'dirichlet/radiation' )  THEN
    754755          outflow_r = .TRUE.
    755        ELSEIF ( bc_lr == 'neumann/dirichlet' )  THEN
     756       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
    756757          inflow_r  = .TRUE.
    757758       ENDIF
    758        uxrp      = 1
     759!       uxrp      = 1
    759760    ENDIF
    760761
    761762    IF ( psouth == MPI_PROC_NULL )  THEN
    762        IF ( bc_ns == 'dirichlet/neumann' )  THEN
     763       IF ( bc_ns == 'dirichlet/radiation' )  THEN
    763764          outflow_s = .TRUE.
    764        ELSEIF ( bc_ns == 'neumann/dirichlet' )  THEN
     765       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
    765766          inflow_s  = .TRUE.
    766767       ENDIF
     
    768769
    769770    IF ( pnorth == MPI_PROC_NULL )  THEN
    770        IF ( bc_ns == 'dirichlet/neumann' )  THEN
     771       IF ( bc_ns == 'dirichlet/radiation' )  THEN
    771772          inflow_n  = .TRUE.
    772        ELSEIF ( bc_ns == 'neumann/dirichlet' )  THEN
     773       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
    773774          outflow_n = .TRUE.
    774775       ENDIF
    775        vynp      = 1
     776!       vynp      = 1
    776777    ENDIF
    777778
     
    787788    ENDIF
    788789#else
    789     IF ( bc_lr == 'dirichlet/neumann' )  THEN
     790    IF ( bc_lr == 'dirichlet/radiation' )  THEN
    790791       inflow_l  = .TRUE.
    791792       outflow_r = .TRUE.
    792        uxrp      = 1
    793     ELSEIF ( bc_lr == 'neumann/dirichlet' )  THEN
     793!       uxrp      = 1
     794    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
    794795       outflow_l = .TRUE.
    795796       inflow_r  = .TRUE.
    796797    ENDIF
    797798
    798     IF ( bc_ns == 'dirichlet/neumann' )  THEN
     799    IF ( bc_ns == 'dirichlet/radiation' )  THEN
    799800       inflow_n  = .TRUE.
    800801       outflow_s = .TRUE.
    801     ELSEIF ( bc_ns == 'neumann/dirichlet' )  THEN
     802    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
    802803       outflow_n = .TRUE.
    803804       inflow_s  = .TRUE.
    804        vynp      = 1
     805!       vynp      = 1
    805806    ENDIF
    806807#endif
  • palm/trunk/SOURCE/modules.f90

    r72 r73  
    66! -----------------
    77! +arrays precipitation_amount, precipitation_rate, precipitation_rate_av,
    8 ! rif_wall, z0_av
     8! rif_wall, z0_av, +arrays u_m_l, u_m_r, etc. for radiation boundary conditions,
    99! +loop_optimization, netcdf_64bit_3d, zu_s_inner, zw_w_inner, id_var_zusi_*,
    1010! id_var_zwwi_*, ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc, pt_reference,
     
    111111          vpt, vpt_m, w, w_m, w_p
    112112
    113     REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall
     113    REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall, u_m_l, u_m_n, u_m_r,   &
     114          u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
    114115                                                   
    115116
  • palm/trunk/SOURCE/pres.f90

    r4 r73  
    396396!--       of non-cyclic boundary conditions). The respective mean velocity
    397397!--       is calculated from this in routine boundary_conds.
    398           IF ( outflow_l  .AND.  i == nxl )  THEN
    399              !$OMP CRITICAL
    400              DO  k = nzb, nzt+1
    401                 uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxl)
    402              ENDDO
    403              !$OMP END CRITICAL
    404           ELSEIF ( outflow_r  .AND.  i == nxr )  THEN
    405              !$OMP CRITICAL
    406              DO  k = nzb, nzt+1
    407                 uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxr)
    408              ENDDO
    409              !$OMP END CRITICAL
    410           ELSEIF ( outflow_s  .AND.  j == nys )  THEN
    411              !$OMP CRITICAL
    412              DO  k = nzb, nzt+1
    413                 uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nys,i)
    414              ENDDO
    415              !$OMP END CRITICAL
    416           ELSEIF ( outflow_n  .AND.  j == nyn )  THEN
    417              !$OMP CRITICAL
    418              DO  k = nzb, nzt+1
    419                 uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nyn,i)
    420              ENDDO
    421              !$OMP END CRITICAL
    422           ENDIF
     398!          IF ( outflow_l  .AND.  i == nxl )  THEN
     399!             !$OMP CRITICAL
     400!             DO  k = nzb, nzt+1
     401!                uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxl)
     402!             ENDDO
     403!             !$OMP END CRITICAL
     404!          ELSEIF ( outflow_r  .AND.  i == nxr )  THEN
     405!             !$OMP CRITICAL
     406!             DO  k = nzb, nzt+1
     407!                uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxr)
     408!             ENDDO
     409!             !$OMP END CRITICAL
     410!          ELSEIF ( outflow_s  .AND.  j == nys )  THEN
     411!             !$OMP CRITICAL
     412!             DO  k = nzb, nzt+1
     413!                uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nys,i)
     414!             ENDDO
     415!             !$OMP END CRITICAL
     416!          ELSEIF ( outflow_n  .AND.  j == nyn )  THEN
     417!             !$OMP CRITICAL
     418!             DO  k = nzb, nzt+1
     419!                uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nyn,i)
     420!             ENDDO
     421!             !$OMP END CRITICAL
     422!          ENDIF
    423423
    424424!
  • palm/trunk/SOURCE/production_e.f90

    r72 r73  
    910910!--                of the eddy diffusivity by km = u* * kappa * zp / phi_m.
    911911          !$OMP PARALLEL DO PRIVATE( ku, kv )
    912           DO  i = nxl, nxr+1
    913              DO  j = nys, nyn+1
     912          DO  i = nxl, nxr
     913             DO  j = nys, nyn
    914914
    915915                ku = nzb_u_inner(j,i)+1
  • palm/trunk/SOURCE/prognostic_equations.f90

    r70 r73  
    44! Actual revisions:
    55! -----------------
     6! checking for negative q and limiting for positive values,
    67! z0 removed from arguments in calls of diffusion_u/v/w,
    78! subroutine names changed to .._noopt, .._cache, and .._vector
     
    481482                                     ) -                                       &
    482483                             tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
     484                IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
    483485             ENDDO
    484486
     
    915917                                        ) -                                    &
    916918                                tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
     919                   IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
    917920                ENDDO
    918921
     
    14191422                                     ) -                                       &
    14201423                             tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
     1424                IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
    14211425             ENDDO
    14221426          ENDDO
  • palm/trunk/SOURCE/read_3d_binary.f90

    r72 r73  
    44! Actual revisions:
    55! -----------------
    6 ! +precipitation_amount, precipitation_rate_av, rif_wall, z0_av
     6! +precipitation_amount, precipitation_rate_av, rif_wall, u_m_l, u_m_r, etc.,
     7! z0_av
    78!
    89! Former revisions:
     
    302303          CASE ( 'u_m' )
    303304             READ ( 13 )  u_m
     305          CASE ( 'u_m_l' )
     306             READ ( 13 )  u_m_l
     307          CASE ( 'u_m_n' )
     308             READ ( 13 )  u_m_n
     309          CASE ( 'u_m_r' )
     310             READ ( 13 )  u_m_r
     311          CASE ( 'u_m_s' )
     312             READ ( 13 )  u_m_s
    304313          CASE ( 'us' )
    305314             READ ( 13 )  us
     
    322331          CASE ( 'v_m' )
    323332             READ (13 )   v_m
     333          CASE ( 'v_m_l' )
     334             READ ( 13 )  v_m_l
     335          CASE ( 'v_m_n' )
     336             READ ( 13 )  v_m_n
     337          CASE ( 'v_m_r' )
     338             READ ( 13 )  v_m_r
     339          CASE ( 'v_m_s' )
     340             READ ( 13 )  v_m_s
    324341          CASE ( 'vpt' )
    325342             READ ( 13 )  vpt
     
    340357          CASE ( 'w_m' )
    341358             READ ( 13 )  w_m
     359          CASE ( 'w_m_l' )
     360             READ ( 13 )  w_m_l
     361          CASE ( 'w_m_n' )
     362             READ ( 13 )  w_m_n
     363          CASE ( 'w_m_r' )
     364             READ ( 13 )  w_m_r
     365          CASE ( 'w_m_s' )
     366             READ ( 13 )  w_m_s
    342367          CASE ( 'z0' )
    343368             READ ( 13 )  z0
  • palm/trunk/SOURCE/time_integration.f90

    r63 r73  
    158158
    159159!
    160 !--       Swap the time levels in preparation for the next time step.
    161           CALL swap_timelevel
    162 
    163 !
    164160!--       Boundary conditions for the prognostic quantities (except of the
    165161!--       velocities at the outflow in case of a non-cyclic lateral wall)
    166162          CALL boundary_conds( 'main' )
     163
     164!
     165!--       Swap the time levels in preparation for the next time step.
     166          CALL swap_timelevel
    167167
    168168!
  • palm/trunk/SOURCE/write_3d_binary.f90

    r72 r73  
    44! Actual revisions:
    55! -----------------
    6 ! +precipitation_amount, precipitation_rate_av, rif_wall, z0_av
     6! +precipitation_amount, precipitation_rate_av, rif_wall, u_m_l, u_m_r, etc.,
     7! z0_av
    78!
    89! Former revisions:
     
    169170    ENDIF
    170171    WRITE ( 14 )  'u_m                 ';  WRITE ( 14 )  u_m
     172    IF ( ALLOCATED( u_m_l ) )  THEN
     173       WRITE ( 14 )  'u_m_l               ';  WRITE ( 14 )  u_m_l
     174    ENDIF
     175    IF ( ALLOCATED( u_m_n ) )  THEN
     176       WRITE ( 14 )  'u_m_n               ';  WRITE ( 14 )  u_m_n
     177    ENDIF
     178    IF ( ALLOCATED( u_m_r ) )  THEN
     179       WRITE ( 14 )  'u_m_r               ';  WRITE ( 14 )  u_m_r
     180    ENDIF
     181    IF ( ALLOCATED( u_m_s ) )  THEN
     182       WRITE ( 14 )  'u_m_s               ';  WRITE ( 14 )  u_m_s
     183    ENDIF
    171184    WRITE ( 14 )  'us                  ';  WRITE ( 14 )  us
    172185    WRITE ( 14 )  'usws                ';  WRITE ( 14 )  usws
     
    184197    ENDIF
    185198    WRITE ( 14 )  'v_m                 ';  WRITE ( 14 )  v_m
     199    IF ( ALLOCATED( v_m_l ) )  THEN
     200       WRITE ( 14 )  'v_m_l               ';  WRITE ( 14 )  v_m_l
     201    ENDIF
     202    IF ( ALLOCATED( v_m_n ) )  THEN
     203       WRITE ( 14 )  'v_m_n               ';  WRITE ( 14 )  v_m_n
     204    ENDIF
     205    IF ( ALLOCATED( v_m_r ) )  THEN
     206       WRITE ( 14 )  'v_m_r               ';  WRITE ( 14 )  v_m_r
     207    ENDIF
     208    IF ( ALLOCATED( v_m_s ) )  THEN
     209       WRITE ( 14 )  'v_m_s               ';  WRITE ( 14 )  v_m_s
     210    ENDIF
    186211    IF ( moisture )  THEN
    187212       WRITE ( 14 )  'vpt                 ';  WRITE ( 14 )  vpt
     
    202227    ENDIF
    203228    WRITE ( 14 )  'w_m                 ';  WRITE ( 14 )  w_m
     229    IF ( ALLOCATED( w_m_l ) )  THEN
     230       WRITE ( 14 )  'w_m_l               ';  WRITE ( 14 )  w_m_l
     231    ENDIF
     232    IF ( ALLOCATED( w_m_n ) )  THEN
     233       WRITE ( 14 )  'w_m_n               ';  WRITE ( 14 )  w_m_n
     234    ENDIF
     235    IF ( ALLOCATED( w_m_r ) )  THEN
     236       WRITE ( 14 )  'w_m_r               ';  WRITE ( 14 )  w_m_r
     237    ENDIF
     238    IF ( ALLOCATED( w_m_s ) )  THEN
     239       WRITE ( 14 )  'w_m_s               ';  WRITE ( 14 )  w_m_s
     240    ENDIF
    204241    WRITE ( 14 )  'z0                  ';  WRITE ( 14 )  z0
    205242    IF ( ALLOCATED( z0_av ) )  THEN
Note: See TracChangeset for help on using the changeset viewer.