File: C:\NOAA\NEMS_11731\src\chem\gocart\src\Components\GOCART_GridComp\CARMA_GridComp\carma_bins.F90

1     ! Colarco, Aug. 21, 2007
2     ! To make the carma size binning more uniform throughout the code/SLOD
3     ! I extract the size binning routine and place here.
4     
5           subroutine carma_bins( NBIN, NGROUP, rhop, &
6                                  rmin, rmrat, rmassmin, &
7                                  r, dr, rmass, rmassup, rlow, rup, &
8                                  dm, vol, rc )
9     
10     !     types
11           use carma_types_mod
12     
13           implicit none
14     
15     !     Input
16           integer :: NBIN, NGROUP
17           real(kind=f), dimension(NBIN,NGROUP) :: rhop
18     
19     !     I/O
20           real(kind=f), dimension(NGROUP) :: rmin, rmrat, rmassmin
21           real(kind=f), dimension(NBIN,NGROUP) :: r, dr, &
22                                                   rmass, rmassup, &
23                                                   rlow, rup, dm, vol
24     
25     !     Output
26           integer :: rc
27     
28     !     Local
29           integer :: igrp, j
30           real(kind=f) :: vrfact, cpi
31     
32           rc = 0
33     
34     #ifdef DEBUG
35            write(*,*) '+ carma_bins'
36     #endif
37     
38     !  Set up the particle bins.
39     !  The following definition applies really to bins set up the CARMA way,
40     !  i.e., specify rmin and rmrat.
41     !  For each particle group, the mass of a particle in
42     !  bin j is <rmrat> times that in bin j-1
43     !
44     !    rmass(NBIN,NGROUP)     =  bin center mass [kg]
45     !    r(NBIN,NGROUP)         =  bin mean (volume-weighted) radius [m]
46     !    vol(NBIN,NGROUP)       =  bin center volume [m^3]
47     !    dr(NBIN,NGROUP)        =  bin width in radius space [m]
48     !    dv(NBIN,NGROUP)        =  bin width in volume space [m^3]
49     !    dm(NBIN,NGROUP)        =  bin width in mass space [kg]
50     !
51           cpi = 4._f/3._f*PI
52     
53     !     We have two methods for constructing the particle size bins
54     !     If rmin provided is negative then must provide valid r, rlow, rup
55     !     If r provided is negative then must provide valid rmin, rmrat
56     !     We only check here to see if both are negative.
57           if(rmin(NGROUP) .lt. 0._f .and. r(NBIN,NGROUP) .lt. 0._f) then
58            rc = 1
59            return
60           endif
61     
62           if(r(NBIN,NGROUP) .lt. 0._f) then
63     
64            do igrp = 1, NGROUP
65     
66             rmassmin(igrp) = cpi*rhop(1,igrp)*rmin(igrp)**3
67     
68             vrfact = ( ( 3._f / 2._f / PI / ( rmrat(igrp) + 1._f ) ) &
69                        ** ( 1._f / 3._f ) ) & 
70                      * ( rmrat(igrp) ** ( 1._f / 3._f ) - 1._f )
71     
72             do j = 1, NBIN
73     
74               rmass(j,igrp)   = rmassmin(igrp) * rmrat(igrp)**(j-1)
75               rmassup(j,igrp) = 2._f*rmrat(igrp)/(rmrat(igrp)+1._f)* &
76                                 rmass(j,igrp)
77               dm(j,igrp)      = 2._f*(rmrat(igrp)-1._f)/(rmrat(igrp)+1._f)* &
78                                 rmass(j,igrp)
79               vol(j,igrp) = rmass(j,igrp) / rhop(1,igrp)
80               r(j,igrp)   = ( rmass(j,igrp)/rhop(1,igrp)/cpi )**(1._f/3._f)
81               rup(j,igrp) = ( rmassup(j,igrp)/rhop(1,igrp)/cpi )**  &
82                                                                     (1._f/3._f)
83               dr(j,igrp)  = vrfact*(rmass(j,igrp)/rhop(1,igrp))**(1._f/3._f)
84               rlow(j,igrp) = rup(j,igrp) - dr(j,igrp)
85     
86             enddo
87            enddo
88     
89           else   ! supplied the radius, rlow, rup
90     
91            rmin(:) = r(1,:)
92            if(NBIN .gt. 1) then
93             rmrat(:) = (r(2,:)/r(1,:)) ** 3
94            else
95             rmrat(:) = 2._f
96            endif
97            dr(:,:) = rup(:,:) - rlow(:,:)
98            rmassmin(:) = cpi*rhop(1,:)*r(1,:)**3
99            rmass(:,:) = cpi*rhop(:,:)*r(:,:)**3
100            vol(:,:) = rmass(:,:) / rhop(:,:)
101            rmassup(:,:) = cpi*rhop(:,:)*rup(:,:)**3
102            dm(:,:) = rmassup(:,:) - cpi*rhop(:,:)*rlow(:,:)**3
103           endif
104     
105           return
106           end subroutine carma_bins
107