Changeset 106 for palm/trunk


Ignore:
Timestamp:
Aug 16, 2007 2:30:26 PM (14 years ago)
Author:
raasch
Message:

preliminary update of bugfixes and extensions for non-cyclic BCs

Location:
palm/trunk/SOURCE
Files:
22 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/CURRENT_MODIFICATIONS

    r105 r106  
    1313inipar parameter top_momentumflux_u|v.
    1414
    15 boundary_conds, check_open, check_parameters, diffusion_u, diffusion_v, flow_statistics, header, init_pegrid, init_3d_model, modules, palm, parin, prognostic_equations, read_var_list, read_3d_binary, swap_timelevel, time_integration, write_var_list, write_3d_binary
     15Non-cyclic boundary conditions can be used along all horizontal directions.
     16
     17Quantities w*p* and w"e can be output as vertical profiles.
     18
     19boundary_conds, check_open, check_parameters, diffusion_u, diffusion_v, flow_statistics, header, init_pegrid, init_3d_model, modules, palm, parin, pres, prognostic_equations, read_var_list, read_3d_binary, swap_timelevel, time_integration, write_var_list, write_3d_binary
    1620
    1721New:
     
    2125Changed:
    2226-------
     27Remaining variables iran changed to iran_part (advec_particles, init_particles).
     28
     29In case that the presure solver is not called for every Runge-Kutta substep
     30(call_psolver_at_all_substeps = .F.), it is called after the first substep
     31instead of the last. In that case, random perturbations are also added to the
     32velocity field after the first substep.
     33
     34advec_particles, init_particles, time_integration
    2335
    2436
    2537Errors:
    2638------
    27 +dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface.f90)
     39Bugs from code parts for non-cyclic boundary conditions are removed: loops for
     40u and v are starting from index nxlu, nysv, respectively. The radiation boundary
     41condition is used for every Runge-Kutta substep. Velocity phase speeds for
     42the radiation boundary conditions are calculated for the first Runge-Kutta
     43substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
     44are defined for this purpose. Several index errors are removed from the
     45radiation boundary condition code parts. Upper bounds for calculating
     46u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
     47values are not available in case of non-cyclic boundary conditions.
     48
     49+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)
     50
     51Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)
    2852
    2953Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).
    3054
    31 check_parameters, user_interface
     55advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, boundary_conds, buoyancy, check_parameters, coriolis, diffusion_u, diffusion_v, init_pegrid, init_3d_model, modules, production_e, prognostic_equations, user_interface
  • palm/trunk/SOURCE/advec_particles.f90

    r98 r106  
    44! Actual revisions:
    55! -----------------
    6 !
     6! remaining variables iran changed to iran_part
    77! TEST: PRINT statements on unit 9 (commented out)
    88!
     
    19551955                   THEN
    19561956                      particles(n)%x = particles(n)%x + &
    1957                                        ( random_function( iran ) - 0.5 ) * &
     1957                                       ( random_function( iran_part ) - 0.5 ) *&
    19581958                                       pdx(particles(n)%group)
    19591959                      IF ( particles(n)%x  <=  ( nxl - 0.5 ) * dx )  THEN
     
    19661966                   THEN
    19671967                      particles(n)%y = particles(n)%y + &
    1968                                        ( random_function( iran ) - 0.5 ) * &
     1968                                       ( random_function( iran_part ) - 0.5 ) *&
    19691969                                       pdy(particles(n)%group)
    19701970                      IF ( particles(n)%y  <=  ( nys - 0.5 ) * dy )  THEN
     
    19771977                   THEN
    19781978                      particles(n)%z = particles(n)%z + &
    1979                                        ( random_function( iran ) - 0.5 ) * &
     1979                                       ( random_function( iran_part ) - 0.5 ) *&
    19801980                                       pdz(particles(n)%group)
    19811981                   ENDIF
  • palm/trunk/SOURCE/advec_u_pw.f90

    r77 r106  
    44! Actual revisions:
    55! -----------------
    6 !
     6! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
    77!
    88! Former revisions:
     
    5858       gu = 2.0 * u_gtrans
    5959       gv = 2.0 * v_gtrans
    60        DO  i = nxl, nxr
     60       DO  i = nxlu, nxr
    6161          DO  j = nys, nyn
    6262             DO  k = nzb_u_inner(j,i)+1, nzt
  • palm/trunk/SOURCE/advec_u_up.f90

    r77 r106  
    44! Actual revisions:
    55! -----------------
    6 !
     6! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
    77!
    88! Former revisions:
     
    5757
    5858
    59        DO  i = nxl, nxr
     59       DO  i = nxlu, nxr
    6060          DO  j = nys, nyn
    6161             DO  k = nzb_u_inner(j,i)+1, nzt
  • palm/trunk/SOURCE/advec_v_pw.f90

    r77 r106  
    44! Actual revisions:
    55! -----------------
    6 !
     6! j loop is starting from nysv (needed for non-cyclic boundary conditions)
    77!
    88! Former revisions:
     
    6060       gv = 2.0 * v_gtrans
    6161       DO  i = nxl, nxr
    62           DO  j = nys, nyn
     62          DO  j = nysv, nyn
    6363             DO  k = nzb_v_inner(j,i)+1, nzt
    6464                tend(k,j,i) = tend(k,j,i) - 0.25 * (                           &
  • palm/trunk/SOURCE/advec_v_up.f90

    r77 r106  
    44! Actual revisions:
    55! -----------------
    6 !
     6! j loop is starting from nysv (needed for non-cyclic boundary conditions)
    77!
    88! Former revisions:
     
    5757
    5858       DO  i = nxl, nxr
    59           DO  j = nys, nyn
     59          DO  j = nysv, nyn
    6060             DO  k = nzb_v_inner(j,i)+1, nzt
    6161!
  • palm/trunk/SOURCE/boundary_conds.f90

    r102 r106  
    44! Actual revisions:
    55! -----------------
    6 ! Boundary conditions for temperature adjusted for coupled runs
     6! Boundary conditions for temperature adjusted for coupled runs,
     7! bugfixes for the radiation boundary conditions at the outflow: radiation
     8! conditions are used for every substep, phase speeds are calculated for the
     9! first Runge-Kutta substep only and then reused, several index values changed
     10!
    711!
    812! Former revisions:
     
    5963    INTEGER ::  i, j, k
    6064
    61     REAL    ::  c_max, c_u, c_v, c_w, denom
     65    REAL    ::  c_max, denom
    6266
    6367
     
    216220!
    217221!-- Radiation boundary condition for the velocities at the respective outflow
    218     IF ( outflow_s  .AND.                                                 &
    219          intermediate_timestep_count == intermediate_timestep_count_max ) &
    220     THEN
     222    IF ( outflow_s )  THEN
    221223
    222224       c_max = dy / dt_3d
     
    226228
    227229!
    228 !--          First calculate the phase speeds for u,v, and w
    229              denom = u_m_s(k,0,i) - u_m_s(k,1,i)
    230 
    231              IF ( denom /= 0.0 )  THEN
    232                 c_u = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom
    233                 IF ( c_u < 0.0 )  THEN
    234                    c_u = 0.0
    235                 ELSEIF ( c_u > c_max )  THEN
    236                    c_u = c_max
    237                 ENDIF
    238              ELSE
    239                 c_u = c_max
     230!--          Calculate the phase speeds for u,v, and w. In case of using a
     231!--          Runge-Kutta scheme, do this for the first substep only and then
     232!--          reuse this values for the further substeps.
     233             IF ( intermediate_timestep_count == 1 )  THEN
     234
     235                denom = u_m_s(k,0,i) - u_m_s(k,1,i)
     236
     237                IF ( denom /= 0.0 )  THEN
     238                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom
     239                   IF ( c_u(k,i) < 0.0 )  THEN
     240                      c_u(k,i) = 0.0
     241                   ELSEIF ( c_u(k,i) > c_max )  THEN
     242                      c_u(k,i) = c_max
     243                   ENDIF
     244                ELSE
     245                   c_u(k,i) = c_max
     246                ENDIF
     247
     248                denom = v_m_s(k,1,i) - v_m_s(k,2,i)
     249
     250                IF ( denom /= 0.0 )  THEN
     251                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / denom
     252                   IF ( c_v(k,i) < 0.0 )  THEN
     253                      c_v(k,i) = 0.0
     254                   ELSEIF ( c_v(k,i) > c_max )  THEN
     255                      c_v(k,i) = c_max
     256                   ENDIF
     257                ELSE
     258                   c_v(k,i) = c_max
     259                ENDIF
     260
     261                denom = w_m_s(k,0,i) - w_m_s(k,1,i)
     262
     263                IF ( denom /= 0.0 )  THEN
     264                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
     265                   IF ( c_w(k,i) < 0.0 )  THEN
     266                      c_w(k,i) = 0.0
     267                   ELSEIF ( c_w(k,i) > c_max )  THEN
     268                      c_w(k,i) = c_max
     269                   ENDIF
     270                ELSE
     271                   c_w(k,i) = c_max
     272                ENDIF
     273
     274!
     275!--             Save old timelevels for the next timestep
     276                u_m_s(k,:,i) = u(k,0:1,i)
     277                v_m_s(k,:,i) = v(k,1:2,i)
     278                w_m_s(k,:,i) = w(k,0:1,i)
     279
    240280             ENDIF
    241              denom = v_m_s(k,0,i) - v_m_s(k,1,i)
    242 
    243              IF ( denom /= 0.0 )  THEN
    244                 c_v = -c_max * ( v(k,0,i) - v_m_s(k,0,i) ) / denom
    245                 IF ( c_v < 0.0 )  THEN
    246                    c_v = 0.0
    247                 ELSEIF ( c_v > c_max )  THEN
    248                    c_v = c_max
    249                 ENDIF
    250              ELSE
    251                 c_v = c_max
    252              ENDIF
    253 
    254              denom = w_m_s(k,0,i) - w_m_s(k,1,i)
    255 
    256              IF ( denom /= 0.0 )  THEN
    257                 c_w = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
    258                 IF ( c_w < 0.0 )  THEN
    259                    c_w = 0.0
    260                 ELSEIF ( c_w > c_max )  THEN
    261                    c_w = c_max
    262                 ENDIF
    263              ELSE
    264                 c_w = c_max
    265              ENDIF
    266281
    267282!
    268283!--          Calculate the new velocities
    269              u_p(k,-1,i) = u(k,-1,i) + dt_3d * c_u * &
     284             u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u(k,i) * &
    270285                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
    271286
    272              v_p(k,-1,i) = v(k,-1,i) + dt_3d * c_v * &
    273                                        ( v(k,-1,i) - v_m_s(k,0,i) ) * ddy
    274 
    275              w_p(k,-1,i) = w(k,-1,i) + dt_3d * c_w * &
     287             v_p(k,-1,i) = v(k,-1,i) - dt_3d * tsc(2) * c_v(k,i) * &
     288                                       ( v(k,0,i) - v(k,1,i) ) * ddy
     289
     290             w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w(k,i) * &
    276291                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
    277 
    278 !
    279 !--          Save old timelevels for the next timestep
    280              u_m_s(k,:,i) = u(k,-1:1,i)
    281              v_m_s(k,:,i) = v(k,-1:1,i)
    282              w_m_s(k,:,i) = w(k,-1:1,i)
    283292
    284293          ENDDO
     
    289298       IF ( ibc_uv_b == 0 )  THEN
    290299          u_p(nzb,-1,:) = -u_p(nzb+1,-1,:)
    291           v_p(nzb,-1,:) = -v_p(nzb+1,-1,:) 
     300          v_p(nzb,0,:)  = -v_p(nzb+1,0,:) 
    292301       ELSE                   
    293302          u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
    294           v_p(nzb,-1,:) =  v_p(nzb+1,-1,:)
    295        ENDIF
    296        w_p(nzb,ny+1,:) = 0.0
     303          v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
     304       ENDIF
     305       w_p(nzb,-1,:) = 0.0
    297306
    298307!
     
    300309       IF ( ibc_uv_t == 0 )  THEN
    301310          u_p(nzt+1,-1,:) = ug(nzt+1)
    302           v_p(nzt+1,-1,:) = vg(nzt+1)
     311          v_p(nzt+1,0,:) = vg(nzt+1)
    303312       ELSE
    304313          u_p(nzt+1,-1,:) = u(nzt,-1,:)
    305           v_p(nzt+1,-1,:) = v(nzt,-1,:)
     314          v_p(nzt+1,0,:)  = v(nzt,0,:)
    306315       ENDIF
    307316       w_p(nzt:nzt+1,-1,:) = 0.0
     
    309318    ENDIF
    310319
    311     IF ( outflow_n  .AND.                                                 &
    312          intermediate_timestep_count == intermediate_timestep_count_max ) &
    313     THEN
     320    IF ( outflow_n )  THEN
    314321
    315322       c_max = dy / dt_3d
     
    319326
    320327!
    321 !--          First calculate the phase speeds for u,v, and w
    322              denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
    323 
    324              IF ( denom /= 0.0 )  THEN
    325                 c_u = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom
    326                 IF ( c_u < 0.0 )  THEN
    327                    c_u = 0.0
    328                 ELSEIF ( c_u > c_max )  THEN
    329                    c_u = c_max
    330                 ENDIF
    331              ELSE
    332                 c_u = c_max
     328!--          Calculate the phase speeds for u,v, and w. In case of using a
     329!--          Runge-Kutta scheme, do this for the first substep only and then
     330!--          reuse this values for the further substeps.
     331             IF ( intermediate_timestep_count == 1 )  THEN
     332
     333                denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
     334
     335                IF ( denom /= 0.0 )  THEN
     336                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom
     337                   IF ( c_u(k,i) < 0.0 )  THEN
     338                      c_u(k,i) = 0.0
     339                   ELSEIF ( c_u(k,i) > c_max )  THEN
     340                      c_u(k,i) = c_max
     341                   ENDIF
     342                ELSE
     343                   c_u(k,i) = c_max
     344                ENDIF
     345
     346                denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
     347
     348                IF ( denom /= 0.0 )  THEN
     349                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
     350                   IF ( c_v(k,i) < 0.0 )  THEN
     351                      c_v(k,i) = 0.0
     352                   ELSEIF ( c_v(k,i) > c_max )  THEN
     353                      c_v(k,i) = c_max
     354                   ENDIF
     355                ELSE
     356                   c_v(k,i) = c_max
     357                ENDIF
     358
     359                denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
     360
     361                IF ( denom /= 0.0 )  THEN
     362                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
     363                   IF ( c_w(k,i) < 0.0 )  THEN
     364                      c_w(k,i) = 0.0
     365                   ELSEIF ( c_w(k,i) > c_max )  THEN
     366                      c_w(k,i) = c_max
     367                   ENDIF
     368                ELSE
     369                   c_w(k,i) = c_max
     370                ENDIF
     371
     372!
     373!--             Swap timelevels for the next timestep
     374                u_m_n(k,:,i) = u(k,ny-1:ny,i)
     375                v_m_n(k,:,i) = v(k,ny-1:ny,i)
     376                w_m_n(k,:,i) = w(k,ny-1:ny,i)
     377
    333378             ENDIF
    334379
    335              denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
    336 
    337              IF ( denom /= 0.0 )  THEN
    338                 c_v = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
    339                 IF ( c_v < 0.0 )  THEN
    340                    c_v = 0.0
    341                 ELSEIF ( c_v > c_max )  THEN
    342                    c_v = c_max
    343                 ENDIF
    344              ELSE
    345                 c_v = c_max
    346              ENDIF
    347 
    348              denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
    349 
    350              IF ( denom /= 0.0 )  THEN
    351                 c_w = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
    352                 IF ( c_w < 0.0 )  THEN
    353                    c_w = 0.0
    354                 ELSEIF ( c_w > c_max )  THEN
    355                    c_w = c_max
    356                 ENDIF
    357              ELSE
    358                 c_w = c_max
    359              ENDIF
    360 
    361380!
    362381!--          Calculate the new velocities
    363              u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * c_u * &
     382             u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u(k,i) * &
    364383                                           ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
    365384
    366              v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * c_v * &
     385             v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v(k,i) * &
    367386                                           ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
    368387
    369              w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * c_w * &
     388             w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w(k,i) * &
    370389                                           ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
    371 
    372 !
    373 !--          Swap timelevels for the next timestep
    374              u_m_n(k,:,i) = u(k,ny-1:ny+1,i)
    375              v_m_n(k,:,i) = v(k,ny-1:ny+1,i)
    376              w_m_n(k,:,i) = w(k,ny-1:ny+1,i)
    377390
    378391          ENDDO
     
    403416    ENDIF
    404417
    405     IF ( outflow_l  .AND.                                                 &
    406          intermediate_timestep_count == intermediate_timestep_count_max ) &
    407     THEN
     418    IF ( outflow_l )  THEN
    408419
    409420       c_max = dx / dt_3d
     
    413424
    414425!
    415 !--          First calculate the phase speeds for u,v, and w
    416              denom = u_m_l(k,j,0) - u_m_l(k,j,1)
    417 
    418              IF ( denom /= 0.0 )  THEN
    419                 c_u = -c_max * ( u(k,j,0) - u_m_r(k,j,0) ) / denom
    420                 IF ( c_u > 0.0 )  THEN
    421                    c_u = 0.0
    422                 ELSEIF ( c_u < -c_max )  THEN
    423                    c_u = -c_max
    424                 ENDIF
    425              ELSE
    426                 c_u = -c_max
     426!--          Calculate the phase speeds for u,v, and w. In case of using a
     427!--          Runge-Kutta scheme, do this for the first substep only and then
     428!--          reuse this values for the further substeps.
     429             IF ( intermediate_timestep_count == 1 )  THEN
     430
     431                denom = u_m_l(k,j,1) - u_m_l(k,j,2)
     432
     433                IF ( denom /= 0.0 )  THEN
     434                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / denom
     435                   IF ( c_u(k,j) > 0.0 )  THEN
     436                      c_u(k,j) = 0.0
     437                   ELSEIF ( c_u(k,j) < -c_max )  THEN
     438                      c_u(k,j) = -c_max
     439                   ENDIF
     440                ELSE
     441                   c_u(k,j) = -c_max
     442                ENDIF
     443
     444                denom = v_m_l(k,j,0) - v_m_l(k,j,1)
     445
     446                IF ( denom /= 0.0 )  THEN
     447                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
     448                   IF ( c_v(k,j) < 0.0 )  THEN
     449                      c_v(k,j) = 0.0
     450                   ELSEIF ( c_v(k,j) > c_max )  THEN
     451                      c_v(k,j) = c_max
     452                   ENDIF
     453                ELSE
     454                   c_v(k,j) = c_max
     455                ENDIF
     456
     457                denom = w_m_l(k,j,0) - w_m_l(k,j,1)
     458
     459                IF ( denom /= 0.0 )  THEN
     460                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
     461                   IF ( c_w(k,j) < 0.0 )  THEN
     462                      c_w(k,j) = 0.0
     463                   ELSEIF ( c_w(k,j) > c_max )  THEN
     464                      c_w(k,j) = c_max
     465                   ENDIF
     466                ELSE
     467                   c_w(k,j) = c_max
     468                ENDIF
     469
     470!
     471!--             Swap timelevels for the next timestep
     472                u_m_l(k,j,:) = u(k,j,1:2)
     473                v_m_l(k,j,:) = v(k,j,0:1)
     474                w_m_l(k,j,:) = w(k,j,0:1)
     475
    427476             ENDIF
    428477
    429              denom = v_m_l(k,j,0) - v_m_l(k,j,1)
    430 
    431              IF ( denom /= 0.0 )  THEN
    432                 c_v = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
    433                 IF ( c_v < 0.0 )  THEN
    434                    c_v = 0.0
    435                 ELSEIF ( c_v > c_max )  THEN
    436                    c_v = c_max
    437                 ENDIF
    438              ELSE
    439                 c_v = c_max
    440              ENDIF
    441 
    442              denom = w_m_l(k,j,0) - w_m_l(k,j,1)
    443 
    444              IF ( denom /= 0.0 )  THEN
    445                 c_w = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
    446                 IF ( c_w < 0.0 )  THEN
    447                    c_w = 0.0
    448                 ELSEIF ( c_w > c_max )  THEN
    449                    c_w = c_max
    450                 ENDIF
    451              ELSE
    452                 c_w = c_max
    453              ENDIF
    454 
    455478!
    456479!--          Calculate the new velocities
    457              u_p(k,j,-1) = u(k,j,-1) + dt_3d * c_u * &
    458                                        ( u(k,j,-1) - u(k,j,0) ) * ddx
    459 
    460              v_p(k,j,-1) = v(k,j,-1) + dt_3d * c_v * &
     480             u_p(k,j,0)  = u(k,j,0)  - dt_3d * tsc(2) * c_u(k,j) * &
     481                                       ( u(k,j,0) - u(k,j,1) ) * ddx
     482
     483             v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v(k,j) * &
    461484                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
    462485
    463              w_p(k,j,-1) = w(k,j,-1) + dt_3d * c_w * &
     486             w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w(k,j) * &
    464487                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
    465 
    466 !
    467 !--          Swap timelevels for the next timestep
    468              u_m_l(k,j,:) = u(k,j,-1:1)
    469              v_m_l(k,j,:) = v(k,j,-1:1)
    470              w_m_l(k,j,:) = w(k,j,-1:1)
    471488
    472489          ENDDO
     
    497514    ENDIF
    498515
    499     IF ( outflow_r  .AND.                                                 &
    500          intermediate_timestep_count == intermediate_timestep_count_max ) &
    501     THEN
     516    IF ( outflow_r )  THEN
    502517
    503518       c_max = dx / dt_3d
     
    507522
    508523!
    509 !--          First calculate the phase speeds for u,v, and w
    510              denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
    511 
    512              IF ( denom /= 0.0 )  THEN
    513                 c_u = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom
    514                 IF ( c_u < 0.0 )  THEN
    515                    c_u = 0.0
    516                 ELSEIF ( c_u > c_max )  THEN
    517                    c_u = c_max
    518                 ENDIF
    519              ELSE
    520                 c_u = c_max
     524!--          Calculate the phase speeds for u,v, and w. In case of using a
     525!--          Runge-Kutta scheme, do this for the first substep only and then
     526!--          reuse this values for the further substeps.
     527             IF ( intermediate_timestep_count == 1 )  THEN
     528
     529                denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
     530
     531                IF ( denom /= 0.0 )  THEN
     532                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom
     533                   IF ( c_u(k,j) < 0.0 )  THEN
     534                      c_u(k,j) = 0.0
     535                   ELSEIF ( c_u(k,j) > c_max )  THEN
     536                      c_u(k,j) = c_max
     537                   ENDIF
     538                ELSE
     539                   c_u(k,j) = c_max
     540                ENDIF
     541
     542                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
     543
     544                IF ( denom /= 0.0 )  THEN
     545                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
     546                   IF ( c_v(k,j) < 0.0 )  THEN
     547                      c_v(k,j) = 0.0
     548                   ELSEIF ( c_v(k,j) > c_max )  THEN
     549                      c_v(k,j) = c_max
     550                   ENDIF
     551                ELSE
     552                   c_v(k,j) = c_max
     553                ENDIF
     554
     555                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
     556
     557                IF ( denom /= 0.0 )  THEN
     558                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
     559                   IF ( c_w(k,j) < 0.0 )  THEN
     560                      c_w(k,j) = 0.0
     561                   ELSEIF ( c_w(k,j) > c_max )  THEN
     562                      c_w(k,j) = c_max
     563                   ENDIF
     564                ELSE
     565                   c_w(k,j) = c_max
     566                ENDIF
     567
     568!
     569!--             Swap timelevels for the next timestep
     570                u_m_r(k,j,:) = u(k,j,nx-1:nx)
     571                v_m_r(k,j,:) = v(k,j,nx-1:nx)
     572                w_m_r(k,j,:) = w(k,j,nx-1:nx)
     573
    521574             ENDIF
    522575
    523              denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
    524 
    525              IF ( denom /= 0.0 )  THEN
    526                 c_v = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
    527                 IF ( c_v < 0.0 )  THEN
    528                    c_v = 0.0
    529                 ELSEIF ( c_v > c_max )  THEN
    530                    c_v = c_max
    531                 ENDIF
    532              ELSE
    533                 c_v = c_max
    534              ENDIF
    535 
    536              denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
    537 
    538              IF ( denom /= 0.0 )  THEN
    539                 c_w = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
    540                 IF ( c_w < 0.0 )  THEN
    541                    c_w = 0.0
    542                 ELSEIF ( c_w > c_max )  THEN
    543                    c_w = c_max
    544                 ENDIF
    545              ELSE
    546                 c_w = c_max
    547              ENDIF
    548 
    549576!
    550577!--          Calculate the new velocities
    551              u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * c_u * &
     578             u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u(k,j) * &
    552579                                           ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
    553580
    554              v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * c_v * &
     581             v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v(k,j) * &
    555582                                           ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
    556583
    557              w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * c_w * &
     584             w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w(k,j) * &
    558585                                           ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
    559 
    560 !
    561 !--          Swap timelevels for the next timestep
    562              u_m_r(k,j,:) = u(k,j,nx-1:nx+1)
    563              v_m_r(k,j,:) = v(k,j,nx-1:nx+1)
    564              w_m_r(k,j,:) = w(k,j,nx-1:nx+1)
    565586
    566587          ENDDO
  • palm/trunk/SOURCE/buoyancy.f90

    r98 r106  
    44! Actual revisions:
    55! -----------------
    6 !
     6! i loop for u-component (sloping surface) is starting from nxlu (needed for
     7! non-cyclic boundary conditions)
    78!
    89! Former revisions:
     
    105106          IF ( wind_component == 1 )  THEN
    106107
    107              DO  i = nxl, nxr
     108             DO  i = nxlu, nxr
    108109                DO  j = nys, nyn
    109110                   DO  k = nzb_s_inner(j,i)+1, nzt-1
  • palm/trunk/SOURCE/check_parameters.f90

    r104 r106  
    55! -----------------
    66! Check coupling_mode and set default (obligatory) values (like boundary
    7 ! conditions for temperature and fluxes) in case of coupled runs
     7! conditions for temperature and fluxes) in case of coupled runs.
     8! +profiles for w*p* and w"e
    89! Bugfix: Error message concerning output of particle concentration (pc)
    910! modified
     
    19681969             dopr_index(i) = 56
    19691970             dopr_unit(i)  = 'm2/s3'
    1970              hom(:,2,56,:) = SPREAD( zu, 2, statistic_regions+1 )
     1971             hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 )
    19711972
    19721973          CASE ( 'w"e/dz' )
     
    20512052                hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
    20522053             ENDIF
     2054
     2055          CASE ( 'w*p*' )
     2056             dopr_index(i) = 68
     2057             dopr_unit(i)  = 'm3/s3'
     2058             hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 )
     2059
     2060          CASE ( 'w"e' )
     2061             dopr_index(i) = 69
     2062             dopr_unit(i)  = 'm3/s3'
     2063             hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 )
    20532064
    20542065
  • palm/trunk/SOURCE/coriolis.f90

    r77 r106  
    44! Actual revisions:
    55! -----------------
    6 !
     6! loops for u and v are starting from index nxlu, nysv, respectively (needed
     7! for non-cyclic boundary conditions)
    78!
    89! Former revisions:
     
    6061!--       u-component
    6162          CASE ( 1 )
    62              DO  i = nxl, nxr
     63             DO  i = nxlu, nxr
    6364                DO  j = nys, nyn
    6465                   DO  k = nzb_u_inner(j,i)+1, nzt
     
    7879          CASE ( 2 )
    7980             DO  i = nxl, nxr
    80                 DO  j = nys, nyn
     81                DO  j = nysv, nyn
    8182                   DO  k = nzb_v_inner(j,i)+1, nzt
    8283                      tend(k,j,i) = tend(k,j,i) - f *     ( 0.25 *          &
  • palm/trunk/SOURCE/data_output_dvrp.f90

    r86 r106  
    151151!
    152152!--       Definition of characteristics of particle material
    153           CALL DVRP_MATERIAL_RGB( m-1, 1, 0.1, 0.7, 0.1, 0.0 )
     153!          CALL DVRP_MATERIAL_RGB( m-1, 1, 0.1, 0.7, 0.1, 0.0 )
     154          CALL DVRP_MATERIAL_RGB( m-1, 1, 0.0, 0.0, 0.0, 0.0 )
    154155
    155156!
     
    318319                   ENDDO
    319320                ENDDO
    320                 DO  k = nzb, nz_do3d
    321                    DO  j = nys+1, nyn
    322                       DO  i = nxl, nxr+1
    323                          local_pf(i,j,k) = 0.25 * local_pf(i,j-1,k) + &
    324                                            0.50 * local_pf(i,j,k)   + &
    325                                            0.25 * local_pf(i,j+1,k)
    326                       ENDDO
    327                    ENDDO
    328                 ENDDO
     321! Averaging for Langmuir circulation
     322!                DO  k = nzb, nz_do3d
     323!                   DO  j = nys+1, nyn
     324!                      DO  i = nxl, nxr+1
     325!                         local_pf(i,j,k) = 0.25 * local_pf(i,j-1,k) + &
     326!                                           0.50 * local_pf(i,j,k)   + &
     327!                                           0.25 * local_pf(i,j+1,k)
     328!                      ENDDO
     329!                   ENDDO
     330!                ENDDO
    329331
    330332             CASE ( 'p', 'p_xy', 'p_xz', 'p_yz' )
     
    475477             CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, &
    476478                             cyclic_dvrp, cyclic_dvrp, cyclic_dvrp )
    477              CALL DVRP_SLICER( m-1, section_mode, slicer_position )
     479!             CALL DVRP_SLICER( m-1, section_mode, slicer_position )
     480             CALL DVRP_SLICER( m-1, 2, 1.0 )
     481             WRITE (9,*) 'nx_dvrp=', nx_dvrp
     482             WRITE (9,*) 'ny_dvrp=', ny_dvrp
     483             WRITE (9,*) 'nz_dvrp=', nz_dvrp
     484             WRITE (9,*) 'section_mode=', section_mode
     485             WRITE (9,*) 'slicer_position=', slicer_position
     486             CALL local_flush( 9 )
     487
    478488             CALL DVRP_VISUALIZE( m-1, 2, dvrp_filecount )
    479489
  • palm/trunk/SOURCE/diffusion_u.f90

    r103 r106  
    44! Actual revisions:
    55! -----------------
    6 ! Momentumflux at top (uswst) included as boundary condition
    7 !
     6! Momentumflux at top (uswst) included as boundary condition,
     7! i loop is starting from nxlu (needed for non-cyclic boundary conditions)
     8!
    89! Former revisions:
    910! -----------------
     
    7980       ENDIF
    8081
    81        DO  i = nxl, nxr
     82       DO  i = nxlu, nxr
    8283          DO  j = nys,nyn
    8384!
  • palm/trunk/SOURCE/diffusion_v.f90

    r103 r106  
    44! Actual revisions:
    55! -----------------
    6 ! Momentumflux at top (vswst) included as boundary condition
     6! Momentumflux at top (vswst) included as boundary condition,
     7! j loop is starting from nysv (needed for non-cyclic boundary conditions)
    78!
    89! Former revisions:
     
    7879
    7980       DO  i = nxl, nxr
    80           DO  j = nys, nyn
     81          DO  j = nysv, nyn
    8182!
    8283!--          Compute horizontal diffusion
  • palm/trunk/SOURCE/flow_statistics.f90

    r102 r106  
    44! Actual revisions:
    55! -----------------
    6 ! Prescribed momentum fluxes at the top surface are used
     6! Prescribed momentum fluxes at the top surface are used,
     7! profiles for w*p* and w"e are calculated
    78!
    89! Former revisions:
     
    620621!
    621622!--    Divergence of vertical flux of resolved scale energy and pressure
    622 !--    fluctuations. First calculate the products, then the divergence.
     623!--    fluctuations as well as flux of pressure fluctuation itself (68).
     624!--    First calculate the products, then the divergence.
    623625!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
    624        IF ( hom(nzb+1,2,55,0) /= 0.0 )  THEN
     626       IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
    625627
    626628          sums_ll = 0.0  ! local array
     
    654656             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
    655657             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
     658             sums_l(k,68,tn) = sums_ll(k,2)
    656659          ENDDO
    657660          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
    658661          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
    659 
    660        ENDIF
    661 
    662 !
    663 !--    Divergence of vertical flux of SGS TKE
    664        IF ( hom(nzb+1,2,57,0) /= 0.0 )  THEN
     662          sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
     663
     664       ENDIF
     665
     666!
     667!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
     668       IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
    665669
    666670          !$OMP DO
     
    669673                DO  k = nzb_s_outer(j,i)+1, nzt
    670674
    671                    sums_l(k,57,tn) = sums_l(k,57,tn) + (                       &
     675                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
    672676                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
    673677                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
    674                                                  ) * ddzw(k)
     678                                                             ) * ddzw(k)
     679
     680                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
     681                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
     682                                                              )
    675683
    676684                ENDDO
     
    678686          ENDDO
    679687          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
     688          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
    680689
    681690       ENDIF
     
    832841       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
    833842       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
    834        hom(:,1,57,sr) = sums(:,57)     ! w"e/dz
     843       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
    835844       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
    836845       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
     
    843852       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
    844853       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
     854       hom(:,1,68,sr) = sums(:,68)     ! w*p*
     855       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
    845856
    846857       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
  • palm/trunk/SOURCE/init_3d_model.f90

    r102 r106  
    77! Actual revisions:
    88! -----------------
    9 ! Flux initialization in case of coupled runs, +momentum fluxes at top boundary
     9! Flux initialization in case of coupled runs, +momentum fluxes at top boundary,
     10! +arrays for phase speed c_u, c_v, c_w, indices for u|v|w_m_l|r changed
    1011!
    1112! Former revisions:
     
    233234
    234235!
    235 !-- Arrays to store velocity data from t-dt needed for radiation boundary
    236 !-- conditions
     236!-- Arrays to store velocity data from t-dt and the phase speeds which
     237!-- are needed for radiation boundary conditions
    237238    IF ( outflow_l )  THEN
    238        ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1), &
    239                  v_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1), &
    240                  w_m_l(nzb:nzt+1,nys-1:nyn+1,-1:1) )
     239       ALLOCATE( u_m_l(nzb:nzt+1,nys-1:nyn+1,1:2), &
     240                 v_m_l(nzb:nzt+1,nys-1:nyn+1,0:1), &
     241                 w_m_l(nzb:nzt+1,nys-1:nyn+1,0:1) )
    241242    ENDIF
    242243    IF ( outflow_r )  THEN
    243        ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1), &
    244                  v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1), &
    245                  w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx+1) )
     244       ALLOCATE( u_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
     245                 v_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx), &
     246                 w_m_r(nzb:nzt+1,nys-1:nyn+1,nx-1:nx) )
     247    ENDIF
     248    IF ( outflow_l  .OR.  outflow_r )  THEN
     249       ALLOCATE( c_u(nzb:nzt+1,nys-1:nyn+1), c_v(nzb:nzt+1,nys-1:nyn+1), &
     250                 c_w(nzb:nzt+1,nys-1:nyn+1) )
    246251    ENDIF
    247252    IF ( outflow_s )  THEN
    248        ALLOCATE( u_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1), &
    249                  v_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1), &
    250                  w_m_s(nzb:nzt+1,-1:1,nxl-1:nxr+1) )
     253       ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1), &
     254                 v_m_s(nzb:nzt+1,1:2,nxl-1:nxr+1), &
     255                 w_m_s(nzb:nzt+1,0:1,nxl-1:nxr+1) )
    251256    ENDIF
    252257    IF ( outflow_n )  THEN
    253        ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1), &
    254                  v_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1), &
    255                  w_m_n(nzb:nzt+1,ny-1:ny+1,nxl-1:nxr+1) )
     258       ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
     259                 v_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1), &
     260                 w_m_n(nzb:nzt+1,ny-1:ny,nxl-1:nxr+1) )
     261    ENDIF
     262    IF ( outflow_s  .OR.  outflow_n )  THEN
     263       ALLOCATE( c_u(nzb:nzt+1,nxl-1:nxr+1), c_v(nzb:nzt+1,nxl-1:nxr+1), &
     264                 c_w(nzb:nzt+1,nxl-1:nxr+1) )
    256265    ENDIF
    257266
     
    808817!--    Initialize old timelevels needed for radiation boundary conditions
    809818       IF ( outflow_l )  THEN
    810           u_m_l(:,:,:) = u(:,:,-1:1)
    811           v_m_l(:,:,:) = v(:,:,-1:1)
    812           w_m_l(:,:,:) = w(:,:,-1:1)
     819          u_m_l(:,:,:) = u(:,:,1:2)
     820          v_m_l(:,:,:) = v(:,:,0:1)
     821          w_m_l(:,:,:) = w(:,:,0:1)
    813822       ENDIF
    814823       IF ( outflow_r )  THEN
    815           u_m_r(:,:,:) = u(:,:,nx-1:nx+1)
    816           v_m_r(:,:,:) = v(:,:,nx-1:nx+1)
    817           w_m_r(:,:,:) = w(:,:,nx-1:nx+1)
     824          u_m_r(:,:,:) = u(:,:,nx-1:nx)
     825          v_m_r(:,:,:) = v(:,:,nx-1:nx)
     826          w_m_r(:,:,:) = w(:,:,nx-1:nx)
    818827       ENDIF
    819828       IF ( outflow_s )  THEN
    820           u_m_s(:,:,:) = u(:,-1:1,:)
    821           v_m_s(:,:,:) = v(:,-1:1,:)
    822           w_m_s(:,:,:) = w(:,-1:1,:)
     829          u_m_s(:,:,:) = u(:,0:1,:)
     830          v_m_s(:,:,:) = v(:,1:2,:)
     831          w_m_s(:,:,:) = w(:,0:1,:)
    823832       ENDIF
    824833       IF ( outflow_n )  THEN
    825           u_m_n(:,:,:) = u(:,ny-1:ny+1,:)
    826           v_m_n(:,:,:) = v(:,ny-1:ny+1,:)
    827           w_m_n(:,:,:) = w(:,ny-1:ny+1,:)
     834          u_m_n(:,:,:) = u(:,ny-1:ny,:)
     835          v_m_n(:,:,:) = v(:,ny-1:ny,:)
     836          w_m_n(:,:,:) = w(:,ny-1:ny,:)
    828837       ENDIF
    829838
  • palm/trunk/SOURCE/init_particles.f90

    r83 r106  
    44! Actual revisions:
    55! -----------------
    6 !
     6! variable iran replaced by iran_part
    77!
    88! Former revisions:
     
    354354       IF ( random_start_position )  THEN
    355355
    356           iran = iran + myid    ! Random positions should be different on
    357                                 ! different PEs
    358 
    359356          DO  n = 1, number_of_initial_particles
    360357             IF ( psl(particles(n)%group) /= psr(particles(n)%group) )  THEN
    361358                particles(n)%x = particles(n)%x + &
    362                                  ( random_function( iran ) - 0.5 ) * &
     359                                 ( random_function( iran_part ) - 0.5 ) * &
    363360                                 pdx(particles(n)%group)
    364361                IF ( particles(n)%x  <=  ( nxl - 0.5 ) * dx )  THEN
     
    370367             IF ( pss(particles(n)%group) /= psn(particles(n)%group) )  THEN
    371368                particles(n)%y = particles(n)%y + &
    372                                  ( random_function( iran ) - 0.5 ) * &
     369                                 ( random_function( iran_part ) - 0.5 ) * &
    373370                                 pdy(particles(n)%group)
    374371                IF ( particles(n)%y  <=  ( nys - 0.5 ) * dy )  THEN
     
    380377             IF ( psb(particles(n)%group) /= pst(particles(n)%group) )  THEN
    381378                particles(n)%z = particles(n)%z + &
    382                                  ( random_function( iran ) - 0.5 ) * &
     379                                 ( random_function( iran_part ) - 0.5 ) * &
    383380                                 pdz(particles(n)%group)
    384381             ENDIF
  • palm/trunk/SOURCE/init_pegrid.f90

    r104 r106  
    55! -----------------
    66! Intercommunicator (comm_inter) and derived data type (type_xy) for
    7 ! coupled model runs created
     7! coupled model runs created,
     8! indices nxlu and nysv are calculated (needed for non-cyclic boundary
     9! conditions)
    810!
    911! Former revisions:
     
    818820!
    819821!-- Setting of flags for inflow/outflow conditions in case of non-cyclic
    820 !-- horizontal boundary conditions. Set variables for extending array u (v)
    821 !-- by one gridpoint on the left/rightmost (northest/southest) processor
     822!-- horizontal boundary conditions.
    822823    IF ( pleft == MPI_PROC_NULL )  THEN
    823824       IF ( bc_lr == 'dirichlet/radiation' )  THEN
     
    869870    ENDIF
    870871#endif
     872!
     873!-- At the outflow, u or v, respectively, have to be calculated for one grid
     874!-- point less.
     875    IF ( outflow_l )  THEN
     876       nxlu = nxl + 1
     877    ELSE
     878       nxlu = nxl
     879    ENDIF
     880    IF ( outflow_s )  THEN
     881       nysv = nys + 1
     882    ELSE
     883       nysv = nys
     884    ENDIF
    871885
    872886    IF ( psolver == 'poisfft_hybrid' )  THEN
  • palm/trunk/SOURCE/modules.f90

    r103 r106  
    66! -----------------
    77! +comm_inter, constant_top_momentumflux, coupling_char, coupling_mode,
    8 ! dt_coupling, ngp_xy, port_name, time_coupling, top_momentumflux_u|v,
    9 ! type_xy, uswst*, vswst*
     8! c_u, c_v, c_w, dt_coupling, ngp_xy, nxlu, nysv, port_name, time_coupling,
     9! top_momentumflux_u|v, type_xy, uswst*, vswst*
    1010!
    1111! Former revisions:
     
    105105
    106106    REAL, DIMENSION(:,:), ALLOCATABLE ::                                       &
    107           dzu_mg, dzw_mg, f1_mg, f2_mg, f3_mg, pt_slope_ref, qs, ts, us, z0
     107          c_u, c_v, c_w, dzu_mg, dzw_mg, f1_mg, f2_mg, f3_mg, pt_slope_ref,    &
     108          qs, ts, us, z0
    108109
    109110    REAL, DIMENSION(:,:), ALLOCATABLE, TARGET ::                               &
     
    565566!------------------------------------------------------------------------------!
    566567
    567     INTEGER ::  ngp_sums, nnx, nx = 0, nxa, nxl, nxr, nxra, nny, ny = 0, nya, &
    568                 nyn, nyna, nys, nnz, nz = 0, nza, nzb, nzb_diff, nzt, nzta,    &
    569                 nzt_diff
     568    INTEGER ::  ngp_sums, nnx, nx = 0, nxa, nxl, nxlu, nxr, nxra, nny, ny = 0, &
     569                nya, nyn, nyna, nys, nysv, nnz, nz = 0, nza, nzb, nzb_diff,    &
     570                nzt, nzta, nzt_diff
    570571
    571572    INTEGER, DIMENSION(:), ALLOCATABLE ::                                      &
  • palm/trunk/SOURCE/pres.f90

    r90 r106  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Volume flow conservation added for the remaining three outflow boundaries
    77!
    88! Former revisions:
     
    7474!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
    7575!-- boundary conditions
    76     IF ( conserve_volume_flow  .AND.  bc_ns == 'radiation/dirichlet')  THEN
     76!-- WARNING: so far, this conservation does not work at the left/south
     77!--          boundary if the topography at the inflow differs from that at the
     78!--          outflow! For this case, volume_flow_area needs adjustment!
     79!
     80!-- Left/right
     81    IF ( conserve_volume_flow  .AND.  ( outflow_l  .OR. outflow_r ) )  THEN
     82
     83       volume_flow(1)   = 0.0
     84       volume_flow_l(1) = 0.0
     85
     86       IF ( outflow_l )  THEN
     87          i = 0
     88       ELSEIF ( outflow_r )  THEN
     89          i = nx+1
     90       ENDIF
     91
     92       DO  j = nys, nyn
     93!
     94!--       Sum up the volume flow through the south/north boundary
     95          DO  k = nzb_2d(j,i) + 1, nzt
     96             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzu(k)
     97          ENDDO
     98       ENDDO
     99
     100#if defined( __parallel )   
     101       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
     102                           MPI_SUM, comm1dy, ierr )   
     103#else
     104       volume_flow = volume_flow_l 
     105#endif
     106       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) )    &
     107                               / volume_flow_area(1)
     108
     109       DO  j = nys, nyn
     110          DO  k = nzb_v_inner(j,i) + 1, nzt
     111             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
     112          ENDDO
     113       ENDDO
     114
     115       CALL exchange_horiz( u )
     116
     117    ENDIF
     118
     119!
     120!-- South/north
     121    IF ( conserve_volume_flow  .AND.  ( outflow_n  .OR. outflow_s ) )  THEN
    77122
    78123       volume_flow(2)   = 0.0
    79124       volume_flow_l(2) = 0.0
    80125
    81        IF ( nyn == ny )  THEN
     126       IF ( outflow_s )  THEN
     127          j = 0
     128       ELSEIF ( outflow_n )  THEN
    82129          j = ny+1
    83           DO  i = nxl, nxr
    84 !
    85 !--          Sum up the volume flow through the north boundary
    86              DO  k = nzb_2d(j,i) + 1, nzt
    87                 volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzu(k)
    88              ENDDO
    89           ENDDO
    90        ENDIF
     130       ENDIF
     131
     132       DO  i = nxl, nxr
     133!
     134!--       Sum up the volume flow through the south/north boundary
     135          DO  k = nzb_2d(j,i) + 1, nzt
     136             volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzu(k)
     137          ENDDO
     138       ENDDO
     139
    91140#if defined( __parallel )   
    92141       CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
     
    96145#endif
    97146       volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) )    &
    98                                / volume_flow_area(2)                         
    99 
    100        IF ( outflow_n )  THEN
    101           j = nyn+1
    102           DO  i = nxl, nxr
    103              DO  k = nzb_v_inner(j,i) + 1, nzt
    104                 v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
    105              ENDDO
    106           ENDDO
    107        ENDIF
     147                               / volume_flow_area(2)
     148
     149       DO  i = nxl, nxr
     150          DO  k = nzb_v_inner(j,i) + 1, nzt
     151             v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
     152          ENDDO
     153       ENDDO
    108154
    109155       CALL exchange_horiz( v )
  • palm/trunk/SOURCE/production_e.f90

    r98 r106  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Bugfix: wrong sign removed from the buoyancy production term in the case
     7! use_reference = .T.,
     8! u_0 and v_0 are calculated for nxr+1, nyn+1 also (otherwise these values are
     9! not available in case of non-cyclic boundary conditions)
    710!
    811! Former revisions:
     
    386389                   DO  j = nys, nyn
    387390                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
    388                          tend(k,j,i) = tend(k,j,i) +                  &
     391                         tend(k,j,i) = tend(k,j,i) -                  &
    389392                                       kh(k,j,i) * g / pt_reference * &
    390393                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
     
    976979!--       (otherwise the timestep may be significantly reduced by large
    977980!--       surface winds).
     981!--       Upper bounds are nxr+1 and nyn+1 because otherwise these values are
     982!--       not available in case of non-cyclic boundary conditions.
    978983!--       WARNING: the exact analytical solution would require the determination
    979984!--                of the eddy diffusivity by km = u* * kappa * zp / phi_m.
    980985          !$OMP PARALLEL DO PRIVATE( ku, kv )
    981           DO  i = nxl, nxr
    982              DO  j = nys, nyn
     986          DO  i = nxl, nxr+1
     987             DO  j = nys, nyn+1
    983988
    984989                ku = nzb_u_inner(j,i)+1
  • palm/trunk/SOURCE/prognostic_equations.f90

    r102 r106  
    44! Actual revisions:
    55! -----------------
    6 ! +uswst, vswst as arguments in calls of diffusion_u|v
     6! +uswst, vswst as arguments in calls of diffusion_u|v,
     7! loops for u and v are starting from index nxlu, nysv, respectively (needed
     8! for non-cyclic boundary conditions)
    79!
    810! Former revisions:
     
    133135!
    134136!-- u-tendency terms with no communication
    135     DO  i = nxl, nxr
     137    DO  i = nxlu, nxr
    136138       DO  j = nys, nyn
    137139!
     
    202204!-- v-tendency terms with no communication
    203205    DO  i = nxl, nxr
    204        DO  j = nys, nyn
     206       DO  j = nysv, nyn
    205207!
    206208!--       Tendency terms
     
    806808!
    807809!--       Tendency terms for u-velocity component
    808           IF ( j < nyn+1 )  THEN
     810          IF ( .NOT. outflow_l  .OR.  i > nxl )  THEN
    809811
    810812             tend(:,j,i) = 0.0
     
    857859!
    858860!--       Tendency terms for v-velocity component
    859           IF ( i < nxr+1 )  THEN
     861          IF ( .NOT. outflow_s  .OR.  j > nys )  THEN
    860862
    861863             tend(:,j,i) = 0.0
     
    906908!
    907909!--       Tendency terms for w-velocity component
    908           IF ( i < nxr+1  .AND.  j < nyn+1 )  THEN
    909 
     910          tend(:,j,i) = 0.0
     911          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
     912             CALL advec_w_pw( i, j )
     913          ELSE
     914             CALL advec_w_up( i, j )
     915          ENDIF
     916          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
     917          THEN
     918             CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x,          &
     919                               km_damp_y, tend, u_m, v_m, w_m )
     920          ELSE
     921             CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
     922                               tend, u, v, w )
     923          ENDIF
     924          CALL coriolis( i, j, 3 )
     925          IF ( ocean )  THEN
     926             CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
     927          ELSE
     928             IF ( .NOT. humidity )  THEN
     929                CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
     930             ELSE
     931                CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
     932             ENDIF
     933          ENDIF
     934          CALL user_actions( i, j, 'w-tendency' )
     935
     936!
     937!--       Prognostic equation for w-velocity component
     938          DO  k = nzb_w_inner(j,i)+1, nzt-1
     939             w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
     940                          dt_3d * (                                         &
     941                                tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
     942                           - tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
     943                                  ) -                                       &
     944                          tsc(5) * rdf(k) * w(k,j,i)
     945          ENDDO
     946
     947!
     948!--       Calculate tendencies for the next Runge-Kutta step
     949          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     950             IF ( intermediate_timestep_count == 1 )  THEN
     951                DO  k = nzb_w_inner(j,i)+1, nzt-1
     952                   tw_m(k,j,i) = tend(k,j,i)
     953                ENDDO
     954             ELSEIF ( intermediate_timestep_count < &
     955                      intermediate_timestep_count_max )  THEN
     956                DO  k = nzb_w_inner(j,i)+1, nzt-1
     957                   tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
     958                ENDDO
     959             ENDIF
     960          ENDIF
     961
     962!
     963!--       Tendency terms for potential temperature
     964          tend(:,j,i) = 0.0
     965          IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
     966             CALL advec_s_pw( i, j, pt )
     967          ELSE
     968             CALL advec_s_up( i, j, pt )
     969          ENDIF
     970          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
     971          THEN
     972             CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
     973                               tswst_m, tend )
     974          ELSE
     975             CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
     976          ENDIF
     977
     978!
     979!--       If required compute heating/cooling due to long wave radiation
     980!--       processes
     981          IF ( radiation )  THEN
     982             CALL calc_radiation( i, j )
     983          ENDIF
     984
     985!
     986!--       If required compute impact of latent heat due to precipitation
     987          IF ( precipitation )  THEN
     988             CALL impact_of_latent_heat( i, j )
     989          ENDIF
     990          CALL user_actions( i, j, 'pt-tendency' )
     991
     992!
     993!--       Prognostic equation for potential temperature
     994          DO  k = nzb_s_inner(j,i)+1, nzt
     995             pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
     996                           dt_3d * (                                        &
     997                               tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
     998                                   ) -                                      &
     999                           tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
     1000          ENDDO
     1001
     1002!
     1003!--       Calculate tendencies for the next Runge-Kutta step
     1004          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1005             IF ( intermediate_timestep_count == 1 )  THEN
     1006                DO  k = nzb_s_inner(j,i)+1, nzt
     1007                   tpt_m(k,j,i) = tend(k,j,i)
     1008                ENDDO
     1009             ELSEIF ( intermediate_timestep_count < &
     1010                      intermediate_timestep_count_max )  THEN
     1011                DO  k = nzb_s_inner(j,i)+1, nzt
     1012                   tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
     1013                                   5.3125 * tpt_m(k,j,i)
     1014                ENDDO
     1015             ENDIF
     1016          ENDIF
     1017
     1018!
     1019!--       If required, compute prognostic equation for salinity
     1020          IF ( ocean )  THEN
     1021
     1022!
     1023!--          Tendency-terms for salinity
    9101024             tend(:,j,i) = 0.0
    911              IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
    912                 CALL advec_w_pw( i, j )
    913              ELSE
    914                 CALL advec_w_up( i, j )
    915              ENDIF
    916              IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
     1025             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
    9171026             THEN
    918                 CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x,          &
    919                                   km_damp_y, tend, u_m, v_m, w_m )
    920              ELSE
    921                 CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
    922                                   tend, u, v, w )
    923              ENDIF
    924              CALL coriolis( i, j, 3 )
    925              IF ( ocean )  THEN
    926                 CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
    927              ELSE
    928                 IF ( .NOT. humidity )  THEN
    929                    CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
    930                 ELSE
    931                    CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
    932                 ENDIF
    933              ENDIF
    934              CALL user_actions( i, j, 'w-tendency' )
    935 
    936 !
    937 !--          Prognostic equation for w-velocity component
    938              DO  k = nzb_w_inner(j,i)+1, nzt-1
    939                 w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
    940                              dt_3d * (                                         &
    941                                    tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
    942                               - tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
    943                                      ) -                                       &
    944                              tsc(5) * rdf(k) * w(k,j,i)
    945              ENDDO
    946 
    947 !
    948 !--          Calculate tendencies for the next Runge-Kutta step
    949              IF ( timestep_scheme(1:5) == 'runge' )  THEN
    950                 IF ( intermediate_timestep_count == 1 )  THEN
    951                    DO  k = nzb_w_inner(j,i)+1, nzt-1
    952                       tw_m(k,j,i) = tend(k,j,i)
    953                    ENDDO
    954                 ELSEIF ( intermediate_timestep_count < &
    955                          intermediate_timestep_count_max )  THEN
    956                    DO  k = nzb_w_inner(j,i)+1, nzt-1
    957                       tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
    958                    ENDDO
    959                 ENDIF
    960              ENDIF
    961 
    962 !
    963 !--          Tendency terms for potential temperature
    964              tend(:,j,i) = 0.0
    965              IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  THEN
    966                 CALL advec_s_pw( i, j, pt )
    967              ELSE
    968                 CALL advec_s_up( i, j, pt )
    969              ENDIF
    970              IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' ) &
    971              THEN
    972                 CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
    973                                   tswst_m, tend )
    974              ELSE
    975                 CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
    976              ENDIF
    977 
    978 !
    979 !--          If required compute heating/cooling due to long wave radiation
    980 !--          processes
    981              IF ( radiation )  THEN
    982                 CALL calc_radiation( i, j )
    983              ENDIF
    984 
    985 !
    986 !--          If required compute impact of latent heat due to precipitation
    987              IF ( precipitation )  THEN
    988                 CALL impact_of_latent_heat( i, j )
    989              ENDIF
    990              CALL user_actions( i, j, 'pt-tendency' )
    991 
    992 !
    993 !--          Prognostic equation for potential temperature
     1027                CALL advec_s_pw( i, j, sa )
     1028             ELSE
     1029                CALL advec_s_up( i, j, sa )
     1030             ENDIF
     1031             CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
     1032                               tend )
     1033
     1034             CALL user_actions( i, j, 'sa-tendency' )
     1035
     1036!
     1037!--          Prognostic equation for salinity
    9941038             DO  k = nzb_s_inner(j,i)+1, nzt
    995                 pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
    996                               dt_3d * (                                        &
    997                                   tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
    998                                       ) -                                      &
    999                               tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
     1039                sa_p(k,j,i) = tsc(1) * sa(k,j,i) +                          &
     1040                              dt_3d * (                                     &
     1041                               tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
     1042                                      ) -                                   &
     1043                             tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
     1044                IF ( sa_p(k,j,i) < 0.0 )  sa_p(k,j,i) = 0.1 * sa(k,j,i)
    10001045             ENDDO
    10011046
     
    10051050                IF ( intermediate_timestep_count == 1 )  THEN
    10061051                   DO  k = nzb_s_inner(j,i)+1, nzt
    1007                       tpt_m(k,j,i) = tend(k,j,i)
     1052                      tsa_m(k,j,i) = tend(k,j,i)
    10081053                   ENDDO
    10091054                ELSEIF ( intermediate_timestep_count < &
    10101055                         intermediate_timestep_count_max )  THEN
    10111056                   DO  k = nzb_s_inner(j,i)+1, nzt
    1012                       tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
    1013                                       5.3125 * tpt_m(k,j,i)
    1014                    ENDDO
    1015                 ENDIF
    1016              ENDIF
    1017 
    1018 !
    1019 !--          If required, compute prognostic equation for salinity
    1020              IF ( ocean )  THEN
    1021 
    1022 !
    1023 !--             Tendency-terms for salinity
    1024                 tend(:,j,i) = 0.0
    1025                 IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
    1026                 THEN
    1027                    CALL advec_s_pw( i, j, sa )
     1057                      tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
     1058                                      5.3125 * tsa_m(k,j,i)
     1059                   ENDDO
     1060                ENDIF
     1061             ENDIF
     1062
     1063!
     1064!--          Calculate density by the equation of state for seawater
     1065             CALL eqn_state_seawater( i, j )
     1066
     1067          ENDIF
     1068
     1069!
     1070!--       If required, compute prognostic equation for total water content /
     1071!--       scalar
     1072          IF ( humidity  .OR.  passive_scalar )  THEN
     1073
     1074!
     1075!--          Tendency-terms for total water content / scalar
     1076             tend(:,j,i) = 0.0
     1077             IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
     1078             THEN
     1079                CALL advec_s_pw( i, j, q )
     1080             ELSE
     1081                CALL advec_s_up( i, j, q )
     1082             ENDIF
     1083             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
     1084             THEN
     1085                CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
     1086                                  qswst_m, tend )
     1087             ELSE
     1088                CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
     1089                                  tend )
     1090             ENDIF
     1091       
     1092!
     1093!--          If required compute decrease of total water content due to
     1094!--          precipitation
     1095             IF ( precipitation )  THEN
     1096                CALL calc_precipitation( i, j )
     1097             ENDIF
     1098             CALL user_actions( i, j, 'q-tendency' )
     1099
     1100!
     1101!--          Prognostic equation for total water content / scalar
     1102             DO  k = nzb_s_inner(j,i)+1, nzt
     1103                q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
     1104                             dt_3d * (                                      &
     1105                                tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
     1106                                     ) -                                    &
     1107                             tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
     1108                IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
     1109             ENDDO
     1110
     1111!
     1112!--          Calculate tendencies for the next Runge-Kutta step
     1113             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1114                IF ( intermediate_timestep_count == 1 )  THEN
     1115                   DO  k = nzb_s_inner(j,i)+1, nzt
     1116                      tq_m(k,j,i) = tend(k,j,i)
     1117                   ENDDO
     1118                ELSEIF ( intermediate_timestep_count < &
     1119                         intermediate_timestep_count_max )  THEN
     1120                   DO  k = nzb_s_inner(j,i)+1, nzt
     1121                      tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
     1122                                     5.3125 * tq_m(k,j,i)
     1123                   ENDDO
     1124                ENDIF
     1125             ENDIF
     1126
     1127          ENDIF
     1128
     1129!
     1130!--       If required, compute prognostic equation for turbulent kinetic
     1131!--       energy (TKE)
     1132          IF ( .NOT. constant_diffusion )  THEN
     1133
     1134!
     1135!--          Tendency-terms for TKE
     1136             tend(:,j,i) = 0.0
     1137             IF ( ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  &
     1138                  .AND.  .NOT. use_upstream_for_tke )  THEN
     1139                CALL advec_s_pw( i, j, e )
     1140             ELSE
     1141                CALL advec_s_up( i, j, e )
     1142             ENDIF
     1143             IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
     1144             THEN
     1145                IF ( .NOT. humidity )  THEN
     1146                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
     1147                                     km_m, l_grid, pt_m, pt_reference,   &
     1148                                     rif_m, tend, zu, zw )
    10281149                ELSE
    1029                    CALL advec_s_up( i, j, sa )
    1030                 ENDIF
    1031                 CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
    1032                                   tend )
    1033        
    1034                 CALL user_actions( i, j, 'sa-tendency' )
    1035 
    1036 !
    1037 !--             Prognostic equation for salinity
    1038                 DO  k = nzb_s_inner(j,i)+1, nzt
    1039                    sa_p(k,j,i) = tsc(1) * sa(k,j,i) +                          &
    1040                                  dt_3d * (                                     &
    1041                                   tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
    1042                                          ) -                                   &
    1043                                 tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
    1044                    IF ( sa_p(k,j,i) < 0.0 )  sa_p(k,j,i) = 0.1 * sa(k,j,i)
    1045                 ENDDO
    1046 
    1047 !
    1048 !--             Calculate tendencies for the next Runge-Kutta step
    1049                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
    1050                    IF ( intermediate_timestep_count == 1 )  THEN
    1051                       DO  k = nzb_s_inner(j,i)+1, nzt
    1052                          tsa_m(k,j,i) = tend(k,j,i)
    1053                       ENDDO
    1054                    ELSEIF ( intermediate_timestep_count < &
    1055                             intermediate_timestep_count_max )  THEN
    1056                       DO  k = nzb_s_inner(j,i)+1, nzt
    1057                          tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
    1058                                          5.3125 * tsa_m(k,j,i)
    1059                       ENDDO
    1060                    ENDIF
    1061                 ENDIF
    1062 
    1063 !
    1064 !--             Calculate density by the equation of state for seawater
    1065                 CALL eqn_state_seawater( i, j )
    1066 
    1067              ENDIF
    1068 
    1069 !
    1070 !--          If required, compute prognostic equation for total water content /
    1071 !--          scalar
    1072              IF ( humidity  .OR.  passive_scalar )  THEN
    1073 
    1074 !
    1075 !--             Tendency-terms for total water content / scalar
    1076                 tend(:,j,i) = 0.0
    1077                 IF ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' ) &
    1078                 THEN
    1079                    CALL advec_s_pw( i, j, q )
    1080                 ELSE
    1081                    CALL advec_s_up( i, j, q )
    1082                 ENDIF
    1083                 IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
    1084                 THEN
    1085                    CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
    1086                                      qswst_m, tend )
    1087                 ELSE
    1088                    CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
    1089                                      tend )
    1090                 ENDIF
    1091        
    1092 !
    1093 !--             If required compute decrease of total water content due to
    1094 !--             precipitation
    1095                 IF ( precipitation )  THEN
    1096                    CALL calc_precipitation( i, j )
    1097                 ENDIF
    1098                 CALL user_actions( i, j, 'q-tendency' )
    1099 
    1100 !
    1101 !--             Prognostic equation for total water content / scalar
    1102                 DO  k = nzb_s_inner(j,i)+1, nzt
    1103                    q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
    1104                                 dt_3d * (                                      &
    1105                                    tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
    1106                                         ) -                                    &
    1107                                 tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
    1108                    IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
    1109                 ENDDO
    1110 
    1111 !
    1112 !--             Calculate tendencies for the next Runge-Kutta step
    1113                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
    1114                    IF ( intermediate_timestep_count == 1 )  THEN
    1115                       DO  k = nzb_s_inner(j,i)+1, nzt
    1116                          tq_m(k,j,i) = tend(k,j,i)
    1117                       ENDDO
    1118                    ELSEIF ( intermediate_timestep_count < &
    1119                             intermediate_timestep_count_max )  THEN
    1120                       DO  k = nzb_s_inner(j,i)+1, nzt
    1121                          tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
    1122                                         5.3125 * tq_m(k,j,i)
    1123                       ENDDO
    1124                    ENDIF
    1125                 ENDIF
    1126 
    1127              ENDIF
    1128 
    1129 !
    1130 !--          If required, compute prognostic equation for turbulent kinetic
    1131 !--          energy (TKE)
    1132              IF ( .NOT. constant_diffusion )  THEN
    1133 
    1134 !
    1135 !--             Tendency-terms for TKE
    1136                 tend(:,j,i) = 0.0
    1137                 IF ( ( tsc(2) == 2.0  .OR.  timestep_scheme(1:5) == 'runge' )  &
    1138                      .AND.  .NOT. use_upstream_for_tke )  THEN
    1139                    CALL advec_s_pw( i, j, e )
    1140                 ELSE
    1141                    CALL advec_s_up( i, j, e )
    1142                 ENDIF
    1143                 IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )&
    1144                 THEN
    1145                    IF ( .NOT. humidity )  THEN
    1146                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
    1147                                         km_m, l_grid, pt_m, pt_reference,   &
    1148                                         rif_m, tend, zu, zw )
     1150                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
     1151                                     km_m, l_grid, vpt_m, pt_reference,  &
     1152                                     rif_m, tend, zu, zw )
     1153                ENDIF
     1154             ELSE
     1155                IF ( .NOT. humidity )  THEN
     1156                   IF ( ocean )  THEN
     1157                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
     1158                                        km, l_grid, rho, prho_reference,  &
     1159                                        rif, tend, zu, zw )
    11491160                   ELSE
    1150                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
    1151                                         km_m, l_grid, vpt_m, pt_reference, &
    1152                                         rif_m, tend, zu, zw )
     1161                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
     1162                                        km, l_grid, pt, pt_reference, rif, &
     1163                                        tend, zu, zw )
    11531164                   ENDIF
    11541165                ELSE
    1155                    IF ( .NOT. humidity )  THEN
    1156                       IF ( ocean )  THEN
    1157                          CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
    1158                                            km, l_grid, rho, prho_reference,  &
    1159                                            rif, tend, zu, zw )
    1160                       ELSE
    1161                          CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
    1162                                            km, l_grid, pt, pt_reference, rif, &
    1163                                            tend, zu, zw )
    1164                       ENDIF
    1165                    ELSE
    1166                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
    1167                                         l_grid, vpt, pt_reference, rif, tend, &
    1168                                         zu, zw )
    1169                    ENDIF
    1170                 ENDIF
    1171                 CALL production_e( i, j )
    1172                 CALL user_actions( i, j, 'e-tendency' )
    1173 
    1174 !
    1175 !--             Prognostic equation for TKE.
    1176 !--             Eliminate negative TKE values, which can occur due to numerical
    1177 !--             reasons in the course of the integration. In such cases the old
    1178 !--             TKE value is reduced by 90%.
    1179                 DO  k = nzb_s_inner(j,i)+1, nzt
    1180                    e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
    1181                                 dt_3d * (                                      &
    1182                                    tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
    1183                                         )
    1184                    IF ( e_p(k,j,i) < 0.0 )  e_p(k,j,i) = 0.1 * e(k,j,i)
    1185                 ENDDO
    1186 
    1187 !
    1188 !--             Calculate tendencies for the next Runge-Kutta step
    1189                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
    1190                    IF ( intermediate_timestep_count == 1 )  THEN
    1191                       DO  k = nzb_s_inner(j,i)+1, nzt
    1192                          te_m(k,j,i) = tend(k,j,i)
    1193                       ENDDO
    1194                    ELSEIF ( intermediate_timestep_count < &
    1195                             intermediate_timestep_count_max )  THEN
    1196                       DO  k = nzb_s_inner(j,i)+1, nzt
    1197                          te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
    1198                                         5.3125 * te_m(k,j,i)
    1199                       ENDDO
    1200                    ENDIF
    1201                 ENDIF
    1202 
    1203              ENDIF   ! TKE equation
    1204 
    1205           ENDIF   ! Gridpoints excluding the non-cyclic wall
     1166                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
     1167                                     l_grid, vpt, pt_reference, rif, tend, &
     1168                                     zu, zw )
     1169                ENDIF
     1170             ENDIF
     1171             CALL production_e( i, j )
     1172             CALL user_actions( i, j, 'e-tendency' )
     1173
     1174!
     1175!--          Prognostic equation for TKE.
     1176!--          Eliminate negative TKE values, which can occur due to numerical
     1177!--          reasons in the course of the integration. In such cases the old
     1178!--          TKE value is reduced by 90%.
     1179             DO  k = nzb_s_inner(j,i)+1, nzt
     1180                e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
     1181                             dt_3d * (                                      &
     1182                                tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
     1183                                     )
     1184                IF ( e_p(k,j,i) < 0.0 )  e_p(k,j,i) = 0.1 * e(k,j,i)
     1185             ENDDO
     1186
     1187!
     1188!--          Calculate tendencies for the next Runge-Kutta step
     1189             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1190                IF ( intermediate_timestep_count == 1 )  THEN
     1191                   DO  k = nzb_s_inner(j,i)+1, nzt
     1192                      te_m(k,j,i) = tend(k,j,i)
     1193                   ENDDO
     1194                ELSEIF ( intermediate_timestep_count < &
     1195                         intermediate_timestep_count_max )  THEN
     1196                   DO  k = nzb_s_inner(j,i)+1, nzt
     1197                      te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
     1198                                     5.3125 * te_m(k,j,i)
     1199                   ENDDO
     1200                ENDIF
     1201             ENDIF
     1202
     1203          ENDIF   ! TKE equation
    12061204
    12071205       ENDDO
     
    12681266!
    12691267!-- Prognostic equation for u-velocity component
    1270     DO  i = nxl, nxr
     1268    DO  i = nxlu, nxr
    12711269       DO  j = nys, nyn
    12721270          DO  k = nzb_u_inner(j,i)+1, nzt
     
    12851283    IF ( timestep_scheme(1:5) == 'runge' )  THEN
    12861284       IF ( intermediate_timestep_count == 1 )  THEN
    1287           DO  i = nxl, nxr
     1285          DO  i = nxlu, nxr
    12881286             DO  j = nys, nyn
    12891287                DO  k = nzb_u_inner(j,i)+1, nzt
     
    12941292       ELSEIF ( intermediate_timestep_count < &
    12951293                intermediate_timestep_count_max )  THEN
    1296           DO  i = nxl, nxr
     1294          DO  i = nxlu, nxr
    12971295             DO  j = nys, nyn
    12981296                DO  k = nzb_u_inner(j,i)+1, nzt
     
    13401338!-- Prognostic equation for v-velocity component
    13411339    DO  i = nxl, nxr
    1342        DO  j = nys, nyn
     1340       DO  j = nysv, nyn
    13431341          DO  k = nzb_v_inner(j,i)+1, nzt
    13441342             v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) +    &
     
    13571355       IF ( intermediate_timestep_count == 1 )  THEN
    13581356          DO  i = nxl, nxr
    1359              DO  j = nys, nyn
     1357             DO  j = nysv, nyn
    13601358                DO  k = nzb_v_inner(j,i)+1, nzt
    13611359                   tv_m(k,j,i) = tend(k,j,i)
     
    13661364                intermediate_timestep_count_max )  THEN
    13671365          DO  i = nxl, nxr
    1368              DO  j = nys, nyn
     1366             DO  j = nysv, nyn
    13691367                DO  k = nzb_v_inner(j,i)+1, nzt
    13701368                   tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
  • palm/trunk/SOURCE/time_integration.f90

    r102 r106  
    44! Actual revisions:
    55! -----------------
    6 ! Call of new routine surface_coupler
     6! Call of new routine surface_coupler,
     7! presure solver is called after the first Runge-Kutta substep instead of the
     8! last in case that call_psolver_at_all_substeps = .F.; for this case, the
     9! random perturbation has to be added to the velocity fields also after the
     10! first substep
    711!
    812! Former revisions:
     
    194198!
    195199!--       Impose a random perturbation on the horizontal velocity field
    196           IF ( create_disturbances  .AND.  &
     200          IF ( create_disturbances  .AND.                                      &
     201               ( call_psolver_at_all_substeps  .AND.                           &
    197202               intermediate_timestep_count == intermediate_timestep_count_max )&
     203          .OR. ( .NOT. call_psolver_at_all_substeps  .AND.                     &
     204               intermediate_timestep_count == 1 ) )                            &
    198205          THEN
    199206             time_disturb = time_disturb + dt_3d
     
    218225!--       Reduce the velocity divergence via the equation for perturbation
    219226!--       pressure.
    220           IF ( intermediate_timestep_count == intermediate_timestep_count_max &
    221                .OR. call_psolver_at_all_substeps )  THEN
     227          IF ( intermediate_timestep_count == 1  .OR. &
     228                call_psolver_at_all_substeps )  THEN
    222229             CALL pres
    223230          ENDIF
    224 
    225 !
    226 !--       In case of a non-cyclic lateral wall, set the boundary conditions for
    227 !--       the velocities at the outflow
    228 !          IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
    229 !             CALL boundary_conds( 'outflow_uvw' )
    230 !          ENDIF
    231231
    232232!
Note: See TracChangeset for help on using the changeset viewer.