Changeset 57


Ignore:
Timestamp:
Mar 9, 2007 12:05:41 PM (15 years ago)
Author:
raasch
Message:

preliminary update of further changes, advec_particles is not running!

Location:
palm/trunk
Files:
21 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/DOC/app/chapter_4.1.html

    r54 r57  
    1140411404
    1140511405
    11406     <tr>
     11406    <tr><td style="vertical-align: top;"><a name="pt_reference"></a><span style="font-weight: bold;">pt_reference</span></td><td style="vertical-align: top;">R</td><td style="vertical-align: top;"><span style="font-style: italic;">use horizontal average as refrence</span></td><td style="vertical-align: top;">Reference temperature to be used in all buoyancy terms (in K).<br><br>By default, the instantaneous horizontal average over the total model domain is used.</td></tr><tr>
    1140711407
    1140811408
  • palm/trunk/DOC/app/chapter_4.6.html

    r54 r57  
    1668116681
    1668216682
    16683     <tr>
     16683    <tr><td><a style="font-weight: bold;" href="chapter_4.1.html#pt_reference">pt_reference</a></td><td>I</td><td>R</td><td><span style="font-style: italic;">use horizontal average as reference</span></td><td>Reference temperature to be used in all buoyancy terms (in K).</td></tr><tr>
    1668416684
    1668516685
  • palm/trunk/SOURCE/CURRENT_MODIFICATIONS

    r54 r57  
    11New:
    22---
     3
     4particle reflection from vertical walls implemented, particle SGS model adjusted to walls
    35
    46Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.
    57
    68new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files
     9
     10new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).
     11
     12new user interface user_advec_particles
    713
    814new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays
     
    1218samples added to the user interface which show how to add user-define time series quantities.
    1319
    14 Makefile, check_open, check_parameters, header, init_3d_model, modules, netcdf, parin, user_interface
     20unit 9 opened for debug output (file DEBUG_<pe#>)
     21
     22Makefile, advec_particles, buoyancy, check_open, check_parameters, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, user_interface, write_var_list
    1523
    1624New: wall_fluxes
     
    1927Changed:
    2028-------
     29+age_m in particle_type
    2130
    2231Move call of user_actions( 'after_integration' ) below increment of times
     
    2736Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.
    2837
    29 check_parameters, data_output_ts, flow_statistics, init_3d_model, modules, parin, time_integration
     38__vtk directives removed from main program.
     39
     40check_parameters, data_output_ts, flow_statistics, init_particles, init_3d_model, modules, palm, parin, time_integration
    3041
    3142
     
    3546Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model
    3647
     48Bugfix in sample for reading user defined data from restart file (user_init)
     49
    3750in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
    3851# case of .mod files
    3952
    4053Makefile
    41 init_1d_model
     54init_1d_model, user_interface
  • palm/trunk/SOURCE/advec_particles.f90

    r39 r57  
    44! Actual revisions:
    55! -----------------
     6! Particle reflection at vertical walls (by Jin Zhang), regarding vertical
     7! walls in the SGS model,
     8! + user interface user_advec_particles
    69! TEST: PRINT statements on unit 9 (commented out)
    710!
     
    3740! sub-timesteps (if sgs velocities are included), +3d-arrays de_dx/y/z,
    3841! izuf renamed iran, output of particle time series
    39 !
    40 ! Revision 1.28  2006/02/23 09:41:21  raasch
    41 ! Only a fraction of the particles may have tails, therefore output for unit
    42 ! 85 and 90 has changed,
    43 ! sorting of particles in gridbox order after each timestep, optimization of
    44 ! the droplet growth by collision,
    45 ! 0.5 replaced by 0.4999999999 in setting the random fluctuations,
    46 ! psl, psr, pdx, etc. are now 1d arrays, i.e. they depend on the particle group,
    47 ! improve particle release at PE boundaries, nt_anz renamed
    48 ! current_timestep_number, idum (particle_type) renamed tail_id,
    49 ! error number argument for handle_netcdf_error,
    50 !
    51 ! Revision 1.27  2005/10/20 14:04:15  raasch
    52 ! droplet growth by collision completed (still has to be optimized for speed),
    53 ! new routine collision_efficiency (see end of this file),
    54 ! 2*r replaced by r in the exponential terms (momentum equation for particles),
    55 ! number of particles really used is additionally output on the netcdf particle
    56 ! data file,
    57 ! test print statement removed
    58 !
    59 ! Revision 1.26  2005/06/26 19:50:14  raasch
    60 ! Cloud droplet features added, radius is used instead of diameter,
    61 ! particle type initial data assignment extended,
    62 ! output on unit 80 for non-parallel runs
    63 !
    64 ! Revision 1.25  2005/05/18 14:53:11  raasch
    65 ! Extensions for NetCDF output
    66 !
    67 ! Revision 1.24  2005/04/23 08:24:58  raasch
    68 ! Revised calculation of output time counters regarding a possible decrease of
    69 ! the output time interval in case of restart runs
    70 !
    71 ! Revision 1.23  2005/03/26 14:20:05  raasch
    72 ! calculation of vertical particle velocity (with inertia) corrected, exp_arg
    73 ! had a wrong sign
    74 !
    75 ! Revision 1.22  2004/04/30 07:48:57  raasch
    76 ! Particle type initial data assignment extended
    77 !
    78 ! Revision 1.21  2004/01/28 15:00:48  raasch
    79 ! Output of particle infos in subroutine allocate_prt_memory on demand only
    80 !
    81 ! Revision 1.20  2003/10/29 08:38:52  raasch
    82 ! Module random_function_mod is used, modifications for new particle group
    83 ! feature, version number is written on binary file
    84 !
    85 ! Revision 1.19  2003/03/16 09:19:32  raasch
    86 ! Two underscores (_) are placed in front of all define-strings
    87 !
    88 ! Revision 1.18  2003/03/14 13:36:31  raasch
    89 ! Error in reflection boundary condition removed: particle speed must also
    90 ! be inverted
    91 !
    92 ! Revision 1.17  2003/03/04 11:23:20  raasch
    93 ! Error in particle inertia part removed (exp_arg must not contain the timestep)
    94 ! Particle velocities are also stored in particles in case of zero density ratio
    95 !
    96 ! Revision 1.16  2002/09/12 12:55:29  raasch
    97 ! Transport of particles with inertia implemented, write density_ratio to
    98 ! restart file
    99 !
    100 ! Revision 1.15  2002/06/11 12:11:56  raasch
    101 ! Declaration of u_int_l, u_int_u, v_int_l, v_int_u added
    102 !
    103 ! Revision 1.14  2002/05/02  18:46:34  18:46:34  raasch (Siegfried Raasch)
    104 ! Horizontal velocity components for particle advection are now interpolated
    105 ! between horizontal grid levels
    106 !
    107 ! Revision 1.13  2002/04/24  18:46:11  18:46:11  raasch (Siegfried Raasch)
    108 ! Temporary particle-arrays trlp, trnp, trrp and trsp are now allocatable arrays
    109 !
    110 ! Revision 1.12  2002/04/16 07:44:07  raasch
    111 ! Additional boundary conditions added, particle data for later analysis can be
    112 ! written on file
    113 !
    114 ! Revision 1.11  2001/11/09 13:38:20  raasch
    115 ! Colour information for particle tails is stored in array
    116 ! particle_tail_coordinates, adding information to the particle tails
    117 ! moved after call of user_particle_attributes
    118 !
    119 ! Revision 1.10  2001/08/21 08:19:48  raasch
    120 ! Storage of particle tails, precision of vertical advection improved,
    121 ! particle attributes (e.g. size and color) can be defined by the user
    122 !
    123 ! Revision 1.9  2001/07/12 12:01:27  raasch
    124 ! Random start positions of initial particles possible
    125 !
    126 ! Revision 1.8  2001/03/29 17:27:10  raasch
    127 ! Translation of remaining German identifiers (variables, subroutines, etc.)
    128 !
    129 ! Revision 1.7  2001/02/13 10:24:00  raasch
    130 ! Particle advection possible in case of galilei transformation
    131 !
    132 ! Revision 1.6  2001/01/25 06:51:04  raasch
    133 ! Module test_variables removed, writing of particle informations is optional
    134 ! now
    135 !
    136 ! Revision 1.5  2001/01/02 17:15:29  raasch
    137 ! Unit 80 is immediately closed after usage every time. Particle positions
    138 ! are written on unit 90 instead of 81 (subroutine write_particles).
    139 !
    140 ! Revision 1.4  2000/12/28 12:43:54  raasch
    141 ! Packing of particles-array is done explicitly (not by PACK-function, which
    142 ! seems to have problems when a large number of PEs is used)
    143 ! Missing translations included,
    144 ! code is used only optionally (cpp-directives are added)
    145 !
    146 ! Revision 1.3  2000/01/20 08:51:37  letzel
    147 ! All comments translated into English
    148 !
    149 ! Revision 1.2  1999/11/25 16:41:48  raasch
    150 ! Indexfehler im nichtparallelen Teil korrigiert
    15142!
    15243! Revision 1.1  1999/11/25 16:16:06  raasch
     
    693584          DO  j = nys, nyn
    694585             DO  k = nzb, nzt+1
    695                 de_dx(k,j,i) = sgs_wfu_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
    696                 de_dy(k,j,i) = sgs_wfv_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
     586
     587                IF ( k <= nzb_s_inner(j,i-1)  .AND.  &
     588                     k  > nzb_s_inner(j,i)    .AND.  &
     589                     k  > nzb_s_inner(j,i+1) )  THEN
     590                   de_dx(k,j,i) = 2.0 * sgs_wfu_part * &
     591                                  ( e(k,j,i+1) - e(k,j,i) ) * ddx
     592                ELSEIF ( k  > nzb_s_inner(j,i-1)  .AND.  &
     593                         k  > nzb_s_inner(j,i)    .AND.  &
     594                         k <= nzb_s_inner(j,i+1) )  THEN
     595                   de_dx(k,j,i) = 2.0 * sgs_wfu_part * &
     596                                  ( e(k,j,i) - e(k,j,i-1) ) * ddx
     597                ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j,i+1) ) &
     598                THEN
     599                   de_dx(k,j,i) = 0.0
     600                ELSEIF ( k < nzb_s_inner(j,i-1)  .AND.  k < nzb_s_inner(j,i) ) &
     601                THEN
     602                   de_dx(k,j,i) = 0.0
     603                ELSE
     604                   de_dx(k,j,i) = sgs_wfu_part * &
     605                                  ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
     606                ENDIF
     607
     608                IF ( k <= nzb_s_inner(j-1,i)  .AND.  &
     609                     k  > nzb_s_inner(j,i)    .AND.  &
     610                     k  > nzb_s_inner(j+1,i) )  THEN
     611                   de_dy(k,j,i) = 2.0 * sgs_wfv_part * &
     612                                  ( e(k,j+1,i) - e(k,j,i) ) * ddy
     613                ELSEIF ( k  > nzb_s_inner(j-1,i)  .AND.  &
     614                         k  > nzb_s_inner(j,i)    .AND.  &
     615                         k <= nzb_s_inner(j+1,i) )  THEN
     616                   de_dy(k,j,i) = 2.0 * sgs_wfv_part * &
     617                                  ( e(k,j,i) - e(k,j-1,i) ) * ddy
     618                ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j+1,i) ) &
     619                THEN
     620                   de_dy(k,j,i) = 0.0
     621                ELSEIF ( k < nzb_s_inner(j-1,i)  .AND.  k < nzb_s_inner(j,i) ) &
     622                THEN
     623                   de_dy(k,j,i) = 0.0
     624                ELSE
     625                   de_dy(k,j,i) = sgs_wfv_part * &
     626                                  ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
     627                ENDIF
     628
    697629             ENDDO
    698630          ENDDO
     
    703635       DO  i = nxl, nxr
    704636          DO  j = nys, nyn
    705              DO  k = nzb+2, nzt-1
     637
     638             DO  k = nzb_s_inner(j,i)+2, nzt-1
    706639                de_dz(k,j,i) = 2.0 * sgs_wfw_part * &
    707640                               ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) )
    708641             ENDDO
    709              de_dz(nzb,j,i)   = 0.0
    710              de_dz(nzb+1,j,i) = 2.0 * sgs_wfw_part *              &
    711                                 ( e(nzb+2,j,i) - e(nzb+1,j,i) ) / &
    712                                 ( zu(nzb+2) - zu(nzb+1) )
     642
     643             k = nzb_s_inner(j,i)
     644             de_dz(nzb:k,j,i)   = 0.0
     645             de_dz(k+1,j,i) = 2.0 * sgs_wfw_part * ( e(k+2,j,i) - e(k+1,j,i) ) &
     646                                                 / ( zu(k+2) - zu(k+1) )
    713647             de_dz(nzt,j,i)   = 0.0
    714648             de_dz(nzt+1,j,i) = 0.0
    715649          ENDDO
    716        ENDDO             
     650       ENDDO
    717651
    718652!
     
    967901             k = ( particles(n)%z + 0.5 * dz ) / dz  ! only exact if eq.dist
    968902
    969              x  = particles(n)%x - i * dx
    970              y  = particles(n)%y - j * dy
    971              aa = x**2          + y**2
    972              bb = ( dx - x )**2 + y**2
    973              cc = x**2          + ( dy - y )**2
    974              dd = ( dx - x )**2 + ( dy - y )**2
    975              gg = aa + bb + cc + dd
    976 
    977              e_int_l = ( ( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1)   &
    978                        + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
    979                        ) / ( 3.0 * gg )
    980 
    981              IF ( k+1 == nzt+1 )  THEN
    982                 e_int = e_int_l
     903             IF ( topography == 'flat' )  THEN       
     904
     905                x  = particles(n)%x - i * dx
     906                y  = particles(n)%y - j * dy
     907                aa = x**2          + y**2
     908                bb = ( dx - x )**2 + y**2
     909                cc = x**2          + ( dy - y )**2
     910                dd = ( dx - x )**2 + ( dy - y )**2
     911                gg = aa + bb + cc + dd
     912
     913                e_int_l = ( ( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1)   &
     914                          + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1) &
     915                          ) / ( 3.0 * gg )
     916
     917                IF ( k+1 == nzt+1 )  THEN
     918                   e_int = e_int_l
     919                ELSE
     920                   e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
     921                               ( gg - bb ) * e(k+1,j,i+1) + &
     922                               ( gg - cc ) * e(k+1,j+1,i) + &
     923                               ( gg - dd ) * e(k+1,j+1,i+1) &
     924                             ) / ( 3.0 * gg )
     925                   e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
     926                                     ( e_int_u - e_int_l )
     927                ENDIF
     928
     929!
     930!--             Interpolate the TKE gradient along x (adopt incides i,j,k and
     931!--             all position variables from above (TKE))
     932                de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
     933                                ( gg - bb ) * de_dx(k,j,i+1) + &
     934                                ( gg - cc ) * de_dx(k,j+1,i) + &
     935                                ( gg - dd ) * de_dx(k,j+1,i+1) &
     936                              ) / ( 3.0 * gg )
     937
     938                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
     939                   de_dx_int = de_dx_int_l
     940                ELSE
     941                   de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
     942                                   ( gg - bb ) * de_dx(k+1,j,i+1) + &
     943                                   ( gg - cc ) * de_dx(k+1,j+1,i) + &
     944                                   ( gg - dd ) * de_dx(k+1,j+1,i+1) &
     945                                 ) / ( 3.0 * gg )
     946                   de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * &
     947                                             ( de_dx_int_u - de_dx_int_l )
     948                ENDIF
     949
     950!
     951!--             Interpolate the TKE gradient along y
     952                de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
     953                                ( gg - bb ) * de_dy(k,j,i+1) + &
     954                                ( gg - cc ) * de_dy(k,j+1,i) + &
     955                                ( gg - dd ) * de_dy(k,j+1,i+1) &
     956                              ) / ( 3.0 * gg )
     957                IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
     958                   de_dy_int = de_dy_int_l
     959                ELSE
     960                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
     961                                   ( gg - bb ) * de_dy(k+1,j,i+1) + &
     962                                   ( gg - cc ) * de_dy(k+1,j+1,i) + &
     963                                   ( gg - dd ) * de_dy(k+1,j+1,i+1) &
     964                                 ) / ( 3.0 * gg )
     965                   de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * &
     966                                             ( de_dy_int_u - de_dy_int_l )
     967                ENDIF
     968
     969!
     970!--             Interpolate the TKE gradient along z
     971                IF ( particles(n)%z < 0.5 * dz ) THEN
     972                   de_dz_int = 0.0
     973                ELSE
     974                   de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
     975                                   ( gg - bb ) * de_dz(k,j,i+1) + &
     976                                   ( gg - cc ) * de_dz(k,j+1,i) + &
     977                                   ( gg - dd ) * de_dz(k,j+1,i+1) &
     978                                 ) / ( 3.0 * gg )
     979
     980                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
     981                      de_dz_int = de_dz_int_l
     982                   ELSE
     983                      de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
     984                                      ( gg - bb ) * de_dz(k+1,j,i+1) + &
     985                                      ( gg - cc ) * de_dz(k+1,j+1,i) + &
     986                                      ( gg - dd ) * de_dz(k+1,j+1,i+1) &
     987                                    ) / ( 3.0 * gg )
     988                      de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * &
     989                                                ( de_dz_int_u - de_dz_int_l )
     990                   ENDIF
     991                ENDIF
     992
     993!
     994!--             Interpolate the dissipation of TKE
     995                diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
     996                               ( gg - bb ) * diss(k,j,i+1) + &
     997                               ( gg - cc ) * diss(k,j+1,i) + &
     998                               ( gg - dd ) * diss(k,j+1,i+1) &
     999                              ) / ( 3.0 * gg )
     1000
     1001                IF ( k+1 == nzt+1 )  THEN
     1002                   diss_int = diss_int_l
     1003                ELSE
     1004                   diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
     1005                                  ( gg - bb ) * diss(k+1,j,i+1) + &
     1006                                  ( gg - cc ) * diss(k+1,j+1,i) + &
     1007                                  ( gg - dd ) * diss(k+1,j+1,i+1) &
     1008                                 ) / ( 3.0 * gg )
     1009                   diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * &
     1010                                           ( diss_int_u - diss_int_l )
     1011                ENDIF
     1012
    9831013             ELSE
    984                 e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
    985                             ( gg - bb ) * e(k+1,j,i+1) + &
    986                             ( gg - cc ) * e(k+1,j+1,i) + &
    987                             ( gg - dd ) * e(k+1,j+1,i+1) &
    988                           ) / ( 3.0 * gg )
    989                 e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
    990                                   ( e_int_u - e_int_l )
    991              ENDIF
    992 
    993 !
    994 !--          Interpolate the TKE gradient along x (adopt incides i,j,k and all
    995 !--          position variables from above (TKE))
    996              de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
    997                              ( gg - bb ) * de_dx(k,j,i+1) + &
    998                              ( gg - cc ) * de_dx(k,j+1,i) + &
    999                              ( gg - dd ) * de_dx(k,j+1,i+1) &
    1000                            ) / ( 3.0 * gg )
    1001 
    1002              IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    1003                 de_dx_int = de_dx_int_l
    1004              ELSE
    1005                 de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
    1006                                 ( gg - bb ) * de_dx(k+1,j,i+1) + &
    1007                                 ( gg - cc ) * de_dx(k+1,j+1,i) + &
    1008                                 ( gg - dd ) * de_dx(k+1,j+1,i+1) &
    1009                               ) / ( 3.0 * gg )
    1010                 de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / dz * &
    1011                                           ( de_dx_int_u - de_dx_int_l )
    1012              ENDIF
    1013 
    1014 !
    1015 !--          Interpolate the TKE gradient along y
    1016              de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
    1017                              ( gg - bb ) * de_dy(k,j,i+1) + &
    1018                              ( gg - cc ) * de_dy(k,j+1,i) + &
    1019                              ( gg - dd ) * de_dy(k,j+1,i+1) &
    1020                            ) / ( 3.0 * gg )
    1021              IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    1022                 de_dy_int = de_dy_int_l
    1023              ELSE
    1024                 de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
    1025                                 ( gg - bb ) * de_dy(k+1,j,i+1) + &
    1026                                 ( gg - cc ) * de_dy(k+1,j+1,i) + &
    1027                                 ( gg - dd ) * de_dy(k+1,j+1,i+1) &
    1028                               ) / ( 3.0 * gg )
    1029                 de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / dz * &
    1030                                           ( de_dy_int_u - de_dy_int_l )
    1031              ENDIF
    1032 
    1033 !
    1034 !--          Interpolate the TKE gradient along z
    1035              de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
    1036                              ( gg - bb ) * de_dz(k,j,i+1) + &
    1037                              ( gg - cc ) * de_dz(k,j+1,i) + &
    1038                              ( gg - dd ) * de_dz(k,j+1,i+1) &
    1039                            ) / ( 3.0 * gg )
    1040 
    1041              IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    1042                 de_dz_int = de_dz_int_l
    1043              ELSE
    1044                 de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
    1045                                 ( gg - bb ) * de_dz(k+1,j,i+1) + &
    1046                                 ( gg - cc ) * de_dz(k+1,j+1,i) + &
    1047                                 ( gg - dd ) * de_dz(k+1,j+1,i+1) &
    1048                               ) / ( 3.0 * gg )
    1049                 de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) / dz * &
    1050                                           ( de_dz_int_u - de_dz_int_l )
    1051              ENDIF
    1052 
    1053 !
    1054 !--          Interpolate the dissipation of TKE
    1055              diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
    1056                             ( gg - bb ) * diss(k,j,i+1) + &
    1057                             ( gg - cc ) * diss(k,j+1,i) + &
    1058                             ( gg - dd ) * diss(k,j+1,i+1) &
    1059                            ) / ( 3.0 * gg )
    1060 
    1061              IF ( k+1 == nzt+1 )  THEN
    1062                 diss_int = diss_int_l
    1063              ELSE
    1064                 diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
    1065                                ( gg - bb ) * diss(k+1,j,i+1) + &
    1066                                ( gg - cc ) * diss(k+1,j+1,i) + &
    1067                                ( gg - dd ) * diss(k+1,j+1,i+1) &
    1068                               ) / ( 3.0 * gg )
    1069                 diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz * &
    1070                                         ( diss_int_u - diss_int_l )
    1071              ENDIF
     1014
     1015!
     1016!--             In case that there are buildings it has to be determined
     1017!--             how many of the gridpoints defining the particle box are
     1018!--             situated within a building
     1019!--             gp_outside_of_building(1): i,j,k
     1020!--             gp_outside_of_building(2): i,j+1,k
     1021!--             gp_outside_of_building(3): i,j,k+1
     1022!--             gp_outside_of_building(4): i,j+1,k+1
     1023!--             gp_outside_of_building(5): i+1,j,k
     1024!--             gp_outside_of_building(6): i+1,j+1,k
     1025!--             gp_outside_of_building(7): i+1,j,k+1
     1026!--             gp_outside_of_building(8): i+1,j+1,k+1
     1027           
     1028                gp_outside_of_building = 0
     1029                location = 0.0
     1030                num_gp = 0
     1031
     1032                IF ( k > nzb_s_inner(j,i)  .OR.  nzb_s_inner(j,i) == 0 )  THEN
     1033                   num_gp = num_gp + 1
     1034                   gp_outside_of_building(1) = 1
     1035                   location(num_gp,1) = i * dx
     1036                   location(num_gp,2) = j * dx
     1037                   location(num_gp,3) = k * dz - 0.5 * dz
     1038                   ei(num_gp)     = e(k,j,i)
     1039                   dissi(num_gp)  = diss(k,j,i)
     1040                   de_dxi(num_gp) = de_dx(k,j,i)
     1041                   de_dyi(num_gp) = de_dy(k,j,i)
     1042                   de_dzi(num_gp) = de_dz(k,j,i)
     1043                ENDIF
     1044
     1045                IF ( k > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
     1046                THEN
     1047                   num_gp = num_gp + 1
     1048                   gp_outside_of_building(2) = 1
     1049                   location(num_gp,1) = i * dx
     1050                   location(num_gp,2) = (j+1) * dx
     1051                   location(num_gp,3) = k * dz - 0.5 * dz
     1052                   ei(num_gp)     = e(k,j+1,i)
     1053                   dissi(num_gp)  = diss(k,j+1,i)       
     1054                   de_dxi(num_gp) = de_dx(k,j+1,i)
     1055                   de_dyi(num_gp) = de_dy(k,j+1,i)
     1056                   de_dzi(num_gp) = de_dz(k,j+1,i)
     1057                ENDIF
     1058
     1059                IF ( k+1 > nzb_s_inner(j,i)  .OR.  nzb_s_inner(j,i) == 0 )  THEN
     1060                   num_gp = num_gp + 1
     1061                   gp_outside_of_building(3) = 1
     1062                   location(num_gp,1) = i * dx
     1063                   location(num_gp,2) = j * dy
     1064                   location(num_gp,3) = (k+1) * dz - 0.5 * dz
     1065                   ei(num_gp)     = e(k+1,j,i)
     1066                   dissi(num_gp)  = diss(k+1,j,i)       
     1067                   de_dxi(num_gp) = de_dx(k+1,j,i)
     1068                   de_dyi(num_gp) = de_dy(k+1,j,i)
     1069                   de_dzi(num_gp) = de_dz(k+1,j,i)
     1070                ENDIF
     1071
     1072                IF ( k+1 > nzb_s_inner(j+1,i)  .OR.  nzb_s_inner(j+1,i) == 0 ) &
     1073                THEN
     1074                   num_gp = num_gp + 1
     1075                   gp_outside_of_building(4) = 1
     1076                   location(num_gp,1) = i * dx
     1077                   location(num_gp,2) = (j+1) * dy
     1078                   location(num_gp,3) = (k+1) * dz - 0.5 * dz
     1079                   ei(num_gp)     = e(k+1,j+1,i)
     1080                   dissi(num_gp)  = diss(k+1,j+1,i)       
     1081                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
     1082                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
     1083                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
     1084                ENDIF
     1085 
     1086                IF ( k > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
     1087                THEN
     1088                   num_gp = num_gp + 1
     1089                   gp_outside_of_building(5) = 1
     1090                   location(num_gp,1) = (i+1) * dx
     1091                   location(num_gp,2) = j * dy
     1092                   location(num_gp,3) = k * dz - 0.5 * dz
     1093                   ei(num_gp)     = e(k,j,i+1)
     1094                   dissi(num_gp)  = diss(k,j,i+1)
     1095                   de_dxi(num_gp) = de_dx(k,j,i+1)
     1096                   de_dyi(num_gp) = de_dy(k,j,i+1)
     1097                   de_dzi(num_gp) = de_dz(k,j,i+1)
     1098                ENDIF
     1099
     1100                IF ( k > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0 ) &
     1101                THEN
     1102                   num_gp = num_gp + 1
     1103                   gp_outside_of_building(6) = 1
     1104                   location(num_gp,1) = (i+1) * dx
     1105                   location(num_gp,2) = (j+1) * dy
     1106                   location(num_gp,3) = k * dz - 0.5 * dz
     1107                   ei(num_gp)     = e(k,j+1,i+1)
     1108                   dissi(num_gp)  = diss(k,j+1,i+1)
     1109                   de_dxi(num_gp) = de_dx(k,j+1,i+1)
     1110                   de_dyi(num_gp) = de_dy(k,j+1,i+1)
     1111                   de_dzi(num_gp) = de_dz(k,j+1,i+1)
     1112                ENDIF
     1113
     1114                IF ( k+1 > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
     1115                THEN
     1116                   num_gp = num_gp + 1
     1117                   gp_outside_of_building(7) = 1
     1118                   location(num_gp,1) = (i+1) * dx
     1119                   location(num_gp,2) = j * dy
     1120                   location(num_gp,3) = (k+1) * dz - 0.5 * dz
     1121                   ei(num_gp)     = e(k+1,j,i+1)
     1122                   dissi(num_gp)  = diss(k+1,j,i+1)
     1123                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
     1124                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
     1125                   de_dzi(num_gp) = de_dz(k+1,j,i+1)
     1126                ENDIF
     1127
     1128                IF ( k+1 > nzb_s_inner(j+1,i+1) .OR. nzb_s_inner(j+1,i+1) == 0)&
     1129                THEN
     1130                   num_gp = num_gp + 1
     1131                   gp_outside_of_building(8) = 1
     1132                   location(num_gp,1) = (i+1) * dx
     1133                   location(num_gp,2) = (j+1) * dy
     1134                   location(num_gp,3) = (k+1) * dz - 0.5 * dz
     1135                   ei(num_gp)     = e(k+1,j+1,i+1)
     1136                   dissi(num_gp)  = diss(k+1,j+1,i+1)
     1137                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
     1138                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
     1139                   de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
     1140                ENDIF
     1141
     1142!
     1143!--             If all gridpoints are situated outside of a building, then the
     1144!--             ordinary interpolation scheme can be used.
     1145                IF ( num_gp == 8 )  THEN
     1146
     1147                   x  = particles(n)%x - i * dx
     1148                   y  = particles(n)%y - j * dy
     1149                   aa = x**2          + y**2
     1150                   bb = ( dx - x )**2 + y**2
     1151                   cc = x**2          + ( dy - y )**2
     1152                   dd = ( dx - x )**2 + ( dy - y )**2
     1153                   gg = aa + bb + cc + dd
     1154
     1155                   e_int_l = (( gg-aa ) * e(k,j,i)   + ( gg-bb ) * e(k,j,i+1) &
     1156                            + ( gg-cc ) * e(k,j+1,i) + ( gg-dd ) * e(k,j+1,i+1)&
     1157                             ) / ( 3.0 * gg )
     1158
     1159                   IF ( k+1 == nzt+1 )  THEN
     1160                      e_int = e_int_l
     1161                   ELSE
     1162                      e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
     1163                                  ( gg - bb ) * e(k+1,j,i+1) + &
     1164                                  ( gg - cc ) * e(k+1,j+1,i) + &
     1165                                  ( gg - dd ) * e(k+1,j+1,i+1) &
     1166                                ) / ( 3.0 * gg )
     1167                      e_int = e_int_l + ( particles(n)%z - zu(k) ) / dz * &
     1168                                        ( e_int_u - e_int_l )
     1169                   ENDIF
     1170
     1171!
     1172!--                Interpolate the TKE gradient along x (adopt incides i,j,k
     1173!--                and all position variables from above (TKE))
     1174                   de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
     1175                                   ( gg - bb ) * de_dx(k,j,i+1) + &
     1176                                   ( gg - cc ) * de_dx(k,j+1,i) + &
     1177                                   ( gg - dd ) * de_dx(k,j+1,i+1) &
     1178                                 ) / ( 3.0 * gg )
     1179
     1180                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
     1181                      de_dx_int = de_dx_int_l
     1182                   ELSE
     1183                      de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
     1184                                      ( gg - bb ) * de_dx(k+1,j,i+1) + &
     1185                                      ( gg - cc ) * de_dx(k+1,j+1,i) + &
     1186                                      ( gg - dd ) * de_dx(k+1,j+1,i+1) &
     1187                                    ) / ( 3.0 * gg )
     1188                      de_dx_int = de_dx_int_l + ( particles(n)%z - zu(k) ) / &
     1189                                              dz * ( de_dx_int_u - de_dx_int_l )
     1190                   ENDIF
     1191
     1192!
     1193!--                Interpolate the TKE gradient along y
     1194                   de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
     1195                                   ( gg - bb ) * de_dy(k,j,i+1) + &
     1196                                   ( gg - cc ) * de_dy(k,j+1,i) + &
     1197                                   ( gg - dd ) * de_dy(k,j+1,i+1) &
     1198                                 ) / ( 3.0 * gg )
     1199                   IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
     1200                      de_dy_int = de_dy_int_l
     1201                   ELSE
     1202                      de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
     1203                                      ( gg - bb ) * de_dy(k+1,j,i+1) + &
     1204                                      ( gg - cc ) * de_dy(k+1,j+1,i) + &
     1205                                      ( gg - dd ) * de_dy(k+1,j+1,i+1) &
     1206                                    ) / ( 3.0 * gg )
     1207                      de_dy_int = de_dy_int_l + ( particles(n)%z - zu(k) ) / &
     1208                                              dz * ( de_dy_int_u - de_dy_int_l )
     1209                   ENDIF
     1210
     1211!
     1212!--                Interpolate the TKE gradient along z
     1213                   IF ( particles(n)%z < 0.5 * dz )  THEN
     1214                      de_dz_int = 0.0
     1215                   ELSE
     1216                      de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
     1217                                      ( gg - bb ) * de_dz(k,j,i+1) + &
     1218                                      ( gg - cc ) * de_dz(k,j+1,i) + &
     1219                                      ( gg - dd ) * de_dz(k,j+1,i+1) &
     1220                                    ) / ( 3.0 * gg )
     1221
     1222                      IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
     1223                         de_dz_int = de_dz_int_l
     1224                      ELSE
     1225                         de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
     1226                                         ( gg - bb ) * de_dz(k+1,j,i+1) + &
     1227                                         ( gg - cc ) * de_dz(k+1,j+1,i) + &
     1228                                         ( gg - dd ) * de_dz(k+1,j+1,i+1) &
     1229                                       ) / ( 3.0 * gg )
     1230                         de_dz_int = de_dz_int_l + ( particles(n)%z - zu(k) ) /&
     1231                                              dz * ( de_dz_int_u - de_dz_int_l )
     1232                      ENDIF
     1233                   ENDIF
     1234
     1235!
     1236!--                Interpolate the dissipation of TKE
     1237                   diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
     1238                                  ( gg - bb ) * diss(k,j,i+1) + &
     1239                                  ( gg - cc ) * diss(k,j+1,i) + &
     1240                                  ( gg - dd ) * diss(k,j+1,i+1) &
     1241                                ) / ( 3.0 * gg )
     1242
     1243                   IF ( k+1 == nzt+1 )  THEN
     1244                      diss_int = diss_int_l
     1245                   ELSE
     1246                      diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
     1247                                     ( gg - bb ) * diss(k+1,j,i+1) + &
     1248                                     ( gg - cc ) * diss(k+1,j+1,i) + &
     1249                                     ( gg - dd ) * diss(k+1,j+1,i+1) &
     1250                                   ) / ( 3.0 * gg )
     1251                      diss_int = diss_int_l + ( particles(n)%z - zu(k) ) / dz *&
     1252                                              ( diss_int_u - diss_int_l )
     1253                   ENDIF
     1254
     1255                ELSE
     1256
     1257!
     1258!--                If wall between gridpoint 1 and gridpoint 5, then
     1259!--                Neumann boundary condition has to be applied
     1260                   IF ( gp_outside_of_building(1) == 1  .AND. &
     1261                        gp_outside_of_building(5) == 0 )  THEN
     1262                      num_gp = num_gp + 1
     1263                      location(num_gp,1) = i * dx + 0.5 * dx
     1264                      location(num_gp,2) = j * dy
     1265                      location(num_gp,3) = k * dz - 0.5 * dz
     1266                      ei(num_gp)     = e(k,j,i)
     1267                      dissi(num_gp)  = diss(k,j,i)
     1268                      de_dxi(num_gp) = 0.0
     1269                      de_dyi(num_gp) = de_dy(k,j,i)
     1270                      de_dzi(num_gp) = de_dz(k,j,i)                     
     1271                   ENDIF
     1272   
     1273                   IF ( gp_outside_of_building(5) == 1  .AND. &
     1274                      gp_outside_of_building(1) == 0 )  THEN
     1275                      num_gp = num_gp + 1
     1276                      location(num_gp,1) = i * dx + 0.5 * dx
     1277                      location(num_gp,2) = j * dy
     1278                      location(num_gp,3) = k * dz - 0.5 * dz
     1279                      ei(num_gp)     = e(k,j,i+1)
     1280                      dissi(num_gp)  = diss(k,j,i+1)
     1281                      de_dxi(num_gp) = 0.0
     1282                      de_dyi(num_gp) = de_dy(k,j,i+1)
     1283                      de_dzi(num_gp) = de_dz(k,j,i+1)
     1284                   ENDIF
     1285
     1286!
     1287!--                If wall between gridpoint 5 and gridpoint 6, then
     1288!--                then Neumann boundary condition has to be applied
     1289                   IF ( gp_outside_of_building(5) == 1  .AND. &
     1290                        gp_outside_of_building(6) == 0 )  THEN
     1291                      num_gp = num_gp + 1
     1292                      location(num_gp,1) = (i+1) * dx
     1293                      location(num_gp,2) = j * dy + 0.5 * dy
     1294                      location(num_gp,3) = k * dz - 0.5 * dz
     1295                      ei(num_gp)     = e(k,j,i+1)
     1296                      dissi(num_gp)  = diss(k,j,i+1)
     1297                      de_dxi(num_gp) = de_dx(k,j,i+1)
     1298                      de_dyi(num_gp) = 0.0
     1299                      de_dzi(num_gp) = de_dz(k,j,i+1)
     1300                   ENDIF
     1301
     1302                   IF ( gp_outside_of_building(6) == 1  .AND. &
     1303                        gp_outside_of_building(5) == 0 )  THEN
     1304                      num_gp = num_gp + 1
     1305                      location(num_gp,1) = (i+1) * dx
     1306                      location(num_gp,2) = j * dy + 0.5 * dy
     1307                      location(num_gp,3) = k * dz - 0.5 * dz
     1308                      ei(num_gp)     = e(k,j+1,i+1)
     1309                      dissi(num_gp)  = diss(k,j+1,i+1)
     1310                      de_dxi(num_gp) = de_dx(k,j+1,i+1)
     1311                      de_dyi(num_gp) = 0.0
     1312                      de_dzi(num_gp) = de_dz(k,j+1,i+1)
     1313                   ENDIF
     1314
     1315!
     1316!--                If wall between gridpoint 2 and gridpoint 6, then
     1317!--                Neumann boundary condition has to be applied
     1318                   IF ( gp_outside_of_building(2) == 1  .AND. &
     1319                        gp_outside_of_building(6) == 0 )  THEN
     1320                      num_gp = num_gp + 1
     1321                      location(num_gp,1) = i * dx + 0.5 * dx
     1322                      location(num_gp,2) = (j+1) * dy
     1323                      location(num_gp,3) = k * dz - 0.5 * dz
     1324                      ei(num_gp)     = e(k,j+1,i)
     1325                      dissi(num_gp)  = diss(k,j+1,i)
     1326                      de_dxi(num_gp) = 0.0
     1327                      de_dyi(num_gp) = de_dy(k,j+1,i)
     1328                      de_dzi(num_gp) = de_dz(k,j+1,i)                     
     1329                   ENDIF
     1330   
     1331                   IF ( gp_outside_of_building(6) == 1  .AND. &
     1332                        gp_outside_of_building(2) == 0 )  THEN
     1333                      num_gp = num_gp + 1
     1334                      location(num_gp,1) = i * dx + 0.5 * dx
     1335                      location(num_gp,2) = (j+1) * dy
     1336                      location(num_gp,3) = k * dz - 0.5 * dz
     1337                      ei(num_gp)     = e(k,j+1,i+1)
     1338                      dissi(num_gp)  = diss(k,j+1,i+1)
     1339                      de_dxi(num_gp) = 0.0
     1340                      de_dyi(num_gp) = de_dy(k,j+1,i+1)
     1341                      de_dzi(num_gp) = de_dz(k,j+1,i+1)
     1342                   ENDIF
     1343
     1344!
     1345!--                If wall between gridpoint 1 and gridpoint 2, then
     1346!--                Neumann boundary condition has to be applied
     1347                   IF ( gp_outside_of_building(1) == 1  .AND. &
     1348                        gp_outside_of_building(2) == 0 )  THEN
     1349                      num_gp = num_gp + 1
     1350                      location(num_gp,1) = i * dx
     1351                      location(num_gp,2) = j * dy + 0.5 * dy
     1352                      location(num_gp,3) = k * dz - 0.5 * dz
     1353                      ei(num_gp)     = e(k,j,i)
     1354                      dissi(num_gp)  = diss(k,j,i)
     1355                      de_dxi(num_gp) = de_dx(k,j,i)
     1356                      de_dyi(num_gp) = 0.0
     1357                      de_dzi(num_gp) = de_dz(k,j,i)
     1358                   ENDIF
     1359
     1360                   IF ( gp_outside_of_building(2) == 1  .AND. &
     1361                        gp_outside_of_building(1) == 0 )  THEN
     1362                      num_gp = num_gp + 1
     1363                      location(num_gp,1) = i * dx
     1364                      location(num_gp,2) = j * dy + 0.5 * dy
     1365                      location(num_gp,3) = k * dz - 0.5 * dz
     1366                      ei(num_gp)     = e(k,j+1,i)
     1367                      dissi(num_gp)  = diss(k,j+1,i)
     1368                      de_dxi(num_gp) = de_dx(k,j+1,i)
     1369                      de_dyi(num_gp) = 0.0
     1370                      de_dzi(num_gp) = de_dz(k,j+1,i)
     1371                   ENDIF
     1372
     1373!
     1374!--                If wall between gridpoint 3 and gridpoint 7, then
     1375!--                Neumann boundary condition has to be applied
     1376                   IF ( gp_outside_of_building(3) == 1  .AND. &
     1377                        gp_outside_of_building(7) == 0 )  THEN
     1378                      num_gp = num_gp + 1
     1379                      location(num_gp,1) = i * dx + 0.5 * dx
     1380                      location(num_gp,2) = j * dy
     1381                      location(num_gp,3) = k * dz + 0.5 * dz
     1382                      ei(num_gp)     = e(k+1,j,i)
     1383                      dissi(num_gp)  = diss(k+1,j,i)
     1384                      de_dxi(num_gp) = 0.0
     1385                      de_dyi(num_gp) = de_dy(k+1,j,i)
     1386                      de_dzi(num_gp) = de_dz(k+1,j,i)                     
     1387                   ENDIF
     1388   
     1389                   IF ( gp_outside_of_building(7) == 1  .AND. &
     1390                        gp_outside_of_building(3) == 0 )  THEN
     1391                      num_gp = num_gp + 1
     1392                      location(num_gp,1) = i * dx + 0.5 * dx
     1393                      location(num_gp,2) = j * dy
     1394                      location(num_gp,3) = k * dz + 0.5 * dz
     1395                      ei(num_gp)     = e(k+1,j,i+1)
     1396                      dissi(num_gp)  = diss(k+1,j,i+1)
     1397                      de_dxi(num_gp) = 0.0
     1398                      de_dyi(num_gp) = de_dy(k+1,j,i+1)
     1399                      de_dzi(num_gp) = de_dz(k+1,j,i+1)
     1400                   ENDIF
     1401
     1402!
     1403!--                If wall between gridpoint 7 and gridpoint 8, then
     1404!--                Neumann boundary condition has to be applied
     1405                   IF ( gp_outside_of_building(7) == 1  .AND. &
     1406                        gp_outside_of_building(8) == 0 )  THEN
     1407                      num_gp = num_gp + 1
     1408                      location(num_gp,1) = (i+1) * dx
     1409                      location(num_gp,2) = j * dy + 0.5 * dy
     1410                      location(num_gp,3) = k * dz + 0.5 * dz
     1411                      ei(num_gp)     = e(k+1,j,i+1)
     1412                      dissi(num_gp)  = diss(k+1,j,i+1)
     1413                      de_dxi(num_gp) = de_dx(k+1,j,i+1)
     1414                      de_dyi(num_gp) = 0.0
     1415                      de_dzi(num_gp) = de_dz(k+1,j,i+1)
     1416                   ENDIF
     1417
     1418                   IF ( gp_outside_of_building(8) == 1  .AND. &
     1419                        gp_outside_of_building(7) == 0 )  THEN
     1420                      num_gp = num_gp + 1
     1421                      location(num_gp,1) = (i+1) * dx
     1422                      location(num_gp,2) = j * dy + 0.5 * dy
     1423                      location(num_gp,3) = k * dz + 0.5 * dz
     1424                      ei(num_gp)     = e(k+1,j+1,i+1)
     1425                      dissi(num_gp)  = diss(k+1,j+1,i+1)
     1426                      de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
     1427                      de_dyi(num_gp) = 0.0
     1428                      de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
     1429                   ENDIF
     1430
     1431!
     1432!--                If wall between gridpoint 4 and gridpoint 8, then
     1433!--                Neumann boundary condition has to be applied
     1434                   IF ( gp_outside_of_building(4) == 1  .AND. &
     1435                        gp_outside_of_building(8) == 0 )  THEN
     1436                      num_gp = num_gp + 1
     1437                      location(num_gp,1) = i * dx + 0.5 * dx
     1438                      location(num_gp,2) = (j+1) * dy
     1439                      location(num_gp,3) = k * dz + 0.5 * dz
     1440                      ei(num_gp)     = e(k+1,j+1,i)
     1441                      dissi(num_gp)  = diss(k+1,j+1,i)
     1442                      de_dxi(num_gp) = 0.0
     1443                      de_dyi(num_gp) = de_dy(k+1,j+1,i)
     1444                      de_dzi(num_gp) = de_dz(k+1,j+1,i)                     
     1445                   ENDIF
     1446   
     1447                   IF ( gp_outside_of_building(8) == 1  .AND. &
     1448                        gp_outside_of_building(4) == 0 )  THEN
     1449                      num_gp = num_gp + 1
     1450                      location(num_gp,1) = i * dx + 0.5 * dx
     1451                      location(num_gp,2) = (j+1) * dy
     1452                      location(num_gp,3) = k * dz + 0.5 * dz
     1453                      ei(num_gp)     = e(k+1,j+1,i+1)
     1454                      dissi(num_gp)  = diss(k+1,j+1,i+1)
     1455                      de_dxi(num_gp) = 0.0
     1456                      de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
     1457                      de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
     1458                   ENDIF
     1459
     1460!
     1461!--                If wall between gridpoint 3 and gridpoint 4, then
     1462!--                Neumann boundary condition has to be applied
     1463                   IF ( gp_outside_of_building(3) == 1  .AND. &
     1464                        gp_outside_of_building(4) == 0 )  THEN
     1465                      num_gp = num_gp + 1
     1466                      location(num_gp,1) = i * dx
     1467                      location(num_gp,2) = j * dy + 0.5 * dy
     1468                      location(num_gp,3) = k * dz + 0.5 * dz
     1469                      ei(num_gp)     = e(k+1,j,i)
     1470                      dissi(num_gp)  = diss(k+1,j,i)
     1471                      de_dxi(num_gp) = de_dx(k+1,j,i)
     1472                      de_dyi(num_gp) = 0.0
     1473                      de_dzi(num_gp) = de_dz(k+1,j,i)
     1474                   ENDIF
     1475
     1476                   IF ( gp_outside_of_building(4) == 1  .AND. &
     1477                        gp_outside_of_building(3) == 0 )  THEN
     1478                      num_gp = num_gp + 1
     1479                      location(num_gp,1) = i * dx
     1480                      location(num_gp,2) = j * dy + 0.5 * dy
     1481                      location(num_gp,3) = k * dz + 0.5 * dz
     1482                      ei(num_gp)     = e(k+1,j+1,i)
     1483                      dissi(num_gp)  = diss(k+1,j+1,i)
     1484                      de_dxi(num_gp) = de_dx(k+1,j+1,i)
     1485                      de_dyi(num_gp) = 0.0
     1486                      de_dzi(num_gp) = de_dz(k+1,j+1,i)
     1487                   ENDIF
     1488
     1489!
     1490!--                If wall between gridpoint 1 and gridpoint 3, then
     1491!--                Neumann boundary condition has to be applied
     1492!--                (only one case as only building beneath is possible)
     1493                   IF ( gp_outside_of_building(1) == 1  .AND. &
     1494                        gp_outside_of_building(3) == 0 )  THEN
     1495                      num_gp = num_gp + 1
     1496                      location(num_gp,1) = i * dx
     1497                      location(num_gp,2) = j * dy
     1498                      location(num_gp,3) = k * dz
     1499                      ei(num_gp)     = e(k+1,j,i)
     1500                      dissi(num_gp)  = diss(k+1,j,i)
     1501                      de_dxi(num_gp) = de_dx(k+1,j,i)
     1502                      de_dyi(num_gp) = de_dy(k+1,j,i)
     1503                      de_dzi(num_gp) = 0.0
     1504                   ENDIF
     1505 
     1506!
     1507!--                If wall between gridpoint 5 and gridpoint 7, then
     1508!--                Neumann boundary condition has to be applied
     1509!--                (only one case as only building beneath is possible)
     1510                   IF ( gp_outside_of_building(5) == 1  .AND. &
     1511                        gp_outside_of_building(7) == 0 )  THEN
     1512                      num_gp = num_gp + 1
     1513                      location(num_gp,1) = (i+1) * dx
     1514                      location(num_gp,2) = j * dy
     1515                      location(num_gp,3) = k * dz
     1516                      ei(num_gp)     = e(k+1,j,i+1)
     1517                      dissi(num_gp)  = diss(k+1,j,i+1)
     1518                      de_dxi(num_gp) = de_dx(k+1,j,i+1)
     1519                      de_dyi(num_gp) = de_dy(k+1,j,i+1)
     1520                      de_dzi(num_gp) = 0.0
     1521                   ENDIF
     1522
     1523!
     1524!--                If wall between gridpoint 2 and gridpoint 4, then
     1525!--                Neumann boundary condition has to be applied
     1526!--                (only one case as only building beneath is possible)
     1527                   IF ( gp_outside_of_building(2) == 1  .AND. &
     1528                        gp_outside_of_building(4) == 0 )  THEN
     1529                      num_gp = num_gp + 1
     1530                      location(num_gp,1) = i * dx
     1531                      location(num_gp,2) = (j+1) * dy
     1532                      location(num_gp,3) = k * dz
     1533                      ei(num_gp)     = e(k+1,j+1,i)
     1534                      dissi(num_gp)  = diss(k+1,j+1,i)
     1535                      de_dxi(num_gp) = de_dx(k+1,j+1,i)
     1536                      de_dyi(num_gp) = de_dy(k+1,j+1,i)
     1537                      de_dzi(num_gp) = 0.0
     1538                   ENDIF
     1539
     1540!
     1541!--                If wall between gridpoint 6 and gridpoint 8, then
     1542!--                Neumann boundary condition has to be applied
     1543!--                (only one case as only building beneath is possible)
     1544                   IF ( gp_outside_of_building(6) == 1  .AND. &
     1545                        gp_outside_of_building(8) == 0 )  THEN
     1546                      num_gp = num_gp + 1
     1547                      location(num_gp,1) = (i+1) * dx
     1548                      location(num_gp,2) = (j+1) * dy
     1549                      location(num_gp,3) = k * dz
     1550                      ei(num_gp)     = e(k+1,j+1,i+1)
     1551                      dissi(num_gp)  = diss(k+1,j+1,i+1)
     1552                      de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
     1553                      de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
     1554                      de_dzi(num_gp) = 0.0
     1555                   ENDIF
     1556
     1557!
     1558!--                Carry out the interpolation
     1559                   IF ( num_gp == 1 )  THEN
     1560!
     1561!--                   If only one of the gridpoints is situated outside of the
     1562!--                   building, it follows that the values at the particle
     1563!--                   location are the same as the gridpoint values
     1564                      e_int = ei(num_gp)
     1565
     1566                   ELSE IF ( num_gp > 1 )  THEN
     1567
     1568                      d_sum = 0.0
     1569!
     1570!--                   Evaluation of the distances between the gridpoints
     1571!--                   contributing to the interpolated values, and the particle
     1572!--                   location
     1573                      DO agp = 1, num_gp
     1574                         d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
     1575                                      + ( particles(n)%y-location(agp,2) )**2  &
     1576                                      + ( particles(n)%z-location(agp,3) )**2
     1577                         d_sum        = d_sum + d_gp_pl(agp)
     1578                      ENDDO
     1579
     1580!
     1581!--                   Finally the interpolation can be carried out
     1582                      e_int     = 0.0
     1583                      diss_int  = 0.0
     1584                      de_dx_int = 0.0
     1585                      de_dy_int = 0.0
     1586                      de_dz_int = 0.0
     1587                      DO agp = 1, num_gp
     1588                         e_int     = e_int     + ( d_sum - d_gp_pl(agp) ) * &
     1589                                                ei(agp) / ( (num_gp-1) * d_sum )
     1590                         diss_int  = diss_int  + ( d_sum - d_gp_pl(agp) ) * &
     1591                                             dissi(agp) / ( (num_gp-1) * d_sum )
     1592                         de_dx_int = de_dx_int + ( d_sum - d_gp_pl(agp) ) * &
     1593                                            de_dxi(agp) / ( (num_gp-1) * d_sum )
     1594                         de_dy_int = de_dy_int + ( d_sum - d_gp_pl(agp) ) * &
     1595                                            de_dyi(agp) / ( (num_gp-1) * d_sum )
     1596                         de_dz_int = de_dz_int + ( d_sum - d_gp_pl(agp) ) * &
     1597                                            de_dzi(agp) / ( (num_gp-1) * d_sum )
     1598                      ENDDO
     1599
     1600                   ENDIF
     1601
     1602                ENDIF
     1603
     1604             ENDIF 
    10721605
    10731606!
     
    11481681
    11491682             ELSE
     1683
     1684!
     1685!--             Restriction of the size of the new timestep: compared to the
     1686!--             previous timestep the increase must not exceed 200%
     1687
     1688                dt_particle_old = particles(n)%age - particles(n)%age_m
     1689                IF ( dt_particle > 2.0 * dt_particle_old )  THEN
     1690                   dt_particle = 2.0 * dt_particle_old
     1691                ENDIF
     1692
    11501693!
    11511694!--             For old particles the SGS components are correlated with the
    11521695!--             values from the previous timestep. Random numbers have also to
    11531696!--             be limited (see above).
     1697!--             As negative values for the subgrid TKE are not allowed, the
     1698!--             change of the subgrid TKE with time cannot be smaller than
     1699!--             -e_int/dt_particle. This value is used as a lower boundary
     1700!--             value for the change of TKE
     1701
     1702                de_dt_min = - e_int / dt_particle
     1703
     1704                de_dt = ( e_int - particles(n)%e_m ) / dt_particle_old
     1705
     1706                IF ( de_dt < de_dt_min )  THEN
     1707                   de_dt = de_dt_min
     1708                ENDIF
     1709
    11541710                particles(n)%speed_x_sgs = particles(n)%speed_x_sgs -          &
    11551711                       fs_int * c_0 * diss_int * particles(n)%speed_x_sgs *    &
    11561712                       dt_particle / ( 4.0 * sgs_wfu_part * e_int ) +          &
    1157                        ( ( 2.0 * sgs_wfu_part * ( e_int - particles(n)%e_m ) / &
    1158                          dt_particle ) * particles(n)%speed_x_sgs /            &
     1713                       ( 2.0 * sgs_wfu_part * de_dt *                          &
     1714                               particles(n)%speed_x_sgs /                      &
    11591715                         ( 2.0 * sgs_wfu_part * e_int ) + de_dx_int            &
    11601716                       ) * dt_particle / 2.0                        +          &
     
    11661722                       fs_int * c_0 * diss_int * particles(n)%speed_y_sgs *    &
    11671723                       dt_particle / ( 4.0 * sgs_wfv_part * e_int ) +          &
    1168                        ( ( 2.0 * sgs_wfv_part * ( e_int - particles(n)%e_m ) / &
    1169                          dt_particle ) * particles(n)%speed_y_sgs /            &
     1724                       ( 2.0 * sgs_wfv_part * de_dt *                          &
     1725                               particles(n)%speed_y_sgs /                      &
    11701726                         ( 2.0 * sgs_wfv_part * e_int ) + de_dy_int            &
    11711727                       ) * dt_particle / 2.0                        +          &
     
    11771733                       fs_int * c_0 * diss_int * particles(n)%speed_z_sgs *    &
    11781734                       dt_particle / ( 4.0 * sgs_wfw_part * e_int ) +          &
    1179                        ( ( 2.0 * sgs_wfw_part * ( e_int - particles(n)%e_m ) / &
    1180                          dt_particle ) * particles(n)%speed_z_sgs /            &
     1735                       ( 2.0 * sgs_wfw_part * de_dt *                          &
     1736                               particles(n)%speed_z_sgs /                      &
    11811737                         ( 2.0 * sgs_wfw_part * e_int ) + de_dz_int            &
    11821738                       ) * dt_particle / 2.0                        +          &
     
    12031759
    12041760          ENDIF
     1761
     1762!
     1763!--       Remember the old age of the particle ( needed to prevent that a
     1764!--       particle crosses several PEs during one timestep and for the
     1765!--       evaluation of the subgrid particle velocity fluctuations )
     1766          particles(n)%age_m = particles(n)%age
     1767
    12051768
    12061769!
     
    12701833
    12711834       ENDDO   ! advection loop
     1835
     1836
     1837!
     1838!--    Particle reflection from walls
     1839       CALL cpu_log( log_point_s(48), 'advec_part_refle', 'start' )
     1840
     1841       DO  n = 1, number_of_particles
     1842
     1843!-------------------------------------------------------
     1844      i2 =  (particles(n)%x + 0.5*dx) * ddx
     1845      j2 =  (particles(n)%y + 0.5*dy) * ddy
     1846      k2 =  particles(n)%z / dz + 1
     1847     
     1848      particles_x = particles(n)%x
     1849      particles_y = particles(n)%y
     1850      particles_z = particles(n)%z
     1851
     1852 
     1853    IF (k2<=nzb_s_inner(j2,i2).and.nzb_s_inner(j2,i2)/=0) THEN
     1854
     1855   old_positions_x = particles(n)%x - particles(n)%speed_x * dt_particle
     1856   old_positions_y = particles(n)%y - particles(n)%speed_y * dt_particle
     1857   old_positions_z = particles(n)%z - particles(n)%speed_z * dt_particle
     1858   i1= (old_positions_x +0.5*dx) * ddx
     1859   j1 =(old_positions_y +0.5*dy) * ddy
     1860   k1 = old_positions_z / dz
     1861   
     1862
     1863!-- CASE 1
     1864
     1865   IF ( particles(n)%x  >  old_positions_x .and. &
     1866         particles(n)%y  >  old_positions_y )THEN
     1867 
     1868     
     1869
     1870
     1871     t_index = 1
     1872   
     1873     Do i = i1, i2
     1874     xline = i*dx+0.5*dx
     1875     
     1876     t(t_index) =(xline-old_positions_x)/(particles(n)%x-old_positions_x)
     1877     t_index = t_index + 1
     1878     ENDDO
     1879
     1880     Do j = j1, j2
     1881     yline = j*dy+0.5*dy
     1882     
     1883     t(t_index) =(yline-old_positions_y)/(particles(n)%y-old_positions_y)
     1884     t_index = t_index + 1
     1885     ENDDO
     1886
     1887     IF ( particles(n)%z < old_positions_z ) THEN
     1888     Do k = k1, k2 , -1
     1889     zline = k*dz
     1890     t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z)
     1891     t_index = t_index + 1
     1892     ENDDO
     1893     ENDIF
     1894     t_index_number = t_index-1
     1895!-- sorting t
     1896                   inc = 1
     1897                   jr  = 1
     1898                   DO WHILE ( inc <= t_index_number )
     1899                      inc = 3 * inc + 1
     1900                   ENDDO
     1901
     1902                   DO WHILE ( inc > 1 )
     1903                      inc = inc / 3
     1904                      DO  ir = inc+1, t_index_number
     1905                         tmp_t = t(ir)
     1906                         jr = ir
     1907                         
     1908                         DO WHILE ( t(jr-inc) > tmp_t )
     1909                            t(jr) = t(jr-inc)
     1910                            jr = jr - inc
     1911                            IF ( jr <= inc )  THEN
     1912                           EXIT
     1913                            ENDIF
     1914                         ENDDO
     1915                        t(jr) = tmp_t
     1916                       
     1917             ENDDO
     1918        ENDDO
     1919   
     1920
     1921   
     1922   
     1923    Do t_index = 1, t_index_number
     1924 
     1925    positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x)
     1926    positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y)
     1927    positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z)
     1928
     1929     
     1930   
     1931     i3 = (positions_x + 0.5*dx) * ddx   
     1932     j3 = (positions_y + 0.5*dy) * ddy
     1933     k3 =  positions_z / dz
     1934
     1935     i5 = positions_x * ddx
     1936     j5 = positions_y * ddy
     1937     k5 = positions_z / dz
     1938
     1939     
     1940 
     1941
     1942
     1943    IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN
     1944             
     1945         IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. &
     1946          k3 == nzb_s_inner(j5,i3) ) THEN
     1947     
     1948          particles(n)%z = 2*positions_z - particles_z
     1949          particles(n)%speed_z = - particles(n)%speed_z
     1950         
     1951           IF ( use_sgs_for_particles )  THEN
     1952             particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
     1953           ENDIF
     1954          GOTO 999
     1955       ENDIF
     1956    ENDIF
     1957   
     1958IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN
     1959             
     1960         IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. &
     1961          k3 == nzb_s_inner(j3,i5) ) THEN
     1962     
     1963          particles(n)%z = 2*positions_z - particles_z
     1964          particles(n)%speed_z = - particles(n)%speed_z
     1965           
     1966            IF ( use_sgs_for_particles )  THEN
     1967              particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
     1968            ENDIF
     1969           
     1970          GOTO 999
     1971       ENDIF
     1972    ENDIF
     1973 
     1974     IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN
     1975
     1976       IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. &
     1977          k3 == nzb_s_inner(j3,i3) ) THEN
     1978     
     1979          particles(n)%z = 2*positions_z - particles_z
     1980          particles(n)%speed_z = - particles(n)%speed_z
     1981           
     1982           IF ( use_sgs_for_particles )  THEN
     1983          particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
     1984           ENDIF
     1985           
     1986          GOTO 999
     1987       ENDIF
     1988   
     1989       IF ( positions_y == j3*dy-0.5*dy .and. &
     1990         positions_z < nzb_s_inner(j3,i3)*dz) THEN
     1991             
     1992          particles(n)%y = 2*positions_y - particles_y
     1993          particles(n)%speed_y = - particles(n)%speed_y
     1994         
     1995          IF ( use_sgs_for_particles )  THEN
     1996            particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
     1997          ENDIF
     1998   
     1999          GOTO 999
     2000       ENDIF
     2001
     2002 
     2003       IF ( positions_x == i3*dx-0.5*dx.and. &
     2004            positions_z < nzb_s_inner(j3,i3)*dz ) THEN
     2005             
     2006          particles(n)%x = 2*positions_x - particles_x
     2007          particles(n)%speed_x = - particles(n)%speed_x
     2008       
     2009        IF ( use_sgs_for_particles )  THEN
     2010          particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
     2011          GOTO 999
     2012        ENDIF 
     2013       
     2014        ENDIF
     2015     
     2016     ENDIF
     2017    ENDDO
     2018   ENDIF
     2019   
     2020!-- CASE 2
     2021 
     2022     IF ( particles(n)%x  >  old_positions_x .and. &
     2023         particles(n)%y  <  old_positions_y )THEN
     2024 
     2025
     2026     t_index = 1
     2027   
     2028     Do i = i1, i2
     2029     xline = i*dx+0.5*dx
     2030     
     2031     t(t_index) =(xline-old_positions_x)/(particles(n)%x-old_positions_x)
     2032     t_index = t_index + 1
     2033     ENDDO
     2034
     2035     Do j = j1, j2,-1
     2036     yline = j*dy-0.5*dy
     2037     
     2038     t(t_index) =(old_positions_y-yline)/(old_positions_y-particles(n)%y)
     2039     t_index = t_index + 1
     2040     ENDDO
     2041
     2042     IF ( particles(n)%z < old_positions_z ) THEN
     2043     Do k = k1, k2 , -1
     2044     zline = k*dz
     2045     t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z)
     2046     t_index = t_index + 1
     2047     ENDDO
     2048     ENDIF
     2049     t_index_number = t_index-1
     2050
     2051!-- sorting t
     2052                   inc = 1
     2053                   jr  = 1
     2054                   DO WHILE ( inc <= t_index_number )
     2055                      inc = 3 * inc + 1
     2056                   ENDDO
     2057
     2058                   DO WHILE ( inc > 1 )
     2059                      inc = inc / 3
     2060                      DO  ir = inc+1, t_index_number
     2061                         tmp_t = t(ir)
     2062                         jr = ir
     2063                         
     2064                         DO WHILE ( t(jr-inc) > tmp_t )
     2065                            t(jr) = t(jr-inc)
     2066                            jr = jr - inc
     2067                            IF ( jr <= inc )  THEN
     2068                           EXIT
     2069                            ENDIF
     2070                         ENDDO
     2071                        t(jr) = tmp_t
     2072                       
     2073             ENDDO
     2074        ENDDO
     2075
     2076     
     2077 
     2078     
     2079   
     2080   
     2081    Do t_index = 1, t_index_number
     2082 
     2083    positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x)
     2084    positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y)
     2085    positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z)
     2086
     2087     
     2088   
     2089     i3 = (positions_x + 0.5*dx) * ddx   
     2090     j3 = (positions_y + 0.5*dy) * ddy
     2091     k3 =  positions_z / dz
     2092
     2093     i5 = positions_x * ddx
     2094     j5 = positions_y * ddy
     2095     k5 = positions_z / dz
     2096   
     2097
     2098 
     2099    IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN
     2100             
     2101         IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. &
     2102          k3 == nzb_s_inner(j3,i5) ) THEN
     2103     
     2104          particles(n)%z = 2*positions_z - particles_z
     2105          particles(n)%speed_z = - particles(n)%speed_z
     2106           
     2107           IF ( use_sgs_for_particles )  THEN
     2108            particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
     2109           ENDIF
     2110       
     2111          GOTO 999
     2112       ENDIF
     2113    ENDIF
     2114     
     2115     IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN
     2116
     2117       IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. &
     2118          k3 == nzb_s_inner(j3,i3) ) THEN
     2119     
     2120          particles(n)%z = 2*positions_z - particles_z
     2121          particles(n)%speed_z = - particles(n)%speed_z
     2122         
     2123         IF ( use_sgs_for_particles )  THEN
     2124          particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
     2125         ENDIF
     2126         
     2127         GOTO 999
     2128       ENDIF
     2129 
     2130       
     2131       IF ( positions_x == i3*dx-0.5*dx.and.&
     2132       positions_z < nzb_s_inner(j3,i3)*dz ) THEN
     2133             
     2134          particles(n)%x = 2*positions_x - particles_x
     2135          particles(n)%speed_x = - particles(n)%speed_x
     2136         
     2137          IF ( use_sgs_for_particles )  THEN
     2138           particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
     2139          ENDIF
     2140         
     2141          GOTO 999
     2142       ENDIF
     2143     
     2144     ENDIF
     2145
     2146     
     2147     IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN
     2148
     2149       IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. &
     2150          k3 == nzb_s_inner(j5,i3) ) THEN
     2151     
     2152          particles(n)%z = 2*positions_z - particles_z
     2153          particles(n)%speed_z = - particles(n)%speed_z
     2154           IF ( use_sgs_for_particles )  THEN
     2155            particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
     2156           ENDIF
     2157
     2158          GOTO 999
     2159       ENDIF
     2160 
     2161       IF ( positions_y == j5 * dy + 0.5*dy.and.&
     2162        positions_z < nzb_s_inner(j5,i3)*dz ) THEN
     2163     
     2164          particles(n)%y = 2*positions_y - particles_y
     2165          particles(n)%speed_y = - particles(n)%speed_y
     2166           
     2167          IF ( use_sgs_for_particles )  THEN
     2168            particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
     2169          ENDIF
     2170          GOTO 999
     2171       ENDIF
     2172
     2173     ENDIF
     2174    ENDDO
     2175   ENDIF
     2176
     2177
     2178!-- CASE 3
     2179
     2180
     2181
     2182        IF ( particles(n)%x < old_positions_x .and. &
     2183         particles(n)%y > old_positions_y )THEN
     2184 
     2185
     2186     t_index = 1
     2187   
     2188     Do i = i1, i2,-1
     2189     xline = i*dx-0.5*dx
     2190     
     2191     t(t_index) =(old_positions_x-xline)/(old_positions_x-particles(n)%x)
     2192     t_index = t_index + 1
     2193     ENDDO
     2194
     2195     Do j = j1, j2
     2196     yline = j*dy+0.5*dy
     2197     
     2198     t(t_index) =(yline-old_positions_y)/(particles(n)%y-old_positions_y)
     2199     t_index = t_index + 1
     2200     ENDDO
     2201     
     2202     IF ( particles(n)%z < old_positions_z ) THEN
     2203     Do k = k1, k2 , -1
     2204     zline = k*dz
     2205     t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z)
     2206     t_index = t_index + 1
     2207     ENDDO
     2208     ENDIF
     2209   
     2210     t_index_number = t_index-1
     2211
     2212!-- sorting t
     2213                   inc = 1
     2214                   jr  = 1
     2215                   DO WHILE ( inc <= t_index_number )
     2216                      inc = 3 * inc + 1
     2217                   ENDDO
     2218
     2219                   DO WHILE ( inc > 1 )
     2220                      inc = inc / 3
     2221                      DO  ir = inc+1, t_index_number
     2222                         tmp_t = t(ir)
     2223                         jr = ir
     2224                         
     2225                         DO WHILE ( t(jr-inc) > tmp_t )
     2226                            t(jr) = t(jr-inc)
     2227                            jr = jr - inc
     2228                            IF ( jr <= inc )  THEN
     2229                           EXIT
     2230                            ENDIF
     2231                         ENDDO
     2232                        t(jr) = tmp_t
     2233                       
     2234             ENDDO
     2235        ENDDO
     2236
     2237   
     2238   
     2239    Do t_index = 1, t_index_number
     2240 
     2241    positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x)
     2242    positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y)
     2243    positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z)
     2244
     2245     
     2246   
     2247     i3 = (positions_x + 0.5*dx) * ddx   
     2248     j3 = (positions_y + 0.5*dy) * ddy
     2249     k3 =  positions_z / dz
     2250
     2251     i5 = positions_x * ddx
     2252     j5 = positions_y * ddy
     2253     k5 = positions_z / dz
     2254
     2255
     2256
     2257     IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0 ) THEN
     2258             
     2259         IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. &
     2260          k3 == nzb_s_inner(j5,i3) ) THEN
     2261     
     2262          particles(n)%z = 2*positions_z - particles_z
     2263          particles(n)%speed_z = - particles(n)%speed_z
     2264         
     2265         IF ( use_sgs_for_particles )  THEN 
     2266          particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
     2267         ENDIF
     2268
     2269          GOTO 999
     2270       ENDIF
     2271    ENDIF
     2272     
     2273
     2274
     2275      IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN
     2276
     2277       
     2278      IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. &
     2279          k3 == nzb_s_inner(j3,i3) ) THEN
     2280     
     2281          particles(n)%z = 2*positions_z - particles_z
     2282          particles(n)%speed_z = - particles(n)%speed_z
     2283         
     2284         IF ( use_sgs_for_particles )  THEN
     2285          particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
     2286         ENDIF
     2287
     2288          GOTO 999
     2289       ENDIF
     2290
     2291          IF ( positions_y == j3*dy-0.5*dy .and. &
     2292            positions_z < nzb_s_inner(j3,i3)*dz) THEN
     2293             
     2294          particles(n)%y = 2*positions_y - particles_y
     2295          particles(n)%speed_y = - particles(n)%speed_y
     2296         
     2297           IF ( use_sgs_for_particles )  THEN
     2298              particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
     2299           ENDIF
     2300
     2301          GOTO 999
     2302       ENDIF 
     2303
     2304     ENDIF
     2305     
     2306       
     2307     IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN
     2308
     2309       IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. &
     2310          k3 == nzb_s_inner(j3,i5) ) THEN
     2311     
     2312          particles(n)%z = 2*positions_z - particles_z
     2313          particles(n)%speed_z = - particles(n)%speed_z
     2314         
     2315         IF ( use_sgs_for_particles )  THEN
     2316          particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
     2317         ENDIF
     2318
     2319
     2320          GOTO 999
     2321       ENDIF
     2322 
     2323       IF ( positions_x == i5 * dx + 0.5*dx .and. &
     2324               positions_z < nzb_s_inner(j3,i5)*dz ) THEN
     2325     
     2326          particles(n)%x= 2*positions_x - particles_x
     2327          particles(n)%speed_x = - particles(n)%speed_x
     2328         
     2329         IF ( use_sgs_for_particles )  THEN
     2330           particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
     2331         ENDIF 
     2332       
     2333          GOTO 999
     2334       ENDIF
     2335
     2336     ENDIF
     2337    ENDDO
     2338   ENDIF
     2339
     2340
     2341!-- CASE 4
     2342     
     2343
     2344       IF ( particles(n)%x < old_positions_x .and. &
     2345         particles(n)%y  < old_positions_y )THEN
     2346 
     2347
     2348     t_index = 1
     2349   
     2350     Do i = i1, i2,-1
     2351     xline = i*dx-0.5*dx
     2352     
     2353     t(t_index) =(old_positions_x-xline)/(old_positions_x-particles(n)%x)
     2354     t_index = t_index + 1
     2355     
     2356     
     2357     
     2358     ENDDO
     2359
     2360     Do j = j1, j2,-1
     2361     yline = j*dy-0.5*dy
     2362     
     2363     t(t_index) =(old_positions_y-yline)/(old_positions_y-particles(n)%y)
     2364     t_index = t_index + 1
     2365     
     2366     
     2367     ENDDO
     2368
     2369     IF ( particles(n)%z < old_positions_z ) THEN
     2370     Do k = k1, k2 , -1
     2371     zline = k*dz
     2372     t(t_index) = (old_positions_z-zline)/(old_positions_z-particles(n)%z)
     2373     t_index = t_index + 1
     2374     
     2375
     2376     ENDDO
     2377     ENDIF
     2378
     2379     t_index_number = t_index-1
     2380
     2381!-- sorting t
     2382                   inc = 1
     2383                   jr  = 1
     2384                   DO WHILE ( inc <= t_index_number )
     2385                      inc = 3 * inc + 1
     2386                   ENDDO
     2387
     2388                   DO WHILE ( inc > 1 )
     2389                      inc = inc / 3
     2390                      DO  ir = inc+1, t_index_number
     2391                         tmp_t = t(ir)
     2392                         jr = ir
     2393                         
     2394                         DO WHILE ( t(jr-inc) > tmp_t )
     2395                            t(jr) = t(jr-inc)
     2396                            jr = jr - inc
     2397                            IF ( jr <= inc )  THEN
     2398                           EXIT
     2399                            ENDIF
     2400                         ENDDO
     2401                        t(jr) = tmp_t
     2402                       
     2403             ENDDO
     2404        ENDDO
     2405   
     2406   
     2407    Do t_index = 1, t_index_number
     2408 
     2409    positions_x = old_positions_x + t(t_index)*(particles_x - old_positions_x)
     2410    positions_y = old_positions_y + t(t_index)*(particles_y - old_positions_y)
     2411    positions_z = old_positions_z + t(t_index)*(particles_z - old_positions_z)
     2412
     2413     
     2414   
     2415     i3 = (positions_x + 0.5*dx) * ddx   
     2416     j3 = (positions_y + 0.5*dy) * ddy
     2417     k3 =  positions_z / dz
     2418
     2419     i5 = positions_x * ddx
     2420     j5 = positions_y * ddy
     2421     k5 = positions_z / dz
     2422
     2423
     2424
     2425
     2426     IF (k3 <= nzb_s_inner(j3,i3).and.nzb_s_inner(j3,i3)/=0 ) THEN
     2427
     2428       IF ( positions_z == nzb_s_inner(j3,i3)*dz.and. &
     2429          k3 == nzb_s_inner(j3,i3) ) THEN
     2430     
     2431          particles(n)%z = 2*positions_z - particles_z
     2432          particles(n)%speed_z = - particles(n)%speed_z
     2433         
     2434          IF ( use_sgs_for_particles )  THEN
     2435           particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
     2436          ENDIF
     2437         
     2438
     2439          GOTO 999
     2440       ENDIF
     2441     ENDIF
     2442           
     2443       
     2444     IF (k5 <= nzb_s_inner(j3,i5).and.nzb_s_inner(j3,i5)/=0 ) THEN
     2445
     2446       IF ( positions_z == nzb_s_inner(j3,i5)*dz.and. &
     2447          k3 == nzb_s_inner(j3,i5) ) THEN
     2448     
     2449          particles(n)%z = 2*positions_z - particles_z
     2450          particles(n)%speed_z = - particles(n)%speed_z
     2451         
     2452         IF ( use_sgs_for_particles )  THEN
     2453          particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
     2454         ENDIF 
     2455       
     2456          GOTO 999
     2457       ENDIF
     2458 
     2459       IF ( positions_x == i5 * dx + 0.5*dx .and. nzb_s_inner(j3,i5)/=0.and.&
     2460                 positions_z < nzb_s_inner(j3,i5)*dz) THEN
     2461     
     2462          particles(n)%x= 2*positions_x - particles_x
     2463          particles(n)%speed_x = - particles(n)%speed_x
     2464         
     2465         IF ( use_sgs_for_particles )  THEN 
     2466           particles(n)%speed_x_sgs = - particles(n)%speed_x_sgs
     2467         ENDIF
     2468          GOTO 999
     2469       ENDIF
     2470     
     2471    ENDIF
     2472   
     2473     IF (k5 <= nzb_s_inner(j5,i3).and.nzb_s_inner(j5,i3)/=0) THEN
     2474       
     2475
     2476       IF ( positions_z == nzb_s_inner(j5,i3)*dz.and. &
     2477          k5 == nzb_s_inner(j5,i3) ) THEN
     2478     
     2479          particles(n)%z = 2*positions_z - particles_z
     2480          particles(n)%speed_z = - particles(n)%speed_z
     2481         
     2482          IF ( use_sgs_for_particles )  THEN
     2483           particles(n)%speed_z_sgs = - particles(n)%speed_z_sgs
     2484          ENDIF
     2485         
     2486          GOTO 999
     2487       ENDIF
     2488
     2489          IF ( positions_y == j5 * dy + 0.5*dy .and. nzb_s_inner(j5,i3)/=0.and.&
     2490                    positions_z < nzb_s_inner(j5,i3)*dz ) THEN
     2491     
     2492          particles(n)%y= 2*positions_y - particles_y
     2493          particles(n)%speed_y = - particles(n)%speed_y
     2494         
     2495         IF ( use_sgs_for_particles )  THEN
     2496           particles(n)%speed_y_sgs = - particles(n)%speed_y_sgs
     2497         ENDIF 
     2498
     2499          GOTO 999
     2500       ENDIF
     2501
     2502
     2503      ENDIF
     2504    ENDDO
     2505   ENDIF
     2506
     2507     999 CONTINUE
     2508   
     2509   ENDIF
     2510
     2511!---------------------------------------------------
     2512
     2513       ENDDO   ! particle reflection from walls
     2514
     2515       CALL cpu_log( log_point_s(48), 'advec_part_refle', 'stop' )
     2516
     2517!
     2518!--    User-defined actions after the evaluation of the new particle position
     2519       CALL user_advec_particles
    12722520
    12732521!
     
    15002748          trlp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    15012749                                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    1502                                 0, 0, 0, 0 )
     2750                                0.0, 0, 0, 0, 0 )
    15032751          trrp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    15042752                                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    1505                                 0, 0, 0, 0 )
     2753                                0.0, 0, 0, 0, 0 )
    15062754
    15072755          IF ( use_particle_tails )  THEN
     
    19333181          trsp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    19343182                                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    1935                                 0, 0, 0, 0 )
     3183                                0.0, 0, 0, 0, 0 )
    19363184          trnp = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    19373185                                0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    1938                                 0, 0, 0, 0 )
     3186                                0.0, 0, 0, 0, 0 )
    19393187
    19403188          IF ( use_particle_tails )  THEN
     
    24603708
    24613709          nn = particles(n)%tail_id
     3710
     3711!
     3712!--       Stop if particles have moved further than the length of one
     3713!--       PE subdomain
     3714          IF ( ABS(particles(n)%speed_x) >                                  &
     3715               ((nxr-nxl+2)*dx)/(particles(n)%age-particles(n)%age_m)  .OR. &
     3716               ABS(particles(n)%speed_y) >                                  &
     3717               ((nyn-nys+2)*dy)/(particles(n)%age-particles(n)%age_m) )  THEN
     3718
     3719             PRINT*, '+++ advec_particles: particle too fast.  n = ', n
     3720#if defined( __parallel )
     3721             CALL MPI_ABORT( comm2d, 9999, ierr )
     3722#else
     3723             CALL local_stop
     3724#endif
     3725          ENDIF
    24623726
    24633727          IF ( particles(n)%age > particle_maximum_age  .AND.  &
  • palm/trunk/SOURCE/buoyancy.f90

    r4 r57  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Reference temperature pt_reference can be used.
    77!
    88! Former revisions:
     
    5959!
    6060!--       Normal case: horizontal surface
    61           DO  i = nxl, nxr
    62              DO  j = nys, nyn
    63                 DO  k = nzb_s_inner(j,i)+1, nzt-1
    64                     tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                    &
     61          IF ( use_pt_reference )  THEN
     62             DO  i = nxl, nxr
     63                DO  j = nys, nyn
     64                   DO  k = nzb_s_inner(j,i)+1, nzt-1
     65                      tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                 &
     66                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / pt_reference + &
     67                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference   &
     68                                                            )
     69                   ENDDO
     70                ENDDO
     71             ENDDO
     72          ELSE
     73             DO  i = nxl, nxr
     74                DO  j = nys, nyn
     75                   DO  k = nzb_s_inner(j,i)+1, nzt-1
     76                      tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                  &
    6577                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
    6678                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
    67                                                           )
    68                 ENDDO
    69              ENDDO
    70           ENDDO
     79                                                            )
     80                   ENDDO
     81                ENDDO
     82             ENDDO
     83          ENDIF
    7184
    7285       ELSE
     
    136149!
    137150!--       Normal case: horizontal surface
    138           DO  k = nzb_s_inner(j,i)+1, nzt-1
    139               tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                          &
     151          IF ( use_pt_reference )  THEN
     152             DO  k = nzb_s_inner(j,i)+1, nzt-1
     153                 tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                      &
     154                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / pt_reference + &
     155                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference   &
     156                                                       )
     157             ENDDO
     158          ELSE
     159             DO  k = nzb_s_inner(j,i)+1, nzt-1
     160                 tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                       &
    140161                        ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
    141162                        ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
    142                                                     )
    143           ENDDO
     163                                                       )
     164             ENDDO
     165          ENDIF
    144166
    145167       ELSE
  • palm/trunk/SOURCE/check_parameters.f90

    r51 r57  
    55! -----------------
    66! "by_user" allowed as initializing action, -data_output_ts,
    7 ! leapfrog with non-flat topography not allowed any more
     7! leapfrog with non-flat topography not allowed any more, pt_reference is
     8! checked
    89!
    910! Former revisions:
     
    574575
    575576!
     577!-- Reference temperature to be used in buoyancy terms
     578    IF ( pt_reference /= 9999999.9 )  use_pt_reference = .TRUE.
     579
     580!
    576581!-- In case of a given slope, compute the relevant quantities
    577582    IF ( alpha_surface /= 0.0 )  THEN
  • palm/trunk/SOURCE/diffusion_e.f90

    r39 r57  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Reference temperature pt_reference can be used in buoyancy term
    77!
    88! Former revisions:
     
    8282                dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
    8383                IF ( dpt_dz > 0.0 ) THEN
    84                    l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    85                                      SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
     84                   IF ( use_pt_reference )  THEN
     85                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
     86                                        SQRT( g / pt_reference * dpt_dz ) + 1E-5
     87                   ELSE
     88                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
     89                                        SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
     90                   ENDIF
    8691                ELSE
    8792                   l_stable = l_grid(k)
     
    195200          dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
    196201          IF ( dpt_dz > 0.0 ) THEN
    197              l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    198                                SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
     202             IF ( use_pt_reference )  THEN
     203                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
     204                                  SQRT( g / pt_reference * dpt_dz ) + 1E-5
     205             ELSE
     206                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
     207                                  SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
     208             ENDIF
    199209          ELSE
    200210             l_stable = l_grid(k)
  • palm/trunk/SOURCE/diffusion_u.f90

    r56 r57  
    44! Actual revisions:
    55! -----------------
    6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes
     6! Wall functions now include diabatic conditions, call of routine wall_fluxes,
     7! z0 removed from argument list
    78!
    89! Former revisions:
     
    5051! Call for all grid points
    5152!------------------------------------------------------------------------------!
    52     SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w, z0 )
     53    SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w )
    5354
    5455       USE control_parameters
     
    6162       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
    6263       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
    63        REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
    6464       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
    6565       REAL, DIMENSION(:,:),   POINTER ::  usws
     
    207207!------------------------------------------------------------------------------!
    208208    SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
    209                                v, w, z0 )
     209                               v, w )
    210210
    211211       USE control_parameters
     
    218218       REAL    ::  kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp
    219219       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
    220        REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
    221220       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
    222221       REAL, DIMENSION(nzb:nzt+1)      ::  usvs
  • palm/trunk/SOURCE/diffusion_v.f90

    r56 r57  
    44! Actual revisions:
    55! -----------------
    6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes
     6! Wall functions now include diabatic conditions, call of routine wall_fluxes,
     7! z0 removed from argument list
    78!
    89! Former revisions:
     
    4849! Call for all grid points
    4950!------------------------------------------------------------------------------!
    50     SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w, z0 )
     51    SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w )
    5152
    5253       USE control_parameters
     
    5960       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
    6061       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
    61        REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
    6262       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
    6363       REAL, DIMENSION(:,:),   POINTER ::  vsws
     
    205205!------------------------------------------------------------------------------!
    206206    SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
    207                                vsws, w, z0 )
     207                               vsws, w )
    208208
    209209       USE control_parameters
     
    216216       REAL    ::  kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp
    217217       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
    218        REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
    219218       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
    220219       REAL, DIMENSION(nzb:nzt+1)      ::  vsus
  • palm/trunk/SOURCE/diffusion_w.f90

    r56 r57  
    44! Actual revisions:
    55! -----------------
    6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes
     6! Wall functions now include diabatic conditions, call of routine wall_fluxes,
     7! z0 removed from argument list
    78!
    89! Former revisions:
     
    4647!------------------------------------------------------------------------------!
    4748    SUBROUTINE diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, &
    48                             w, z0 )
     49                            w )
    4950
    5051       USE control_parameters
     
    5960       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
    6061                   km_damp_y(nys-1:nyn+1)
    61        REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
    6262       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
    6363       REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
     
    194194!------------------------------------------------------------------------------!
    195195    SUBROUTINE diffusion_w_ij( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
    196                                tend, u, v, w, z0 )
     196                               tend, u, v, w )
    197197
    198198       USE control_parameters
     
    207207       REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1), &
    208208                   km_damp_y(nys-1:nyn+1)
    209        REAL    ::  z0(nys-1:nyn+1,nxl-1:nxr+1)
    210209       REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
    211210       REAL, DIMENSION(nzb:nzt+1)      ::  wsus, wsvs
  • palm/trunk/SOURCE/diffusivities.f90

    r4 r57  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Reference temperature pt_reference can be used in buoyancy term
    77!
    88! Former revisions:
     
    9090             dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
    9191             IF ( dpt_dz > 0.0 ) THEN
    92                 l_stable = 0.76 * sqrt_e(k) / &
    93                                   SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
     92                IF ( use_pt_reference )  THEN
     93                   l_stable = 0.76 * sqrt_e(k) / &
     94                                     SQRT( g / pt_reference * dpt_dz ) + 1E-5
     95                ELSE
     96                   l_stable = 0.76 * sqrt_e(k) / &
     97                                     SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
     98                ENDIF
    9499             ELSE
    95100                l_stable = l_grid(k)
  • palm/trunk/SOURCE/header.f90

    r46 r57  
    861861!-- Other quantities
    862862    WRITE ( io, 411 )  g
     863    IF ( use_pt_reference )  WRITE ( io, 412 )  pt_reference
    863864
    864865!
    865866!-- Cloud physics parameters
    866867    IF ( cloud_physics ) THEN
    867        WRITE ( io, 412 )
    868        WRITE ( io, 413 ) surface_pressure, r_d, rho_surface, cp, l_v
     868       WRITE ( io, 415 )
     869       WRITE ( io, 416 ) surface_pressure, r_d, rho_surface, cp, l_v
    869870    ENDIF
    870871
     
    13261327            '                            f*    = ',F9.6,' 1/s')
    13271328411 FORMAT (/'    Gravity             :   g     = ',F4.1,' m/s**2')
    1328 412 FORMAT (/'    Cloud physics parameters:'/ &
     1329412 FORMAT (/'    Reference temperature in buoyancy terms: ',F8.4,' K')
     1330415 FORMAT (/'    Cloud physics parameters:'/ &
    13291331             '    ------------------------'/)
    1330 413 FORMAT ('        Surface pressure   :   p_0   = ',F7.2,' hPa'/      &
     1332416 FORMAT ('        Surface pressure   :   p_0   = ',F7.2,' hPa'/      &
    13311333            '        Gas constant       :   R     = ',F5.1,' J/(kg K)'/ &
    13321334            '        Density of air     :   rho_0 = ',F5.3,' kg/m**3'/  &
  • palm/trunk/SOURCE/init_particles.f90

    r39 r57  
    44! Actual revisions:
    55! -----------------
    6 !
     6! displacements for mpi_particle_type changed, age_m initialized
    77!
    88! Former revisions:
     
    6161    blocklengths(1)  = 18; blocklengths(2)  = 4; blocklengths(3)  =  1
    6262#if defined( __t3eb )
    63     displacements(1) = 0; displacements(2) = 144; displacements(3) = 176
     63    displacements(1) = 0; displacements(2) = 152; displacements(3) = 184
    6464#else
    65     displacements(1) = 0; displacements(2) = 144; displacements(3) = 160
     65    displacements(1) = 0; displacements(2) = 152; displacements(3) = 168
    6666#endif
    6767    types(1) = MPI_REAL
     
    192192       particles = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    193193                                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    194                                   0, 0, 0, 0 )
     194                                  0.0, 0, 0, 0, 0 )
    195195       particle_groups = particle_groups_type( 0.0, 0.0, 0.0, 0.0 )
    196196
     
    267267                            particles(n)%z             = pos_z
    268268                            particles(n)%age           = 0.0
     269                            particles(n)%age_m         = 0.0
    269270                            particles(n)%dt_sum        = 0.0
    270271                            particles(n)%dvrp_psize    = dvrp_psize
  • palm/trunk/SOURCE/modules.f90

    r51 r57  
    77! +array rif_wall
    88! +netcdf_64bit_3d, zu_s_inner, zw_w_inner, id_var_zusi_*, id_var_zwwi_*,
    9 ! ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc
     9! ts_value, u_nzb_p1_for_vfc, v_nzb_p1_for_vfc, pt_reference, use_pt_reference,
     10! +age_m in particle_type
    1011! -data_output_ts, dots_n
    1112! arrays dots_label and dots_unit now dimensioned with dots_max
     
    296297                sloping_surface = .FALSE., stop_dt = .FALSE., &
    297298                terminate_run = .FALSE., use_prior_plot1d_parameters = .FALSE.,&
    298                 use_surface_fluxes = .FALSE., use_top_fluxes = .FALSE., &
    299                 use_ug_for_galilei_tr = .TRUE., &
     299                use_pt_reference = .FALSE., use_surface_fluxes = .FALSE., &
     300                use_top_fluxes = .FALSE., use_ug_for_galilei_tr = .TRUE., &
    300301                use_upstream_for_tke = .FALSE., wall_adjustment = .TRUE.
    301302
     
    330331             overshoot_limit_u = 0.0, overshoot_limit_v = 0.0, &
    331332             overshoot_limit_w = 0.0, particle_maximum_age = 9999999.9, &
    332              phi = 55.0, prandtl_number = 1.0, pt_slope_offset = 0.0, &
    333              pt_surface = 300.0, pt_surface_initial_change = 0.0, &
    334              q_surface = 0.0, q_surface_initial_change = 0.0, &
    335              rayleigh_damping_factor = -1.0, rayleigh_damping_height = -1.0, &
    336              residual_limit = 1.0E-4, restart_time = 9999999.9, rif_max = 1.0, &
    337              rif_min = -5.0, roughness_length = 0.1, simulated_time = 0.0, &
     333             phi = 55.0, prandtl_number = 1.0, pt_reference = 9999999.9, &
     334             pt_slope_offset = 0.0, pt_surface = 300.0, &
     335             pt_surface_initial_change = 0.0, q_surface = 0.0, &
     336             q_surface_initial_change = 0.0, rayleigh_damping_factor = -1.0, &
     337             rayleigh_damping_height = -1.0, residual_limit = 1.0E-4, &
     338             restart_time = 9999999.9, rif_max = 1.0, rif_min = -5.0, &
     339             roughness_length = 0.1, simulated_time = 0.0, &
    338340             simulated_time_at_begin, sin_alpha_surface, &
    339341             skip_time_data_output = 0.0, skip_time_data_output_av = 9999999.9,&
     
    828830    TYPE particle_type
    829831       SEQUENCE
    830        REAL    ::  age, dt_sum, dvrp_psize, e_m, origin_x, origin_y, origin_z, &
    831                    radius, speed_x, speed_x_sgs, speed_y, speed_y_sgs,         &
    832                    speed_z, speed_z_sgs, weight_factor, x, y, z
     832       REAL    ::  age, age_m, dt_sum, dvrp_psize, e_m, origin_x, origin_y, &
     833                   origin_z, radius, speed_x, speed_x_sgs, speed_y,         &
     834                   speed_y_sgs, speed_z, speed_z_sgs, weight_factor, x, y, z
    833835       INTEGER ::  color, group, tailpoints, tail_id
    834836    END TYPE particle_type
  • palm/trunk/SOURCE/palm.f90

    r4 r57  
    1 #if defined( __vtk )
    2  SUBROUTINE palm
    3 #else
    41 PROGRAM palm
    5 #endif
    62
    73!------------------------------------------------------------------------------!
    84! Actual revisions:
    95! -----------------
    10 ! TEST: open(9)
     6! __vtk directives removed, open unit 9 for debug output
    117!
    128! Former revisions:
     
    9995
    10096!
     97!-- Open a file for debug output
     98    OPEN( 9, FILE='DEBUG'//myid_char, FORM='FORMATTED' )
     99
     100!
    101101!-- Generate grid parameters
    102102    CALL init_grid
     
    106106    CALL check_parameters
    107107
    108     OPEN( 9, FILE='DEBUG'//myid_char, FORM='FORMATTED' )
    109108!
    110109!-- Initialize all necessary variables
     
    169168#endif
    170169
    171 #if defined( __vtk )
    172  END SUBROUTINE palm
    173 #else
    174170 END PROGRAM palm
    175 #endif
    176171
    177 
  • palm/trunk/SOURCE/parin.f90

    r48 r57  
    44! Actual revisions:
    55! -----------------
    6 ! +netcdf_64bit_3d in d3par, -data_output_ts
     6! +netcdf_64bit_3d in d3par, +pt_reference in inipar, -data_output_ts
    77!
    88! Former revisions:
     
    6565                       overshoot_limit_pt, overshoot_limit_u, &
    6666                       overshoot_limit_v, overshoot_limit_w, passive_scalar, &
    67                        phi, prandtl_layer, precipitation, pt_surface, &
    68                        pt_surface_initial_change, pt_vertical_gradient, &
    69                        pt_vertical_gradient_level, q_surface, &
    70                        q_surface_initial_change, q_vertical_gradient, &
    71                        q_vertical_gradient_level, radiation, random_generator, &
    72                        random_heatflux, rif_max, rif_min, roughness_length, &
    73                        scalar_advec, statistic_regions, surface_heatflux, &
    74                        surface_pressure, surface_scalarflux, surface_waterflux,&
    75                        s_surface, s_surface_initial_change, &
    76                        s_vertical_gradient, s_vertical_gradient_level, &
    77                        top_heatflux, timestep_scheme, topography, ug_surface, &
     67                       phi, prandtl_layer, precipitation, pt_reference, &
     68                       pt_surface, pt_surface_initial_change, &
     69                       pt_vertical_gradient, pt_vertical_gradient_level, &
     70                       q_surface, q_surface_initial_change, &
     71                       q_vertical_gradient, q_vertical_gradient_level, &
     72                       radiation, random_generator, random_heatflux, rif_max, &
     73                       rif_min, roughness_length, scalar_advec, &
     74                       statistic_regions, surface_heatflux, surface_pressure, &
     75                       surface_scalarflux, surface_waterflux, s_surface, &
     76                       s_surface_initial_change, s_vertical_gradient, &
     77                       s_vertical_gradient_level, top_heatflux, &
     78                       timestep_scheme, topography, ug_surface, &
    7879                       ug_vertical_gradient, ug_vertical_gradient_level, &
    7980                       ups_limit_e, ups_limit_pt, ups_limit_u, ups_limit_v, &
  • palm/trunk/SOURCE/production_e.f90

    r56 r57  
    44! Actual revisions:
    55! -----------------
    6 ! Wall functions now include diabatic conditions, call of routine wall_fluxes_e
     6! Wall functions now include diabatic conditions, call of routine wall_fluxes_e,
     7! reference temperature pt_reference can be used in buoyancy term
    78!
    89! Former revisions:
     
    359360          IF ( .NOT. moisture )  THEN
    360361
    361              DO  j = nys, nyn
    362 
    363                 DO  k = nzb_diff_s_inner(j,i), nzt_diff
    364                    tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
    365                                     ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
     362             IF ( use_pt_reference )  THEN
     363
     364                DO  j = nys, nyn
     365                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
     366                      tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g/pt_reference * &
     367                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
     368                   ENDDO
     369
     370                   IF ( use_surface_fluxes )  THEN
     371                      k = nzb_diff_s_inner(j,i)-1
     372                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
     373                   ENDIF
     374
     375                   IF ( use_top_fluxes )  THEN
     376                      k = nzt
     377                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
     378                   ENDIF
    366379                ENDDO
    367380
    368                 IF ( use_surface_fluxes )  THEN
    369                    k = nzb_diff_s_inner(j,i)-1
    370                    tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
    371                 ENDIF
    372 
    373                 IF ( use_top_fluxes )  THEN
    374                    k = nzt
    375                    tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
    376                 ENDIF
    377 
    378              ENDDO
     381             ELSE
     382
     383                DO  j = nys, nyn
     384                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
     385                      tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
     386                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
     387                   ENDDO
     388
     389                   IF ( use_surface_fluxes )  THEN
     390                      k = nzb_diff_s_inner(j,i)-1
     391                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
     392                   ENDIF
     393
     394                   IF ( use_top_fluxes )  THEN
     395                      k = nzt
     396                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
     397                   ENDIF
     398                ENDDO
     399
     400             ENDIF
    379401
    380402          ELSE
     
    736758       IF ( .NOT. moisture )  THEN
    737759
    738           DO  k = nzb_diff_s_inner(j,i), nzt_diff
    739              tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
    740                                     ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
    741           ENDDO
    742 
    743           IF ( use_surface_fluxes )  THEN
    744              k = nzb_diff_s_inner(j,i)-1
    745              tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
    746           ENDIF
    747 
    748           IF ( use_top_fluxes )  THEN
    749              k = nzt
    750              tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
    751           ENDIF
     760          IF ( use_pt_reference )  THEN
     761
     762             DO  k = nzb_diff_s_inner(j,i), nzt_diff
     763                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
     764                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
     765             ENDDO
     766
     767             IF ( use_surface_fluxes )  THEN
     768                k = nzb_diff_s_inner(j,i)-1
     769                tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
     770             ENDIF
     771
     772             IF ( use_top_fluxes )  THEN
     773                k = nzt
     774                tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
     775             ENDIF
     776
     777          ELSE
     778
     779             DO  k = nzb_diff_s_inner(j,i), nzt_diff
     780                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
     781                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
     782             ENDDO
     783
     784             IF ( use_surface_fluxes )  THEN
     785                k = nzb_diff_s_inner(j,i)-1
     786                tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
     787             ENDIF
     788
     789             IF ( use_top_fluxes )  THEN
     790                k = nzt
     791                tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
     792             ENDIF
     793
     794         ENDIF
    752795
    753796       ELSE
  • palm/trunk/SOURCE/prognostic_equations.f90

    r39 r57  
    44! Actual revisions:
    55! -----------------
    6 !
     6! z0 removed from arguments in calls of diffusion_u/v/w
    77!
    88! Former revisions:
     
    130130          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    131131             CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, u_m,   &
    132                                usws_m, v_m, w_m, z0 )
     132                               usws_m, v_m, w_m )
    133133          ELSE
    134134             CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
    135                                v, w, z0 )
     135                               v, w )
    136136          ENDIF
    137137          CALL coriolis( i, j, 1 )
     
    198198          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    199199             CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, u_m, &
    200                                v_m, vsws_m, w_m, z0 )
     200                               v_m, vsws_m, w_m )
    201201          ELSE
    202202             CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v,  &
    203                                vsws, w, z0 )
     203                               vsws, w )
    204204          ENDIF
    205205          CALL coriolis( i, j, 2 )
     
    265265          IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    266266             CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, km_damp_y,  &
    267                                tend, u_m, v_m, w_m, z0 )
     267                               tend, u_m, v_m, w_m )
    268268          ELSE
    269269             CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y,    &
    270                                tend, u, v, w, z0 )
     270                               tend, u, v, w )
    271271          ENDIF
    272272          CALL coriolis( i, j, 3 )
     
    668668             THEN
    669669                CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend,  &
    670                                   u_m, usws_m, v_m, w_m, z0 )
     670                                  u_m, usws_m, v_m, w_m )
    671671             ELSE
    672672                CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
    673                                   usws, v, w, z0 )
     673                                  usws, v, w )
    674674             ENDIF
    675675             CALL coriolis( i, j, 1 )
     
    718718             THEN
    719719                CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend,     &
    720                                   u_m, v_m, vsws_m, w_m, z0 )
     720                                  u_m, v_m, vsws_m, w_m )
    721721             ELSE
    722722                CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
    723                                   vsws, w, z0 )
     723                                  vsws, w )
    724724             ENDIF
    725725             CALL coriolis( i, j, 2 )
     
    767767             THEN
    768768                CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x,          &
    769                                   km_damp_y, tend, u_m, v_m, w_m, z0 )
     769                                  km_damp_y, tend, u_m, v_m, w_m )
    770770             ELSE
    771771                CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
    772                                   tend, u, v, w, z0 )
     772                                  tend, u, v, w )
    773773             ENDIF
    774774             CALL coriolis( i, j, 3 )
     
    10411041    IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    10421042       CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, v_m, &
    1043                          w_m, z0 )
     1043                         w_m )
    10441044    ELSE
    1045        CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w, z0 )
     1045       CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w )
    10461046    ENDIF
    10471047    CALL coriolis( 1 )
     
    11131113    IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    11141114       CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
    1115                          w_m, z0 )
     1115                         w_m )
    11161116    ELSE
    1117        CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w, z0 )
     1117       CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w )
    11181118    ENDIF
    11191119    CALL coriolis( 2 )
     
    11841184    IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
    11851185       CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, &
    1186                          v_m, w_m, z0 )
     1186                         v_m, w_m )
    11871187    ELSE
    1188        CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w, &
    1189                          z0 )
     1188       CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w )
    11901189    ENDIF
    11911190    CALL coriolis( 3 )
  • palm/trunk/SOURCE/read_var_list.f90

    r39 r57  
    44! Actual revisions:
    55! -----------------
    6 !
     6! +pt_reference
    77!
    88! Former revisions:
     
    252252          CASE ( 'pt_init' )
    253253             READ ( 13 )  pt_init
     254          CASE ( 'pt_reference' )
     255             READ ( 13 )  pt_reference
    254256          CASE ( 'pt_surface' )
    255257             READ ( 13 )  pt_surface
  • palm/trunk/SOURCE/user_interface.f90

    r48 r57  
    88! routine user_statistics now has one argument (sr),
    99! sample for generating time series quantities added
     10! Bugfix in sample for reading user defined data from restart file (user_init)
    1011!
    1112! Former revisions:
     
    174175!
    175176!          END SELECT
     177!
     178!          READ ( 13 )  field_chr
     179!
    176180!       ENDDO
    177181!    ENDIF
  • palm/trunk/SOURCE/write_var_list.f90

    r39 r57  
    44! Actual revisions:
    55! -----------------
    6 !
     6! +pt_refrence
    77!
    88! Former revisions:
     
    218218    WRITE ( 14 )  'pt_init                       '
    219219    WRITE ( 14 )  pt_init
     220    WRITE ( 14 )  'pt_reference                  '
     221    WRITE ( 14 )  pt_reference
    220222    WRITE ( 14 )  'pt_surface                    '
    221223    WRITE ( 14 )  pt_surface
Note: See TracChangeset for help on using the changeset viewer.