Changeset 420


Ignore:
Timestamp:
Jan 13, 2010 3:10:53 PM (14 years ago)
Author:
franke
Message:

collision of cloud droplets has changed in advec_particles
bugfixes in advec_particles and collision_efficiency
change in disturb_field to make runs reproducible on HLRN

Location:
palm/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/DOC/app/chapter_3.8.html

    r349 r420  
    128128<P STYLE="line-height: 100%">It's also possible to perform precursor
    129129runs (one atmospheric and one oceanic) followed by a coupled restart
    130 run. In order to achieve this, the parameter <A HREF="chapter_4.1.html#coupling_start_time">coupling_time_start</A>
     130run. In order to achieve this, the parameter <A HREF="chapter_4.1.html#coupling_start_time">coupling_start_time</A>
    131131must be set according to the <A HREF="../misc/precursor_run_control.pdf">documentation.</A></P>
    132132<HR>
  • palm/trunk/SOURCE/CURRENT_MODIFICATIONS

    r412 r420  
    11New:
    22---
    3 humidity=.T. is now usable for runs with topography. wall_humidityflux and 
     3humidity=.T. is now usable for runs with topography. wall_humidityflux and
    44wall_scalarflux are the corresponding new parin arrays.
    55(check_parameters, init_3d_model, parin)
    66
    7 Large scale vertical motion (subsidence/ascent) can be applied to the 
    8 prognostic equation for the potential temperature. (check_parameters, header, 
     7Large scale vertical motion (subsidence/ascent) can be applied to the
     8prognostic equation for the potential temperature. (check_parameters, header,
    99Makefile, modules, parin, prognostic_equations, read_var_list, subsidence,
    1010write_var_list)
     
    1818var_ts is replaced by dots_max (modules,init_3d_model)
    1919
    20 init_3d_model, modules
     20Every cloud droplet has now an own weighting factor and can be deleted due to
     21collisions. Condensation and collision of cloud droplets are adjusted
     22accordingly. (advec_particles)
     23
     24Collision efficiency for large cloud droplets has changed according to table of
     25Rogers and Yau. (collision_efficiency)
     26
     27advec_particles, collision_efficiency, init_3d_model, modules
    2128
    2229
     
    2532Dimension of array stat in cascade change to prevent type problems with
    2633mpi2 libraries (poisfft_hybrid)
     34
     35A loop has been splitted to make runs reproducible on HLRN
     36systems. (disturb_field)
    2737
    2838Bugfix: exchange of ghost points for prho included (time_integration)
     
    3545messages from gfortran compiler (modules)
    3646
    37 calc_precipitation, modules, poisfft_hybrid, sum_up_3d_data, time_integration
     47Bugfix: calculation of cloud droplet velocity (advec_particles)
     48
     49Bugfix: transfer of particles at south/left edge (advec_particles)
     50
     51Bugfix: calculation of collision_efficiency (collision_efficiency)
     52
     53advec_particles, calc_precipitation, collision_efficiency, disturb_field, modules, poisfft_hybrid, sum_up_3d_data, time_integration
    3854
    3955
  • palm/trunk/SOURCE/advec_particles.f90

    r392 r420  
    55! -----------------
    66! TEST: PRINT statements on unit 9 (commented out)
     7! Own weighting factor for every cloud droplet is implemented and condensation
     8! and collision of cloud droplets are adjusted accordingly. +delta_v, -s_r3,
     9! -s_r4
     10! Initialization of variables for the (sub-) timestep is moved to the beginning
     11! in order to enable deletion of cloud droplets due to collision processes.
     12! Collision efficiency for large cloud droplets has changed according to table
     13! of Rogers and Yau. (collision_efficiency)
     14! Bugfix: calculation of cloud droplet velocity
     15! Bugfix: transfer of particles at south/left edge when new position
     16!         y>=(ny+0.5)*dy-1.e-12 or x>=(nx+0.5)*dx-1.e-12, very rare
     17! Bugfix: calculation of y (collision_efficiency)
    718!
    819! Former revisions:
     
    117128    LOGICAL ::  dt_3d_reached, dt_3d_reached_l, prt_position
    118129
    119     REAL    ::  aa, arg, bb, cc, dd, delta_r, dens_ratio, de_dt, de_dt_min,    &
    120                 de_dx_int, de_dx_int_l, de_dx_int_u, de_dy_int, de_dy_int_l,   &
    121                 de_dy_int_u, de_dz_int, de_dz_int_l, de_dz_int_u, diss_int,    &
    122                 diss_int_l, diss_int_u, distance, dt_gap, dt_particle,         &
    123                 dt_particle_m, d_radius, d_sum, e_a, e_int, e_int_l, e_int_u,  &
    124                 e_mean_int, e_s, exp_arg, exp_term, fs_int, gg,                &
    125                 lagr_timescale, mean_r, new_r, p_int, pt_int, pt_int_l,        &
     130    REAL    ::  aa, arg, bb, cc, dd, delta_r, delta_v, dens_ratio, de_dt,      &
     131                de_dt_min, de_dx_int, de_dx_int_l, de_dx_int_u, de_dy_int,     &
     132                de_dy_int_l, de_dy_int_u, de_dz_int, de_dz_int_l, de_dz_int_u, &
     133                diss_int, diss_int_l, diss_int_u, distance, dt_gap,            &
     134                dt_particle, dt_particle_m, d_radius, d_sum, e_a, e_int,       &
     135                e_int_l, e_int_u, e_mean_int, e_s, exp_arg, exp_term, fs_int,  &
     136                gg,lagr_timescale, mean_r, new_r, p_int, pt_int, pt_int_l,     &
    126137                pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, ql_int_u, &
    127                 random_gauss, sl_r3, sl_r4, s_r3, s_r4, t_int, u_int, u_int_l, &
    128                 u_int_u, vv_int, v_int, v_int_l, v_int_u, w_int, w_int_l,      &
    129                 w_int_u, x, y
     138                random_gauss, sl_r3, sl_r4, t_int, u_int, u_int_l, u_int_u,    &
     139                vv_int, v_int, v_int_l, v_int_u, w_int, w_int_l, w_int_u,      &
     140                x, y
    130141
    131142    REAL, DIMENSION(1:30) ::  de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
     
    180191       CALL cpu_log( log_point_s(40), 'advec_part_io', 'stop' )
    181192    ENDIF
     193
     194!
     195!--    Initialize variables for the (sub-) timestep, i.e. for
     196!--    marking those particles to be deleted after the timestep
     197       particle_mask     = .TRUE.
     198       deleted_particles = 0
     199       trlp_count_recv   = 0
     200       trnp_count_recv   = 0
     201       trrp_count_recv   = 0
     202       trsp_count_recv   = 0
     203       trlpt_count_recv  = 0
     204       trnpt_count_recv  = 0
     205       trrpt_count_recv  = 0
     206       trspt_count_recv  = 0
     207       IF ( use_particle_tails )  THEN
     208          tail_mask     = .TRUE.
     209       ENDIF
     210       deleted_tails = 0
    182211
    183212!
     
    215244!--    Reset summation arrays
    216245       ql_c = 0.0;  ql_v = 0.0;  ql_vp = 0.0
     246       delta_r=0.0;  delta_v=0.0
    217247
    218248!
     
    343373!--       Sum up the total volume of liquid water (needed below for
    344374!--       re-calculating the weighting factors)
    345           ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor * &
     375          ql_v(k,j,i) = ql_v(k,j,i) + particles(n)%weight_factor *             &
    346376                                      particles(n)%radius**3
    347377       ENDDO
     
    351381!--    Particle growth by collision
    352382       CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' )
    353        
     383
    354384       DO  i = nxl, nxr
    355385          DO  j = nys, nyn
     
    374404                         tmp_particle = particles(psi+is)
    375405                         js = is
    376                          DO WHILE ( particles(psi+js-inc)%radius > &
     406                         DO WHILE ( particles(psi+js-inc)%radius >             &
    377407                                    tmp_particle%radius )
    378408                            particles(psi+js) = particles(psi+js-inc)
     
    386416!
    387417!--                Calculate the mean radius of all those particles which
    388 !--                are of smaller or equal size than the current particle
     418!--                are of smaller size than the current particle
    389419!--                and use this radius for calculating the collision efficiency
    390420                   psi = prt_start_index(k,j,i)
    391                    s_r3 = 0.0
    392                    s_r4 = 0.0
    393                    DO  n = psi, psi+prt_count(k,j,i)-1
    394 !
    395 !--                   There may be some particles of size equal to the
    396 !--                   current particle but with larger index
     421
     422                   DO  n = psi+prt_count(k,j,i)-1, psi+1, -1
    397423                      sl_r3 = 0.0
    398424                      sl_r4 = 0.0
    399                       DO is = n, psi+prt_count(k,j,i)-2
    400                          IF ( particles(is+1)%radius == &
    401                               particles(is)%radius ) THEN
    402                             sl_r3 = sl_r3 + particles(is+1)%radius**3
    403                             sl_r4 = sl_r4 + particles(is+1)%radius**4
    404                          ELSE
    405                             EXIT
    406                          ENDIF
    407                       ENDDO
    408 
    409                       IF ( ( s_r3 + sl_r3 ) > 0.0 )  THEN
    410 
    411                          mean_r = ( s_r4 + sl_r4 ) / ( s_r3 + sl_r3 )
    412 
    413                          CALL collision_efficiency( mean_r,              &
    414                                                     particles(n)%radius, &
     425
     426                      DO is = n-1, psi, -1
     427                         IF ( particles(is)%radius < particles(n)%radius ) THEN
     428                            sl_r3 = sl_r3 + particles(is)%weight_factor        &
     429                                          *( particles(is)%radius**3 )
     430                            sl_r4 = sl_r4 + particles(is)%weight_factor        &
     431                                          *( particles(is)%radius**4 )
     432                         ENDIF
     433                      ENDDO
     434
     435                      IF ( ( sl_r3 ) > 0.0 )  THEN
     436                         mean_r = ( sl_r4 ) / ( sl_r3 )
     437
     438                         CALL collision_efficiency( mean_r,                    &
     439                                                    particles(n)%radius,       &
    415440                                                    effective_coll_efficiency )
    416441
     
    418443                         effective_coll_efficiency = 0.0
    419444                      ENDIF
    420        
    421 !
    422 !--                   Contribution of the current particle to the next one
    423                       s_r3 = s_r3 + particles(n)%radius**3
    424                       s_r4 = s_r4 + particles(n)%radius**4
    425 
    426                       IF ( effective_coll_efficiency > 1.0  .OR. &
    427                            effective_coll_efficiency < 0.0 )     &
     445
     446                      IF ( effective_coll_efficiency > 1.0  .OR.               &
     447                           effective_coll_efficiency < 0.0 )                   &
    428448                      THEN
    429449                         WRITE( message_string, * )  'collision_efficiency ' , &
     
    431451                         CALL message( 'advec_particles', 'PA0145', 2, 2, -1,  &
    432452                                       6, 1 )
    433                       END IF
     453                      ENDIF
    434454
    435455!
     
    459479                                 ) / ( 3.0 * gg )
    460480
    461                       ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz * &
     481                      ql_int = ql_int_l + ( particles(n)%z - zu(kk) ) / dz *   &
    462482                                          ( ql_int_u - ql_int_l )
    463483
     
    478498                      gg = aa + bb + cc + dd
    479499
    480                       u_int_l = ( ( gg-aa ) * u(kk,jj,ii)   + ( gg-bb ) *     &
    481                                                               u(kk,jj,ii+1)   &
    482                                 + ( gg-cc ) * u(kk,jj+1,ii) + ( gg-dd ) *     &
    483                                                               u(kk,jj+1,ii+1) &
     500                      u_int_l = ( ( gg-aa ) * u(kk,jj,ii)   + ( gg-bb ) *      &
     501                                                              u(kk,jj,ii+1)    &
     502                                + ( gg-cc ) * u(kk,jj+1,ii) + ( gg-dd ) *      &
     503                                                              u(kk,jj+1,ii+1)  &
    484504                                ) / ( 3.0 * gg ) - u_gtrans
    485505                      IF ( kk+1 == nzt+1 )  THEN
     
    491511                                                             u(kk+1,jj+1,ii+1) &
    492512                                   ) / ( 3.0 * gg ) - u_gtrans
    493                          u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz * &
     513                         u_int = u_int_l + ( particles(n)%z - zu(kk) ) / dz *  &
    494514                                           ( u_int_u - u_int_l )
    495515                      ENDIF
     
    558578
    559579!
    560 !--                   Change in radius due to collision
    561                       delta_r = effective_coll_efficiency * &
    562                                 ql_int * rho_surface / ( 1.0 - ql_int ) *   &
    563                                 0.25 / rho_l *                              &
    564                                 SQRT( ( u_int - particles(n)%speed_x )**2 + &
    565                                       ( v_int - particles(n)%speed_y )**2 + &
    566                                       ( w_int - particles(n)%speed_z )**2   &
    567                                     ) * dt_3d
     580!--                Change in radius due to collision
     581                      delta_r = effective_coll_efficiency / 3.0                &
     582                                * pi * sl_r3 * ddx * ddy / dz                  &
     583                                * SQRT( ( u_int - particles(n)%speed_x )**2    &
     584                                      + ( v_int - particles(n)%speed_y )**2    &
     585                                      + ( w_int - particles(n)%speed_z )**2    &
     586                                      ) * dt_3d
     587!
     588!--                Change in volume due to collision
     589                      delta_v = particles(n)%weight_factor                     &
     590                                * ( ( particles(n)%radius + delta_r )**3       &
     591                                    - particles(n)%radius**3 )
     592
     593!
     594!--                Check if collected particles provide enough LWC for volume
     595!--                change of collector particle
     596                      IF ( delta_v >= sl_r3 .and. sl_r3 > 0.0 ) THEN
     597
     598                         delta_r = ( ( sl_r3/particles(n)%weight_factor )      &
     599                                     + particles(n)%radius**3 )**( 1./3. )     &
     600                                     - particles(n)%radius
     601
     602                         DO is = n-1, psi, -1
     603                            IF ( particles(is)%radius < particles(n)%radius )  &
     604                            THEN
     605                                 particles(is)%weight_factor = 0.0
     606                                 particle_mask(is)  = .FALSE.
     607                                 deleted_particles = deleted_particles + 1
     608                            ENDIF
     609                         ENDDO
     610
     611                      ELSE IF ( delta_v < sl_r3 .and. sl_r3 > 0.0 ) THEN
     612
     613                         DO is = n-1, psi, -1
     614                            IF ( particles(is)%radius < particles(n)%radius    &
     615                                 .and. sl_r3 > 0.0 ) THEN
     616                                 particles(is)%weight_factor =                 &
     617                                               ( ( particles(is)%weight_factor &
     618                                               * ( particles(is)%radius**3 ) ) &
     619                                               - ( delta_v                     &
     620                                               * particles(is)%weight_factor   &
     621                                               * ( particles(is)%radius**3 )   &
     622                                               / sl_r3 ) )                     &
     623                                               / ( particles(is)%radius**3 )
     624
     625                               IF (particles(is)%weight_factor < 0.0) THEN
     626                                  WRITE( message_string, * ) 'negative ',      &
     627                                                     'weighting factor: ',     &
     628                                                     particles(is)%weight_factor
     629                                  CALL message( 'advec_particles', '', 2, 2,   &
     630                                                -1, 6, 1 )
     631                               ENDIF
     632                            ENDIF
     633                         ENDDO
     634                      ENDIF
    568635
    569636                      particles(n)%radius = particles(n)%radius + delta_r
    570 
    571                       ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%radius**3
    572 
     637                      ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%weight_factor &
     638                                      *( particles(n)%radius**3 )
    573639                   ENDDO
    574640
    575                 ENDIF
    576 
    577 !
    578 !--             Re-calculate the weighting factor (total liquid water content
    579 !--             must be conserved during collision)
    580                 IF ( ql_vp(k,j,i) /= 0.0 )  THEN
    581 
    582                    ql_vp(k,j,i) = ql_v(k,j,i) / ql_vp(k,j,i)
    583 !
    584 !--                Re-assign this weighting factor to the particles of the
    585 !--                current gridbox
     641                   ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
     642                                  *( particles(psi)%radius**3 )
     643
     644                ELSE IF ( prt_count(k,j,i) == 1 )  THEN
    586645                   psi = prt_start_index(k,j,i)
    587                    DO  n = psi, psi + prt_count(k,j,i)-1
    588                       particles(n)%weight_factor = ql_vp(k,j,i)
    589                    ENDDO
    590 
     646                   ql_vp(k,j,i) =  particles(psi)%weight_factor                &
     647                                  *( particles(psi)%radius**3 )
     648                ENDIF
     649!
     650!--            Check if condensation of LWC was conserved
     651!--            during collision process
     652                IF (ql_v(k,j,i) /= 0.0 ) THEN
     653                   IF (ql_vp(k,j,i)/ql_v(k,j,i) >= 1.00001 .or.                &
     654                       ql_vp(k,j,i)/ql_v(k,j,i) <= 0.99999  ) THEN
     655                      WRITE( message_string, * ) 'LWC is not conserved during',&
     656                                                 ' collision! ',               &
     657                                                 'LWC after condensation: ',   &
     658                                                 ql_v(k,j,i),                  &
     659                                                 ' LWC after collision: ',     &
     660                                                 ql_vp(k,j,i)
     661                      CALL message( 'advec_particles', '', 2, 2, -1, 6, 1 )
     662                   ENDIF
    591663                ENDIF
    592664
     
    726798!
    727799!--       Final values are obtained by division by the total number of grid
    728 !--       points used for the summation. 
     800!--       points used for the summation.
    729801          hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0)   ! u
    730802          hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0)   ! v
     
    761833          CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, &
    762834                              MPI_REAL, MPI_SUM, comm2d, ierr )
    763                  
     835
    764836#else
    765837          sums(:,8)  = sums_l(:,8,0)
     
    779851       ENDIF
    780852
    781     ENDIF 
     853    ENDIF
    782854
    783855!
     
    815887       dt_3d_reached_l = .TRUE.
    816888
    817 !
    818 !--    Initialize variables for the (sub-) timestep, i.e. for marking those
    819 !--    particles to be deleted after the timestep
    820        particle_mask     = .TRUE.
    821        deleted_particles = 0
    822        trlp_count_recv   = 0
    823        trnp_count_recv   = 0
    824        trrp_count_recv   = 0
    825        trsp_count_recv   = 0
    826        trlpt_count_recv  = 0
    827        trnpt_count_recv  = 0
    828        trrpt_count_recv  = 0
    829        trspt_count_recv  = 0
    830        IF ( use_particle_tails )  THEN
    831           tail_mask     = .TRUE.
    832        ENDIF
    833        deleted_tails = 0
    834 
    835 
    836889       DO  n = 1, number_of_particles
    837890!
     
    847900          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
    848901              + offset_ocean_nzt                     ! only exact if equidistant
    849                                                
     902
    850903!
    851904!--       Interpolation of the velocity components in the xy-plane
     
    942995                 + offset_ocean_nzt                      ! only exact if eq.dist
    943996
    944              IF ( topography == 'flat' )  THEN        
     997             IF ( topography == 'flat' )  THEN
    945998
    946999                x  = particles(n)%x - i * dx
     
    10661119!--             gp_outside_of_building(7): i+1,j,k+1
    10671120!--             gp_outside_of_building(8): i+1,j+1,k+1
    1068            
     1121
    10691122                gp_outside_of_building = 0
    10701123                location = 0.0
     
    10921145                   location(num_gp,3) = k * dz - 0.5 * dz
    10931146                   ei(num_gp)     = e(k,j+1,i)
    1094                    dissi(num_gp)  = diss(k,j+1,i)       
     1147                   dissi(num_gp)  = diss(k,j+1,i)
    10951148                   de_dxi(num_gp) = de_dx(k,j+1,i)
    10961149                   de_dyi(num_gp) = de_dy(k,j+1,i)
     
    11051158                   location(num_gp,3) = (k+1) * dz - 0.5 * dz
    11061159                   ei(num_gp)     = e(k+1,j,i)
    1107                    dissi(num_gp)  = diss(k+1,j,i)       
     1160                   dissi(num_gp)  = diss(k+1,j,i)
    11081161                   de_dxi(num_gp) = de_dx(k+1,j,i)
    11091162                   de_dyi(num_gp) = de_dy(k+1,j,i)
     
    11191172                   location(num_gp,3) = (k+1) * dz - 0.5 * dz
    11201173                   ei(num_gp)     = e(k+1,j+1,i)
    1121                    dissi(num_gp)  = diss(k+1,j+1,i)       
     1174                   dissi(num_gp)  = diss(k+1,j+1,i)
    11221175                   de_dxi(num_gp) = de_dx(k+1,j+1,i)
    11231176                   de_dyi(num_gp) = de_dy(k+1,j+1,i)
     
    11331186                   location(num_gp,3) = k * dz - 0.5 * dz
    11341187                   ei(num_gp)     = e(k,j,i+1)
    1135                    dissi(num_gp)  = diss(k,j,i+1) 
     1188                   dissi(num_gp)  = diss(k,j,i+1)
    11361189                   de_dxi(num_gp) = de_dx(k,j,i+1)
    11371190                   de_dyi(num_gp) = de_dy(k,j,i+1)
     
    11471200                   location(num_gp,3) = k * dz - 0.5 * dz
    11481201                   ei(num_gp)     = e(k,j+1,i+1)
    1149                    dissi(num_gp)  = diss(k,j+1,i+1) 
     1202                   dissi(num_gp)  = diss(k,j+1,i+1)
    11501203                   de_dxi(num_gp) = de_dx(k,j+1,i+1)
    11511204                   de_dyi(num_gp) = de_dy(k,j+1,i+1)
     
    11611214                   location(num_gp,3) = (k+1) * dz - 0.5 * dz
    11621215                   ei(num_gp)     = e(k+1,j,i+1)
    1163                    dissi(num_gp)  = diss(k+1,j,i+1) 
     1216                   dissi(num_gp)  = diss(k+1,j,i+1)
    11641217                   de_dxi(num_gp) = de_dx(k+1,j,i+1)
    11651218                   de_dyi(num_gp) = de_dy(k+1,j,i+1)
     
    11751228                   location(num_gp,3) = (k+1) * dz - 0.5 * dz
    11761229                   ei(num_gp)     = e(k+1,j+1,i+1)
    1177                    dissi(num_gp)  = diss(k+1,j+1,i+1) 
     1230                   dissi(num_gp)  = diss(k+1,j+1,i+1)
    11781231                   de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
    11791232                   de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
     
    11821235
    11831236!
    1184 !--             If all gridpoints are situated outside of a building, then the 
     1237!--             If all gridpoints are situated outside of a building, then the
    11851238!--             ordinary interpolation scheme can be used.
    11861239                IF ( num_gp == 8 )  THEN
     
    13091362                      de_dxi(num_gp) = 0.0
    13101363                      de_dyi(num_gp) = de_dy(k,j,i)
    1311                       de_dzi(num_gp) = de_dz(k,j,i)                     
    1312                    ENDIF 
    1313    
     1364                      de_dzi(num_gp) = de_dz(k,j,i)
     1365                   ENDIF
     1366
    13141367                   IF ( gp_outside_of_building(5) == 1  .AND. &
    13151368                      gp_outside_of_building(1) == 0 )  THEN
     
    13671420                      de_dxi(num_gp) = 0.0
    13681421                      de_dyi(num_gp) = de_dy(k,j+1,i)
    1369                       de_dzi(num_gp) = de_dz(k,j+1,i)                     
    1370                    ENDIF 
    1371    
     1422                      de_dzi(num_gp) = de_dz(k,j+1,i)
     1423                   ENDIF
     1424
    13721425                   IF ( gp_outside_of_building(6) == 1  .AND. &
    13731426                        gp_outside_of_building(2) == 0 )  THEN
     
    14251478                      de_dxi(num_gp) = 0.0
    14261479                      de_dyi(num_gp) = de_dy(k+1,j,i)
    1427                       de_dzi(num_gp) = de_dz(k+1,j,i)                    
    1428                    ENDIF 
    1429    
     1480                      de_dzi(num_gp) = de_dz(k+1,j,i)
     1481                   ENDIF
     1482
    14301483                   IF ( gp_outside_of_building(7) == 1  .AND. &
    14311484                        gp_outside_of_building(3) == 0 )  THEN
     
    14831536                      de_dxi(num_gp) = 0.0
    14841537                      de_dyi(num_gp) = de_dy(k+1,j+1,i)
    1485                       de_dzi(num_gp) = de_dz(k+1,j+1,i)                     
    1486                    ENDIF 
    1487    
     1538                      de_dzi(num_gp) = de_dz(k+1,j+1,i)
     1539                   ENDIF
     1540
    14881541                   IF ( gp_outside_of_building(8) == 1  .AND. &
    14891542                        gp_outside_of_building(4) == 0 )  THEN
     
    16001653                   IF ( num_gp == 1 )  THEN
    16011654!
    1602 !--                   If only one of the gridpoints is situated outside of the 
    1603 !--                   building, it follows that the values at the particle 
     1655!--                   If only one of the gridpoints is situated outside of the
     1656!--                   building, it follows that the values at the particle
    16041657!--                   location are the same as the gridpoint values
    16051658                      e_int     = ei(num_gp)
     
    16121665                      d_sum = 0.0
    16131666!
    1614 !--                   Evaluation of the distances between the gridpoints 
    1615 !--                   contributing to the interpolated values, and the particle 
     1667!--                   Evaluation of the distances between the gridpoints
     1668!--                   contributing to the interpolated values, and the particle
    16161669!--                   location
    16171670                      DO agp = 1, num_gp
     
    18391892             IF ( cloud_droplets )  THEN
    18401893                exp_arg  = 4.5 * dens_ratio * molecular_viscosity /        &
    1841                            ( particles(n)%radius )**2 /                    &
     1894                           ( particles(n)%radius )**2 *                    &
    18421895                           ( 1.0 + 0.15 * ( 2.0 * particles(n)%radius *    &
    18431896                             SQRT( ( u_int - particles(n)%speed_x )**2 +   &
     
    21782231                      deleted_particles = deleted_particles + 1
    21792232
     2233                      IF ( trlp(trlp_count)%x >= (nx + 0.5)* dx - 1.e-12 ) THEN
     2234                         trlp(trlp_count)%x = trlp(trlp_count)%x - 1.e-10
     2235                         trlp(trlp_count)%origin_x = trlp(trlp_count)%origin_x &
     2236                                                     - 1
     2237                      ENDIF
     2238
    21802239                      IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    21812240                         trlpt_count = trlpt_count + 1
     
    26172676                         deleted_particles = deleted_particles + 1
    26182677
     2678                         IF ( trsp(trsp_count)%y >= (ny+0.5)* dy - 1.e-12 ) THEN
     2679                            trsp(trsp_count)%y = trsp(trsp_count)%y - 1.e-10
     2680                            trsp(trsp_count)%origin_y =                        &
     2681                                                    trsp(trsp_count)%origin_y - 1
     2682                         ENDIF
     2683
    26192684                         IF ( use_particle_tails  .AND.  nn /= 0 )  THEN
    26202685                            trspt_count = trspt_count + 1
     
    32963361       IF ( dt_3d_reached )  EXIT
    32973362
     3363!
     3364!--    Initialize variables for the next (sub-) timestep, i.e. for marking those
     3365!--    particles to be deleted after the timestep
     3366       particle_mask     = .TRUE.
     3367       deleted_particles = 0
     3368       trlp_count_recv   = 0
     3369       trnp_count_recv   = 0
     3370       trrp_count_recv   = 0
     3371       trsp_count_recv   = 0
     3372       trlpt_count_recv  = 0
     3373       trnpt_count_recv  = 0
     3374       trrpt_count_recv  = 0
     3375       trspt_count_recv  = 0
     3376       IF ( use_particle_tails )  THEN
     3377          tail_mask     = .TRUE.
     3378       ENDIF
     3379       deleted_tails = 0
     3380
    32983381    ENDDO   ! timestep loop
    3299 
    33003382
    33013383!
     
    33083390    ENDIF
    33093391
    3310 
    3311 !
    3312 !-- Re-evaluate the weighting factors. After advection, particles within a
    3313 !-- grid box may have different weighting factors if some have been advected
    3314 !-- from a neighbouring box. The factors are re-evaluated so that they are
    3315 !-- the same for all particles of one box. This procedure must conserve the
    3316 !-- liquid water content within one box.
    33173392    IF ( cloud_droplets )  THEN
    33183393
     
    33223397
    33233398!
    3324 !--    Re-calculate the weighting factors and calculate the liquid water content
     3399!--    Calculate the liquid water content
    33253400       DO  i = nxl, nxr
    33263401          DO  j = nys, nyn
     
    33283403
    33293404!
    3330 !--             Calculate the total volume of particles in the boxes (ql_vp) as
    3331 !--             well as the real volume (ql_v, weighting factor has to be
    3332 !--             included)
     3405!--             Calculate the total volume in the boxes (ql_v, weighting factor
     3406!--             has to beincluded)
    33333407                psi = prt_start_index(k,j,i)
    33343408                DO  n = psi, psi+prt_count(k,j,i)-1
    3335                    ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%radius**3
    3336 
    3337                    ql_v(k,j,i)  = ql_v(k,j,i)  + particles(n)%weight_factor * &
     3409                   ql_v(k,j,i)  = ql_v(k,j,i)  + particles(n)%weight_factor *  &
    33383410                                                 particles(n)%radius**3
    33393411                ENDDO
    33403412
    33413413!
    3342 !--             Re-calculate the weighting factors and calculate the liquid
    3343 !--             water content
    3344                 IF ( ql_vp(k,j,i) /= 0.0 )  THEN
    3345                    ql_vp(k,j,i) = ql_v(k,j,i) / ql_vp(k,j,i)
    3346                    ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333 * pi * &
    3347                                            ql_v(k,j,i) /             &
     3414!--             Calculate the liquid water content
     3415                IF ( ql_v(k,j,i) /= 0.0 )  THEN
     3416                   ql(k,j,i) = ql(k,j,i) + rho_l * 1.33333333 * pi *           &
     3417                                           ql_v(k,j,i) /                       &
    33483418                                           ( rho_surface * dx * dy * dz )
     3419
     3420                   IF ( ql(k,j,i) < 0.0 ) THEN
     3421                      WRITE( message_string, * )  'LWC out of range: ' , &
     3422                                                  ql(k,j,i)
     3423                         CALL message( 'advec_particles', '', 2, 2, -1, 6, 1 )
     3424                   ENDIF
     3425
    33493426                ELSE
    33503427                   ql(k,j,i) = 0.0
    33513428                ENDIF
    3352 
    3353 !
    3354 !--             Re-assign the weighting factor to the particles
    3355                 DO  n = psi, psi+prt_count(k,j,i)-1
    3356                    particles(n)%weight_factor = ql_vp(k,j,i)
    3357                 ENDDO
    33583429
    33593430             ENDDO
     
    38333904       e = 0.001
    38343905    ELSEIF ( mean_rm >= 25.0 )  THEN
    3835        IF( j <= 3 )  e = 0.55
     3906       IF( j <= 2 )  e = 0.0
     3907       IF( j == 3 )  e = 0.47
    38363908       IF( j == 4 )  e = 0.8
    38373909       IF( j == 5 )  e = 0.9
    38383910       IF( j >=6  )  e = 1.0
    38393911    ELSEIF ( rm >= 3000.0 )  THEN
    3840        e = 1.0
     3912       IF( i == 1 )  e = 0.02
     3913       IF( i == 2 )  e = 0.16
     3914       IF( i == 3 )  e = 0.33
     3915       IF( i == 4 )  e = 0.55
     3916       IF( i == 5 )  e = 0.71
     3917       IF( i == 6 )  e = 0.81
     3918       IF( i == 7 )  e = 0.90
     3919       IF( i >= 8 )  e = 0.94
    38413920    ELSE
    38423921       x  = mean_rm - collected_r(i)
    3843        y  = rm - collected_r(j)
    3844        dx = collected_r(i+1) - collected_r(i) 
    3845        dy = collector_r(j+1) - collector_r(j) 
     3922       y  = rm - collector_r(j)
     3923       dx = collected_r(i+1) - collected_r(i)
     3924       dy = collector_r(j+1) - collector_r(j)
    38463925       aa = x**2 + y**2
    38473926       bb = ( dx - x )**2 + y**2
     
    38543933    ENDIF
    38553934
    3856  END SUBROUTINE collision_efficiency 
     3935 END SUBROUTINE collision_efficiency
    38573936
    38583937
     
    39023981          WRITE( message_string, * ) ' particle out of range: i=', i, ' j=', &
    39033982                          j, ' k=', k,                                       &
    3904                           ' &nxl=', nxl, ' nxr=', nxr,                       &
     3983                          ' nxl=', nxl, ' nxr=', nxr,                        &
    39053984                          ' nys=', nys, ' nyn=', nyn,                        &
    39063985                          ' nzb=', nzb, ' nzt=', nzt
  • palm/trunk/SOURCE/disturb_field.f90

    r77 r420  
    44! Actual revisions:
    55! -----------------
    6 !
     6! A loop has been splitted to make runs reproducible on HLRN systems
    77!
    88! Former revisions:
     
    9898!
    9999!-- Exchange of ghost points for the random perturbation
     100
    100101    CALL exchange_horiz( dist1 )
    101102
     
    104105!-- Neighboured grid points in all three directions are used for the
    105106!-- filter operation.
    106     DO  i = nxl, nxr
    107        DO  j = nys, nyn
     107!-- Loop has been splitted to make runs reproducible on HLRN systems using
     108!-- compiler option -O3
     109     DO  i = nxl, nxr
     110        DO  j = nys, nyn
    108111          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
    109              dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1) + dist1(k,j-1,i) &
    110                             + dist1(k,j+1,i) + dist1(k+1,j,i) + dist1(k-1,j,i) &
     112             dist2(k,j,i) = ( dist1(k,j,i-1) + dist1(k,j,i+1) &
     113                            + dist1(k,j+1,i) + dist1(k+1,j,i) &
     114                            ) / 12.0
     115          ENDDO
     116          DO  k = disturbance_level_ind_b-1, disturbance_level_ind_t+1
     117              dist2(k,j,i) = dist2(k,j,i) + ( dist1(k,j-1,i) + dist1(k-1,j,i)  &
    111118                            + 6.0 * dist1(k,j,i)                               &
    112119                            ) / 12.0
    113120          ENDDO
    114        ENDDO
    115     ENDDO
     121        ENDDO
     122     ENDDO
    116123
    117124!
     
    119126!-- Afterwards, filter operation and exchange of ghost points are repeated.
    120127    CALL exchange_horiz( dist2 )
     128
    121129    DO  i = nxl, nxr
    122130       DO  j = nys, nyn
     
    129137       ENDDO
    130138    ENDDO
     139
    131140    CALL exchange_horiz( dist1 )
    132141
Note: See TracChangeset for help on using the changeset viewer.