Changeset 1918 for palm/trunk/SOURCE


Ignore:
Timestamp:
May 27, 2016 2:35:57 PM (9 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

Location:
palm/trunk/SOURCE
Files:
6 edited

Legend:

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

    r1917 r1918  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! setting of a fixed reference state ('single_value') for ocean runs removed
    2222!
    2323! Former revisions:
     
    15261526
    15271527!
    1528 !-- Ocean runs always use reference values in the buoyancy term
    1529     IF ( ocean )  THEN
    1530        reference_state = 'single_value'
    1531        use_single_reference_value = .TRUE.
    1532     ENDIF
    1533 
    1534 !
    15351528!-- Sign of buoyancy/stability terms
    15361529    IF ( ocean )  atmos_ocean_sign = -1.0_wp
  • palm/trunk/SOURCE/flow_statistics.f90

    r1854 r1918  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! in case of Wicker-Skamarock scheme, calculate disturbance kinetic energy here,
     22! if flow_statistics is called before the first initial time step
    2223!
    2324! Former revisions:
     
    630631!--    Computation of statistics when ws-scheme is not used. Else these
    631632!--    quantities are evaluated in the advection routines.
    632        IF ( .NOT. ws_scheme_mom .OR. sr /= 0 )  THEN
     633       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
     634       THEN
    633635          !$OMP DO
    634636          DO  i = nxl, nxr
     
    16231625
    16241626!
    1625 !-- If required, sum up horizontal averages for subsequent time averaging
    1626     IF ( do_sum )  THEN
     1627!-- If required, sum up horizontal averages for subsequent time averaging.
     1628!-- Do not sum, if flow statistics is called before the first initial time step.
     1629    IF ( do_sum  .AND.  simulated_time /= 0.0_wp )  THEN
    16271630       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
    16281631       hom_sum = hom_sum + hom(:,1,:,:)
     
    22302233!--    Computation of statistics when ws-scheme is not used. Else these
    22312234!--    quantities are evaluated in the advection routines.
    2232        IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0 )  THEN
     2235       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
     2236       THEN
    22332237
    22342238          !$OMP DO
     
    35663570!
    35673571!-- If required, sum up horizontal averages for subsequent time averaging
    3568     IF ( do_sum )  THEN
     3572!-- Do not sum, if flow statistics is called before the first initial time step.
     3573    IF ( do_sum  .AND.  simulated_time /= 0.0_wp )  THEN
    35693574       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
    35703575       hom_sum = hom_sum + hom(:,1,:,:)
  • palm/trunk/SOURCE/init_3d_model.f90

    r1917 r1918  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! intermediate_timestep_count is set 0 instead 1 for first call of pres,
     22! bugfix: initialization of local sum arrays are moved to the beginning of the
     23!         routine because otherwise results from pres are overwritten
    2224!
    2325! Former revisions:
     
    709711    weight_substep = 1.0_wp
    710712    weight_pres    = 1.0_wp
    711     intermediate_timestep_count = 1  ! needed when simulated_time = 0.0
     713    intermediate_timestep_count = 0  ! needed when simulated_time = 0.0
    712714       
    713715    CALL location_message( 'finished', .TRUE. )
     716
     717!
     718!-- Initialize local summation arrays for routine flow_statistics.
     719!-- This is necessary because they may not yet have been initialized when they
     720!-- are called from flow_statistics (or - depending on the chosen model run -
     721!-- are never initialized)
     722    sums_divnew_l      = 0.0_wp
     723    sums_divold_l      = 0.0_wp
     724    sums_l_l           = 0.0_wp
     725    sums_up_fraction_l = 0.0_wp
     726    sums_wsts_bc_l     = 0.0_wp
     727
     728
    714729!
    715730!-- Initialize model variables
     
    18211836       ENDDO
    18221837    ENDIF
    1823 
    1824 !
    1825 !-- Initialize local summation arrays for routine flow_statistics.
    1826 !-- This is necessary because they may not yet have been initialized when they
    1827 !-- are called from flow_statistics (or - depending on the chosen model run -
    1828 !-- are never initialized)
    1829     sums_divnew_l      = 0.0_wp
    1830     sums_divold_l      = 0.0_wp
    1831     sums_l_l           = 0.0_wp
    1832     sums_up_fraction_l = 0.0_wp
    1833     sums_wsts_bc_l     = 0.0_wp
    18341838
    18351839!
  • palm/trunk/SOURCE/modules.f90

    r1907 r1918  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! default timestep switched from -1.0 to +1.0 in order to avoid wrong sign of
     22! initially calculated divergence
    2223!
    2324! Former revisions:
     
    751752                 dt_max = 20.0_wp, &
    752753                 dt_restart = 9999999.9_wp, &
    753                  dt_run_control = 60.0_wp, dt_3d = -1.0_wp, dz = -1.0_wp, &
     754                 dt_run_control = 60.0_wp, dt_3d = 1.0_wp, dz = -1.0_wp, &
    754755                 dz_max = 9999999.9_wp, dz_stretch_factor = 1.08_wp, &
    755756                 dz_stretch_level = 100000.0_wp, e_init = 0.0_wp, e_min = 0.0_wp, &
  • 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
  • palm/trunk/SOURCE/time_integration.f90

    r1917 r1918  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! determination of time step moved to the end of the time step loop,
     22! the first time step is now always calculated before the time step loop (i.e.
     23! also in case of restart runs)
    2224!
    2325! Former revisions:
     
    260262    USE indices,                                                               &
    261263        ONLY:  i_left, i_right, j_north, j_south, nbgp, nx, nxl, nxlg, nxr,    &
    262                nxrg, nyn, nys, nzb, nzt, nzb_u_inner, nzb_v_inner
     264               nxrg, nyn, nyng, nys, nysg, nzb, nzt, nzb_u_inner, nzb_v_inner
    263265
    264266    USE interaction_droplets_ptq_mod,                                          &
    265267        ONLY:  interaction_droplets_ptq
     268
     269    USE interfaces
    266270
    267271    USE kinds
     
    311315
    312316    USE statistics,                                                            &
    313         ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l
     317        ONLY:  flow_statistics_called, hom, pr_palm, sums_ls_l, u_max,         &
     318               u_max_ijk, v_max, v_max_ijk, w_max, w_max_ijk
    314319
    315320    USE surface_layer_fluxes_mod,                                              &
     
    326331    CHARACTER (LEN=9) ::  time_to_string          !<
    327332
    328 !
    329 !-- At the beginning of a simulation determine the time step as well as
    330 !-- determine and print out the run control parameters
    331     IF ( simulated_time == 0.0_wp )  CALL timestep
     333    REAL(wp) ::  dt_3d_old  !< temporary storage of timestep to be used for
     334                            !< steering of run control output interval
     335
     336!
     337!-- At beginning determine the first time step
     338    CALL timestep
    332339
    333340!
     
    340347    ENDIF
    341348
     349!
     350!-- Determine and print out the run control quantities before the first time
     351!-- step of this run. For the initial run, some statistics (e.g. divergence)
     352!-- need to be determined first.
     353    IF ( simulated_time == 0.0_wp )  CALL flow_statistics
    342354    CALL run_control
    343355
     
    349361!--    In case of model termination initiated by the local model the coupler
    350362!--    must not be called because this would again cause an MPI hang.
    351        DO WHILE ( time_coupling >= dt_coupling .AND. terminate_coupled == 0 )
     363       DO WHILE ( time_coupling >= dt_coupling  .AND. terminate_coupled == 0 )
    352364          CALL surface_coupler
    353365          time_coupling = time_coupling - dt_coupling
    354366       ENDDO
    355        IF (time_coupling == 0.0_wp .AND. time_since_reference_point < dt_coupling)&
     367       IF (time_coupling == 0.0_wp  .AND.                                      &
     368           time_since_reference_point < dt_coupling )                          &
    356369       THEN
    357370          time_coupling = time_since_reference_point
     
    372385
    373386       CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
    374 
    375 !
    376 !--    Determine size of next time step
    377        IF ( simulated_time /= 0.0_wp )  THEN
    378           CALL timestep
    379 !
    380 !--       Synchronize the timestep in case of nested run.
    381           IF ( nested_run )  THEN
    382 !
    383 !--          Synchronize by unifying the time step.
    384 !--          Global minimum of all time-steps is used for all.     
    385              CALL pmci_synchronize
    386           ENDIF
    387        ENDIF
    388387
    389388!
     
    10261025
    10271026!
    1028 !--    Computation and output of run control parameters.
    1029 !--    This is also done whenever perturbations have been imposed
    1030        IF ( time_run_control >= dt_run_control  .OR.                     &
    1031             timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created ) &
    1032        THEN
    1033           CALL run_control
    1034           IF ( time_run_control >= dt_run_control )  THEN
    1035              time_run_control = MOD( time_run_control, &
    1036                                      MAX( dt_run_control, dt_3d ) )
    1037           ENDIF
    1038        ENDIF
    1039 
    1040 !
    10411027!--    Profile output (ASCII) on file
    10421028       IF ( time_dopr_listing >= dt_dopr_listing )  THEN
     
    11491135
    11501136!
     1137!--    Determine size of next time step. Save timestep dt_3d because it is
     1138!--    newly calculated in routine timestep, but required further below for
     1139!--    steering the run control output interval
     1140       dt_3d_old = dt_3d
     1141       CALL timestep
     1142
     1143!
     1144!--    Computation and output of run control parameters.
     1145!--    This is also done whenever perturbations have been imposed
     1146       IF ( time_run_control >= dt_run_control  .OR.                     &
     1147            timestep_scheme(1:5) /= 'runge'  .OR.  disturbance_created ) &
     1148       THEN
     1149!          IF ( current_timestep_number == 1 )  THEN
     1150!             IF ( nxl < 7 .AND.  nxr > 7 .AND. nys < 7 .AND. nyn > 7 )  THEN
     1151!                u(10,7,7) = 0.55
     1152!             ENDIF
     1153!             PRINT*, 'calculating minmax'
     1154!             CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, u,       &
     1155!                                  'abs', 0.0_wp, u_max, u_max_ijk )
     1156!             CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, v,       &
     1157!                                  'abs', 0.0_wp, v_max, v_max_ijk )
     1158!             CALL global_min_max( nzb, nzt+1, nysg, nyng, nxlg, nxrg, w,       &
     1159!                                  'abs', 0.0_wp, w_max, w_max_ijk )
     1160!             PRINT*, 'calculated u_max = ', u_max, '  myid = ', myid
     1161!          ENDIF
     1162          CALL run_control
     1163          IF ( time_run_control >= dt_run_control )  THEN
     1164             time_run_control = MOD( time_run_control, &
     1165                                     MAX( dt_run_control, dt_3d_old ) )
     1166          ENDIF
     1167       ENDIF
     1168
     1169!
     1170!--    Synchronize the timestep in case of nested run.
     1171       IF ( nested_run )  THEN
     1172!
     1173!--       Synchronize by unifying the time step.
     1174!--       Global minimum of all time-steps is used for all.
     1175          CALL pmci_synchronize
     1176       ENDIF
     1177
     1178!
    11511179!--    Output elapsed simulated time in form of a progress bar on stdout
    11521180       IF ( myid == 0 )  CALL output_progress_bar
Note: See TracChangeset for help on using the changeset viewer.