Line data    Source code

       1             : ! Include shortname defintions, so that the pre-existing F77 code does not have to be modified to
       2             : ! reference the CARMA structure.
       3             : #include "carma_globaer.h"
       4             : 
       5             : !! This routine calculates aerosol hygroscopicity
       6             : !!
       7             : !! @author  Pengfei Yu, Mike Mills
       8             : !! @version Oct 2020
       9     1050624 : subroutine hygroscopicity(carma, cstate, rc)
      10             :   ! types
      11             :   use carma_precision_mod
      12             :   use carma_enums_mod
      13             :   use carma_constants_mod
      14             :   use carma_types_mod
      15             :   use carmastate_mod
      16             :   use carma_mod
      17             :   use carmaelement_mod
      18             : 
      19             :   implicit none
      20             : 
      21             :   type(carma_type), intent(in)         :: carma   !! the carma object
      22             :   type(carmastate_type), intent(inout) :: cstate  !! the carma state object
      23             :   integer, intent(inout) :: rc      !! return code, negative indicates failure
      24             : 
      25             :   ! Local declarations
      26             :   integer         :: igroup         !! group index
      27             :   integer         :: ibin           !! bin index
      28             :   integer         :: iepart         !! element in group containing the particle concentration
      29             :   integer         :: icore, i, z
      30             :   real(kind=f)    :: coremass, shellmass
      31             :   real(kind=f), parameter :: thresh  = 1e-14_f     ! relative threshold for core mass and kappa checking
      32             :   real(kind=f), parameter :: thresh1 = 1._f+thresh
      33             :   real(kind=f), parameter :: thresh0 = 0._f-thresh
      34             : 
      35             :   1 format('hygroscopicity::ibin=',i4,',core mass=',e10.3,',shell mass=',e10.3,',hygroscopicity=',f9.4)
      36             : 
      37  1389975552 :   kappahygro(:NZ,:NBIN,:NGROUP) = -huge(1._f) ! default
      38             : 
      39             :   ! loop through all bins, groups, and elements to calculate hygro for each group:
      40     3151872 :   do igroup = 1,NGROUP
      41             :     ! Only calculate hygro for groups that use it
      42     3151872 :     if (irhswell(igroup) == I_PETTERS) then
      43     1050624 :       iepart = ienconc(igroup)     ! element of particle number concentration
      44    22063104 :       do ibin = 1, NBIN
      45   694462464 :         do z = 1, NZ
      46             :         
      47   672399360 :           kappahygro(z,ibin,igroup) = 0._f
      48             :           
      49   672399360 :           if (pc(z, ibin, iepart).gt.0._f) then
      50             :           
      51             :             ! Weight hygro by mass of each core
      52   672346747 :             coremass = 0._f
      53  6723467470 :             do i = 1, ncore(igroup)
      54  6051120723 :               icore = icorelem(i, igroup)
      55  6051120723 :               coremass = coremass + pc(z, ibin, icore)
      56  6723467470 :               kappahygro(z,ibin,igroup) = kappahygro(z,ibin,igroup) + pc(z,ibin,icore) * kappaelem(icore)
      57             :             end do ! i = 1, ncore(igroup)
      58             : 
      59             :             ! Add shell mass to hygro weighting
      60             :             !
      61             :             ! NOTE: Check for coremass being to big for PC and adjust
      62             :             ! accordingly.
      63   672346747 :             shellmass = (pc(z, ibin, iepart) * rmass(ibin, igroup)) - coremass
      64   672346747 :             if (shellmass < 0._f) then
      65      219210 :               shellmass = 0._f
      66      219210 :               pc(z, ibin, iepart) = coremass / rmass(ibin, igroup)
      67             :             else 
      68   672127537 :               kappahygro(z,ibin,igroup) = kappahygro(z,ibin,igroup) + shellmass * kappaelem(iepart)
      69             :             end if
      70             : 
      71             :             ! Divide by total mass of all particles in the bin to normalize:
      72   672346747 :             kappahygro(z,ibin,igroup) = kappahygro(z,ibin,igroup) / pc(z, ibin, iepart) / rmass(ibin, igroup)
      73             :           end if
      74             :           
      75   672399360 :           if (kappahygro(z,ibin,igroup).gt.thresh1.or.kappahygro(z,ibin,igroup).lt.thresh0) then
      76           0 :             write(LUNOPRT,*) "hygro77: z,ibin,kappahygro,pc,rmass,shellmass,coremass", &
      77           0 :               z,ibin,kappahygro(z,ibin,igroup),pc(z, ibin, iepart),rmass(ibin, igroup),shellmass,coremass
      78           0 :             rc=RC_ERROR
      79           0 :             return
      80             :           end if
      81             :           
      82   672399360 :           kappahygro(z,ibin,igroup) = min(kappahygro(z,ibin,igroup),1._f)
      83   693411840 :           kappahygro(z,ibin,igroup) = max(kappahygro(z,ibin,igroup),0._f)
      84             :         end do ! z = 1, NZ
      85             :       end do ! ibin = 1, NBIN
      86             :     end if ! irhswell(igroup) == I_PETTERS
      87             :   end do ! igroup = 1,NGROUP
      88             : 
      89     1050624 :   rc = RC_OK
      90             : 
      91     1050624 :   return
      92     1050624 : end subroutine !hygroscopicity