Home


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

Changeset 792


Ignore:
Timestamp:
Dec 1, 2011 12:23:23 AM (12 years ago)
Author:
raasch
Message:

further adjustments for speedup of particle code

Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/advec_particles.f90

    r760 r792  
    44! Current revisions:
    55! ------------------
    6 !
     6! particle arrays (particles, particles_temp) implemented as pointers in
     7! order to speed up sorting (see routine sort_particles)
    78!
    89! Former revisions:
    … …  
    39953996    INTEGER ::  i, ilow, j, k, n
    39963997
    3997     TYPE(particle_type), DIMENSION(1:number_of_particles) ::  particles_temp
     3998    TYPE(particle_type), DIMENSION(:), POINTER ::  particles_temp
    39983999
    39994000
    … …  
    40014002
    40024003!
    4003 !-- Initialize the array used for counting and indexing the particles
    4004     prt_count = 0
     4004!-- Initialize counters and set pointer of the temporary array into which
     4005!-- particles are sorted to free memory
     4006    prt_count  = 0
     4007    sort_count = sort_count +1
     4008
     4009    SELECT CASE ( MOD( sort_count, 2 ) )
     4010
     4011       CASE ( 0 )
     4012
     4013          particles_temp => part_1
     4014          part_1 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
     4015                                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
     4016                                  0.0, 0, 0, 0, 0 )
     4017
     4018       CASE ( 1 )
     4019
     4020          particles_temp => part_2
     4021          part_2 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
     4022                                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
     4023                                  0.0, 0, 0, 0, 0 )
     4024
     4025    END SELECT
    40054026
    40064027!
    … …  
    40554076    ENDDO
    40564077
    4057     particles(1:number_of_particles) = particles_temp
     4078!
     4079!-- Redirect the particles pointer to the sorted array
     4080    SELECT CASE ( MOD( sort_count, 2 ) )
     4081
     4082       CASE ( 0 )
     4083
     4084          particles => part_1
     4085
     4086       CASE ( 1 )
     4087
     4088          particles => part_2
     4089
     4090    END SELECT
    40584091
    40594092!

  • palm/trunk/SOURCE/init_particles.f90

    r668 r792  
    44! Current revisions:
    55! -----------------
     6! array particles implemented as pointer
    67!
    78! Former revisions:
    … …  
    190191       ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
    191192                 prt_start_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    192                  particle_mask(maximum_number_of_particles),         &
    193                  particles(maximum_number_of_particles) )
     193                 particle_mask(maximum_number_of_particles),     &
     194                 part_1(maximum_number_of_particles),            &
     195                 part_2(maximum_number_of_particles) )
     196
     197       part_1 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
     198                               0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
     199                               0.0, 0, 0, 0, 0 )
     200
     201       part_2 = particle_type( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
     202                               0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
     203                               0.0, 0, 0, 0, 0 )
     204
     205       sort_count = 0
     206
     207       particles => part_1
    194208
    195209       READ ( 90 )  prt_count, prt_start_index
    … …  
    214228       ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg),       &
    215229                 prt_start_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    216                  particle_mask(maximum_number_of_particles),         &
    217                  particles(maximum_number_of_particles) )
     230                 particle_mask(maximum_number_of_particles),     &
     231                 part_1(maximum_number_of_particles),            &
     232                 part_2(maximum_number_of_particles)  )
     233
     234       particles => part_1
     235
     236       sort_count = 0
    218237
    219238!

  • palm/trunk/SOURCE/modules.f90

    r791 r792  
    55! Current revisions:
    66! -----------------
    7 !
     7! particle arrays (particles, parrticles_temp) implemented as pointers,
     8! +particles_1, particles_2, sort_count
    89!
    910! Former revisions:
    … …  
    11661167                offset_ocean_nzt_m1 = 0, particles_per_point = 1,              &
    11671168                particle_file_count = 0, skip_particles_for_tail = 100,        &
    1168                 total_number_of_particles, total_number_of_tails = 0
     1169                sort_count = 0, total_number_of_particles,                     &
     1170                total_number_of_tails = 0
    11691171
    11701172    INTEGER, PARAMETER ::  max_number_of_particle_groups = 10
    … …  
    12121214    END TYPE particle_type
    12131215
    1214     TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  initial_particles, &
    1215                                                        particles
     1216    TYPE(particle_type), DIMENSION(:), ALLOCATABLE ::  initial_particles
     1217    TYPE(particle_type), DIMENSION(:), ALLOCATABLE, TARGET ::  part_1, part_2
     1218    TYPE(particle_type), DIMENSION(:), POINTER  ::  particles
    12161219
    12171220    TYPE particle_groups_type

  • palm/trunk/SOURCE/wang_kernel.f90

    r791 r792  
    44! Current revisions:
    55! -----------------
    6 !
     6! speed optimizations + formatting
    77!
    88! Former revisions:
    … …  
    2929    PRIVATE
    3030
    31     PUBLIC TurbSD, turb_enhance_eff, effic, fallg, PHI, ZHI, colker
    32 
    33     INTEGER, SAVE :: ip, jp, kp, pstart, pend, psum
    34     REAL, SAVE :: epsilon, urms
    35     REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: gck, ec, ecf
    36     REAL, DIMENSION(:),   ALLOCATABLE, SAVE :: winf
    37     REAL, SAVE :: eps2, urms2
     31    PUBLIC  colker, effic, fallg, phi, turbsd, turb_enhance_eff, zhi
     32
     33    INTEGER, SAVE ::  ip, jp, kp, pend, pstart, psum
     34
     35    REAL, SAVE :: epsilon, eps2, urms, urms2
     36
     37    REAL, DIMENSION(:),   ALLOCATABLE, SAVE ::  winf
     38    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  ec, ecf, gck
    3839
    3940!
    4041!-- Public interfaces
    41     INTERFACE TurbSD
    42        MODULE PROCEDURE TurbSD
    43     END INTERFACE TurbSD
    44 
    45     INTERFACE PHI
    46        MODULE PROCEDURE PHI
    47     END INTERFACE PHI
    48 
    49     INTERFACE ZHI
    50        MODULE PROCEDURE ZHI
    51     END INTERFACE ZHI
     42    INTERFACE colker
     43       MODULE PROCEDURE colker
     44    END INTERFACE colker
     45
     46    INTERFACE effic
     47       MODULE PROCEDURE effic
     48    END INTERFACE effic
     49
     50    INTERFACE fallg
     51       MODULE PROCEDURE fallg
     52    END INTERFACE fallg
     53
     54    INTERFACE phi
     55       MODULE PROCEDURE phi
     56    END INTERFACE phi
     57
     58    INTERFACE turbsd
     59       MODULE PROCEDURE turbsd
     60    END INTERFACE turbsd
    5261
    5362    INTERFACE turb_enhance_eff
    … …  
    5564    END INTERFACE turb_enhance_eff
    5665
    57     INTERFACE effic
    58        MODULE PROCEDURE effic
    59     END INTERFACE effic
    60 
    61     INTERFACE fallg
    62        MODULE PROCEDURE fallg
    63     END INTERFACE fallg
    64 
    65     INTERFACE colker
    66        MODULE PROCEDURE colker
    67     END INTERFACE colker
     66    INTERFACE zhi
     67       MODULE PROCEDURE zhi
     68    END INTERFACE zhi
    6869
    6970    CONTAINS
    … …  
    7273! SUBROUTINE for calculation of collision kernel
    7374!------------------------------------------------------------------------------!
    74     SUBROUTINE colker(i1,j1,k1,kernel)
     75    SUBROUTINE colker( i1, j1, k1, kernel )
    7576
    7677       USE arrays_3d
    7778       USE cloud_parameters
    7879       USE constants
     80       USE cpulog
    7981       USE indices
     82       USE interfaces
    8083       USE particle_attributes
    8184
    8285       IMPLICIT NONE
    8386
    84        INTEGER :: i,i1,j,j1,k1
     87       INTEGER ::  i, i1, j, j1, k1
    8588
    8689       REAL, DIMENSION(prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1,prt_start_index(k1,j1,i1):prt_start_index(k1,j1,i1)+prt_count(k1,j1,i1)-1) :: kernel
     90
     91!       CALL cpu_log( log_point_s(46), 'colker', 'start' )
    8792
    8893       ip = i1
    … …  
    9499       psum   = prt_count(kp,jp,ip)
    95100
    96        ALLOCATE(winf(pstart:pend), ec(pstart:pend,pstart:pend))
    97 
    98        IF ( turbulence_effects_on_collision ) THEN
    99 
    100           ALLOCATE(gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend))
    101 
    102           epsilon = diss(kp,jp,ip) * 10000.0  !dissipation rate in cm**2/s**-3
    103           urms     = 202.0 * (epsilon/400.0)**(1.0/3.0)
    104 
    105           IF (epsilon <= 0.001) THEN
     101       ALLOCATE( ec(pstart:pend,pstart:pend), winf(pstart:pend) )
     102
     103       IF ( turbulence_effects_on_collision )  THEN
     104
     105          ALLOCATE( gck(pstart:pend,pstart:pend), ecf(pstart:pend,pstart:pend) )
     106
     107          epsilon = diss(kp,jp,ip) * 1.0E5  !dissipation rate in cm**2/s**-3
     108          urms     = 202.0 * ( epsilon/ 400.0 )**( 1.0 / 3.0 )
     109
     110          IF ( epsilon <= 0.001 ) THEN
    106111
    107112             CALL fallg
    108113             CALL effic
    109114
    110              DO j = pstart, pend, 1
    111                 DO i =  pstart, pend, 1
    112                    kernel(i,j) = pi * (particles(j)%radius*100.0 + particles(i)%radius*100.0)**2.0 * ec(i,j) * abs(winf(j)-winf(i))
     115             DO  j = pstart, pend
     116                DO  i =  pstart, pend
     117                   kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
     118                                        particles(i)%radius * 100.0 )**2 &
     119                                    * ec(i,j) * ABS( winf(j) - winf(i) )
    113120                ENDDO
    114121             ENDDO
    … …  
    116123          ELSE
    117124
    118              CALL TurbSD
     125             CALL turbsd
    119126             CALL turb_enhance_eff
    120127             CALL effic
    121128
    122              DO j = pstart, pend, 1
    123                 DO i =  pstart, pend, 1
     129             DO  j = pstart, pend, 1
     130                DO  i =  pstart, pend, 1
    124131                   kernel(i,j) = ec(i,j) * gck(i,j) * ecf(i,j)
    125132                ENDDO
    126133             ENDDO
     134
    127135          ENDIF
    128136
    … …  
    131139       ELSE
    132140
     141!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'start' )
    133142          CALL fallg
     143!          CALL cpu_log( log_point_s(50), 'colker_fallg', 'stop' )
     144!          CALL cpu_log( log_point_s(51), 'colker_effic', 'start' )
    134145          CALL effic
    135 
    136           DO j = pstart, pend, 1
    137              DO i =  pstart, pend, 1
    138                 kernel(i,j) = pi * (particles(j)%radius*100.0 + particles(i)%radius*100.0)**2.0 * ec(i,j) * abs(winf(j)-winf(i))
     146!          CALL cpu_log( log_point_s(51), 'colker_effic', 'stop' )
     147
     148          DO  j = pstart, pend
     149             DO  i =  pstart, pend
     150                kernel(i,j) = pi * ( particles(j)%radius * 100.0 +    &
     151                                     particles(i)%radius * 100.0 )**2 &
     152                                 * ec(i,j) * ABS( winf(j) - winf(i) )
    139153             ENDDO
    140154          ENDDO
    … …  
    142156       ENDIF
    143157
    144        DEALLOCATE(winf, ec)
    145 
    146        RETURN
     158       DEALLOCATE( ec, winf )
     159
     160!       CALL cpu_log( log_point_s(46), 'colker', 'stop' )
    147161
    148162    END SUBROUTINE colker
    … …  
    151165! SUBROUTINE for calculation of w, g and gck
    152166!------------------------------------------------------------------------------!
    153     SUBROUTINE TurbSD
    154 !from Aayala 2008b, page 37ff, necessary input parameter water density, radii of droplets, air density, air viscosity, turbulent dissipation rate, taylor microscale reynolds number, gravitational acceleration
     167    SUBROUTINE turbsd
     168! from Aayala 2008b, page 37ff, necessary input parameter water density, radii
     169! of droplets, air density, air viscosity, turbulent dissipation rate,
     170! taylor microscale reynolds number, gravitational acceleration
    155171       USE constants
    156172       USE cloud_parameters
    … …  
    301317       ENDDO
    302318
    303 
    304        RETURN
    305319    END SUBROUTINE TurbSD
    306320
    … …  
    317331
    318332       PHI = 1./aa1 - vsett/2.0/b/aa1**2.0  !in s
    319        RETURN
     333
    320334    END FUNCTION PHI
    321335
    … …  
    339353             + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2.0 + vsett2/aa1**2.0)   &
    340354                                                           * 1./2./b/aa3             ! in s**2
    341         RETURN
    342355    END FUNCTION ZHI
    343356
    … …  
    430443         ENDIF
    431444      ENDDO
    432       RETURN
     445
    433446   END SUBROUTINE fallg
    434447
    … …  
    439452   SUBROUTINE effic
    440453
    441        USE constants
    442        USE cloud_parameters
    443        USE particle_attributes
    444        USE arrays_3d
    445 
    446 !collision efficiencies of hall kernel
     454      USE arrays_3d
     455      USE constants
     456      USE cloud_parameters
     457      USE particle_attributes
     458
    447459      IMPLICIT NONE
    448460
    449461      INTEGER :: i, ir, iq, j, k, kk
    450       REAL    :: ek, rq, pp, qq
    451 
    452       REAL, DIMENSION(1:21) :: rat
    453       REAL, DIMENSION(1:15) :: r0
    454       REAL, DIMENSION(1:15,1:21) :: ecoll
    455 
    456       r0 = (/6.,8.,10.,15.,20.,25.,30.,40.,50., &
    457               60.,70.,100.,150.,200.,300./)
    458       rat = (/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5, &
    459                0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0/)
    460 
    461       ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001/)
    462       ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100,0.050,0.200,0.500,0.770,0.870,0.970/)
    463       ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400,0.430,0.580,0.790,0.930,0.960,1.000/)
    464       ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600,0.640,0.750,0.910,0.970,0.980,1.000/)
    465       ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700,0.770,0.840,0.950,0.970,1.000,1.000/)
    466       ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780,0.840,0.880,0.950,1.000,1.000,1.000/)
    467       ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830,0.870,0.900,0.950,1.000,1.000,1.000/)
    468       ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860,0.890,0.920,1.000,1.000,1.000,1.000/)
    469       ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880,0.900,0.940,1.000,1.000,1.000,1.000/)
    470       ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900,0.910,0.950,1.000,1.000,1.000,1.000/)
    471       ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900,0.910,0.950,1.000,1.000,1.000,1.000/)
    472       ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900,0.910,0.950,1.000,1.000,1.000,1.000/)
    473       ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900,0.910,0.950,1.000,1.000,1.000,1.000/)
    474       ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890,0.910,0.950,1.000,1.000,1.000,1.000/)
    475       ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880,0.920,0.950,1.000,1.000,1.000,1.000/)
    476       ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880,0.930,0.970,1.000,1.000,1.000,1.000/)
    477       ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890,0.950,1.000,1.000,1.000,1.000,1.000/)
    478       ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920,1.000,1.020,1.020,1.020,1.020,1.020/)
    479       ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010,1.030,1.040,1.040,1.040,1.040,1.040/)
    480       ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300,1.700,2.300,2.300,2.300,2.300,2.300/)
    481       ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300,3.000,4.000,4.000,4.000,4.000,4.000/)
    482 
    483 ! two-dimensional linear interpolation of the collision efficiency
    484 ! radius has to be in µm
    485 
    486       DO j = pstart, pend
    487          DO i = pstart, j
    488 
    489             DO k = 2, 15
    490                IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN
    491                   ir = k
    492                ELSEIF ((particles(j)%radius*1.0E06).gt.r0(15)) THEN
    493                   ir = 16
    494                ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN
    495                   ir = 1
    496                ENDIF
    497             ENDDO
    498 
    499             rq = particles(i)%radius*1.0E06/(particles(j)%radius*1.0E06)
    500 
    501             DO kk = 2, 21
    502                IF (rq.le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk
    503             ENDDO
    504 
    505             IF (ir.lt.16) THEN
    506                IF (ir.ge.2) THEN
    507                   pp =((particles(j)%radius * 1.0E06) - r0(ir-1)) / (r0(ir)-r0(ir-1))
    508                   qq =(rq-rat(iq-1))/(rat(iq)-rat(iq-1))
    509                   ec(j,i) = (1.-pp) * (1.-qq) * ecoll(ir-1,iq-1) +  &
    510                           pp * (1.-qq) * ecoll(ir,iq-1) +         &
    511                           qq * (1.-pp) * ecoll(ir-1,iq) +         &
    512                           pp * qq * ecoll(ir,iq)
     462
     463      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
     464
     465      LOGICAL, SAVE ::  first = .TRUE.
     466
     467      REAL ::  ek, particle_radius, pp, qq, rq
     468
     469      REAL, DIMENSION(1:21), SAVE :: rat
     470      REAL, DIMENSION(1:15), SAVE :: r0
     471      REAL, DIMENSION(1:15,1:21), SAVE :: ecoll
     472
     473!
     474!--   Initial assignment of constants
     475      IF ( first )  THEN
     476
     477         first = .FALSE.
     478         r0  = (/6.,8.,10.,15.,20.,25.,30.,40.,50., 60.,70.,100.,150.,200., &
     479                 300./)
     480         rat = (/0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,  &
     481                 0.65, 0.7,0.75,0.8,0.85,0.9,0.95,1.0/)
     482
     483         ecoll(:,1) = (/0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001, &
     484                        0.001,0.001,0.001,0.001,0.001,0.001/)
     485         ecoll(:,2) = (/0.003,0.003,0.003,0.004,0.005,0.005,0.005,0.010,0.100, &
     486                        0.050,0.200,0.500,0.770,0.870,0.970/)
     487         ecoll(:,3) = (/0.007,0.007,0.007,0.008,0.009,0.010,0.010,0.070,0.400, &
     488                        0.430,0.580,0.790,0.930,0.960,1.000/)
     489         ecoll(:,4) = (/0.009,0.009,0.009,0.012,0.015,0.010,0.020,0.280,0.600, &
     490                        0.640,0.750,0.910,0.970,0.980,1.000/)
     491         ecoll(:,5) = (/0.014,0.014,0.014,0.015,0.016,0.030,0.060,0.500,0.700, &
     492                        0.770,0.840,0.950,0.970,1.000,1.000/)
     493         ecoll(:,6) = (/0.017,0.017,0.017,0.020,0.022,0.060,0.100,0.620,0.780, &
     494                        0.840,0.880,0.950,1.000,1.000,1.000/)
     495         ecoll(:,7) = (/0.030,0.030,0.024,0.022,0.032,0.062,0.200,0.680,0.830, &
     496                        0.870,0.900,0.950,1.000,1.000,1.000/)
     497         ecoll(:,8) = (/0.025,0.025,0.025,0.036,0.043,0.130,0.270,0.740,0.860, &
     498                        0.890,0.920,1.000,1.000,1.000,1.000/)
     499         ecoll(:,9) = (/0.027,0.027,0.027,0.040,0.052,0.200,0.400,0.780,0.880, &
     500                        0.900,0.940,1.000,1.000,1.000,1.000/)
     501         ecoll(:,10)= (/0.030,0.030,0.030,0.047,0.064,0.250,0.500,0.800,0.900, &
     502                        0.910,0.950,1.000,1.000,1.000,1.000/)
     503         ecoll(:,11)= (/0.040,0.040,0.033,0.037,0.068,0.240,0.550,0.800,0.900, &
     504                        0.910,0.950,1.000,1.000,1.000,1.000/)
     505         ecoll(:,12)= (/0.035,0.035,0.035,0.055,0.079,0.290,0.580,0.800,0.900, &
     506                        0.910,0.950,1.000,1.000,1.000,1.000/)
     507         ecoll(:,13)= (/0.037,0.037,0.037,0.062,0.082,0.290,0.590,0.780,0.900, &
     508                        0.910,0.950,1.000,1.000,1.000,1.000/)
     509         ecoll(:,14)= (/0.037,0.037,0.037,0.060,0.080,0.290,0.580,0.770,0.890, &
     510                        0.910,0.950,1.000,1.000,1.000,1.000/)
     511         ecoll(:,15)= (/0.037,0.037,0.037,0.041,0.075,0.250,0.540,0.760,0.880, &
     512                        0.920,0.950,1.000,1.000,1.000,1.000/)
     513         ecoll(:,16)= (/0.037,0.037,0.037,0.052,0.067,0.250,0.510,0.770,0.880, &
     514                        0.930,0.970,1.000,1.000,1.000,1.000/)
     515         ecoll(:,17)= (/0.037,0.037,0.037,0.047,0.057,0.250,0.490,0.770,0.890, &
     516                        0.950,1.000,1.000,1.000,1.000,1.000/)
     517         ecoll(:,18)= (/0.036,0.036,0.036,0.042,0.048,0.230,0.470,0.780,0.920, &
     518                        1.000,1.020,1.020,1.020,1.020,1.020/)
     519         ecoll(:,19)= (/0.040,0.040,0.035,0.033,0.040,0.112,0.450,0.790,1.010, &
     520                        1.030,1.040,1.040,1.040,1.040,1.040/)
     521         ecoll(:,20)= (/0.033,0.033,0.033,0.033,0.033,0.119,0.470,0.950,1.300, &
     522                        1.700,2.300,2.300,2.300,2.300,2.300/)
     523         ecoll(:,21)= (/0.027,0.027,0.027,0.027,0.027,0.125,0.520,1.400,2.300, &
     524                        3.000,4.000,4.000,4.000,4.000,4.000/)
     525      ENDIF
     526
     527!
     528!--   Calculate the radius class index of particles with respect to array r
     529      ALLOCATE( ira(pstart:pend) )
     530      DO  j = pstart, pend
     531         particle_radius = particles(j)%radius * 1.0E6
     532         DO  k = 1, 15
     533            IF ( particle_radius < r0(k) )  THEN
     534               ira(j) = k
     535               EXIT
     536            ENDIF
     537         ENDDO
     538         IF ( particle_radius >= r0(15) )  ira(j) = 16
     539      ENDDO
     540
     541!
     542!--   Two-dimensional linear interpolation of the collision efficiency.
     543!--   Radius has to be in µm
     544      DO  j = pstart, pend
     545         DO  i = pstart, j
     546
     547!            DO k = 2, 15
     548!               IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN
     549!                  ir = k
     550!               ELSEIF ((particles(j)%radius*1.0E06).gt.r0(15)) THEN
     551!                  ir = 16
     552!               ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN
     553!                  ir = 1
     554!               ENDIF
     555!            ENDDO
     556
     557            ir = ira(j)
     558
     559            rq = particles(i)%radius / particles(j)%radius
     560
     561!            DO kk = 2, 21
     562!               IF ( rq .le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk
     563!            ENDDO
     564
     565            iq = INT( rq * 20 ) + 1
     566
     567            IF ( ir < 16 )  THEN
     568
     569               IF ( ir >= 2 )  THEN
     570                  pp = ( ( particles(j)%radius * 1.0E06 ) - r0(ir-1) ) / &
     571                       ( r0(ir) - r0(ir-1) )
     572                  qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     573                  ec(j,i) = ( 1.0-pp ) * ( 1.0-qq ) * ecoll(ir-1,iq-1)   &
     574                            + pp * ( 1.0-qq ) * ecoll(ir,iq-1)           &
     575                            + qq * ( 1.0-pp ) * ecoll(ir-1,iq)           &
     576                            + pp * qq * ecoll(ir,iq)
    513577               ELSE
    514578                  qq = (rq-rat(iq-1))/(rat(iq)-rat(iq-1))
    515579                  ec(j,i) = (1.-qq) * ecoll(1,iq-1) + qq * ecoll(1,iq)
    516580               ENDIF
     581
    517582            ELSE
    518                qq = (rq - rat(iq-1))/(rat(iq)-rat(iq-1))
    519                ek = (1.-qq)* ecoll(15,iq-1) + qq * ecoll(15,iq)
    520                ec(j,i) = min(ek,1.0)
     583               qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     584               ek = ( 1.0 - qq ) * ecoll(15,iq-1) + qq * ecoll(15,iq)
     585               ec(j,i) = MIN( ek, 1.0 )
    521586            ENDIF
     587
    522588            ec(i,j) = ec(j,i)
    523             IF (ec(i,j).lt.1.e-20) stop 99
     589            IF ( ec(i,j) < 1.0E-20 )  STOP 99
     590
    524591         ENDDO
    525592      ENDDO
    526       RETURN
     593
     594      DEALLOCATE( ira )
     595
    527596   END SUBROUTINE effic
    528597
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |