Changeset 799 for palm


Ignore:
Timestamp:
Dec 21, 2011 5:48:03 PM (12 years ago)
Author:
franke
Message:

Implementation of Wang collision kernel and bugfix for calculation of vpt, pt_p, and ec in case of cloud droplets

Location:
palm/trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/subjob

    r789 r799  
    140140    # VARIABLENVEREINBARUNGEN + DEFAULTWERTE
    141141 delete_dayfile=false
    142  email_notification=none 
     142 email_notification=none
    143143 group_number=none
    144144 locat=normal
     
    250250     (rte*)                  local_addres=133.5.185.60;   local_host=lcrte;;
    251251     (scirocco)              local_addres=172.20.25.41;   local_host=lcmuk;;
     252     (solano)                local_addres=130.75.105.110; local_host=lcmuk;;
    252253     (sun1|sun2)             local_addres=130.75.6.1;     local_host=unics;;
    253254     (sx-*)                  local_addres=172.16.1.131;   local_host=necriam;;
     
    885886       feature_directive=""
    886887    fi
    887    
     888
    888889    if [[ $queue = dataq || $queue = permq ]]
    889890    then
     
    952953    then
    953954       cat > $job_to_send << %%END%%
    954 #!/bin/ksh 
     955#!/bin/ksh
    955956#PBS -S /bin/ksh
    956957#PBS -N $job_name
     
    994995
    995996       cat > $job_to_send << %%END%%
    996 #!/bin/ksh 
     997#!/bin/ksh
    997998# @\$-o $remote_dayfile
    998999# @\$-eo -oi
    9991000# @\$-lP 16
    1000 # @\$-lp 1   
     1001# @\$-lp 1
    10011002# @\$-lm 28gb  -llm unlimited -ls unlimited
    10021003# @\$-q $queue
     
    10211022    then
    10221023       cat > $job_to_send << %%END%%
    1023 #!/bin/ksh 
     1024#!/bin/ksh
    10241025#PBS -S /bin/ksh
    10251026#PBS -N $job_name
     
    14421443          if [[ $queue = default ]]
    14431444          then
    1444              eval  $submcom  $job_on_remhost 
     1445             eval  $submcom  $job_on_remhost
    14451446          else
    1446              eval  $submcom  -q $queue  $job_on_remhost 
     1447             eval  $submcom  -q $queue  $job_on_remhost
    14471448          fi
    14481449       else
  • palm/trunk/SOURCE/advec_particles.f90

    r792 r799  
    44! Current revisions:
    55! ------------------
    6 ! particle arrays (particles, particles_temp) implemented as pointers in
    7 ! order to speed up sorting (see routine sort_particles)
     6! Implementation of Wang collision kernel and corresponding new parameter
     7! wang_collision_kernel
    88!
    99! Former revisions:
     
    1111! $Id$
    1212!
     13! 792 2011-12-01 raasch
     14! particle arrays (particles, particles_temp) implemented as pointers in
     15! order to speed up sorting (see routine sort_particles)
     16!
    1317! 759 2011-09-15 13:58:31Z raasch
    1418! Splitting of parallel I/O (routine write_particles)
     
    1721! Declaration of de_dx, de_dy, de_dz adapted to additional
    1822! ghost points. Furthermore the calls of exchange_horiz were modified.
    19 !
    2023!
    2124! 622 2010-12-10 08:08:13Z raasch
     
    127130    USE random_function_mod
    128131    USE statistics
     132    USE wang_kernel_mod
    129133
    130134    IMPLICIT NONE
    131135
    132136    INTEGER ::  agp, deleted_particles, deleted_tails, i, ie, ii, inc, is, j,  &
    133                 jj, js, k, kk, kw, m, n, nc, nn, num_gp, psi, tlength,         &
    134                 trlp_count, trlp_count_sum, trlp_count_recv,                   &
     137                jj, js, k, kk, kw, m, n, nc, nd, nn, num_gp, pse, psi,         &
     138                tlength, trlp_count, trlp_count_sum, trlp_count_recv,          &
    135139                trlp_count_recv_sum, trlpt_count, trlpt_count_recv,            &
    136140                trnp_count, trnp_count_sum, trnp_count_recv,                   &
     
    139143                trrp_count_recv_sum, trrpt_count, trrpt_count_recv,            &
    140144                trsp_count, trsp_count_sum, trsp_count_recv,                   &
    141                 trsp_count_recv_sum, trspt_count, trspt_count_recv, nd
     145                trsp_count_recv_sum, trspt_count, trspt_count_recv 
    142146
    143147    INTEGER ::  gp_outside_of_building(1:8)
     
    151155                dt_particle, dt_particle_m, d_radius, d_sum, e_a, e_int,       &
    152156                e_int_l, e_int_u, e_mean_int, e_s, exp_arg, exp_term, fs_int,  &
    153                 gg,lagr_timescale, mean_r, new_r, p_int, pt_int, pt_int_l,     &
    154                 pt_int_u, q_int, q_int_l, q_int_u, ql_int, ql_int_l, ql_int_u, &
    155                 random_gauss, sl_r3, sl_r4, t_int, u_int, u_int_l, u_int_u,    &
    156                 vv_int, v_int, v_int_l, v_int_u, w_int, w_int_l, w_int_u,      &
    157                 x, y
     157                gg, integral, lagr_timescale, lw_max, mean_r, new_r, p_int,    &
     158                pt_int, pt_int_l, pt_int_u, q_int, q_int_l, q_int_u, ql_int,  &
     159                ql_int_l, ql_int_u, random_gauss, sl_r3, sl_r4, t_int, u_int,  &
     160                u_int_l, u_int_u,vv_int, v_int, v_int_l, v_int_u, w_int,       &
     161                w_int_l, w_int_u, x, y
    158162
    159163    REAL, DIMENSION(1:30) ::  de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
     164    REAL, DIMENSION(:,:), ALLOCATABLE ::  kern
    160165
    161166    REAL    ::  location(1:30,1:3)
     
    399404       CALL cpu_log( log_point_s(43), 'advec_part_coll', 'start' )
    400405
     406       IF ( wang_collision_kernel ) THEN
     407
     408       DO  i = nxl, nxr
     409          DO  j = nys, nyn
     410             DO  k = nzb+1, nzt
     411!
     412!--             Collision requires at least two particles in the box
     413                IF ( prt_count(k,j,i) > 1 )  THEN
     414!
     415!--                First, sort particles within the gridbox by their size,
     416!--                using Shell's method (see Numerical Recipes)
     417!--                NOTE: In case of using particle tails, the re-sorting of
     418!--                ----  tails would have to be included here!
     419                   psi = prt_start_index(k,j,i) - 1
     420                   inc = 1
     421                   DO WHILE ( inc <= prt_count(k,j,i) )
     422                      inc = 3 * inc + 1
     423                   ENDDO
     424
     425                   DO WHILE ( inc > 1 )
     426                      inc = inc / 3
     427                      DO  is = inc+1, prt_count(k,j,i)
     428                         tmp_particle = particles(psi+is)
     429                         js = is
     430                         DO WHILE ( particles(psi+js-inc)%radius >             &
     431                                    tmp_particle%radius )
     432                            particles(psi+js) = particles(psi+js-inc)
     433                            js = js - inc
     434                            IF ( js <= inc )  EXIT
     435                         ENDDO
     436                         particles(psi+js) = tmp_particle
     437                      ENDDO
     438                   ENDDO
     439
     440                   psi = prt_start_index(k,j,i)
     441                   pse = psi + prt_count(k,j,i)-1
     442
     443                   ALLOCATE( kern(psi:pse,psi:pse) )
     444
     445!
     446!--                Calculate collision kernel for all particles in grid box
     447                   CALL colker( i, j, k, kern )
     448!
     449!--                collison kernel is calculated in cm**3/s but needed in m**3/s
     450                   kern = kern * 1.0E-6
     451
     452                   DO  n = pse, psi+1, -1
     453
     454                      integral = 0.0
     455                      lw_max   = 0.0
     456
     457!
     458!--                   Calculate growth of collector particle radius using all
     459!--                   droplets smaller than current droplet
     460                      DO  is = psi, n-1
     461
     462                         integral = integral +                               &
     463                                    ( particles(is)%radius**3 * kern(n,is) * &
     464                                      particles(is)%weight_factor )
     465!
     466!--                      Calculate volume of liquid water of the collected
     467!--                      droplets which is the maximum liquid water available
     468!--                      for droplet growth   
     469                         lw_max = lw_max + ( particles(is)%radius**3 * &
     470                                             particles(is)%weight_factor )
     471
     472                      ENDDO
     473
     474!
     475!--                   Change in radius of collector droplet due to collision
     476                      delta_r = 1.0 / ( 3.0 * particles(n)%radius**2 ) * &
     477                                integral * dt_3d * ddx * ddy / dz
     478
     479!
     480!--                   Change in volume of collector droplet due to collision
     481                      delta_v = particles(n)%weight_factor                     &
     482                                * ( ( particles(n)%radius + delta_r )**3       &
     483                                    - particles(n)%radius**3 )
     484
     485                      IF ( lw_max < delta_v  .AND.  delta_v > 0.0 )  THEN
     486
     487                         PRINT*, 'Particle has grown to much because the       &
     488                                  change of volume of particle is larger than  &
     489                                  liquid water available!'
     490
     491                      ELSEIF ( lw_max == delta_v  .AND.  delta_v > 0.0 )  THEN
     492
     493                         DO  is = psi, n-1
     494
     495                            particles(is)%weight_factor = 0.0
     496                            particle_mask(is)  = .FALSE.
     497                            deleted_particles = deleted_particles + 1
     498
     499                         ENDDO
     500
     501                      ELSEIF ( lw_max > delta_v  .AND.  delta_v > 0.0 )  THEN
     502!
     503!--                      Calculate new weighting factor of collected droplets
     504                         DO  is = psi, n-1
     505
     506                            particles(is)%weight_factor = &
     507                               particles(is)%weight_factor - &
     508                               ( ( kern(n,is) * particles(is)%weight_factor / &
     509                                   integral ) * delta_v )
     510
     511                            IF ( particles(is)%weight_factor < 0.0 )  THEN
     512                                  WRITE( message_string, * ) 'negative ',      &
     513                                                     'weighting factor: ',     &
     514                                                     particles(is)%weight_factor
     515                                  CALL message( 'advec_particles', '', 2, 2,   &
     516                                                -1, 6, 1 )
     517
     518                            ELSEIF ( particles(is)%weight_factor == 0.0 )  THEN
     519
     520                                particles(is)%weight_factor = 0.0
     521                                particle_mask(is)  = .FALSE.
     522                                deleted_particles = deleted_particles + 1
     523
     524                            ENDIF
     525
     526                         ENDDO
     527
     528                      ENDIF
     529
     530                      particles(n)%radius = particles(n)%radius + delta_r
     531                      ql_vp(k,j,i) = ql_vp(k,j,i) + particles(n)%weight_factor &
     532                                                    * particles(n)%radius**3
     533
     534                   ENDDO
     535
     536                   ql_vp(k,j,i) = ql_vp(k,j,i) + particles(psi)%weight_factor  &
     537                                                 * particles(psi)%radius**3
     538
     539                   DEALLOCATE( kern )
     540
     541                ELSE IF ( prt_count(k,j,i) == 1 )  THEN
     542
     543                   psi = prt_start_index(k,j,i)
     544                   ql_vp(k,j,i) =  particles(psi)%weight_factor *              &
     545                                   particles(psi)%radius**3
     546                ENDIF
     547
     548!
     549!--             Check if condensation of LWC was conserved during collision
     550!--             process
     551                IF ( ql_v(k,j,i) /= 0.0 )  THEN
     552                   IF ( ql_vp(k,j,i) / ql_v(k,j,i) >= 1.0001  .OR.             &
     553                        ql_vp(k,j,i) / ql_v(k,j,i) <= 0.9999 )  THEN
     554                      WRITE( message_string, * ) 'LWC is not conserved during',&
     555                                                 ' collision! ',               &
     556                                                 'LWC after condensation: ',   &
     557                                                 ql_v(k,j,i),                  &
     558                                                 ' LWC after collision: ',     &
     559                                                 ql_vp(k,j,i)
     560                      CALL message( 'advec_particles', '', 2, 2, -1, 6,  &
     561                                    1 )
     562                   ENDIF
     563                ENDIF
     564
     565             ENDDO
     566          ENDDO
     567       ENDDO
     568
     569       ELSE
     570
    401571       DO  i = nxl, nxr
    402572          DO  j = nys, nyn
     
    683853          ENDDO
    684854       ENDDO
     855
     856       ENDIF
    685857
    686858       CALL cpu_log( log_point_s(43), 'advec_part_coll', 'stop' )
     
    12001372                   de_dzi(num_gp) = de_dz(k+1,j+1,i)
    12011373                ENDIF
    1202  
     1374
    12031375                IF ( k > nzb_s_inner(j,i+1)  .OR.  nzb_s_inner(j,i+1) == 0 ) &
    12041376                THEN
     
    16201792                      de_dzi(num_gp) = 0.0
    16211793                   ENDIF
    1622  
     1794
    16231795!
    16241796!--                If wall between gridpoint 5 and gridpoint 7, then
     
    19712143#else
    19722144       dt_3d_reached = dt_3d_reached_l
    1973 #endif     
     2145#endif
    19742146
    19752147       CALL cpu_log( log_point_s(44), 'advec_part_advec', 'stop' )
     
    21622334!--       data
    21632335          DO  n = 1, number_of_particles
    2164              i = ( particles(n)%x + 0.5 * dx ) * ddx
     2336             i = ( particles(n)%x + 0.5 * dx ) * ddx
    21652337!
    21662338!--          Above calculation does not work for indices less than zero
    2167              IF ( particles(n)%x < -0.5 * dx )  i = -1
     2339             IF ( particles(n)%x < -0.5 * dx )  i = -1
    21682340
    21692341             IF ( i < nxl )  THEN
     
    22252397          nn = particles(n)%tail_id
    22262398
    2227           i = ( particles(n)%x + 0.5 * dx ) * ddx
     2399          i = ( particles(n)%x + 0.5 * dx ) * ddx
    22282400!
    22292401!--       Above calculation does not work for indices less than zero
    2230           IF ( particles(n)%x < - 0.5 * dx )  i = -1
     2402          IF ( particles(n)%x < - 0.5 * dx )  i = -1
    22312403
    22322404          IF ( i <  nxl )  THEN
     
    26022774          DO  n = 1, number_of_particles
    26032775             IF ( particle_mask(n) )  THEN
    2604                 j = ( particles(n)%y + 0.5 * dy ) * ddy
     2776                j = ( particles(n)%y + 0.5 * dy ) * ddy
    26052777!
    26062778!--             Above calculation does not work for indices less than zero
    2607                 IF ( particles(n)%y < -0.5 * dy )  j = -1
     2779                IF ( particles(n)%y < -0.5 * dy )  j = -1
    26082780
    26092781                IF ( j < nys )  THEN
     
    26692841!--       moved.
    26702842          IF ( particle_mask(n) )  THEN
    2671              j = ( particles(n)%y + 0.5 * dy ) * ddy
     2843             j = ( particles(n)%y + 0.5 * dy ) * ddy
    26722844!
    26732845!--          Above calculation does not work for indices less than zero
    2674              IF ( particles(n)%y < -0.5 * dy )  j = -1
     2846             IF ( particles(n)%y < -0.5 * dy )  j = -1
    26752847
    26762848             IF ( j < nys )  THEN
     
    38153987
    38163988    CHARACTER (LEN=10) ::  particle_binary_version
    3817 
    38183989    INTEGER ::  i
    38193990
     
    38794050
    38804051 END SUBROUTINE write_particles
     4052
    38814053
    38824054
  • palm/trunk/SOURCE/compute_vpt.f90

    r484 r799  
    44! Current revisions:
    55! -----------------
    6 !
     6! Bugfix: ql is now included in calculation of vpt in case of
     7!         cloud droplets
    78!
    89! Former revisions:
     
    3233    INTEGER :: k
    3334
    34     IF ( .NOT. cloud_physics ) THEN
     35    IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
    3536       vpt = pt * ( 1.0 + 0.61 * q )
    36     ELSE
     37    ELSE IF (cloud_physics) THEN
    3738       DO  k = nzb, nzt+1
    3839          vpt(k,:,:) = ( pt(k,:,:) + pt_d_t(k) * l_d_cp * ql(k,:,:) ) * &
    3940                       ( 1.0 + 0.61 * q(k,:,:) - 1.61 * ql(k,:,:) )
    4041       ENDDO
     42    ELSE
     43       vpt = pt * ( 1.0 + 0.61 * q - 1.61 * ql )
    4144    ENDIF
    4245
  • palm/trunk/SOURCE/interaction_droplets_ptq.f90

    r484 r799  
    44! Current revisions:
    55! -----------------
    6 !
     6! Bugfix: pt_d_t(k) was missing in calculation of pt_p
    77!
    88! Former revisions:
     
    5353             DO  k = nzb_2d(j,i)+1, nzt
    5454                q_p(k,j,i)  = q_p(k,j,i)  - ql_c(k,j,i)
    55                 pt_p(k,j,i) = pt_p(k,j,i) + l_d_cp * ql_c(k,j,i)
     55                pt_p(k,j,i) = pt_p(k,j,i) + l_d_cp * ql_c(k,j,i) * pt_d_t(k)
    5656             ENDDO
    5757          ENDDO
     
    8080       DO  k = nzb_2d(j,i)+1, nzt
    8181          q_p(k,j,i)  = q_p(k,j,i)  - ql_c(k,j,i)
    82           pt_p(k,j,i) = pt_p(k,j,i) + l_d_cp * ql_c(k,j,i)
     82          pt_p(k,j,i) = pt_p(k,j,i) + l_d_cp * ql_c(k,j,i) * pt_d_t(k)
    8383       ENDDO
    8484
  • palm/trunk/SOURCE/wang_kernel.f90

    r792 r799  
    44! Current revisions:
    55! -----------------
    6 ! speed optimizations + formatting
     6! speed optimizations and formatting
     7! Bugfix: iq=1 is not allowed (routine effic)
     8! Bugfix: replaced stop by ec=0.0 in case of very small ec (routine effic) 
    79!
    810! Former revisions:
     
    169171! of droplets, air density, air viscosity, turbulent dissipation rate,
    170172! taylor microscale reynolds number, gravitational acceleration
     173
    171174       USE constants
    172175       USE cloud_parameters
     
    176179       IMPLICIT NONE
    177180
    178        REAL :: airvisc, airdens, waterdens, gravity, anu, Relamda, &
     181       REAL :: Relamda, &
    179182               Tl, Lf, tauk, eta, vk, ao, lambda, tt, z, be,       &
    180183               bbb, d1, e1, d2, e2, ccc, b1, c1, b2, c2, v1xysq,  &
     
    182185               wrFIN, SSt, XX, YY, c1_gr, ao_gr, fao_gr, rc, grFIN, v1, t1, v2, t2, rrp
    183186
    184        REAL, DIMENSION(pstart:pend) :: vST, tau, St
     187       REAL, SAVE :: airvisc, airdens, anu, gravity, waterdens
     188
     189       REAL, DIMENSION(pstart:pend) :: St, tau
     190
     191       LOGICAL, SAVE ::  first = .TRUE.
    185192
    186193       INTEGER :: i, j
    187194
    188        airvisc   = 0.1818   !dynamic viscosity in mg/cm*s
    189        airdens   = 1.2250   !air density in mg/cm**3
    190        waterdens = 1000.0   !water density in mg/cm**3
    191        gravity   = 980.6650 !in cm/s**2
    192 
    193        anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
    194 
    195        Relamda = urms**2.0*sqrt(15.0/epsilon/anu)
    196 
    197        Tl = urms**2.0/epsilon       !in s
    198        Lf = 0.5 * (urms**3.0)/epsilon !in cm
     195!
     196!--   Initial assignment of constants
     197       IF ( first )  THEN
     198
     199          first = .FALSE.
     200          airvisc   = 0.1818   !dynamic viscosity in mg/cm*s
     201          airdens   = 1.2250   !air density in mg/cm**3
     202          waterdens = 1000.0   !water density in mg/cm**3
     203          gravity   = 980.6650 !in cm/s**2
     204          anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
     205
     206       ENDIF   
     207
     208       Relamda = urms**2*sqrt(15.0/epsilon/anu)
     209
     210       Tl = urms**2/epsilon       !in s
     211       Lf = 0.5 * (urms**3)/epsilon !in cm
    199212
    200213       tauk = (anu/epsilon)**0.5       !in s
    201        eta  = (anu**3.0/epsilon)**0.25  !in cm
     214       eta  = (anu**3/epsilon)**0.25  !in cm
    202215       vk   = eta/tauk                   !in cm/s
    203216
     
    211224
    212225       DO i = pstart, pend
    213           vST(i) = winf(i)           !in cm/s
    214           tau(i) = vST(i)/gravity    !in s
     226          tau(i) = winf(i)/gravity    !in s
    215227          St(i)  = tau(i)/tauk
    216228       ENDDO
     
    222234       be = sqrt(2.0)*lambda/Lf
    223235
    224        bbb = sqrt(1.0-2.0*be**2.0)
     236       bbb = sqrt(1.0-2.0*be**2)
    225237       d1 = (1.+bbb)/2.0/bbb
    226238       e1 = Lf*(1.0+bbb)/2.0   !in cm
     
    228240       e2 = Lf*(1.0-bbb)/2.0   !in cm
    229241
    230        ccc = sqrt(1.0-2.0*z**2.0)
     242       ccc = sqrt(1.0-2.0*z**2)
    231243       b1 = (1.+ccc)/2./ccc
    232244       c1 = Tl*(1.+ccc)/2.   !in s
     
    235247
    236248       DO i = pstart, pend
    237           v1 = vST(i)         !in cm/s
     249          v1 = winf(i)         !in cm/s
    238250          t1 = tau(i)         !in s
    239251
    240252          DO j = pstart,i
    241253             rrp = (particles(i)%radius + particles(j)%radius) * 100.0 !radius in cm
    242              v2  = vST(j)         !in cm/s
     254             v2  = winf(j)         !in cm/s
    243255             t2  = tau(j)         !in s
    244256
     
    247259                     - b2*d1*PHI(c2,e1,v1,t1) &
    248260                     + b2*d2*PHI(c2,e2,v1,t1)
    249              v1xysq = v1xysq * urms**2.0/t1    !in cm**2/s**2
     261             v1xysq = v1xysq * urms**2/t1    !in cm**2/s**2
    250262
    251263             vrms1xy= sqrt(v1xysq)             !in cm/s
     
    255267                     - b2*d1*PHI(c2,e1,v2,t2) &
    256268                     + b2*d2*PHI(c2,e2,v2,t2)
    257              v2xysq = v2xysq * urms**2.0/t2    !in cm**2/s**2
     269             v2xysq = v2xysq * urms**2/t2    !in cm**2/s**2
    258270
    259271              vrms2xy= sqrt(v2xysq)             !in cm/s
    260272
    261              IF(vST(i).ge.vST(j)) THEN
    262                 v1 = vST(i)
     273             IF(winf(i).ge.winf(j)) THEN
     274                v1 = winf(i)
    263275                t1 = tau(i)
    264                 v2 = vST(j)
     276                v2 = winf(j)
    265277                t2 = tau(j)
    266278             ELSE
    267                 v1 = vST(j)
     279                v1 = winf(j)
    268280                t1 = tau(j)
    269                 v2 = vST(i)
     281                v2 = winf(i)
    270282                t2 = tau(i)
    271283             ENDIF
     
    276288                     + b2*d2*ZHI(c2,e2,v1,t1,v2,t2)
    277289             fR = d1 * exp(-rrp/e1) - d2 * exp(-rrp/e2)
    278              v1v2xy = v1v2xy * fR * urms**2.0/tau(i)/tau(j)    !in cm**2/s**2
    279 
    280              wrtur2xy=vrms1xy**2.0 + vrms2xy**2.0 - 2.*v1v2xy  !in cm**2/s**2
     290             v1v2xy = v1v2xy * fR * urms**2/tau(i)/tau(j)    !in cm**2/s**2
     291
     292             wrtur2xy=vrms1xy**2 + vrms2xy**2 - 2.*v1v2xy  !in cm**2/s**2
    281293
    282294              IF (wrtur2xy.le.0.0) wrtur2xy=0.0
    283295
    284               wrgrav2=pi/8.*(vST(j)-vST(i))**2.0
     296              wrgrav2=pi/8.*(winf(j)-winf(i))**2
    285297
    286298              wrFIN = sqrt((2.0/pi)*(wrtur2xy+wrgrav2))    !in cm/s
     
    295307             ENDIF
    296308
    297              XX = -0.1988*SSt**4.0 + 1.5275*SSt**3.0 &
    298                   -4.2942*SSt**2.0 + 5.3406*SSt
     309             XX = -0.1988*SSt**4 + 1.5275*SSt**3 &
     310                  -4.2942*SSt**2 + 5.3406*SSt
    299311
    300312             IF(XX.le.0.0) XX = 0.0
     
    304316             c1_gr = XX/(gravity/(vk/tauk))**YY
    305317
    306              ao_gr  = ao + (pi/8.)*(gravity/(vk/tauk))**2.0
     318             ao_gr  = ao + (pi/8.)*(gravity/(vk/tauk))**2
    307319             fao_gr = 20.115 * (ao_gr/Relamda)**0.5
    308320             rc     = sqrt( fao_gr * abs(St(j)-St(i)) ) * eta   !in cm
    309321
    310              grFIN = ((eta**2.0+rc**2.0)/(rrp**2.0+rc**2.0))**(c1_gr/2.)
     322             grFIN = ((eta**2+rc**2)/(rrp**2+rc**2))**(c1_gr/2.)
    311323             IF (grFIN.lt.1.0) grFIN = 1.0
    312324
    313              gck(i,j) = 2. * pi * rrp**2.0 * wrFIN * grFIN       ! in cm**3/s
     325             gck(i,j) = 2. * pi * rrp**2 * wrFIN * grFIN       ! in cm**3/s
    314326             gck(j,i) = gck(i,j)
    315327
     
    330342       aa1 = 1./tau0 + 1./a + vsett/b
    331343
    332        PHI = 1./aa1 - vsett/2.0/b/aa1**2.0  !in s
     344       PHI = 1./aa1 - vsett/2.0/b/aa1**2  !in s
    333345
    334346    END FUNCTION PHI
     
    346358        aa2 = vsett1/b + 1./tau1 + 1./a
    347359        aa3 = (vsett1-vsett2)/b + 1./tau1 + 1./tau2
    348         aa4 = (vsett2/b)**2.0 - (1./tau2 + 1./a)**2.0
     360        aa4 = (vsett2/b)**2 - (1./tau2 + 1./a)**2
    349361        aa5 = vsett2/b + 1./tau2 + 1./a
    350362        aa6 = 1./tau1 - 1./a + (1./tau2 + 1./a) * vsett1/vsett2
    351         ZHI =  (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2.0          &
    352              + (4./aa4 - 1./aa5**2.0 - 1./aa1**2.0) * vsett2/2./b/aa6      &
    353              + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2.0 + vsett2/aa1**2.0)   &
     363        ZHI =  (1./aa1 - 1./aa2) * (vsett1-vsett2)/2./b/aa3**2          &
     364             + (4./aa4 - 1./aa5**2 - 1./aa1**2) * vsett2/2./b/aa6      &
     365             + (2.*b/aa2 - 2.*b/aa1 - vsett1/aa2**2 + vsett2/aa1**2)   &
    354366                                                           * 1./2./b/aa3             ! in s**2
     367
    355368    END FUNCTION ZHI
    356369
    357 
    358370!------------------------------------------------------------------------------!
    359371! SUBROUTINE for calculation of terminal velocity winf
    360372!------------------------------------------------------------------------------!
    361 
    362373   SUBROUTINE fallg
    363374
     
    371382      INTEGER :: i, j
    372383
    373       REAL :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py,  &
    374               x, y, xrey, bond
    375 
    376       REAL, DIMENSION(1:7)  :: b
    377       REAL, DIMENSION(1:6)  :: c
    378       REAL, DIMENSION(1:20) :: rat
    379       REAL, DIMENSION(1:15) :: r0
    380       REAL, DIMENSION(1:15,1:20) :: ecoll
    381 
    382       b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, &
     384      LOGICAL, SAVE ::  first = .TRUE.
     385
     386      REAL, SAVE :: eta, xlamb, rhoa, rhow, grav, cunh, t0, sigma, stok, stb, phy, py
     387
     388      REAL :: bond, x, xrey, y
     389
     390      REAL, DIMENSION(1:7), SAVE  :: b
     391      REAL, DIMENSION(1:6), SAVE  :: c
     392
     393!
     394!--   Initial assignment of constants
     395      IF ( first )  THEN
     396
     397         first = .FALSE.
     398         b = (/-0.318657e1,0.992696e0,-0.153193e-2,-0.987059e-3, &
    383399             -0.578878e-3,0.855176e-4,-0.327815e-5/)
    384       c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0,    &
     400         c = (/-0.500015e1,0.523778e1,-0.204914e1,0.475294e0,    &
    385401             -0.542819e-1, 0.238449e-2/)
    386402
    387       eta   = 1.818e-4     !in poise = g/(cm s)
    388       xlamb = 6.62e-6    !in cm
    389 
    390       rhow = 1.0       !in g/cm**3
    391       rhoa = 1.225e-3    !in g/cm**3
    392 
    393       grav = 980.665     !in cm/s**2
    394       cunh = 1.257 * xlamb !in cm
    395       t0   = 273.15        !in K
    396       sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2
    397       stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s)
    398       stb  = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3
    399       phy  = (sigma**3.0) * (rhoa**2.0)/((eta**4.0) * grav * (rhow-rhoa))
    400       py   = phy**(1.0/6.0)
     403         eta   = 1.818e-4     !in poise = g/(cm s)
     404         xlamb = 6.62e-6    !in cm
     405
     406         rhow = 1.0       !in g/cm**3
     407         rhoa = 1.225e-3    !in g/cm**3
     408
     409         grav = 980.665     !in cm/s**2
     410         cunh = 1.257 * xlamb !in cm
     411         t0   = 273.15        !in K
     412         sigma= 76.1 - 0.155 * (293.15 - t0) !in N/m = g/s**2
     413         stok = 2.0 * grav * (rhow - rhoa)/(9.0 * eta) ! in 1/(cm s)
     414         stb  = 32.0 * rhoa * (rhow-rhoa) * grav/(3.0 * eta * eta) ! in 1/cm**3
     415         phy  = (sigma**3) * (rhoa**2)/((eta**4) * grav * (rhow-rhoa))
     416         py   = phy**(1.0/6.0)
     417
     418      ENDIF
    401419
    402420!particle radius has to be in cm
     
    405423         IF (particles(j)%radius*100.0 .le. 1.e-3) THEN
    406424
    407             winf(j)=stok*((particles(j)%radius*100.0)**2.+cunh* particles(j)%radius*100.0) !in cm/s
     425            winf(j)=stok*((particles(j)%radius*100.0)**2+cunh* particles(j)%radius*100.0) !in cm/s
    408426
    409427         ELSEIF (particles(j)%radius*100.0.gt.1.e-3.and.particles(j)%radius*100.0.le.5.35e-2) THEN
    410428
    411             x = log(stb*(particles(j)%radius*100.0)**3.0)
     429            x = log(stb*(particles(j)%radius*100.0)**3)
    412430            y = 0.0
    413431
     
    421439         ELSEIF (particles(j)%radius*100.0.gt.5.35e-2) THEN
    422440
    423             bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2.0/sigma
    424 
    425441            IF (particles(j)%radius*100.0.gt.0.35)  THEN
    426442               bond = grav*(rhow-rhoa) * 0.35 * 0.35/sigma
     443            ELSE
     444               bond = grav*(rhow-rhoa)*(particles(j)%radius*100.0)**2/sigma
    427445            ENDIF
    428446
     
    435453
    436454            xrey = py*exp(y)
    437             winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
    438455
    439456            IF (particles(j)%radius*100.0 .gt.0.35)  THEN
    440457              winf(j) = xrey * eta/(2.0 * rhoa * 0.35) !in cm/s
     458            ELSE
     459               winf(j) = xrey*eta/(2.0*rhoa*particles(j)%radius*100.0)  !in cm/s
    441460            ENDIF
    442461
    443462         ENDIF
    444463      ENDDO
    445 
     464      RETURN
    446465   END SUBROUTINE fallg
    447466
     
    457476      USE particle_attributes
    458477
     478!collision efficiencies of hall kernel
    459479      IMPLICIT NONE
    460480
     
    545565         DO  i = pstart, j
    546566
    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
     567            ir = ira(j)
     568
     569            rq = particles(i)%radius / particles(j)%radius
     570
     571!            DO kk = 2, 21
     572!               IF ( rq <= rat(kk) )  THEN
     573!                  iq = kk
     574!                  EXIT
    554575!               ENDIF
    555576!            ENDDO
    556577
    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 
    565578            iq = INT( rq * 20 ) + 1
     579            iq = MAX(iq , 2)
    566580
    567581            IF ( ir < 16 )  THEN
     
    587601
    588602            ec(i,j) = ec(j,i)
    589             IF ( ec(i,j) < 1.0E-20 )  STOP 99
     603            IF ( ec(i,j) < 1.0E-20 )  ec(i,j) = 0.0
    590604
    591605         ENDDO
     
    599613! SUBROUTINE for calculation of enhancement factor collision efficencies
    600614!------------------------------------------------------------------------------!
    601 
    602615   SUBROUTINE turb_enhance_eff
    603616
     
    607620       USE arrays_3d
    608621
    609 !collision efficiencies of hall kernel
    610622      IMPLICIT NONE
    611623
    612624      INTEGER :: i, ik, ir, iq, j, k, kk
    613       REAL    :: rq, y1, pp, qq, y2, y3, x1, x2, x3
    614 
    615       REAL, DIMENSION(1:11) :: rat
    616       REAL, DIMENSION(1:7)  :: r0
    617       REAL, DIMENSION(1:7,1:11) :: ecoll_100, ecoll_400
    618 
    619       r0  = (/10., 20., 30.,40., 50., 60.,100./)
    620       rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/)
     625
     626      INTEGER, DIMENSION(:), ALLOCATABLE ::  ira
     627
     628      REAL    :: rq, y1, particle_radius, pp, qq, y2, y3, x1, x2, x3
     629
     630      LOGICAL, SAVE ::  first = .TRUE.
     631
     632      REAL, DIMENSION(1:11), SAVE :: rat
     633      REAL, DIMENSION(1:7), SAVE  :: r0
     634      REAL, DIMENSION(1:7,1:11), SAVE :: ecoll_100, ecoll_400
     635
     636!
     637!--   Initial assignment of constants
     638      IF ( first )  THEN
     639
     640         first = .FALSE.
     641
     642         r0  = (/10., 20., 30.,40., 50., 60.,100./)
     643         rat = (/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/)
    621644
    622645! 100 cm^2/s^3
    623       ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
    624       ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
    625       ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
    626       ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
    627       ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
    628       ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
    629       ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
    630       ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
    631       ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
    632       ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
    633       ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
     646         ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
     647         ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
     648         ecoll_100(:,3) = (/1.32,  1.32,  1.245, 1.123, 1.000,  1.000,  1.0 /)
     649         ecoll_100(:,4) = (/1.250, 1.250, 1.148, 1.087, 1.025,  1.025,  1.0 /)
     650         ecoll_100(:,5) = (/1.186, 1.186, 1.066, 1.060, 1.056,  1.056,  1.0 /)
     651         ecoll_100(:,6) = (/1.045, 1.045, 1.000, 1.014, 1.028,  1.028,  1.0 /)
     652         ecoll_100(:,7) = (/1.070, 1.070, 1.030, 1.038, 1.046,  1.046,  1.0 /)
     653         ecoll_100(:,8) = (/1.000, 1.000, 1.054, 1.042, 1.029,  1.029,  1.0 /)
     654         ecoll_100(:,9) = (/1.223, 1.223, 1.117, 1.069, 1.021,  1.021,  1.0 /)
     655         ecoll_100(:,10)= (/1.570, 1.570, 1.244, 1.166, 1.088,  1.088,  1.0 /)
     656         ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
    634657
    635658! 400 cm^2/s^3
    636       ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
    637       ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
    638       ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
    639       ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
    640       ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
    641       ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
    642       ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
    643       ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
    644       ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
    645       ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
    646       ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
     659         ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
     660         ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
     661         ecoll_400(:,3) = (/1.988, 1.988,  1.475, 1.313, 1.150,  1.150,  1.0 /)
     662         ecoll_400(:,4) = (/1.490, 1.490,  1.187, 1.156, 1.126,  1.126,  1.0 /)
     663         ecoll_400(:,5) = (/1.249, 1.249,  1.088, 1.090, 1.092,  1.092,  1.0 /)
     664         ecoll_400(:,6) = (/1.139, 1.139,  1.130, 1.091, 1.051,  1.051,  1.0 /)
     665         ecoll_400(:,7) = (/1.220, 1.220,  1.190, 1.138, 1.086,  1.086,  1.0 /)
     666         ecoll_400(:,8) = (/1.325, 1.325,  1.267, 1.165, 1.063,  1.063,  1.0 /)
     667         ecoll_400(:,9) = (/1.716, 1.716,  1.345, 1.223, 1.100,  1.100,  1.0 /)
     668         ecoll_400(:,10)= (/3.788, 3.788,  1.501, 1.311, 1.120,  1.120,  1.0 /)
     669         ecoll_400(:,11)= (/36.52, 36.52,  19.16, 22.80,  26.0,   26.0,  1.0 /)
     670
     671      ENDIF
     672
     673!
     674!--   Calculate the radius class index of particles with respect to array r
     675      ALLOCATE( ira(pstart:pend) )
     676
     677      DO  j = pstart, pend
     678         particle_radius = particles(j)%radius * 1.0E6
     679         DO  k = 1, 7
     680            IF ( particle_radius < r0(k) )  THEN
     681               ira(j) = k
     682               EXIT
     683            ENDIF
     684         ENDDO
     685         IF ( particle_radius >= r0(7) )  ira(j) = 8
     686      ENDDO
    647687
    648688! two-dimensional linear interpolation of the collision efficiency
     
    650690         DO i = pstart, j
    651691
    652             DO k = 2, 7
    653                IF ((particles(j)%radius*1.0E06).le.r0(k).and.(particles(j)%radius*1.0E06).ge.r0(k-1)) THEN
    654                   ir = k
    655                ELSEIF ((particles(j)%radius*1.0E06).gt.r0(7)) THEN
    656                   ir = 8
    657                ELSEIF ((particles(j)%radius*1.0E06).lt.r0(1)) THEN
    658                   ir = 1
     692            ir = ira(j)
     693
     694            rq = particles(i)%radius/particles(j)%radius
     695
     696            DO kk = 2, 11
     697               IF ( rq <= rat(kk) )  THEN
     698                  iq = kk
     699                  EXIT
    659700               ENDIF
    660             ENDDO
    661 
    662             rq = particles(i)%radius*1.0E06/(particles(j)%radius*1.0E06)
    663 
    664             DO kk = 2, 11
    665                IF (rq.le.rat(kk).and.rq.gt.rat(kk-1)) iq = kk
    666701            ENDDO
    667702
Note: See TracChangeset for help on using the changeset viewer.