Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Dec 23, 2010 12:06:00 PM (13 years ago)
Author:
suehring
Message:

summary:

Gryschka:

  • Coupling with different resolution and different numbers of PEs in ocean and atmosphere is available
  • Exchange of u and v from ocean surface to atmosphere surface
  • Mirror boundary condition for u and v at the bottom are replaced by dirichlet boundary conditions
  • Inflow turbulence is now defined by flucuations around spanwise mean
  • Bugfixes for cyclic_fill and constant_volume_flow

Suehring:

  • New advection added ( Wicker and Skamarock 5th order ), therefore:
    • New module advec_ws.f90
    • Modified exchange of ghost boundaries.
    • Modified evaluation of turbulent fluxes
    • New index bounds nxlg, nxrg, nysg, nyng

advec_ws.f90

Advection scheme for scalars and momentum using the flux formulation of
Wicker and Skamarock 5th order.
Additionally the module contains of a routine using for initialisation and
steering of the statical evaluation. The computation of turbulent fluxes takes
place inside the advection routines.
In case of vector architectures Dirichlet and Radiation boundary conditions are
outstanding and not available. Furthermore simulations within topography are
not possible so far. A further routine local_diss_ij is available and is used
if a control of dissipative fluxes is desired.

check_parameters.f90

Exchange of parameters between ocean and atmosphere via PE0
Check for illegal combination of ws-scheme and timestep scheme.
Check for topography and ws-scheme.
Check for not cyclic boundary conditions in combination with ws-scheme and
loop_optimization = 'vector'.
Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.

Different processor/grid topology in atmosphere and ocean is now allowed!
Bugfixes in checking for conserve_volume_flow_mode.

exchange_horiz.f90

Dynamic exchange of ghost points with nbgp_local to ensure that no useless
ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0) used for
normal grid, the remaining types used for the several grid levels.
Exchange is done via MPI-Vectors with a dynamic value of ghost points which
depend on the advection scheme. Exchange of left and right PEs is 10% faster
with MPI-Vectors than without.

flow_statistics.f90

When advection is computed with ws-scheme, turbulent fluxes are already
computed in the respective advection routines and buffered in arrays
sums_xxxx_ws_l(). This is due to a consistent treatment of statistics
with the numerics and to avoid unphysical kinks near the surface. So some if-
requests has to be done to dicern between fluxes from ws-scheme other advection
schemes. Furthermore the computation of z_i is only done if the heat flux
exceeds a minimum value. This affects only simulations of a neutral boundary
layer and is due to reasons of computations in the advection scheme.

inflow_turbulence.f90

Using nbgp recycling planes for a better resolution of the turbulent flow near
the inflow.

init_grid.f90

Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
Furthermore the allocation of arrays and steering of loops is done with these
parameters. Call of exchange_horiz are modified.
In case of dirichlet bounday condition at the bottom zu(0)=0.0
dzu_mg has to be set explicitly for a equally spaced grid near bottom.
ddzu_pres added to use a equally spaced grid near bottom.

init_pegrid.f90

Moved determination of target_id's from init_coupling
Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)
with different grid/processor-topology in ocean and atmosphere

Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
This distinction is due to reasons of data exchange and performance for the
normal grid and grids in poismg.
The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
New MPI-Vectors for data exchange between left and right boundaries added.
This is due to reasons of performance (10% faster).

ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values

parin.f90

Steering parameter dissipation_control added in inipar.

Makefile

Module advec_ws added.

Modules

Removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc

For coupling with different resolution in ocean and atmophere:
+nx_a, +nx_o, ny_a, +ny_o, ngp_a, ngp_o, +total_2d_o, +total_2d_a,
+coupling_topology

Buffer arrays for the left sided advective fluxes added in arrays_3d.
+flux_s_u, +flux_s_v, +flux_s_w, +diss_s_u, +diss_s_v, +diss_s_w,
+flux_s_pt, +diss_s_pt, +flux_s_e, +diss_s_e, +flux_s_q, +diss_s_q,
+flux_s_sa, +diss_s_sa
3d arrays for dissipation control added. (only necessary for vector arch.)
+var_x, +var_y, +var_z, +gamma_x, +gamma_y, +gamma_z
Default of momentum_advec and scalar_advec changed to 'ws-scheme' .
+exchange_mg added in control_parameters to steer the data exchange.
Parameters +nbgp, +nxlg, +nxrg, +nysg, +nyng added in indices.
flag array +boundary_flags added in indices to steer the degradation of order
of the advective fluxes when non-cyclic boundaries are used.
MPI-datatypes +type_y, +type_y_int and +type_yz for data_exchange added in
pegrid.
+sums_wsus_ws_l, +sums_wsvs_ws_l, +sums_us2_ws_l, +sums_vs2_ws_l,
+sums_ws2_ws_l, +sums_wspts_ws_l, +sums_wssas_ws_l, +sums_wsqs_ws_l
and +weight_substep added in statistics to steer the statistical evaluation
of turbulent fluxes in the advection routines.
LOGICALS +ws_scheme_sca and +ws_scheme_mom added to get a better performance
in prognostic_equations.
LOGICAL +dissipation_control control added to steer numerical dissipation
in ws-scheme.

Changed length of string run_description_header

pres.f90

New allocation of tend when ws-scheme and multigrid is used. This is due to
reasons of perforance of the data_exchange. The same is done with p after
poismg is called.
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
multigrid is used. Calls of exchange_horiz are modified.

bugfix: After pressure correction no volume flow correction in case of
non-cyclic boundary conditions
(has to be done only before pressure correction)

Call of SOR routine is referenced with ddzu_pres.

prognostic_equations.f90

Calls of the advection routines with WS5 added.
Calls of ws_statistics added to set the statistical arrays to zero after each
time step.

advec_particles.f90

Declaration of de_dx, de_dy, de_dz adapted to additional ghost points.
Furthermore the calls of exchange_horiz were modified.

asselin_filter.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

average_3d_data.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

boundary_conds.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
Removed mirror boundary conditions for u and v at the bottom in case of
ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
in init_3d_model

calc_liquid_water_content.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

calc_spectra.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for
allocation of tend.

check_open.f90

Output of total array size was adapted to nbgp.

data_output_2d.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays local_2d and total_2d.
Calls of exchange_horiz are modified.

data_output_2d.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Skip-value skip_do_avs changed to a dynamic adaption of ghost points.

data_output_mask.f90

Calls of exchange_horiz are modified.

diffusion_e.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_s.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_u.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_v.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_w.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

exchange_horiz_2d.f90

Dynamic exchange of ghost points with nbgp, which depends on the advection
scheme. Exchange between left and right PEs is now done with MPI-vectors.

global_min_max.f90

Adapting of the index arrays, because MINLOC assumes lowerbound
at 1 and not at nbgp.

init_3d_model.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Call ws_init to initialize arrays needed for statistical evaluation and
optimization when ws-scheme is used.
Initial volume flow is now calculated by using the variable hom_sum.
Therefore the correction of initial volume flow for non-flat topography
removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
representative for the height z0

Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)

init_coupling.f90

determination of target_id's moved to init_pegrid

init_pt_anomaly.f90

Call of exchange_horiz are modified.

init_rankine.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

init_slope.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

header.f90

Output of advection scheme.

poismg.f90

Calls of exchange_horiz are modified.

prandtl_fluxes.f90

Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

production_e.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

read_3d_binary.f90

+/- 1 replaced with +/- nbgp when swapping and allocating variables.

sor.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Call of exchange_horiz are modified.
bug removed in declaration of ddzw(), nz replaced by nzt+1

subsidence.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

sum_up_3d_data.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

surface_coupler.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in
MPI_SEND() and MPI_RECV.
additional case for nonequivalent processor and grid topopolgy in ocean and
atmosphere added (coupling_topology = 1)

Added exchange of u and v from Ocean to Atmosphere

time_integration.f90

Calls of exchange_horiz are modified.
Adaption to slooping surface.

timestep.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_3d_data_averaging.f90, user_data_output_2d.f90, user_data_output_3d.f90,
user_actions.f90, user_init.f90, user_init_plant_canopy.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_read_restart_data.f90

Allocation with nbgp.

wall_fluxes.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

write_compressed.f90

Array bounds and nx, ny adapted with nbgp.

sor.f90

