Ignore:
Timestamp:
Jan 18, 2011 4:19:48 PM (11 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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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.
Note: See TracChangeset for help on using the changeset viewer.