Changeset 673 for palm/trunk/SOURCE


Ignore:
Timestamp:
Jan 18, 2011 4:19:48 PM (13 years ago)
Author:
suehring
Message:

Right computation of the pressure using Runge-Kutta weighting coefficients. Consideration of the pressure gradient during the time integration removed. Removed bugfix concerning velocity variances.

Location:
palm/trunk/SOURCE
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/advec_ws.f90

    r668 r673  
    55! Current revisions:
    66! -----------------
    7 !
     7! Allocation of weight_substep moved to init_3d_model.
     8! Declaration of ws_scheme_sca and ws_scheme_mom moved to check_parameters.
     9! Setting bc for the horizontal velocity variances added (moved from
     10! flow_statistics)
    811!
    912! Former revisions:
     
    8184       USE statistics
    8285
    83 !       
    84 !--    Set the LOGICALS to enhance the performance.
    85        IF ( momentum_advec == 'ws-scheme' )  ws_scheme_mom = .TRUE.
    86        IF ( scalar_advec == 'ws-scheme' )  ws_scheme_sca = .TRUE.
    8786!
    8887!--    Allocate arrays needed for dissipation control.
     
    101100       adv_mom_5 = 0.0083333333333333
    102101       adv_mom_3 = 0.041666666666666
    103 !
    104 !--    Allocate arrays needed for statistical evaluation of fluxes when
    105 !--    ws-scheme is used.
    106 !--    The following array contains the weighting factors of the RK3 average
    107 !--    of tendecies.
    108        ALLOCATE (weight_substep(1:intermediate_timestep_count_max) )
    109        IF ( intermediate_timestep_count_max == 3)  THEN   !RK3
    110           weight_substep(1) = 0.166666666666666
    111           weight_substep(2) = 0.3
    112           weight_substep(3) = 0.533333333333333
    113        ENDIF
     102
    114103!         
    115104!--    Arrays needed for statical evaluation of fluxes.
     
    197186       
    198187    END SUBROUTINE ws_init
     188   
     189
    199190   
    200191    SUBROUTINE ws_statistics
     
    821812             / (abs(u_comp(k) - gu) + 1.0E-20)         )                      &
    822813             * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
    823        ENDDO   
     814       ENDDO
     815       sums_us2_ws_l(nzb_u_inner(j,i),:) =                                   &
     816                                           sums_us2_ws_l(nzb_u_inner(j,i)+1,:)
     817
    824818!
    825819!--    Vertical advection, degradation of order near surface and top.
     
    11531147
    11541148       ENDDO
     1149       sums_vs2_ws_l(nzb_v_inner(j,i),:) =                                   & 
     1150                                           sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)
    11551151!
    11561152!--    Vertical advection, degradation of order near surface and top.
     
    22082204                 * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
    22092205            ENDDO
     2206            sums_us2_ws_l(nzb_u_inner(j,i),:) =                               &
     2207                                           sums_us2_ws_l(nzb_u_inner(j,i)+1,:)
    22102208         ENDDO
    22112209       ENDDO
     
    25702568                  * weight_substep(intermediate_timestep_count) * rmask(j,i,:)
    25712569            ENDDO
     2570            sums_vs2_ws_l(nzb_v_inner(j,i),:) =                              & 
     2571                                           sums_vs2_ws_l(nzb_v_inner(j,i)+1,:)
    25722572         ENDDO
    25732573       ENDDO
  • palm/trunk/SOURCE/check_parameters.f90

    r668 r673  
    44! Current revisions:
    55! -----------------
    6 !
     6! Declaration of ws_scheme_sca and ws_scheme_mom added (moved from
     7! advec_ws.f90).
    78!
    89! Former revisions:
     
    570571!
    571572!-- Advection schemes:
     573!       
     574!-- Set the LOGICALS to enhance the performance.
     575    IF ( momentum_advec == 'ws-scheme' )    ws_scheme_mom = .TRUE.
     576    IF ( scalar_advec   == 'ws-scheme'   )  ws_scheme_sca = .TRUE.
     577   
    572578    IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' .AND. &
    573579         momentum_advec /= 'ups-scheme' ) THEN
  • palm/trunk/SOURCE/data_output_2d.f90

    r668 r673  
    44! Current revisions:
    55! -----------------
    6 !
     6! When using Multigrid or SOR solver an additional CALL exchange_horiz is
     7! is needed for pressure output.
    78! Former revisions:
    89! -----------------
     
    239240
    240241             CASE ( 'p_xy', 'p_xz', 'p_yz' )
     242                IF ( psolver == 'multigrid' .OR. psolver == 'sor' )           &
     243                   CALL exchange_horiz( p, nbgp )
    241244                IF ( av == 0 )  THEN
    242245                   to_be_resorted => p
  • palm/trunk/SOURCE/data_output_3d.f90

    r668 r673  
    44! Current revisions:
    55! -----------------
    6 !
     6! When using Multigrid or SOR solver an additional CALL exchange_horiz is
     7! is needed for pressure output.
    78! Former revisions:
    89! -----------------
     
    152153
    153154          CASE ( 'p' )
     155             IF ( psolver == 'multigrid' .OR. psolver == 'sor' )              &
     156                CALL exchange_horiz( p, nbgp )
    154157             IF ( av == 0 )  THEN
    155158                to_be_resorted => p
  • palm/trunk/SOURCE/flow_statistics.f90

    r668 r673  
    44! Current revisions:
    55! -----------------
    6 ! When advection is computed with ws-scheme, turbulent fluxes are already
    7 ! computed in the respective advection routines and buffered in arrays
    8 ! sums_xx_ws_l(). This is due to a consistent treatment of statistics with the
    9 ! numerics and to avoid unphysical kinks near the surface.
    10 ! So some if requests has to be done to dicern between fluxes from ws-scheme
    11 ! other advection schemes.
    12 ! Furthermore the computation of z_i is only done if the heat flux exceeds a
    13 ! minimum value. This affects only simulations of a neutral boundary layer and
    14 ! is due to reasons of computations in the advection scheme.
    15 !
     6! Top bc for the horizontal velocity variances added for ocean runs.
     7! Setting the corresponding bottom bc moved to advec_ws.
    168!
    179! Former revisions:
     
    157149       IF ( ws_scheme_mom )  THEN
    158150!       
    159 !--       Boundary condition for u'u' and v'v', because below the surface no
    160 !--       computation for these quantities is done.
    161           DO  i = nxl, nxr
    162              DO  j =  nys, nyn
    163                 sums_us2_ws_l(nzb_u_inner(j,i),sr) =                          &
    164                     sums_us2_ws_l(nzb_u_inner(j,i)+1,sr)
    165                 sums_vs2_ws_l(nzb_v_inner(j,i),sr) =                          & 
    166                     sums_vs2_ws_l(nzb_v_inner(j,i)+1,sr)
    167              ENDDO
    168           ENDDO
     151!--       According to the Neumann bc for the horizontal velocity components,
     152!--       the corresponding fluxes has to satisfiy the same bc.
     153          IF ( ocean )  THEN
     154             sums_us2_ws_l(nzt+1,sr) = sums_us2_ws_l(nzt,sr)
     155             sums_vs2_ws_l(nzt+1,sr) = sums_vs2_ws_l(nzt,sr)   
     156          ENDIF
    169157!         
    170158!--       Swap the turbulent quantities evaluated in advec_ws.
     
    184172       ENDIF
    185173       IF ( ws_scheme_sca )  THEN
    186           sums_l(:,17,0) = sums_wspts_ws_l(:,sr)      ! w*pts* from advec_s_ws
     174          sums_l(:,17,0) = sums_wspts_ws_l(:,sr)      ! w*pt* from advec_s_ws
    187175          IF ( ocean ) sums_l(:,66,0) = sums_wssas_ws_l(:,sr) ! w*sa*
    188176          IF ( humidity  .OR.  passive_scalar ) sums_l(:,49,0) =              &
  • palm/trunk/SOURCE/init_3d_model.f90

    r668 r673  
    77! Current revisions:
    88! -----------------
    9 !
     9! weight_substep (moved from advec_ws) and weight_pres added.
     10! Allocate p_sub when using Multigrid or SOR solver.
     11! Call of ws_init moved behind the if requests.
    1012! Former revisions:
    1113! -----------------
     
    208210              w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
    209211              w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     212!
     213!-- Following array is required to buffer the perturbation pressure during
     214!-- Runge-Kutta 3rd order time integration.
     215    IF ( psolver == 'multigrid' .OR. psolver == 'sor' )  THEN
     216       ALLOCATE( p_sub(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     217    ENDIF
    210218
    211219    IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     
    425433
    426434    ENDIF
    427 
     435   
     436!
     437!-- Allocate arrays containing the RK coefficient for right evaluation of
     438!-- perturbation pressure and turbulent fluxes. At this point it is needed
     439!-- for right pressure correction during initialization. Further below
     440!-- the real values will be set.
     441    ALLOCATE (weight_substep(1:intermediate_timestep_count_max),           &
     442              weight_pres(1:intermediate_timestep_count_max) )
     443    weight_substep = 1.
     444    weight_pres = 1.
     445    intermediate_timestep_count = 1 ! needed for simulated_time=0
     446       
    428447!
    429448!-- Initialize model variables
     
    823842          IF ( precipitation )  precipitation_amount = 0.0
    824843       ENDIF
     844       
     845!
     846!--    Initialize quantities for special advections schemes
     847       CALL init_advec
    825848
    826849!
     
    13041327
    13051328!
    1306 !-- Initialize quantities for special advections schemes
    1307     CALL init_advec
    1308     IF ( momentum_advec == 'ws-scheme' .OR.  &
    1309          scalar_advec == 'ws-scheme' ) CALL ws_init
     1329!-- Initialize the ws-scheme.   
     1330    IF ( ws_scheme_sca .OR. ws_scheme_mom )  CALL ws_init       
     1331
     1332!
     1333!-- Setting weighting factors for right evaluation of perturbation pressure
     1334!-- and turbulent quantities during the RK substeps.               
     1335    IF ( TRIM(timestep_scheme) == 'runge-kutta-3' )  THEN    ! RK3       
     1336       weight_substep(1) = 0.166666666666666
     1337       weight_substep(2) = 0.3
     1338       weight_substep(3) = 0.533333333333333
     1339         
     1340       weight_pres(1) = 0.333333333333333
     1341       weight_pres(2) = 0.416666666666666
     1342       weight_pres(3) = 0.25         
     1343    ELSEIF ( TRIM(timestep_scheme) == 'runge-kutta-2' )  THEN  ! RK2       
     1344       weight_substep(1) = 0.5
     1345       weight_substep(2) = 0.5
     1346         
     1347       weight_pres(1) = 0.5
     1348       weight_pres(2) = 0.5         
     1349    ELSE                                                  ! Euler and Leapfrog       
     1350       weight_substep(1) = 1.0     
     1351       weight_pres(1) = 1.0                   
     1352    ENDIF
    13101353
    13111354!
  • palm/trunk/SOURCE/modules.f90

    r668 r673  
    55! Current revisions:
    66! -----------------
    7 !
     7! +weight_pres to weight the respective contribution of the Runge-Kutta
     8! substeps. +p_sub to buffer the intermediate contributions for Multigrid and
     9! SOR.
    810! Former revisions:
    911! -----------------
     
    252254    REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
    253255          canopy_heat_flux, cdc, d, diss, lad_s, lad_u, lad_v, lad_w, lai,     &
    254           l_wall, sec, sls, tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n,    &
    255           v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
     256          l_wall, p_sub, sec, sls, tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l,    &
     257          v_m_n, v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
    256258
    257259    REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
     
    13821384    REAL ::     u_max, v_max, w_max
    13831385    REAL, DIMENSION(:), ALLOCATABLE       ::  sums_divnew_l, sums_divold_l, &
    1384                                               weight_substep
     1386                                              weight_substep, weight_pres
    13851387    REAL, DIMENSION(:,:), ALLOCATABLE     ::  sums, sums_wsts_bc_l,        &
    13861388                                              sums_wsus_ws_l, sums_wsvs_ws_l,&
  • palm/trunk/SOURCE/pres.f90

    r668 r673  
    44! Current revisions:
    55! -----------------
    6 !
     6! Weighting coefficients added for right computation of the pressure during
     7! Runge-Kutta substeps.
    78! Former revisions:
    89! -----------------
     
    7273    INTEGER ::  i, j, k, sr
    7374
    74     REAL    ::  ddt_3d, localsum, threadsum
     75    REAL    ::  ddt_3d, localsum, threadsum, d_weight_pres
    7576
    7677    REAL, DIMENSION(1:2) ::  volume_flow_l, volume_flow_offset
     
    8283
    8384    ddt_3d = 1.0 / dt_3d
     85    d_weight_pres = 1. / weight_pres(intermediate_timestep_count)
    8486
    8587!
     
    253255             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
    254256                          ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
    255                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d
     257                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
     258                                                                * d_weight_pres
    256259          ENDDO
    257260!
     
    268271                                       - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
    269272                                         sums(k+1,4)                         &
    270                                        ) * ddzw(k+1) * ddt_3d
     273                                       ) * ddzw(k+1) * ddt_3d * d_weight_pres
    271274          ENDIF
    272275
     
    291294             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
    292295                          ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
    293                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d
     296                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
     297                                                                * d_weight_pres
    294298             ENDDO
    295299          ENDDO
     
    307311                           - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
    308312                             sums(k+1,4)                         &
    309                                         ) * ddzw(k+1) * ddt_3d
     313                                        ) * ddzw(k+1) * ddt_3d   &
     314                                          * d_weight_pres
    310315          ENDDO
    311316       ENDDO
     
    321326                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
    322327                          ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
    323                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d
     328                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
     329                                                                * d_weight_pres
    324330             ENDDO
    325331          ENDDO
     
    340346       ENDDO
    341347    ENDDO
    342     localsum = ( localsum + threadsum ) * dt_3d
     348    localsum = ( localsum + threadsum ) * dt_3d                               &
     349               * weight_pres(intermediate_timestep_count)
    343350    !$OMP END PARALLEL
    344351#endif
     
    496503          DEALLOCATE( p )
    497504          ALLOCATE( p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    498           DO  i = nxl, nxr
    499              DO  j = nys, nyn
    500                 DO  k = nzb_s_inner(j,i), nzt
    501                    p(k,j,i) = tend(k,j,i)
    502                 ENDDO
    503              ENDDO
    504           ENDDO
     505         
    505506       ENDIF
    506507
     
    510511!-- Store perturbation pressure on array p, used in the momentum equations
    511512    IF ( psolver(1:7) == 'poisfft' )  THEN
    512 !
    513 !--    Here, only the values from the left and right boundaries are copied
    514 !--    The remaining values are copied in the following loop due to speed
    515 !--    optimization
    516        !$OMP PARALLEL DO
    517        DO  j = nysg, nyng
    518           DO  k = nzb, nzt+1
    519              p(k,j,nxlg:nxl-1) = tend(k,j,nxlg:nxl-1)
    520              p(k,j,nxr+1:nxrg) = tend(k,j,nxr+1:nxrg)
    521           ENDDO
    522        ENDDO
     513
     514       IF ( intermediate_timestep_count == 1 )  THEN
     515          !$OMP PARALLEL PRIVATE (i,j,k)
     516          !$OMP DO
     517          DO  i = nxlg, nxrg
     518             DO  j = nysg, nyng
     519                DO  k = nzb, nzt+1
     520                   p(k,j,i) = tend(k,j,i)                                     &
     521                            * weight_substep(intermediate_timestep_count)
     522                   p(k,j,i) = tend(k,j,i)                                     &
     523                            * weight_substep(intermediate_timestep_count)
     524                ENDDO
     525             ENDDO
     526          ENDDO
     527          !$OMP END PARALLEL
     528         
     529       ELSE
     530          !$OMP PARALLEL PRIVATE (i,j,k)
     531          !$OMP DO
     532          DO  i = nxlg, nxrg
     533             DO  j = nysg, nyng
     534                DO  k = nzb, nzt+1
     535                   p(k,j,i) = p(k,j,i) + tend(k,j,i)                          &
     536                            * weight_substep(intermediate_timestep_count)
     537                   p(k,j,i) = p(k,j,i) + tend(k,j,i)                          &
     538                            * weight_substep(intermediate_timestep_count)
     539                ENDDO
     540             ENDDO
     541          ENDDO
     542          !$OMP END PARALLEL
     543         
     544       ENDIF
     545       
    523546    ENDIF
    524547
     
    533556    !$OMP PARALLEL PRIVATE (i,j,k)
    534557    !$OMP DO
    535     DO  i = nxl, nxr
    536        IF ( psolver(1:7) == 'poisfft' )  THEN
    537           DO  j = nysg, nyng
    538              DO  k = nzb, nzt+1
    539                 p(k,j,i) = tend(k,j,i)
    540              ENDDO
    541           ENDDO
    542        ENDIF
     558    DO  i = nxl, nxr   
    543559       DO  j = nys, nyn
    544560          DO  k = nzb_w_inner(j,i)+1, nzt
    545              w(k,j,i) = w(k,j,i) - dt_3d * &
    546                                    ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1)
     561             w(k,j,i) = w(k,j,i) - dt_3d *                                    &
     562                        ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1)           &
     563                        * weight_pres(intermediate_timestep_count)
    547564          ENDDO
    548565          DO  k = nzb_u_inner(j,i)+1, nzt
    549              u(k,j,i) = u(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j,i-1) ) * ddx
     566             u(k,j,i) = u(k,j,i) - dt_3d *                                 &
     567                        ( tend(k,j,i) - tend(k,j,i-1) ) * ddx              &
     568                        * weight_pres(intermediate_timestep_count)
    550569          ENDDO
    551570          DO  k = nzb_v_inner(j,i)+1, nzt
    552              v(k,j,i) = v(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j-1,i) ) * ddy
    553           ENDDO
    554 
     571             v(k,j,i) = v(k,j,i) - dt_3d *                                 &
     572                        ( tend(k,j,i) - tend(k,j-1,i) ) * ddy              &
     573                        * weight_pres(intermediate_timestep_count)
     574          ENDDO                                                         
    555575!
    556576!--       Sum up the volume flow through the right and north boundary
     
    575595    ENDDO
    576596    !$OMP END PARALLEL
    577 
     597   
     598    IF ( psolver == 'multigrid' .OR. psolver == 'sor' )  THEN       
     599       IF ( intermediate_timestep_count == 1 .OR. simulated_time == 0)  THEN
     600          !$OMP PARALLEL PRIVATE (i,j,k)
     601          !$OMP DO
     602          DO  i = nxl, nxr
     603             DO  j = nys, nyn
     604                DO  k = nzb, nzt+1
     605                   p_sub(k,j,i) = tend(k,j,i)                                 &
     606                                * weight_substep(intermediate_timestep_count)
     607                ENDDO
     608             ENDDO
     609          ENDDO
     610          !$OMP END PARALLEL
     611       ELSE
     612          !$OMP PARALLEL PRIVATE (i,j,k)
     613          !$OMP DO
     614          DO  i = nxl, nxr
     615             DO  j = nys, nyn
     616                DO  k = nzb, nzt+1
     617                   p_sub(k,j,i) = p_sub(k,j,i) + tend(k,j,i)                  &
     618                                * weight_substep(intermediate_timestep_count)
     619                ENDDO
     620             ENDDO
     621          ENDDO
     622          !$OMP END PARALLEL
     623       ENDIF
     624         
     625       IF ( intermediate_timestep_count == intermediate_timestep_count_max )  &
     626          THEN
     627          !$OMP PARALLEL PRIVATE (i,j,k)
     628          !$OMP DO
     629          DO  i = nxl, nxr
     630             DO  j = nys, nyn
     631                DO  k = nzb, nzt+1
     632                   p(k,j,i) = p_sub(k,j,i)
     633                ENDDO
     634             ENDDO
     635          ENDDO
     636          !$OMP END PARALLEL
     637       ENDIF
     638    ENDIF
     639 
    578640!
    579641!-- Resize tend to its normal size in case of multigrid and ws-scheme.
  • palm/trunk/SOURCE/prognostic_equations.f90

    r669 r673  
    44! Current revisions:
    55! -----------------
    6 !
     6! Consideration of the pressure gradient (steered by tsc(4)) during the time
     7! integration removed.
    78! Former revisions:
    89! -----------------
     
    212213                          dt_3d * (                                            &
    213214                                   tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
    214                                  - tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx    &
    215215                                  ) -                                          &
    216216                           tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
     
    297297                          dt_3d * (                                            &
    298298                                   tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
    299                                  - tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy    &
    300299                                  ) -                                          &
    301300                          tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
     
    383382                          dt_3d * (                                            &
    384383                                   tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
    385                               - tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
    386384                                  ) -                                          &
    387385                          tsc(5) * rdf(k) * w(k,j,i)
     
    899897!------------------------------------------------------------------------------!
    900898! Version with one optimized loop over all equations. It is only allowed to
    901 ! be called for the standard Piascek-Williams advection scheme.
     899! be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
    902900!
    903901! Here the calls of most subroutines are embedded in two DO loops over i and j,
     
    983981                             dt_3d * (                                         &
    984982                                   tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
    985                                  - tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx    &
    986983                                     ) -                                       &
    987984                             tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
     
    10501047                             dt_3d * (                                         &
    10511048                                   tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
    1052                                  - tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy    &
    10531049                                     ) -                                       &
    10541050                             tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
     
    11161112                          dt_3d * (                                         &
    11171113                                tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
    1118                            - tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
    11191114                                  ) -                                       &
    11201115                          tsc(5) * rdf(k) * w(k,j,i)
     
    15251520                          dt_3d * (                                            &
    15261521                                   tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
    1527                                  - tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx    &
    15281522                                  ) -                                          &
    15291523                          tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
     
    16171611                          dt_3d * (                                            &
    16181612                                   tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
    1619                                  - tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy    &
    16201613                                  ) -                                          &
    16211614                          tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
     
    17061699                          dt_3d * (                                            &
    17071700                                   tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
    1708                               - tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
    17091701                                  ) -                                          &
    17101702                          tsc(5) * rdf(k) * w(k,j,i)
  • palm/trunk/SOURCE/timestep_scheme_steering.f90

    r484 r673  
    44! Current revisions:
    55! -----------------
    6 !
     6! No pressure term during time integration (tsc(4)=0.0).
    77!
    88! Former revisions:
     
    6060             IF ( timestep_scheme == 'leapfrog+euler'  .OR. &
    6161                  timestep_scheme == 'euler' .OR. simulated_time == 0.0 )  THEN
    62                 tsc(1:5) = (/ 1.0, 1.0, 0.0, 1.0, 1.0 /)
     62                tsc(1:5) = (/ 1.0, 1.0, 0.0, 0.0, 1.0 /)
    6363             ELSE
    64                 tsc(1:5) = (/ 0.0, 2.0, 0.0, 1.0, 2.0 /)
     64                tsc(1:5) = (/ 0.0, 2.0, 0.0, 0.0, 2.0 /)
    6565             ENDIF
    6666          ELSE
     
    6969!--          user.
    7070             IF ( timestep_scheme == 'euler' )  THEN
    71                 tsc(1:5) = (/ 1.0, 1.0, 0.0, 1.0, 1.0 /)
     71                tsc(1:5) = (/ 1.0, 1.0, 0.0, 0.0, 1.0 /)
    7272             ELSE
    73                 tsc(1:5) = (/ 0.0, 2.0, 0.0, 1.0, 2.0 /)
     73                tsc(1:5) = (/ 0.0, 2.0, 0.0, 0.0, 2.0 /)
    7474             ENDIF
    7575          ENDIF
     
    8484          IF ( simulated_time == 0.0 )  THEN
    8585             dt_changed = .TRUE.
    86              tsc(1:5) = (/ 1.0, 1.0, 0.0, 1.0, 1.0 /)
     86             tsc(1:5) = (/ 1.0, 1.0, 0.0, 0.0, 1.0 /)
    8787          ELSE
    8888             dt_changed = .FALSE.
    8989             IF ( timestep_scheme == 'euler' )  THEN
    90                 tsc(1:5) = (/ 1.0, 1.0, 0.0, 1.0, 1.0 /)
     90                tsc(1:5) = (/ 1.0, 1.0, 0.0, 0.0, 1.0 /)
    9191             ELSE
    92                 tsc(1:5) = (/ 0.0, 2.0, 0.0, 1.0, 2.0 /)
     92                tsc(1:5) = (/ 0.0, 2.0, 0.0, 0.0, 2.0 /)
    9393             ENDIF
    9494          ENDIF
Note: See TracChangeset for help on using the changeset viewer.