Line data    Source code

       1             : ! Include shortname defintions, so that the F77 code does not have to be modified to
       2             : ! reference the CARMA structure.
       3             : #include "carma_globaer.h"
       4             : 
       5             : !!  This routine evaluates the derived mapping arrays and sets up
       6             : !!  the particle size bins.
       7             : !!
       8             : !!  @author Eric Jensen
       9             : !!  @ version Oct-1995
      10        1536 : subroutine setupbins(carma, rc)
      11             : 
      12             :   ! types
      13             :   use carma_precision_mod
      14             :   use carma_enums_mod
      15             :   use carma_constants_mod
      16             :   use carma_types_mod
      17             :   use carma_mod
      18             : 
      19             :   implicit none
      20             : 
      21             :   type(carma_type), intent(inout) :: carma   !! the carma object
      22             :   integer, intent(inout)          :: rc      !! return code, negative indicates failure
      23             : 
      24             :   ! Local declarations
      25             :   integer :: ielem, ibin, i, j, ix, iy, iz, ie, ig, ip, igrp, jgrp
      26        3072 :   real(kind=f)  :: tmp_rhop(NBIN, NGROUP)
      27             :   real(kind=f)  :: vrfact
      28             :   real(kind=f)  :: cpi
      29             :   ! Local declarations needed for creation of fractal bin structure
      30             :   real(kind=f)  :: rf, rp
      31             :   real(kind=f)  :: vpor, upor, gamma, happel, perm, brinkman, epsil, omega
      32             : 
      33             :   !  Define formats
      34             :   !
      35             :   1 format(a,':  ',12i6)
      36             :   2 format(a,':  ',i6)
      37             :   3 format(a,':  ',f12.2)
      38             :   4 format(a,':  ',12f12.2)
      39             :   5 format(/,'Particle grid structure (setupbins):')
      40             :   6 format(a,':  ',1p12e12.3)
      41             :   7 format(a,':  ',12l6)
      42             : 
      43             : 
      44             :   !  Determine which elements are particle number concentrations
      45             :   !  <ienconc(igroup)> is the element corresponding to particle number 
      46             :   !  concentration in group <igroup>
      47             :   !
      48        1536 :   igrp = 0
      49       12288 :   do ielem = 1, NELEM
      50       10752 :     if( itype(ielem) .eq. I_INVOLATILE .or. &
      51        1536 :         itype(ielem) .eq. I_VOLATILE )then
      52             : 
      53        3072 :       igrp = igrp + 1
      54        3072 :       ienconc(igrp) = ielem
      55             :     endif
      56             :   enddo
      57             :   
      58        1536 :   if( igrp .gt. NGROUP )then
      59           0 :     if (do_print) write(LUNOPRT,'(/,a)') 'CARMA_setupbin:: ERROR - bad itype array'
      60           0 :     rc = -1
      61           0 :     return
      62             :   endif
      63             : 
      64             :   !  Determine which group each element belongs to
      65             :   !  i.e., <igelem(ielem)> is the group to which element <ielem> belongs!
      66        1536 :   igrp = 0
      67       12288 :   do ielem = 1, NELEM
      68       10752 :     if( itype(ielem) .eq. I_INVOLATILE .or. &
      69             :        itype(ielem) .eq. I_VOLATILE )then
      70        3072 :       igrp = igrp + 1
      71             :     endif
      72       12288 :     igelem(ielem) = igrp
      73             :   enddo
      74             : 
      75             :   !  Determine how many cores are in each group <ncore>.
      76             :   !  The core elements in a group are given by <icorelem(1:ncore,igroup)>.
      77             :   !
      78             :   !  Also evaluate whether or not second moment is used <if_sec_mom> for each group.
      79        1536 :   ielem = 0
      80             :   
      81        4608 :   do igrp = 1, NGROUP
      82             :   
      83        3072 :     ncore(igrp) = 0
      84        3072 :     if_sec_mom(igrp) = .false.
      85        3072 :     imomelem(igrp) = 0
      86             :   
      87       15360 :     do j = 1, nelemg(igrp)
      88             :   
      89       10752 :       ielem = ielem + 1
      90             :   
      91       10752 :       if( itype(ielem) .eq. I_COREMASS .or. &
      92        3072 :           itype(ielem) .eq. I_VOLCORE )then
      93             :   
      94        7680 :         ncore(igrp) = ncore(igrp) + 1
      95        7680 :         icorelem(ncore(igrp),igrp) = ielem
      96             :   
      97        3072 :       elseif( itype(ielem) .eq. I_CORE2MOM )then
      98             :   
      99           0 :         if_sec_mom(igrp) = .true.
     100           0 :         imomelem(igrp) = ielem
     101             :   
     102             :       endif
     103             :   
     104             :     enddo
     105             :   enddo
     106             : 
     107             :   !  Particle mass densities (NBIN for each group) -- the user might want
     108             :   !  to modify this (this code segment does not appear in setupaer subroutine
     109             :   !  because <igelem> is not defined until this subroutine).
     110        4608 :   do ig = 1,NGROUP
     111        3072 :     ie = ienconc(ig)
     112       66048 :     do ibin = 1,NBIN
     113       64512 :       tmp_rhop(ibin, ig) = rhoelem(ibin, ie)
     114             : 
     115             :         !  Set initial density of all hydrometeor groups to 1 such that nucleation
     116             :         !  mapping arrays are calculated correctly.
     117             :         !  or not
     118             : !           if( itype(ie) .ne. I_INVOLATILE ) then
     119             : !             rhop3(ixyz,ibin,ig) = 1.
     120             : !           endif
     121             :     enddo
     122             :   enddo
     123             :   
     124             :   !  Set up the particle bins.
     125             :   !  For each particle group, the mass of a particle in
     126             :   !  bin j is <rmrat> times that in bin j-1
     127             :   !
     128             :   !    rmass(NBIN,NGROUP)     =  bin center mass [g]
     129             :   !    r(NBIN,NGROUP)         =  bin mean (volume-weighted) radius [cm]
     130             :   !    vol(NBIN,NGROUP)       =  bin center volume [cm^3]
     131             :   !    dr(NBIN,NGROUP)        =  bin width in radius space [cm]
     132             :   !    dv(NBIN,NGROUP)        =  bin width in volume space [cm^3]
     133             :   !    dm(NBIN,NGROUP)        =  bin width in mass space [g]
     134        1536 :   cpi = 4._f/3._f*PI
     135             : 
     136        4608 :   do igrp = 1, NGROUP
     137             : 
     138           0 :     vrfact = ( (3._f/2._f/PI/(rmrat(igrp)+1._f))**(ONE/3._f) )* &
     139        3072 :             ( rmrat(igrp)**(ONE/3._f) - 1._f )
     140             : 
     141             :     ! If rmassmin wasn't specified, then use rmin to determine the mass
     142             :     ! of the first bin.
     143        3072 :     if (rmassmin(igrp) == 0._f) then
     144        3072 :       rmassmin(igrp) = cpi*tmp_rhop(1,igrp)*rmin(igrp)**3
     145             :     else
     146             :       
     147             :       ! Just for internal consistency, recalculate rmin based on the rmass
     148             :       ! that is being used.
     149             :       rmin(igrp) = (rmassmin(igrp) / cpi / tmp_rhop(1,igrp)) ** (1._f / 3._f)
     150           0 :     end if
     151             :     
     152             :     do j = 1, NBIN
     153       66048 :       rmass(j,igrp)   = rmassmin(igrp) * rmrat(igrp)**(j-1)
     154       61440 :       rmassup(j,igrp) = 2._f*rmrat(igrp)/(rmrat(igrp)+1._f)*rmass(j,igrp)
     155       61440 :       dm(j,igrp)      = 2._f*(rmrat(igrp)-1._f)/(rmrat(igrp)+1._f)*rmass(j,igrp)
     156       61440 :       vol(j,igrp) = rmass(j,igrp) / tmp_rhop(j,igrp)
     157       61440 :       r(j,igrp)   = ( rmass(j,igrp)/tmp_rhop(j,igrp)/cpi )**(ONE/3._f)
     158       61440 :       rup(j,igrp) = ( rmassup(j,igrp)/tmp_rhop(j,igrp)/cpi )**(ONE/3._f)
     159       61440 :       dr(j,igrp)  = vrfact*(rmass(j,igrp)/tmp_rhop(j,igrp))**(ONE/3._f)
     160       61440 :       rlow(j,igrp) = rup(j,igrp) - dr(j,igrp)
     161       61440 :  
     162             :       if (is_grp_fractal(igrp)) then
     163       64512 :       ! fractal flag is true
     164             : 
     165             :         if (r(j,igrp) .le. rmon(igrp)) then   ! if the bin radius is less than the monomer size
     166           0 :                                      
     167             :           nmon(j,igrp) = 1.0_f
     168           0 :           rrat(j,igrp) = 1.0_f
     169           0 :           arat(j,igrp) = 1.0_f
     170           0 :           rprat(j,igrp) = 1.0_f
     171           0 :           df(j,igrp) = 3.0_f  ! Reset fractal dimension to 3 (this is a formality)
     172           0 : 
     173             :         else   ! if bin radius is greater than the monomer size
     174             : 
     175             :           rf = (1.0_f/falpha(igrp))**(1.0_f/df(j,igrp))*r(j,igrp)**(3.0_f/df(j,igrp))*rmon(igrp)**(1.0_f-3.0_f/df(j,igrp))
     176           0 :           nmon(j,igrp) = falpha(igrp)*(rf/rmon(igrp))**df(j,igrp)
     177           0 : 
     178             :           rrat(j,igrp) = rf/r(j,igrp)           
     179           0 :                                                                                          
     180             :           ! Calculate mobility radius for permeable aggregates
     181             :           ! using Vainshtein (2003) formulation     
     182             :           vpor = 1.0_f - (nmon(j,igrp))**(1.0_f-3.0_f/df(j,igrp))           ! Volume average porosity (eq. 3.2)
     183           0 :           upor = 1.0_f-(1.0_f - vpor)*sqrt(df(j,igrp)/3.0_f)                ! Uniform poroisty (eq. 3.10)
     184           0 :           gamma = (1.0_f - upor)**(1.0_f/3.0_f)
     185           0 :           happel = 2.0_f/(9.0_f*(1.0_f-upor))*   &                          ! Happel permeability model
     186             :                   (3.0_f-4.5_f*gamma+4.5_f*gamma**5.0_f-3.0_f*gamma**6.0_f)/  &
     187             :                   (3.0_f+2.0_f*gamma**5.0_f)    
     188           0 :           perm = happel*rmon(igrp)**2.0_f                                   ! Permeability (eq. 3.3)
     189             :           brinkman = nmon(j,igrp)**(1.0_f/df(j,igrp))*1.0_f/sqrt(happel)    ! Brinkman parameter (eq. 3.9) 
     190           0 :           epsil = 1.0_f - brinkman**(-1._f)*tanh(brinkman)                    !
     191           0 :           omega = 2.0_f/3.0_f*epsil/(2.0_f/3.0_f+epsil/brinkman**2.0_f)     ! drag coefficient (eq. 2.7)
     192           0 :           rp = rf * omega
     193           0 :           rprat(j,igrp) = rp/r(j,igrp)
     194           0 : 
     195             :           arat(j,igrp) = (rprat(j,igrp) / rrat(j, igrp))**2.0_f
     196           0 :         endif
     197             :       else
     198             :          ! Not a fractal.
     199             :          nmon(j,igrp) = 1.0_f
     200       61440 :          rprat(j,igrp) = 1.0_f
     201       61440 :          df(j,igrp) = 3.0_f
     202       61440 :       endif
     203             :    enddo
     204             :   enddo
     205             :   
     206             :   !  Evaluate differences between valuse of <rmass> in different bins.
     207             :   do igrp = 1, NGROUP
     208        4608 :    do jgrp = 1, NGROUP
     209       10752 :     do i = 1, NBIN
     210      132096 :      do j = 1, NBIN
     211     2586624 :        diffmass(i,igrp,j,jgrp) = rmass(i,igrp) - rmass(j,jgrp)
     212     2580480 :      enddo
     213             :     enddo
     214             :    enddo
     215             :   enddo
     216             :   
     217             :   !  Report some initialization values
     218             :   if (do_print_init) then
     219        1536 :     write(LUNOPRT,5)
     220           2 :     write(LUNOPRT,2) 'NGROUP ',NGROUP
     221           2 :     write(LUNOPRT,2) 'NELEM  ',NELEM
     222           2 :     write(LUNOPRT,2) 'NBIN   ',NBIN
     223           2 :     write(LUNOPRT,6) 'Massmin',(rmassmin(i),i=1,NGROUP)
     224           6 :     write(LUNOPRT,4) 'Mrat   ',(rmrat(i),i=1,NGROUP)
     225           6 :     write(LUNOPRT,1) 'nelemg ',(nelemg(i),i=1,NGROUP)
     226           6 :     write(LUNOPRT,1) 'itype  ',(itype(i),i=1,NELEM)
     227          16 :     write(LUNOPRT,1) 'ienconc',(ienconc(i),i=1,NGROUP)
     228           6 :     write(LUNOPRT,1) 'igelem ',(igelem(i),i=1,NELEM)
     229          16 :     write(LUNOPRT,1) 'ncore  ',(ncore(i),i=1,NGROUP)
     230           6 :     write(LUNOPRT,7) 'fractal',(is_grp_fractal(i),i=1,NGROUP)
     231           6 :   end if
     232             :  
     233             :   !  Return to caller with particle grid initialized
     234             :   return
     235             : end