Home


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

Ignore:
Timestamp:
Sep 19, 2012 2:30:36 PM (12 years ago)
Author:
franke
Message:

Bugfixes

Missing calculation of mean particle weighting factor for output added. (data_output_2d, data_output_3d, data_output_mask, sum_up_3d_data)
Calculation of mean particle radius for output now considers the weighting factor. (data_output_mask)
Calculation of sugrid-scale buoyancy flux for humidity and cloud droplets corrected. (flow_statistics)
Factor in calculation of enhancement factor for collision efficencies corrected. (lpm_collision_kernels)
Calculation of buoyancy production now considers the liquid water mixing ratio in case of cloud droplets. (production_e)

Changes

Calculation of buoyancy flux for humidity in case of WS-scheme is now using turbulent fluxes of WS-scheme. (flow_statistics)
Calculation of the collision kernels now in SI units. (lpm_collision_kernels)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/lpm_collision_kernels.f90

    r850 r1007  
    44! Current revisions:
    55! -----------------
    6 !
     6! converted all units to SI units and replaced some parameters by corresponding
     7! PALM parameters
     8! Bugfix: factor in calculation of enhancement factor for collision efficencies
     9! changed from 10. to 1.0
    710!
    811! Former revisions:
    … …  
    6467
    6568    PUBLIC  ckernel, collision_efficiency_rogers, init_kernels, &
    66             rclass_lbound, rclass_ubound, recalculate_kernel           
     69            rclass_lbound, rclass_ubound, recalculate_kernel
    6770
    6871    REAL ::  epsilon, eps2, rclass_lbound, rclass_ubound, urms, urms2
    … …  
    126129!             ENDIF
    127130          ENDDO
    128 !
    129 !--       Collision routines expect radius to be in cm
    130           radclass = radclass * 100.0
    131 
    132 !
    133 !--       Set the class bounds for dissipation in interval [0, 1000] cm**2/s**3
     131
     132!
     133!--       Set the class bounds for dissipation in interval [0.0, 0.1] m**2/s**3
    134134          DO  i = 1, dissipation_classes
    135              epsclass(i) = 1000.0 * REAL( i ) / dissipation_classes
     135             epsclass(i) = 0.1 * REAL( i ) / dissipation_classes
    136136!             IF ( myid == 0 )  THEN
    137137!                PRINT*, 'i=', i, ' eps = ', epsclass(i)
    … …  
    148148
    149149             epsilon = epsclass(k)
    150              urms    = 202.0 * ( epsilon / 400.0 )**( 1.0 / 3.0 )
     150             urms    = 2.02 * ( epsilon / 0.04 )**( 1.0 / 3.0 )
    151151
    152152             CALL turbsd
    … …  
    183183
    184184             PRINT*, '*** Hall kernel'
    185              WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
     185             WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6, &
     186                                              i = 1,radius_classes )
    186187             DO  j = 1, radius_classes
    187                 WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j), ( hkernel(i,j), i = 1,radius_classes )
     188                WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j),  &
     189                                          ( hkernel(i,j), i = 1,radius_classes )
    188190             ENDDO
    189191
    … …  
    200202
    201203                PRINT*, '*** epsilon = ', epsclass(k)
    202                 WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
     204                WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6, &
     205                                                 i = 1,radius_classes )
    203206                DO  j = 1, radius_classes
    204 !                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( ckernel(i,j,k), i = 1,radius_classes )
    205                    WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( hwratio(i,j), i = 1,radius_classes )
     207!                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E6, &
     208!                                       ( ckernel(i,j,k), i = 1,radius_classes )
     209                   WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j)*1.0E6, &
     210                                          ( hwratio(i,j), i = 1,radius_classes )
    206211                ENDDO
    207212             ENDDO
    … …  
    210215
    211216          DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf )
    212 
    213           ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
    214217
    215218       ELSEIF( collision_kernel == 'hall'  .OR.  collision_kernel == 'wang' ) &
    … …  
    249252
    250253!
    251 !--    Store particle radii on the radclass array. Collision routines
    252 !--    expect radii to be in cm.
    253        radclass(1:radius_classes) = particles(pstart:pend)%radius * 100.0
     254!--    Store particle radii on the radclass array
     255       radclass(1:radius_classes) = particles(pstart:pend)%radius
    254256
    255257       IF ( wang_kernel )  THEN
    256           epsilon = diss(k1,j1,i1) * 1.0E4   ! dissipation rate in cm**2/s**-3
     258          epsilon = diss(k1,j1,i1)   ! dissipation rate in m**2/s**3
    257259       ELSE
    258260          epsilon = 0.0
    259261       ENDIF
    260        urms    = 202.0 * ( epsilon / 400.0 )**( 0.33333333333 )
    261 
    262        IF ( wang_kernel  .AND.  epsilon > 0.001 )  THEN
     262       urms    = 2.02 * ( epsilon / 0.04 )**( 0.33333333333 )
     263
     264       IF ( wang_kernel  .AND.  epsilon > 1.0E-7 )  THEN
    263265!
    264266!--       Call routines to calculate efficiencies for the Wang kernel
    … …  
    294296       ENDIF
    295297
    296        ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
    297 
    298298       DEALLOCATE( ec, radclass, winf )
    299299
    … …  
    314314       USE particle_attributes
    315315       USE arrays_3d
     316       USE control_parameters
    316317
    317318       IMPLICIT NONE
    … …  
    326327                v1xysq, v2, v2xysq, wrfin, wrgrav2, wrtur2xy, xx, yy, z
    327328
    328        REAL, SAVE ::  airdens, airvisc, anu, gravity, waterdens
    329 
    330329       REAL, DIMENSION(1:radius_classes) ::  st, tau
    331330
    … …  
    336335
    337336          first = .FALSE.
    338           airvisc   = 0.1818     ! dynamic viscosity in mg/cm*s
    339           airdens   = 1.2250     ! air density in mg/cm**3
    340           waterdens = 1000.0     ! water density in mg/cm**3
    341           gravity   = 980.6650   ! in cm/s**2
    342           anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
    343 
    344        ENDIF   
    345 
    346        lambda    = urms * SQRT( 15.0 * anu / epsilon )     ! in cm
    347        lambda_re = urms**2 * SQRT( 15.0 / epsilon / anu )
     337
     338       ENDIF
     339
     340       lambda    = urms * SQRT( 15.0 * molecular_viscosity / epsilon )    ! in m
     341       lambda_re = urms**2 * SQRT( 15.0 / epsilon / molecular_viscosity )
    348342       tl        = urms**2 / epsilon                       ! in s
    349        lf        = 0.5 * urms**3 / epsilon                 ! in cm
    350        tauk      = SQRT( anu / epsilon )                   ! in s
    351        eta       = ( anu**3 / epsilon )**0.25              ! in cm
    352        vk        = eta / tauk                              ! in cm/s
     343       lf        = 0.5 * urms**3 / epsilon                 ! in m
     344       tauk      = SQRT( molecular_viscosity / epsilon )                  ! in s
     345       eta       = ( molecular_viscosity**3 / epsilon )**0.25             ! in m
     346       vk        = eta / tauk
    353347
    354348       ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re )
    355349       tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk   ! in s
    356350
    357        CALL fallg    ! gives winf in cm/s
     351       CALL fallg    ! gives winf in m/s
    358352
    359353       DO  i = 1, radius_classes
    360           tau(i) = winf(i) / gravity    ! in s
     354          tau(i) = winf(i) / g    ! in s
    361355          st(i)  = tau(i) / tauk
    362356       ENDDO
    … …  
    368362       bbb = SQRT( 1.0 - 2.0 * be**2 )
    369363       d1  = ( 1.0 + bbb ) / ( 2.0 * bbb )
    370        e1  = lf * ( 1.0 + bbb ) * 0.5   ! in cm
     364       e1  = lf * ( 1.0 + bbb ) * 0.5   ! in m
    371365       d2  = ( 1.0 - bbb ) * 0.5 / bbb
    372        e2  = lf * ( 1.0 - bbb ) * 0.5   ! in cm
     366       e2  = lf * ( 1.0 - bbb ) * 0.5   ! in m
    373367       ccc = SQRT( 1.0 - 2.0 * z**2 )
    374368       b1  = ( 1.0 + ccc ) * 0.5 / ccc
    … …  
    379373       DO  i = 1, radius_classes
    380374
    381           v1 = winf(i)        ! in cm/s
     375          v1 = winf(i)        ! in m/s
    382376          t1 = tau(i)         ! in s
    383377
    384378          DO  j = 1, i
    385              rrp = radclass(i) + radclass(j)               ! radius in cm
    386              v2  = winf(j)                                 ! in cm/s
     379             rrp = radclass(i) + radclass(j)
     380             v2  = winf(j)                                 ! in m/s
    387381             t2  = tau(j)                                  ! in s
    388382
    389              v1xysq  = b1 * d1 * phi(c1,e1,v1,t1) - b1 * d2 * phi(c1,e2,v1,t1) &
    390                      - b2 * d1 * phi(c2,e1,v1,t1) + b2 * d2 * phi(c2,e2,v1,t1)
    391              v1xysq  = v1xysq * urms**2 / t1                ! in cm**2/s**2
    392              vrms1xy = SQRT( v1xysq )                       ! in cm/s
    393 
    394              v2xysq  = b1 * d1 * phi(c1,e1,v2,t2) - b1 * d2 * phi(c1,e2,v2,t2) &
    395                      - b2 * d1 * phi(c2,e1,v2,t2) + b2 * d2 * phi(c2,e2,v2,t2)
    396              v2xysq  = v2xysq * urms**2 / t2                ! in cm**2/s**2
    397              vrms2xy = SQRT( v2xysq )                       ! in cm/s
     383             v1xysq  = b1 * d1 * phi_w(c1,e1,v1,t1) - b1 * d2 * phi_w(c1,e2,v1,t1) &
     384                     - b2 * d1 * phi_w(c2,e1,v1,t1) + b2 * d2 * phi_w(c2,e2,v1,t1)
     385             v1xysq  = v1xysq * urms**2 / t1                ! in m**2/s**2
     386             vrms1xy = SQRT( v1xysq )                       ! in m/s
     387
     388             v2xysq  = b1 * d1 * phi_w(c1,e1,v2,t2) - b1 * d2 * phi_w(c1,e2,v2,t2) &
     389                     - b2 * d1 * phi_w(c2,e1,v2,t2) + b2 * d2 * phi_w(c2,e2,v2,t2)
     390             v2xysq  = v2xysq * urms**2 / t2                ! in m**2/s**2
     391             vrms2xy = SQRT( v2xysq )                       ! in m/s
    398392
    399393             IF ( winf(i) >= winf(j) )  THEN
    … …  
    414408                         b2 * d2* zhi(c2,e2,v1,t1,v2,t2)
    415409             fr       = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 )
    416              v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)  ! in cm**2/s**2
    417              wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy   ! in cm**2/s**2
     410             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)  ! in m**2/s**2
     411             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy   ! in m**2/s**2
    418412             IF ( wrtur2xy < 0.0 )  wrtur2xy = 0.0
    419413             wrgrav2  = pi / 8.0 * ( winf(j) - winf(i) )**2
    420              wrfin    = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) )   ! in cm/s
     414             wrfin    = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) )   ! in m/s
    421415
    422416!
    … …  
    433427             yy = 0.1886 * EXP( 20.306 / lambda_re )
    434428
    435              c1_gr  =  xx / ( gravity / vk * tauk )**yy
    436 
    437              ao_gr  = ao + ( pi / 8.0) * ( gravity / vk * tauk )**2
     429             c1_gr  =  xx / ( g / vk * tauk )**yy
     430
     431             ao_gr  = ao + ( pi / 8.0) * ( g / vk * tauk )**2
    438432             fao_gr = 20.115 * SQRT( ao_gr / lambda_re )
    439433             rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta   ! in cm
    … …  
    452446
    453447!------------------------------------------------------------------------------!
    454 ! phi as a function
    455 !------------------------------------------------------------------------------!
    456     REAL FUNCTION phi( a, b, vsett, tau0 )
     448! phi_w as a function
     449!------------------------------------------------------------------------------!
     450    REAL FUNCTION phi_w( a, b, vsett, tau0 )
    457451
    458452       IMPLICIT NONE
    … …  
    461455
    462456       aa1 = 1.0 / tau0 + 1.0 / a + vsett / b
    463        phi = 1.0 / aa1  - 0.5 * vsett / b / aa1**2  ! in s
    464 
    465     END FUNCTION phi
    466 
    467 
    468 !------------------------------------------------------------------------------!
    469 ! zeta as a function
     457       phi_w = 1.0 / aa1  - 0.5 * vsett / b / aa1**2  ! in s
     458
     459    END FUNCTION phi_w
     460
     461
     462!------------------------------------------------------------------------------!
     463! zhi as a function
    470464!------------------------------------------------------------------------------!
    471465    REAL FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 )
    … …  
    490484
    491485!------------------------------------------------------------------------------!
    492 ! Calculation of terminal velocity winf
     486! Calculation of terminal velocity winf following Equations 10-138 to 10-145
     487! from (Pruppacher and Klett, 1997)
    493488!------------------------------------------------------------------------------!
    494489    SUBROUTINE fallg
    … …  
    498493       USE particle_attributes
    499494       USE arrays_3d
     495       USE control_parameters
    500496
    501497       IMPLICIT NONE
    … …  
    505501       LOGICAL, SAVE ::  first = .TRUE.
    506502
    507        REAL, SAVE ::  cunh, eta, grav, phy, py, rhoa, rhow, sigma, stb, stok, &
     503       REAL, SAVE ::  cunh, eta, phy, py, rho_a, sigma, stb, stok, &
    508504                      t0, xlamb
    509505
    … …  
    523519                 -0.542819E-1, 0.238449E-2 /)
    524520
    525           eta   = 1.818E-4          ! in poise = g/(cm s)
    526           xlamb = 6.62E-6           ! in cm
    527           rhow  = 1.0               ! in g/cm**3
    528           rhoa  = 1.225E-3          ! in g/cm**3
    529           grav  = 980.665           ! in cm/s**2
    530           cunh  = 1.257 * xlamb     ! in cm
    531           t0    = 273.15            ! in K
    532           sigma = 76.1 - 0.155 * ( 293.15 - t0 )              ! in N/m = g/s**2
    533           stok  = 2.0  * grav * ( rhow - rhoa ) / ( 9.0 * eta ) ! in 1/(cm s)
    534           stb   = 32.0 * rhoa * ( rhow - rhoa) * grav / (3.0 * eta * eta)
    535                                                                 ! in 1/cm**3
    536           phy   = sigma**3 * rhoa**2 / ( eta**4 * grav * ( rhow - rhoa ) )
     521!
     522!--       Parameter values for p = 1013,25 hPa and T = 293,15 K
     523          eta   = 1.818E-5         ! in kg/(m s)
     524          xlamb = 6.6E-8           ! in m
     525          rho_a = 1.204            ! in kg/m**3
     526          cunh  = 1.26 * xlamb     ! in m
     527          sigma = 0.07363          ! in kg/s**2
     528          stok  = 2.0  * g * ( rho_l - rho_a ) / ( 9.0 * eta ) ! in 1/(m s)
     529          stb   = 32.0 * rho_a * ( rho_l - rho_a) * g / (3.0 * eta * eta)
     530          phy   = sigma**3 * rho_a**2 / ( eta**4 * g * ( rho_l - rho_a ) )
    537531          py    = phy**( 1.0 / 6.0 )
    538532
    … …  
    541535       DO  j = 1, radius_classes
    542536
    543           IF ( radclass(j) <= 1.0E-3 ) THEN
    544 
    545              winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) ) ! in cm/s
    546 
    547           ELSEIF ( radclass(j) > 1.0E-3  .AND.  radclass(j) <= 5.35E-2 )  THEN
     537          IF ( radclass(j) <= 1.0E-5 ) THEN
     538
     539             winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) )
     540
     541          ELSEIF ( radclass(j) > 1.0E-5  .AND.  radclass(j) <= 5.35E-4 )  THEN
    548542
    549543             x = LOG( stb * radclass(j)**3 )
    … …  
    553547                y = y + b(i) * x**(i-1)
    554548             ENDDO
    555 
    556              xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y )
    557              winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )      ! in cm/s
    558 
    559           ELSEIF ( radclass(j) > 5.35E-2 )  THEN
    560 
    561              IF ( radclass(j) > 0.35 )  THEN
    562                 bond = grav * ( rhow - rhoa ) * 0.35**2 / sigma
     549!
     550!--          Note: this Eq. is wrong in (Pruppacher and Klett, 1997, p. 418)
     551!--          for correct version see (Beard, 1976)
     552             xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y )
     553
     554             winf(j) = xrey * eta / ( 2.0 * rho_a * radclass(j) )
     555
     556          ELSEIF ( radclass(j) > 5.35E-4 )  THEN
     557
     558             IF ( radclass(j) > 0.0035 )  THEN
     559                bond = g * ( rho_l - rho_a ) * 0.0035**2 / sigma
    563560             ELSE
    564                 bond = grav * ( rhow - rhoa ) * radclass(j)**2 / sigma
     561               bond = g * ( rho_l - rho_a ) * radclass(j)**2 / sigma
    565562             ENDIF
    566563
    … …  
    574571             xrey = py * EXP( y )
    575572
    576              IF ( radclass(j) > 0.35 )  THEN
    577                 winf(j) = xrey * eta / ( 2.0 * rhoa * 0.35 )           ! in cm/s
     573             IF ( radclass(j) > 0.0035 )  THEN
     574                winf(j) = xrey * eta / ( 2.0 * rho_a * 0.0035 )
    578575             ELSE
    579                 winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )    ! in cm/s
     576                winf(j) = xrey * eta / ( 2.0 * rho_a * radclass(j) )
    580577             ENDIF
    581578
    … …  
    668665!
    669666!--    Calculate the radius class index of particles with respect to array r
     667!--    Radius has to be in µm
    670668       ALLOCATE( ira(1:radius_classes) )
    671669       DO  j = 1, radius_classes
    672           particle_radius = radclass(j) * 1.0E4   ! in microm
     670          particle_radius = radclass(j) * 1.0E6
    673671          DO  k = 1, 15
    674672             IF ( particle_radius < r0(k) )  THEN
    … …  
    693691             IF ( ir < 16 )  THEN
    694692                IF ( ir >= 2 )  THEN
    695                    pp = ( ( radclass(j) * 1.0E04 ) - r0(ir-1) ) / &
     693                   pp = ( ( radclass(j) * 1.0E06 ) - r0(ir-1) ) / &
    696694                        ( r0(ir) - r0(ir-1) )
    697695                   qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    … …  
    754752          rat = (/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /)
    755753!
    756 !--       In 100 cm^2/s^3
     754!--       for 100 cm**2/s**3
    757755          ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
    758756          ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
    … …  
    767765          ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
    768766!
    769 !--       In 400 cm^2/s^3
     767!--       for 400 cm**2/s**3
    770768          ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
    771769          ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
    … …  
    784782!
    785783!--    Calculate the radius class index of particles with respect to array r0
     784!--    Radius has to be in µm
    786785       ALLOCATE( ira(1:radius_classes) )
    787786
    788787       DO  j = 1, radius_classes
    789           particle_radius = radclass(j) * 1.0E4    ! in microm
     788          particle_radius = radclass(j) * 1.0E6
    790789          DO  k = 1, 7
    791790             IF ( particle_radius < r0(k) )  THEN
    … …  
    799798!
    800799!--    Two-dimensional linear interpolation of the collision efficiencies
     800!--    Radius has to be in µm
    801801       DO  j =  1, radius_classes
    802802          DO  i = 1, j
    … …  
    812812             ENDDO
    813813
    814              y1 = 1.0      ! 0  cm2/s3
    815 !
    816 !--          100 cm2/s3, 400 cm2/s3
     814             y1 = 0.0001      ! for 0 m**2/s**3
     815
    817816             IF ( ir < 8 )  THEN
    818817                IF ( ir >= 2 )  THEN
    819                    pp = ( radclass(j)*1.0E4 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
     818                   pp = ( radclass(j)*1.0E6 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
    820819                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    821820                   y2 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_100(ir-1,iq-1) +  &
    822821                                pp * ( 1.0-qq ) * ecoll_100(ir,iq-1)   +  &
    823                                 qq * ( 10.-pp ) * ecoll_100(ir-1,iq)   +  &
     822                                qq * ( 1.0-pp ) * ecoll_100(ir-1,iq)   +  &
    824823                                pp * qq         * ecoll_100(ir,iq)
    825824                   y3 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_400(ir-1,iq-1) +  &
    … …  
    838837             ENDIF
    839838!
    840 !--          Linear interpolation of dissipation rate in cm2/s3
    841              IF ( epsilon <= 100.0 )  THEN
    842                 ecf(j,i) = ( epsilon - 100.0 ) / (   0.0 - 100.0 ) * y1 &
    843                          + ( epsilon -   0.0 ) / ( 100.0 -   0.0 ) * y2
    844              ELSEIF ( epsilon <= 600.0 )  THEN
    845                 ecf(j,i) = ( epsilon - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
    846                          + ( epsilon - 100.0 ) / ( 400.0 - 100.0 ) * y3
     839!--          Linear interpolation of dissipation rate in m**2/s**3
     840             IF ( epsilon <= 0.01 )  THEN
     841                ecf(j,i) = ( epsilon - 0.01 ) / (   0.0 - 0.01 ) * y1 &
     842                         + ( epsilon -   0.0 ) / ( 0.01 -   0.0 ) * y2
     843             ELSEIF ( epsilon <= 0.06 )  THEN
     844                ecf(j,i) = ( epsilon - 0.04 ) / ( 0.01 - 0.04 ) * y2 &
     845                         + ( epsilon - 0.01 ) / ( 0.04 - 0.01 ) * y3
    847846             ELSE
    848                 ecf(j,i) = (   600.0 - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
    849                          + (   600.0 - 100.0 ) / ( 400.0 - 100.0 ) * y3
     847                ecf(j,i) = (   0.06 - 0.04 ) / ( 0.01 - 0.04 ) * y2 &
     848                         + (   0.06 - 0.01 ) / ( 0.04 - 0.01 ) * y3
    850849             ENDIF
    851850
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |