Line data    Source code

       1             : !===============================================================================
       2             : ! used to compute sea salt surface emissions for modal and sectional aerosol models
       3             : !===============================================================================
       4             : module sslt_sections
       5             :   use shr_kind_mod, only: r8 => shr_kind_r8
       6             : 
       7             :   implicit none
       8             : 
       9             :   private
      10             : 
      11             :   public :: sslt_sections_init
      12             :   public :: fluxes
      13             :   public :: nsections
      14             :   public :: Dg
      15             :   public :: rdry
      16             : 
      17             :   integer,parameter :: nsections = 31
      18             : 
      19             :   ! only use up to ~20um
      20             :   real(r8),parameter :: Dg(nsections) = (/  &
      21             :        0.0020e-5_r8, 0.0025e-5_r8, 0.0032e-5_r8,  &
      22             :        0.0040e-5_r8, 0.0051e-5_r8, 0.0065e-5_r8,  &
      23             :        0.0082e-5_r8, 0.0104e-5_r8, 0.0132e-5_r8,  &
      24             :        0.0167e-5_r8, 0.0211e-5_r8, 0.0267e-5_r8,  &
      25             :        0.0338e-5_r8, 0.0428e-5_r8, 0.0541e-5_r8,  &
      26             :        0.0685e-5_r8, 0.0867e-5_r8, 0.1098e-5_r8,  &
      27             :        0.1389e-5_r8, 0.1759e-5_r8, 0.2226e-5_r8,  &
      28             :        0.2818e-5_r8, 0.3571e-5_r8, 0.4526e-5_r8,  &
      29             :        0.5735e-5_r8, 0.7267e-5_r8, 0.9208e-5_r8,  &
      30             :        1.1668e-5_r8, 1.4786e-5_r8, 1.8736e-5_r8,  &
      31             :        2.3742e-5_r8 /)
      32             : 
      33             :   real(r8), dimension(nsections) :: bm, rdry, rm
      34             :   real(r8), dimension(4,nsections) :: consta, constb  !constants for calculating emission polynomial
      35             : 
      36             : contains
      37             :   
      38             :   !===========================================================================
      39             :   !===========================================================================
      40        1536 :   subroutine sslt_sections_init()
      41             : 
      42             :     integer :: m
      43             : 
      44             :     ! use Ekman's ss
      45        1536 :     rdry(:)=Dg(:)/2._r8   ! meter
      46             :     ! multiply rm with 1.814 because it should be RH=80% and not dry particles
      47             :     ! for the parameterization
      48       49152 :     rm(:)=1.814_r8*rdry(:)*1.e6_r8   ! um
      49       49152 :     bm(:)=(0.380_r8-log10(rm(:)))/0.65_r8  ! use in Manahan
      50             : 
      51             :     ! calculate constants form emission polynomials
      52       49152 :     do m=1,nsections
      53       49152 :        if ((m).le.9)then
      54       13824 :           consta(1,m) = (-2.576_r8)*10._r8**35*Dg(m)**4+5.932_r8*10._r8**28  &
      55             :                * Dg(m)**3+(-2.867_r8)*10._r8**21*Dg(m)**2+(-3.003_r8)  &
      56       13824 :                * 10._r8**13*Dg(m) + (-2.881_r8)*10._r8**6
      57             :           constb(1,m) = 7.188_r8*10._r8**37  &
      58             :                * Dg(m)**4+(-1.616_r8)*10._r8**31*Dg(m)**3+6.791_r8*10._r8**23  &
      59       13824 :                * Dg(m)**2+1.829_r8*10._r8**16*Dg(m)+7.609_r8*10._r8**8
      60       33792 :        elseif ((m).ge.10.and.(m).le.13)then
      61        6144 :           consta(2,m) = (-2.452_r8)*10._r8**33*Dg(m)**4+2.404_r8*10._r8**27  &
      62             :                * Dg(m)**3+(-8.148_r8)*10._r8**20*Dg(m)**2+(1.183_r8)*10._r8**14  &
      63        6144 :                * Dg(m)+(-6.743_r8)*10._r8**6
      64             :           constb(2,m) = 7.368_r8*10._r8**35  &
      65             :                * Dg(m)**4+(-7.310_r8)*10._r8**29*Dg(m)**3+ 2.528_r8*10._r8**23  &
      66        6144 :                * Dg(m)**2+(-3.787_r8)*10._r8**16*Dg(m)+ 2.279_r8*10._r8**9
      67       27648 :        elseif ((m).ge.14.and.(m).lt.22)then
      68       12288 :           consta(3,m) = (1.085_r8)*10._r8**29*Dg(m)**4+(-9.841_r8)*10._r8**23  &
      69             :                * Dg(m)**3+(3.132_r8)*10._r8**18*Dg(m)**2+(-4.165_r8)*10._r8**12  &
      70       12288 :                * Dg(m)+(2.181_r8)*10._r8**6
      71             :           constb(3,m) = (-2.859_r8)*10._r8**31  &
      72             :                * Dg(m)**4+(2.601_r8)*10._r8**26*Dg(m)**3+(-8.297_r8)*10._r8**20  &
      73       12288 :                * Dg(m)**2+(1.105_r8)*10._r8**15*Dg(m)+(-5.800_r8)*10._r8**8
      74       15360 :        elseif (m.ge.22.and.m.le.40)then
      75             :           ! use monahan
      76       15360 :           consta(4,m) = (1.373_r8*rm(m)**(-3)*(1+0.057_r8*rm(m)**1.05_r8)  &
      77             :                * 10**(1.19_r8*exp(-bm(m)**2)))  &
      78       30720 :                * (rm(m)-rm(m-1))
      79             :        endif
      80             :     enddo
      81        1536 :   end subroutine sslt_sections_init
      82             : 
      83             :   !===========================================================================
      84             :   !===========================================================================
      85     1489176 :   function fluxes ( sst, u10cubed, ncol ) result(fi)
      86             : 
      87             :     real (r8),intent(in) :: sst(:)
      88             :     real (r8),intent(in) :: u10cubed(:)
      89             :     integer  ,intent(in) :: ncol
      90             : 
      91             :     real (r8) :: fi(ncol,nsections)
      92             : 
      93             :     integer :: m
      94     2978352 :     real (r8) :: W(ncol)
      95             : 
      96             :     ! Calculations of source strength and size distribution
      97             :     ! NB the 0.1 is the dlogDp we have to multiplie with to get the flux, but the value dependence
      98             :     ! of course on what dlogDp you have. You will also have to change the sections of Dg if you use
      99             :     ! a different number of size bins with different intervals.
     100             : 
     101    26354952 :     W(:ncol)=3.84e-6_r8*u10cubed(:ncol)*0.1_r8 ! whitecap area
     102             : 
     103             :     ! calculate number flux fi (#/m2/s)
     104   772328232 :     fi(:,:)=0._r8
     105    47653632 :     do m=1,nsections
     106    47653632 :        if (m.le.9)then
     107   223791984 :           fi(:ncol,m)=W(:ncol)*((sst(:ncol))*consta(1,m)+constb(1,m))
     108    32761872 :        elseif (m.ge.10.and.m.le.13)then
     109    99463104 :           fi(:ncol,m)=W(:ncol)*((sst(:ncol))*consta(2,m)+constb(2,m))
     110    26805168 :        elseif (m.ge.14.and.m.lt.22)then
     111   198926208 :           fi(:ncol,m)=W(:ncol)*((sst(:ncol))*consta(3,m)+constb(3,m))
     112    14891760 :        elseif (m.ge.22.and.m.le.40)then
     113             :           ! use Monahan
     114   248657760 :           fi(:ncol,m)=consta(4,m)*u10cubed(:ncol)
     115             :        endif
     116             :     enddo
     117             : 
     118     1489176 :   end function fluxes
     119             : 
     120             : end module sslt_sections