Changeset 1322 for palm/trunk/SOURCE/lpm_collision_kernels.f90

Timestamp:

Mar 20, 2014 4:38:49 PM (10 years ago)

Author:

raasch

Message:

REAL functions and a lot of REAL constants provided with KIND-attribute,
some small bugfixes

File:

• palm/trunk/SOURCE/lpm_collision_kernels.f90 (modified) (16 diffs)

palm/trunk/SOURCE/lpm_collision_kernels.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
+! REAL constants defined as wp_kind
 !
 ! Former revisions:
@@ -167 +168 @@
 !--       Calculate the radius class bounds with logarithmic distances
 !--       in the interval [1.0E-6, 2.0E-4] m
-          rclass_lbound = LOG( 1.0E-6 )
-          rclass_ubound = LOG( 2.0E-4 )
+          rclass_lbound = LOG( 1.0E-6_wp )
+          rclass_ubound = LOG( 2.0E-4_wp )
           radclass(1)   = 1.0E-6
           DO  i = 2, radius_classes
@@ -182 +183 @@
 !--       Set the class bounds for dissipation in interval [0.0, 0.1] m**2/s**3
           DO  i = 1, dissipation_classes
-             epsclass(i) = 0.1 * REAL( i ) / dissipation_classes
+             epsclass(i) = 0.1 * REAL( i, KIND=wp ) / dissipation_classes
 !             IF ( myid == 0 )  THEN
 !                PRINT*, 'i=', i, ' eps = ', epsclass(i)
@@ -197 +198 @@
 
              epsilon = epsclass(k)
-             urms    = 2.02 * ( epsilon / 0.04 )**( 1.0 / 3.0 )
+             urms    = 2.02 * ( epsilon / 0.04_wp )**( 1.0_wp / 3.0_wp )
 
              CALL turbsd
@@ -313 +314 @@
           epsilon = 0.0
        ENDIF
-       urms    = 2.02 * ( epsilon / 0.04 )**( 0.33333333333 )
+       urms    = 2.02 * ( epsilon / 0.04_wp )**( 0.33333333333_wp )
 
        IF ( wang_kernel  .AND.  epsilon > 1.0E-7 )  THEN
@@ -431 +432 @@
        ENDIF
 
-       lambda    = urms * SQRT( 15.0 * molecular_viscosity / epsilon )    ! in m
-       lambda_re = urms**2 * SQRT( 15.0 / epsilon / molecular_viscosity )
+       lambda    = urms * SQRT( 15.0_wp * molecular_viscosity / epsilon ) ! in m
+       lambda_re = urms**2 * SQRT( 15.0_wp / epsilon / molecular_viscosity )
        tl        = urms**2 / epsilon                       ! in s
        lf        = 0.5 * urms**3 / epsilon                 ! in m
        tauk      = SQRT( molecular_viscosity / epsilon )                  ! in s
-       eta       = ( molecular_viscosity**3 / epsilon )**0.25             ! in m
+       eta       = ( molecular_viscosity**3 / epsilon )**0.25_wp          ! in m
        vk        = eta / tauk
 
        ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re )
-       tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk   ! in s
+       tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0_wp ) * ao ) ) * tauk   ! in s
 
        CALL fallg    ! gives winf in m/s
@@ -452 +453 @@
 !--    Calculate wr (from Aayala 2008b, page 38f)
        z   = tt / tl
-       be  = SQRT( 2.0 ) * lambda / lf
+       be  = SQRT( 2.0_wp ) * lambda / lf
        bbb = SQRT( 1.0 - 2.0 * be**2 )
        d1  = ( 1.0 + bbb ) / ( 2.0 * bbb )
@@ -504 +505 @@
              wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy   ! in m**2/s**2
              IF ( wrtur2xy < 0.0 )  wrtur2xy = 0.0
-             wrgrav2  = pi / 8.0 * ( winf(j) - winf(i) )**2
-             wrfin    = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) )   ! in m/s
+             wrgrav2  = pi / 8.0_wp * ( winf(j) - winf(i) )**2
+             wrfin    = SQRT( ( 2.0_wp / pi ) * ( wrtur2xy + wrgrav2) ) ! in m/s
 
 !
@@ -518 +519 @@
                    5.3406 * sst
              IF ( xx < 0.0 )  xx = 0.0
-             yy = 0.1886 * EXP( 20.306 / lambda_re )
+             yy = 0.1886 * EXP( 20.306_wp / lambda_re )
 
              c1_gr  =  xx / ( g / vk * tauk )**yy
 
-             ao_gr  = ao + ( pi / 8.0) * ( g / vk * tauk )**2
+             ao_gr  = ao + ( pi / 8.0_wp) * ( g / vk * tauk )**2
              fao_gr = 20.115 * SQRT( ao_gr / lambda_re )
              rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta   ! in cm
@@ -652 +653 @@
           stb   = 32.0 * rho_a * ( rho_l - rho_a) * g / (3.0 * eta * eta)
           phy   = sigma**3 * rho_a**2 / ( eta**4 * g * ( rho_l - rho_a ) )
-          py    = phy**( 1.0 / 6.0 )
+          py    = phy**( 1.0_wp / 6.0_wp )
 
        ENDIF
@@ -685 +686 @@
              ENDIF
 
-             x = LOG( 16.0 * bond * py / 3.0 )
+             x = LOG( 16.0 * bond * py / 3.0_wp )
              y = 0.0
 
@@ -695 +696 @@
 
              IF ( radclass(j) > 0.0035 )  THEN
-                winf(j) = xrey * eta / ( 2.0 * rho_a * 0.0035 )
+                winf(j) = xrey * eta / ( 2.0 * rho_a * 0.0035_wp )
              ELSE
                 winf(j) = xrey * eta / ( 2.0 * rho_a * radclass(j) )
@@ -799 +800 @@
        ALLOCATE( ira(1:radius_classes) )
        DO  j = 1, radius_classes
-          particle_radius = radclass(j) * 1.0E6
+          particle_radius = radclass(j) * 1.0E6_wp
           DO  k = 1, 15
              IF ( particle_radius < r0(k) )  THEN
@@ -822 +823 @@
              IF ( ir < 16 )  THEN
                 IF ( ir >= 2 )  THEN
-                   pp = ( ( radclass(j) * 1.0E06 ) - r0(ir-1) ) / &
+                   pp = ( ( radclass(j) * 1.0E06_wp ) - r0(ir-1) ) / &
                         ( r0(ir) - r0(ir-1) )
                    qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
@@ -930 +931 @@
 
        DO  j = 1, radius_classes
-          particle_radius = radclass(j) * 1.0E6
+          particle_radius = radclass(j) * 1.0E6_wp
           DO  k = 1, 7
              IF ( particle_radius < r0(k) )  THEN
@@ -1039 +1040 @@
        REAL(wp), DIMENSION(1:9,1:19), SAVE ::  ef = 0.0          !:
 
-       mean_rm = mean_r * 1.0E06
-       rm      = r      * 1.0E06
+       mean_rm = mean_r * 1.0E06_wp
+       rm      = r      * 1.0E06_wp
 
        IF ( first )  THEN