Ignore:
Timestamp:
May 27, 2016 2:35:57 PM (8 years ago)
Author:
raasch
Message:

bugfixes for calculating run control quantities, bugfix for calculating pressure with fft-method in case of Neumann conditions both at bottom and top, steering of pres modified, ocean mode now uses initial density profile as reference in the buoyancy term

File:
1 edited

Legend:

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

    r1909 r1918  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! sum of divergence is also calculated when pres is called before the initial
     22! first time step,
     23! stearing is modified by using intermediate_timestep_count = 0 in order to
     24! determine if pres is called before the first initial timestep,
     25! bugfix: in case of Neumann conditions for pressure both at bottom and top,
     26!         mean vertical velocity is also removed for the first time step
     27! bugfix for calculating divergences
    2228!
    2329! Former revisions:
     
    114120               nest_domain, nest_bound_l, nest_bound_n, nest_bound_r,          &
    115121               nest_bound_s, on_device, outflow_l, outflow_n, outflow_r,       &
    116                outflow_s, psolver, simulated_time, subdomain_size, topography, &
    117                volume_flow, volume_flow_area, volume_flow_initial
     122               outflow_s, psolver, subdomain_size, topography, volume_flow,    &
     123               volume_flow_area, volume_flow_initial
    118124
    119125    USE cpulog,                                                                &
     
    149155
    150156    REAL(wp)     ::  ddt_3d         !<
     157    REAL(wp)     ::  d_weight_pres  !<
    151158    REAL(wp)     ::  localsum       !<
    152159    REAL(wp)     ::  threadsum      !<
    153     REAL(wp)     ::  d_weight_pres  !<
     160    REAL(wp)     ::  weight_pres_l  !<
    154161
    155162    REAL(wp), DIMENSION(1:3)   ::  volume_flow_l       !<
     
    162169
    163170
     171!
     172!-- Calculate quantities to be used locally
    164173    ddt_3d = 1.0_wp / dt_3d
    165     d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
     174    IF ( intermediate_timestep_count == 0 )  THEN
     175!
     176!--    If pres is called before initial time step
     177       weight_pres_l = 1.0_wp
     178       d_weight_pres = 1.0_wp
     179    ELSE
     180       weight_pres_l = weight_pres(intermediate_timestep_count)
     181       d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
     182    ENDIF
    166183
    167184!
     
    187204       ENDIF
    188205       
    189     ELSEIF ( psolver == 'sor'  .AND.  intermediate_timestep_count == 1 )  THEN
     206    ELSEIF ( psolver == 'sor'  .AND.  intermediate_timestep_count <= 1 )  THEN
    190207
    191208!
     
    283300!
    284301!-- Remove mean vertical velocity in case that Neumann conditions are
    285 !-- used both at bottom and top boundary, and if not a nested domain
    286     IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.  .NOT. nest_domain)  THEN
    287        IF ( simulated_time > 0.0_wp )  THEN ! otherwise nzb_w_inner not yet known
    288           w_l = 0.0_wp;  w_l_l = 0.0_wp
    289           DO  i = nxl, nxr
    290              DO  j = nys, nyn
    291                 DO  k = nzb_w_inner(j,i)+1, nzt
    292                    w_l_l(k) = w_l_l(k) + w(k,j,i)
    293                 ENDDO
    294              ENDDO
    295           ENDDO
     302!-- used both at bottom and top boundary, and if not a nested domain.
     303!-- This cannot be done before the first initial time step because ngp_2dh_outer
     304!-- is not yet known then.
     305    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.  .NOT. nest_domain  .AND.    &
     306         intermediate_timestep_count /= 0 )                                    &
     307    THEN
     308       w_l = 0.0_wp;  w_l_l = 0.0_wp
     309       DO  i = nxl, nxr
     310          DO  j = nys, nyn
     311             DO  k = nzb_w_inner(j,i)+1, nzt
     312                w_l_l(k) = w_l_l(k) + w(k,j,i)
     313             ENDDO
     314          ENDDO
     315       ENDDO
    296316#if defined( __parallel )   
    297           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    298           CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
    299                               comm2d, ierr )
     317       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
     318       CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
     319                           comm2d, ierr )
    300320#else
    301           w_l = w_l_l 
     321       w_l = w_l_l
    302322#endif
    303           DO  k = 1, nzt
    304              w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
    305           ENDDO
    306           DO  i = nxlg, nxrg
    307              DO  j = nysg, nyng
    308                 DO  k = nzb_w_inner(j,i)+1, nzt
    309                    w(k,j,i) = w(k,j,i) - w_l(k)
    310                 ENDDO
    311              ENDDO
    312           ENDDO
    313        ENDIF
     323       DO  k = 1, nzt
     324          w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
     325       ENDDO
     326       DO  i = nxlg, nxrg
     327          DO  j = nysg, nyng
     328             DO  k = nzb_w_inner(j,i)+1, nzt
     329                w(k,j,i) = w(k,j,i) - w_l(k)
     330             ENDDO
     331          ENDDO
     332       ENDDO
    314333    ENDIF
    315334
     
    363382    ENDDO
    364383
    365     IF ( intermediate_timestep_count == intermediate_timestep_count_max )  THEN
    366        localsum = localsum + threadsum * dt_3d * &
    367                              weight_pres(intermediate_timestep_count)
     384    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
     385         intermediate_timestep_count == 0 )  THEN
     386       localsum = localsum + threadsum * dt_3d * weight_pres_l
    368387    ENDIF
    369388    !$OMP END PARALLEL
     
    390409!-- Compute possible PE-sum of divergences for flow_statistics. Carry out
    391410!-- computation only at last Runge-Kutta substep.
    392     IF ( intermediate_timestep_count == intermediate_timestep_count_max )  THEN
     411    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
     412         intermediate_timestep_count == 0 )  THEN
    393413       !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
    394414       !$OMP DO SCHEDULE( STATIC )
     
    402422       ENDDO
    403423       !$acc end parallel loop
    404        localsum = localsum + threadsum * dt_3d * &
    405                              weight_pres(intermediate_timestep_count)
     424       localsum = localsum + threadsum * dt_3d * weight_pres_l
    406425       !$OMP END PARALLEL
    407426    ENDIF
     
    411430!-- For completeness, set the divergence sum of all statistic regions to those
    412431!-- of the total domain
    413     IF ( intermediate_timestep_count == intermediate_timestep_count_max )      &
     432    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
     433         intermediate_timestep_count == 0 )  THEN
    414434       sums_divold_l(0:statistic_regions) = localsum
     435    ENDIF
    415436
    416437    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
     
    564585!-- Store perturbation pressure on array p, used for pressure data output.
    565586!-- Ghost layers are added in the output routines (except sor-method: see below)
    566     IF ( intermediate_timestep_count == 1 )  THEN
     587    IF ( intermediate_timestep_count <= 1 )  THEN
    567588       !$OMP PARALLEL PRIVATE (i,j,k)
    568589       !$OMP DO
     
    616637    !$OMP PARALLEL PRIVATE (i,j,k)
    617638    !$OMP DO
    618     !$acc kernels present( ddzu, nzb_u_inner, nzb_v_inner, nzb_w_inner, tend, u, v, w, weight_pres )
     639    !$acc kernels present( ddzu, nzb_u_inner, nzb_v_inner, nzb_w_inner, tend, u, v, w )
    619640    !$acc loop independent
    620641    DO  i = nxl, nxr   
     
    626647                w(k,j,i) = w(k,j,i) - dt_3d *                                 &
    627648                           ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) *      &
    628                            weight_pres(intermediate_timestep_count)
     649                           weight_pres_l
    629650             ENDIF
    630651          ENDDO
     
    634655                u(k,j,i) = u(k,j,i) - dt_3d *                                 &
    635656                           ( tend(k,j,i) - tend(k,j,i-1) ) * ddx *            &
    636                            weight_pres(intermediate_timestep_count)
     657                           weight_pres_l
    637658             ENDIF
    638659          ENDDO
     
    642663                v(k,j,i) = v(k,j,i) - dt_3d *                                 &
    643664                           ( tend(k,j,i) - tend(k,j-1,i) ) * ddy *            &
    644                            weight_pres(intermediate_timestep_count)
     665                           weight_pres_l
    645666             ENDIF
    646667          ENDDO                                                         
     
    737758!-- A possible PE-sum is computed in flow_statistics. Carry out computation
    738759!-- only at last Runge-Kutta step.
    739     IF ( intermediate_timestep_count == intermediate_timestep_count_max )  THEN
     760    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
     761         intermediate_timestep_count == 0 )  THEN
    740762       CALL cpu_log( log_point_s(1), 'divergence', 'start' )
    741763       sums_divnew_l = 0.0_wp
Note: See TracChangeset for help on using the changeset viewer.