bug removed in declaration of ddzw(), nz replaced by nzt+1

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE

    • Property svn:mergeinfo set to (toggle deleted branches)
      /palm/branches/suehring423-666
      /palm/branches/letzel/masked_output/SOURCE296-409

  • palm/trunk/SOURCE/flow_statistics.f90

    r625 r667  
    44! Current revisions:
    55! -----------------
    6 !
     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!
    716!
    817! Former revisions:
    … …  
    102111    REAL    ::  sums_ll(nzb:nzt+1,2)
    103112
    104 
    105113    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
    106114
    … …  
    133141       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
    134142
     143!
     144!--    Copy the turbulent quantities, evaluated in the advection routines to
     145!--    the local array sums_l() for further computations
     146       IF ( ws_scheme_mom )  THEN
     147!       
     148!--       Boundary condition for u'u' and v'v', because below the surface no
     149!--       computation for these quantities is done.
     150          DO  i = nxl, nxr
     151             DO  j =  nys, nyn
     152                sums_us2_ws_l(nzb_u_inner(j,i),sr) =                          &
     153                    sums_us2_ws_l(nzb_u_inner(j,i)+1,sr)
     154                sums_vs2_ws_l(nzb_v_inner(j,i),sr) =                          & 
     155                    sums_vs2_ws_l(nzb_v_inner(j,i)+1,sr)
     156             ENDDO
     157          ENDDO
     158!         
     159!--       Swap the turbulent quantities evaluated in advec_ws.
     160          sums_l(:,13,0) = sums_wsus_ws_l(:,sr)       ! w*u*
     161          sums_l(:,15,0) = sums_wsvs_ws_l(:,sr)       ! w*v*
     162          sums_l(:,30,0) = sums_us2_ws_l(:,sr)        ! u*2
     163          sums_l(:,31,0) = sums_vs2_ws_l(:,sr)        ! v*2
     164          sums_l(:,32,0) = sums_ws2_ws_l(:,sr)        ! w*2
     165          sums_l(:,34,0) = sums_l(:,34,0) + 0.5 *                             &
     166                (sums_us2_ws_l(:,sr) + sums_vs2_ws_l(:,sr)                    &
     167                + sums_ws2_ws_l(:,sr))                      ! e*
     168          DO  k = nzb, nzt
     169             sums_l(nzb+5,pr_palm,0) = sums_l(nzb+5,pr_palm,0) + 0.5 * (      &
     170                sums_us2_ws_l(k,sr) + sums_vs2_ws_l(k,sr) +                   &
     171                sums_ws2_ws_l(k,sr))
     172          ENDDO
     173       ENDIF
     174       IF ( ws_scheme_sca )  THEN
     175          sums_l(:,17,0) = sums_wspts_ws_l(:,sr)      ! w*pts* from advec_s_ws
     176          IF ( ocean ) sums_l(:,66,0) = sums_wssas_ws_l(:,sr) ! w*sa*
     177          IF ( humidity  .OR.  passive_scalar ) sums_l(:,49,0) =              &
     178                                                   sums_wsqs_ws_l(:,sr) !w*q*
     179       ENDIF
    135180!
    136181!--    Horizontally averaged profiles of horizontal velocities and temperature.
    … …  
    138183!--    for other horizontal averages.
    139184       tn = 0
     185
    140186       !$OMP PARALLEL PRIVATE( i, j, k, tn )
    141187#if defined( __intel_openmp_bug )
    … …  
    215261       ENDIF
    216262       !$OMP END PARALLEL
    217 
    218263!
    219264!--    Summation of thread sums
    … …  
    304349       hom(:,1,2,sr) = sums(:,2)             ! v
    305350       hom(:,1,4,sr) = sums(:,4)             ! pt
     351
    306352
    307353!
    … …  
    354400          DO  j =  nys, nyn
    355401             sums_l_etot = 0.0
    356              sums_l_eper = 0.0
    357402             DO  k = nzb_s_inner(j,i), nzt+1
    358                 u2   = u(k,j,i)**2
    359                 v2   = v(k,j,i)**2
    360                 w2   = w(k,j,i)**2
    361                 ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
    362                 vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
    363403!
    364404!--             Prognostic and diagnostic variables
    … …  
    369409                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
    370410
    371 !
    372 !--             Variances
    373                 sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
    374                 sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
    375                 sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
    376411                sums_l(k,33,tn) = sums_l(k,33,tn) + &
    377412                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
    … …  
    381416                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
    382417                ENDIF
    383 !
    384 !--             Higher moments
    385 !--             (Computation of the skewness of w further below)
    386                 sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
    387                                                     rmask(j,i,sr)
    388 !
    389 !--             Perturbation energy
    390                 sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 * ( ust2 + vst2 + w2 ) &
    391                                                     * rmask(j,i,sr)
     418
    392419                sums_l_etot  = sums_l_etot + &
    393                                         0.5 * ( u2 + v2 + w2 ) * rmask(j,i,sr)
    394                 sums_l_eper  = sums_l_eper + &
    395                                         0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
     420                                        0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 +    &
     421                                        w(k,j,i)**2 ) * rmask(j,i,sr)
    396422             ENDDO
    397423!
    … …  
    401427!--          allow vectorization of that loop.
    402428             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
    403              sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) + sums_l_eper
    404429!
    405430!--          2D-arrays (being collected in the last column of sums_l)
    … …  
    420445
    421446!
     447!--    Computation of statistics when ws-scheme is not used. Else these
     448!--    quantities are evaluated in the advection routines.
     449       IF ( .NOT. ws_scheme_mom )  THEN
     450          !$OMP DO
     451          DO  i = nxl, nxr
     452             DO  j =  nys, nyn
     453                sums_l_eper = 0.0
     454                DO  k = nzb_s_inner(j,i), nzt+1
     455                   u2   = u(k,j,i)**2
     456                   v2   = v(k,j,i)**2
     457                   w2   = w(k,j,i)**2
     458                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
     459                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
     460
     461                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
     462                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
     463                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
     464!
     465!--   Higher moments
     466!--  (Computation of the skewness of w further below)
     467                   sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
     468                                                    rmask(j,i,sr)
     469!
     470!--             Perturbation energy
     471
     472                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 *       &
     473                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
     474                   sums_l_eper  = sums_l_eper + &
     475                                        0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
     476
     477                ENDDO
     478                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)   &
     479                     + sums_l_eper
     480             ENDDO
     481          ENDDO
     482       ELSE
     483          !$OMP DO
     484          DO  i = nxl, nxr
     485             DO  j =  nys, nyn
     486                DO  k = nzb_s_inner(j,i), nzt + 1
     487                   w2   = w(k,j,i)**2
     488!
     489!--                Higher moments
     490!--                (Computation of the skewness of w further below)
     491                   sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
     492                                                    rmask(j,i,sr)
     493                ENDDO
     494             ENDDO
     495          ENDDO
     496       ENDIF
     497
     498!
    422499!--    Horizontally averaged profiles of the vertical fluxes
     500
    423501       !$OMP DO
    424502       DO  i = nxl, nxr
    … …  
    587665                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
    588666                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
    589 !
    590 !--             Momentum flux w*u*
    591                 sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                     &
    592                                                     ( w(k,j,i-1) + w(k,j,i) ) &
    593                                                     * ust * rmask(j,i,sr)
    594 !
    595 !--             Momentum flux w*v*
    596                 sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                     &
    597                                                     ( w(k,j-1,i) + w(k,j,i) ) &
    598                                                     * vst * rmask(j,i,sr)
    599 !
    600 !--             Heat flux w*pt*
    601 !--             The following formula (comment line, not executed) does not
    602 !--             work if applied to subregions
    603 !                sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5 *                     &
    604 !                                                    ( pt(k,j,i)+pt(k+1,j,i) ) &
    605 !                                                    * w(k,j,i) * rmask(j,i,sr)
    606                 sums_l(k,17,tn) = sums_l(k,17,tn) + pts * w(k,j,i) * &
    607                                                     rmask(j,i,sr)
    608 !
     667
    609668!--             Higher moments
    610669                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * &
    … …  
    617676!--             but so far there is no other suitable place to calculate)
    618677                IF ( ocean )  THEN
    619                    pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
     678                   IF( .NOT. ws_scheme_sca )  THEN
     679                      pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
    620680                                 sa(k+1,j,i) - hom(k+1,1,23,sr) )
    621                    sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
     681                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
    622682                                                       rmask(j,i,sr)
     683                   ENDIF
    623684                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
    624685                                                       rmask(j,i,sr)
    … …  
    635696                   sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
    636697                                                       rmask(j,i,sr)
    637                    pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
    638                                  q(k+1,j,i) - hom(k+1,1,41,sr) )
    639                    sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
    640                                                        rmask(j,i,sr)
     698
    641699                   IF ( cloud_physics  .OR.  cloud_droplets )  THEN
    642700                      pts = 0.5 *                                            &
    … …  
    652710!
    653711!--             Passive scalar flux
    654                 IF ( passive_scalar )  THEN
     712                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca ))  THEN
    655713                   pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
    656714                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
    … …  
    661719!
    662720!--             Energy flux w*e*
    663                 sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *           &
    664                                               ( ust**2 + vst**2 + w(k,j,i)**2 )&
    665                                               * rmask(j,i,sr)
    666          
     721!--             has to be adjusted
     722                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *          &
     723                                             ( ust**2 + vst**2 + w(k,j,i)**2 )&
     724                                             * rmask(j,i,sr)
    667725             ENDDO
    668726          ENDDO
    669727       ENDDO
     728!-     for reasons of speed optimization the loop is splitted, to avoid if-else
     729!-     statements inside the loop
     730!-     Fluxes which have been computed in part directly inside the advection routines
     731!-     treated seperatly.
     732!-     First treat the momentum fluxes
     733       IF ( .NOT. ws_scheme_mom )  THEN
     734         !$OMP DO
     735         DO  i = nxl, nxr
     736            DO  j = nys, nyn
     737               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
     738                  ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
     739                              u(k+1,j,i) - hom(k+1,1,1,sr) )
     740                  vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
     741                              v(k+1,j,i) - hom(k+1,1,2,sr) )
     742!                             
     743!--               Momentum flux w*u*
     744                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                   &
     745                                                    ( w(k,j,i-1) + w(k,j,i) ) &
     746                                                    * ust * rmask(j,i,sr)
     747!
     748!--               Momentum flux w*v*
     749                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                   &
     750                                                    ( w(k,j-1,i) + w(k,j,i) ) &
     751                                                    * vst * rmask(j,i,sr)
     752               ENDDO
     753            ENDDO
     754         ENDDO
     755
     756       ENDIF
     757       IF ( .NOT. ws_scheme_sca )  THEN
     758         !$OMP DO
     759         DO  i = nxl, nxr
     760            DO  j = nys, nyn
     761               DO  k = nzb_diff_s_inner(j,i) - 1, nzt_diff
     762!-                vertical heat flux
     763                  sums_l(k,17,tn) = sums_l(k,17,tn) +  0.5 * &
     764                           ( pt(k,j,i)   - hom(k,1,4,sr) + &
     765                           pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
     766                           * w(k,j,i) * rmask(j,i,sr)
     767                  IF ( humidity )  THEN
     768                     pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
     769                                q(k+1,j,i) - hom(k+1,1,41,sr) )
     770                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
     771                                                      rmask(j,i,sr)
     772                  ENDIF
     773               ENDDO
     774            ENDDO
     775         ENDDO
     776
     777       ENDIF
     778
    670779
    671780!
    … …  
    814923
    815924#if defined( __parallel )
     925
    816926!
    817927!--    Compute total sum from local sums
    … …  
    843953          sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
    844954       ENDDO
     955
    845956!--    Upstream-parts
    846957       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
    … …  
    868979          ENDDO
    869980       ENDIF
    870 
    871981!
    872982!--    Collect horizontal average in hom.
    … …  
    9341044                                       ! upstream-parts u_x, u_y, u_z, v_x,
    9351045                                       ! v_y, usw. (in last but one profile)
    936        hom(:,1,pr_palm,sr) =   sums(:,pr_palm) 
     1046       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
    9371047                                       ! u*, w'u', w'v', t* (in last profile)
    9381048
    … …  
    9511061       z_i(1) = 0.0
    9521062       first = .TRUE.
     1063
    9531064       IF ( ocean )  THEN
    9541065          DO  k = nzt, nzb+1, -1
    955              IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
     1066             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
     1067                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
    9561068                first = .FALSE.
    9571069                height = zw(k)
    9581070             ENDIF
    9591071             IF ( hom(k,1,18,sr) < 0.0  .AND. &
     1072                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
    9601073                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
    9611074                IF ( zw(k) < 1.5 * height )  THEN
    … …  
    9691082       ELSE
    9701083          DO  k = nzb, nzt-1
    971              IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
     1084             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
     1085                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
    9721086                first = .FALSE.
    9731087                height = zw(k)
    9741088             ENDIF
    975              IF ( hom(k,1,18,sr) < 0.0  .AND. &
     1089             IF ( hom(k,1,18,sr) < 0.0  .AND. &
     1090                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
    9761091                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
    9771092                IF ( zw(k) < 1.5 * height )  THEN
    … …  
    10261141!--    the characteristic convective boundary layer temperature.
    10271142!--    The horizontal average at nzb+1 is input for the average temperature.
    1028        IF ( hom(nzb,1,18,sr) > 0.0  .AND.  z_i(1) /= 0.0 )  THEN
     1143       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
     1144           .AND.  z_i(1) /= 0.0 )  THEN
    10291145          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
    10301146                                       hom(nzb,1,18,sr) *      &
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |