
!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!


C RCS file, release, date & time of last delta, author, state, [and locker]
C $Header: /home/neumannd/data/M3HOME/CMAQ/std_5.0.1_seasalt/CMAQv5.0.1/models/CCTM/src/emis/emis/SSEMIS.F,v 1.1.1.1 2012/05/10 17:23:39 sjr Exp $

C what(1) key, module and SID; SCCS file; date and time of last delta:
C %W% %P% %G% %U%
 
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      MODULE SSEMIS
C-----------------------------------------------------------------------
C Description:

C Revision History:
C   04 Apr 11 S.Howard: generalization to accomodate different sea-salt speciation
C                       in aero5 and aero6
C   07 Apr 11 J.Young: additional mods
C   31 Oct 14 D. Neumann: External sea salt emissions can be read in. The env
C                         variables CTM_READSSEMIS [Y|N] and SSEMIS_1 
C                         [DIR/SSEMIS_INPUT_FILE] need to be specified for this.
C                         Emissions can be in I, J, and K mode. Num, srf, and
C                         specified mass emissions are read. Only mass emissions
C                         which are specified in AEROSPC are read from the file.
C                         When mass emis in one mode exist it is not checked
C                         whether num and srf emis exist as well. We trust the
C                         creator of the file :-) .
C                         
C-----------------------------------------------------------------------

      USE HGRD_DEFN
      USE AERO_DATA, ONLY: N_MODE

      IMPLICIT NONE

C Number of chemical species in fresh sea-salt aerosol
      INTEGER :: NSSSPC ! REMARK NEUMANND: ASO4J and ASO4K are ONE species

C Sea-Salt Emissions Rates
      REAL, ALLOCATABLE, SAVE :: SSOUTM( :,:,:,: )  ! mass emission rates [ug/m3/s]
      REAL, ALLOCATABLE, SAVE :: SSOUTN( :,:,: )    ! number emission rates [1/m3/s]
      REAL, ALLOCATABLE, SAVE :: SSOUTS( :,:,: )    ! surface-area emisrates [m2/m3/s]

      PUBLIC NSSSPC, SSOUTM, SSOUTN, SSOUTS, SSSPC_MAP, SSEMIS_INIT, GET_SSEMIS
      PRIVATE

C Position of H2O in sea-salt emission arrays
      INTEGER :: KH2O

C Variables for handling sea-salt emissions
      REAL,    ALLOCATABLE, SAVE :: OCEAN( :,: )   ! fractional sea water cover
      REAL,    ALLOCATABLE, SAVE :: SZONE( :,: )   ! fractional surf-zone cover
      INTEGER, ALLOCATABLE       :: SSSPC_MAP( : ) ! map to aerospc array

C Variables for the sea-salt diagnostic file
      INTEGER :: NSSDIAG                ! number of species in sea-salt
                                        ! diagnostic emission file
      REAL,    ALLOCATABLE       :: SSOUTD( : )   ! emission rates
      LOGICAL, SAVE              :: SSEMDIAG      ! flag for SSEMIS diagnostic file
      LOGICAL, SAVE              :: READSSEMIS    ! TEST WHETHER VARIABLES SHOULD BE READ IN FROM AN EXTERNAL FILE
      REAL,    ALLOCATABLE, SAVE :: SSBF( :,:,: ) ! seasalt emiss accumulator
      REAL,    ALLOCATABLE, SAVE :: WRSS( :,: )   ! seasalt emiss write buffer
      INTEGER, SAVE              :: SDATE, STIME  ! scenario start date & time

C Species names in the speciated sea-salt-emissions
      CHARACTER( 16 ), ALLOCATABLE :: WRSS_SPC( : )  ! species names

C Variables interpolated from the meteorological input files
      REAL,    ALLOCATABLE, SAVE :: PRES( :,: )    ! atmospheric pressure [Pa]
      REAL,    ALLOCATABLE, SAVE :: TA  ( :,: )    ! air temperature [K]
      REAL,    ALLOCATABLE, SAVE :: QV  ( :,: )    ! H2O mass mixing ratio [kg/kg]
      REAL,    ALLOCATABLE, SAVE :: WSPD10( :,: )  ! wind speed at 10m [m/s]

C Lognormal parameters fit to the open-ocean flux function
      REAL,    SAVE :: DGNJ( 136 )     ! geom mean diam of accum-mode [um]
      REAL,    SAVE :: DGNK( 136 )     ! geom mean diam of coarse-mode [um]
      REAL,    SAVE :: SIGJ( 136 )     ! geom std deviation of accum-mode flux
      REAL,    SAVE :: SIGK( 136 )     ! geom std deviation of coarse-mode flux
      REAL,    SAVE :: FNJ ( 136 )     ! magnitude of accum-mode flux [1/m2/s]
      REAL,    SAVE :: FNK ( 136 )     ! magnitude of coarse-mode flux [1/m2/s]

C Open-ocean flux numerically integrated over size range of interest
      REAL,    SAVE :: VFLUX( 136 )    ! total particle volume flux [m3/m2/s]

C Polycoeffs from Zhang et al. (2005)
      REAL( 8 ), SAVE :: C0_COEFF( 6 )   ! Eq 8
      REAL( 8 ), SAVE :: X_COEFF( 6 )    ! Eq 1

      INTEGER I

C NEW VARIABLES BY NEUMANND
      INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: SS_MAP_N, SS_MAP_S 
                        ! mapping from mode
                        ! id of particle number (*_N) and particle surface (*_S)
                        ! variables to variable id in the source file
      INTEGER, DIMENSION(:,:), ALLOCATABLE, SAVE :: SS_MAP_M, AERO_MAP_M
                        ! mapping from mode
                        ! id and aerospc/ssspc id to variable id in the source 
                        ! file
      
      

C-------------------------- Data Statements ----------------------------

C RH-dependent values calculated using MATLAB

C Geometric mean diameter of accumulation mode [um]
      DATA ( DGNJ( I ), I = 1, 136 ) /
     &      0.2651, 0.2677, 0.2696, 0.2713, 0.2736, 0.2752, 0.2768, 0.2784, 
     &      0.2806, 0.2823, 0.2840, 0.2858, 0.2883, 0.2903, 0.2924, 0.2951, 
     &      0.2974, 0.2998, 0.3028, 0.3054, 0.3080, 0.3113, 0.3140, 0.3174, 
     &      0.3203, 0.3232, 0.3267, 0.3297, 0.3332, 0.3363, 0.3394, 0.3430, 
     &      0.3462, 0.3499, 0.3533, 0.3573, 0.3611, 0.3652, 0.3700, 0.3752, 
     &      0.3806, 0.3871, 0.3945, 0.4030, 0.4129, 0.4246, 0.4259, 0.4271, 
     &      0.4284, 0.4298, 0.4311, 0.4328, 0.4342, 0.4356, 0.4371, 0.4385, 
     &      0.4401, 0.4416, 0.4434, 0.4451, 0.4467, 0.4484, 0.4501, 0.4518, 
     &      0.4538, 0.4556, 0.4575, 0.4594, 0.4613, 0.4633, 0.4655, 0.4676, 
     &      0.4697, 0.4718, 0.4740, 0.4764, 0.4787, 0.4811, 0.4835, 0.4859, 
     &      0.4886, 0.4912, 0.4938, 0.4965, 0.4994, 0.5022, 0.5051, 0.5081, 
     &      0.5112, 0.5143, 0.5175, 0.5208, 0.5242, 0.5276, 0.5312, 0.5348, 
     &      0.5386, 0.5424, 0.5463, 0.5503, 0.5544, 0.5587, 0.5631, 0.5675, 
     &      0.5722, 0.5769, 0.5818, 0.5867, 0.5919, 0.5972, 0.6026, 0.6083, 
     &      0.6140, 0.6201, 0.6263, 0.6325, 0.6392, 0.6458, 0.6529, 0.6600, 
     &      0.6676, 0.6751, 0.6833, 0.6914, 0.7001, 0.7088, 0.7177, 0.7276, 
     &      0.7372, 0.7478, 0.7582, 0.7690, 0.7810, 0.7928, 0.8050, 0.8187 / 

C Geometric mean diameter of coarse mode [um]
      DATA ( DGNK( I ), I = 1, 136 ) /
     &      1.7618, 1.8130, 1.8257, 1.8375, 1.8861, 1.8973, 1.9084, 1.9194, 
     &      1.9678, 1.9796, 1.9919, 2.0046, 2.0548, 2.0690, 2.0840, 2.1359, 
     &      2.1525, 2.1698, 2.2236, 2.2423, 2.2617, 2.3171, 2.3378, 2.3936, 
     &      2.4154, 2.4376, 2.4940, 2.5169, 2.5732, 2.5969, 2.6209, 2.6776, 
     &      2.7027, 2.7599, 2.7868, 2.8453, 2.8754, 2.9076, 2.9723, 3.0398, 
     &      3.0835, 3.1604, 3.2439, 3.3357, 3.4386, 3.5555, 3.5660, 3.5767, 
     &      3.5876, 3.5986, 3.6099, 3.6421, 3.6539, 3.6659, 3.6782, 3.6906, 
     &      3.7034, 3.7163, 3.7487, 3.7623, 3.7761, 3.7903, 3.8047, 3.8194, 
     &      3.8519, 3.8673, 3.8830, 3.8990, 3.9154, 3.9321, 3.9651, 3.9826, 
     &      4.0005, 4.0188, 4.0374, 4.0710, 4.0906, 4.1106, 4.1311, 4.1520, 
     &      4.1860, 4.2080, 4.2305, 4.2535, 4.2878, 4.3120, 4.3367, 4.3621, 
     &      4.3969, 4.4236, 4.4510, 4.4790, 4.5146, 4.5442, 4.5745, 4.6104, 
     &      4.6424, 4.6752, 4.7089, 4.7463, 4.7816, 4.8182, 4.8565, 4.8952, 
     &      4.9350, 4.9746, 5.0169, 5.0563, 5.1010, 5.1470, 5.1884, 5.2369, 
     &      5.2789, 5.3308, 5.3844, 5.4289, 5.4856, 5.5313, 5.5922, 5.6394, 
     &      5.7042, 5.7532, 5.8228, 5.8741, 5.9486, 6.0022, 6.0564, 6.1392, 
     &      6.1968, 6.2860, 6.3476, 6.4099, 6.5098, 6.5772, 6.6466, 6.7590 /

C Geometric standard deviation of accumulation mode
      DATA ( SIGJ( I ), I = 1, 136 ) /
     &      1.7572, 1.7612, 1.7612, 1.7612, 1.7651, 1.7651, 1.7651, 1.7651, 
     &      1.7689, 1.7689, 1.7689, 1.7689, 1.7727, 1.7727, 1.7727, 1.7764, 
     &      1.7764, 1.7764, 1.7801, 1.7801, 1.7801, 1.7837, 1.7837, 1.7873, 
     &      1.7873, 1.7873, 1.7907, 1.7907, 1.7941, 1.7941, 1.7941, 1.7974, 
     &      1.7974, 1.8006, 1.8006, 1.8037, 1.8037, 1.8037, 1.8066, 1.8095, 
     &      1.8095, 1.8122, 1.8147, 1.8172, 1.8195, 1.8217, 1.8217, 1.8217, 
     &      1.8217, 1.8217, 1.8217, 1.8237, 1.8237, 1.8237, 1.8237, 1.8237, 
     &      1.8237, 1.8237, 1.8255, 1.8255, 1.8255, 1.8255, 1.8255, 1.8255, 
     &      1.8272, 1.8272, 1.8272, 1.8272, 1.8272, 1.8272, 1.8287, 1.8287, 
     &      1.8287, 1.8287, 1.8287, 1.8300, 1.8300, 1.8300, 1.8300, 1.8300, 
     &      1.8312, 1.8312, 1.8312, 1.8312, 1.8322, 1.8322, 1.8322, 1.8322, 
     &      1.8330, 1.8330, 1.8330, 1.8330, 1.8336, 1.8336, 1.8336, 1.8340, 
     &      1.8340, 1.8340, 1.8340, 1.8342, 1.8342, 1.8342, 1.8343, 1.8343, 
     &      1.8343, 1.8342, 1.8342, 1.8339, 1.8339, 1.8339, 1.8334, 1.8334, 
     &      1.8327, 1.8327, 1.8327, 1.8318, 1.8318, 1.8307, 1.8307, 1.8295, 
     &      1.8295, 1.8281, 1.8281, 1.8266, 1.8266, 1.8248, 1.8229, 1.8230, 
     &      1.8209, 1.8209, 1.8188, 1.8164, 1.8164, 1.8140, 1.8115, 1.8115 /
              
C Geometric standard deviation of coarse mode
      DATA ( SIGK( I ), I = 1, 136 ) /
     &      2.1039, 2.0834, 2.0834, 2.0834, 2.0639, 2.0639, 2.0639, 2.0639, 
     &      2.0452, 2.0452, 2.0452, 2.0452, 2.0273, 2.0273, 2.0273, 2.0101, 
     &      2.0101, 2.0101, 1.9937, 1.9937, 1.9937, 1.9780, 1.9780, 1.9631, 
     &      1.9631, 1.9631, 1.9488, 1.9488, 1.9351, 1.9351, 1.9351, 1.9222, 
     &      1.9222, 1.9098, 1.9098, 1.8981, 1.8981, 1.8981, 1.8870, 1.8765, 
     &      1.8765, 1.8666, 1.8573, 1.8486, 1.8405, 1.8329, 1.8329, 1.8329, 
     &      1.8329, 1.8329, 1.8329, 1.8259, 1.8259, 1.8259, 1.8259, 1.8259, 
     &      1.8259, 1.8259, 1.8195, 1.8195, 1.8195, 1.8195, 1.8195, 1.8195, 
     &      1.8137, 1.8137, 1.8137, 1.8137, 1.8137, 1.8137, 1.8085, 1.8085, 
     &      1.8085, 1.8085, 1.8085, 1.8038, 1.8038, 1.8038, 1.8038, 1.8038, 
     &      1.7997, 1.7997, 1.7997, 1.7997, 1.7963, 1.7963, 1.7963, 1.7963, 
     &      1.7934, 1.7934, 1.7934, 1.7934, 1.7913, 1.7913, 1.7913, 1.7897, 
     &      1.7897, 1.7897, 1.7897, 1.7888, 1.7889, 1.7889, 1.7886, 1.7886, 
     &      1.7886, 1.7891, 1.7891, 1.7904, 1.7903, 1.7903, 1.7923, 1.7924, 
     &      1.7952, 1.7952, 1.7952, 1.7988, 1.7989, 1.8034, 1.8034, 1.8089, 
     &      1.8090, 1.8155, 1.8154, 1.8230, 1.8230, 1.8318, 1.8417, 1.8417, 
     &      1.8529, 1.8529, 1.8656, 1.8799, 1.8799, 1.8960, 1.9138, 1.9138 /
              
C Number flux in accumulation mode at 1 m/s [1/m2/s]
            DATA ( FNJ( I ), I = 1, 136 ) /
     &      89.326, 89.524, 89.524, 89.524, 89.709, 89.709, 89.709, 89.709, 
     &      89.881, 89.881, 89.881, 89.881, 90.040, 90.040, 90.040, 90.188, 
     &      90.188, 90.188, 90.325, 90.324, 90.324, 90.451, 90.451, 90.569, 
     &      90.569, 90.569, 90.678, 90.678, 90.779, 90.779, 90.779, 90.873, 
     &      90.873, 90.959, 90.959, 91.039, 91.039, 91.039, 91.113, 91.181, 
     &      91.181, 91.243, 91.301, 91.353, 91.401, 91.445, 91.445, 91.445, 
     &      91.445, 91.445, 91.445, 91.484, 91.484, 91.484, 91.484, 91.484, 
     &      91.484, 91.484, 91.520, 91.520, 91.520, 91.520, 91.520, 91.520, 
     &      91.551, 91.551, 91.551, 91.551, 91.551, 91.551, 91.579, 91.579, 
     &      91.579, 91.579, 91.579, 91.604, 91.604, 91.604, 91.604, 91.604, 
     &      91.625, 91.625, 91.625, 91.625, 91.642, 91.642, 91.642, 91.642, 
     &      91.656, 91.656, 91.656, 91.656, 91.667, 91.667, 91.667, 91.675, 
     &      91.675, 91.675, 91.675, 91.679, 91.679, 91.679, 91.680, 91.680, 
     &      91.680, 91.678, 91.678, 91.672, 91.672, 91.672, 91.663, 91.663, 
     &      91.650, 91.650, 91.650, 91.633, 91.633, 91.612, 91.612, 91.587, 
     &      91.587, 91.558, 91.558, 91.524, 91.524, 91.485, 91.441, 91.441, 
     &      91.392, 91.392, 91.337, 91.274, 91.274, 91.205, 91.129, 91.129 /
              
C Number flux in coarse mode at 1 m/s [1/m2/s]
      DATA ( FNK( I ), I = 1, 136 ) /
     &      11.315, 11.110, 11.110, 11.109, 10.919, 10.919, 10.919, 10.919, 
     &      10.740, 10.740, 10.740, 10.740, 10.573, 10.573, 10.573, 10.417, 
     &      10.417, 10.417, 10.270, 10.271, 10.271, 10.134, 10.134, 10.006, 
     &      10.006, 10.006, 9.886, 9.886, 9.774, 9.774, 9.774, 9.669, 
     &      9.669, 9.571, 9.571, 9.480, 9.480, 9.480, 9.395, 9.315, 
     &      9.315, 9.243, 9.174, 9.111, 9.053, 8.999, 8.999, 8.999, 
     &      9.000, 9.000, 9.000, 8.951, 8.951, 8.951, 8.951, 8.951, 
     &      8.951, 8.951, 8.907, 8.907, 8.907, 8.907, 8.907, 8.907, 
     &      8.867, 8.867, 8.867, 8.867, 8.867, 8.867, 8.832, 8.832, 
     &      8.832, 8.832, 8.832, 8.800, 8.800, 8.800, 8.800, 8.800, 
     &      8.773, 8.773, 8.773, 8.773, 8.751, 8.751, 8.751, 8.751, 
     &      8.732, 8.732, 8.732, 8.732, 8.718, 8.718, 8.718, 8.708, 
     &      8.708, 8.708, 8.708, 8.702, 8.703, 8.703, 8.701, 8.701, 
     &      8.701, 8.704, 8.704, 8.712, 8.712, 8.711, 8.724, 8.725, 
     &      8.742, 8.742, 8.742, 8.765, 8.765, 8.794, 8.794, 8.828, 
     &      8.828, 8.869, 8.869, 8.915, 8.915, 8.970, 9.031, 9.031, 
     &      9.101, 9.101, 9.179, 9.268, 9.268, 9.366, 9.476, 9.476 /
              
C Volume flux at 1 m/s numerically integrated over full size range [m3/m2/s]
      DATA ( VFLUX( I ), I = 1, 136 ) /
     &      3.3799E-16, 3.4012E-16, 3.4728E-16, 3.5407E-16, 3.5489E-16,  
     &      3.6124E-16, 3.6758E-16, 3.7402E-16, 3.7473E-16, 3.8152E-16,  
     &      3.8864E-16, 3.9617E-16, 3.9799E-16, 4.0632E-16, 4.1518E-16,  
     &      4.1825E-16, 4.2807E-16, 4.3846E-16, 4.4284E-16, 4.5420E-16,  
     &      4.6612E-16, 4.7169E-16, 4.8449E-16, 4.9073E-16, 5.0429E-16,  
     &      5.1829E-16, 5.2529E-16, 5.3992E-16, 5.4734E-16, 5.6260E-16,  
     &      5.7831E-16, 5.8648E-16, 6.0311E-16, 6.1215E-16, 6.3024E-16,  
     &      6.4087E-16, 6.6140E-16, 6.8386E-16, 6.9948E-16, 7.1775E-16,  
     &      7.4926E-16, 7.7574E-16, 8.0807E-16, 8.4803E-16, 8.9805E-16,  
     &      9.6150E-16, 9.7003E-16, 9.7878E-16, 9.8774E-16, 9.9692E-16,  
     &      1.0063E-15, 1.0026E-15, 1.0124E-15, 1.0224E-15, 1.0327E-15,  
     &      1.0432E-15, 1.0541E-15, 1.0652E-15, 1.0624E-15, 1.0739E-15,  
     &      1.0858E-15, 1.0980E-15, 1.1106E-15, 1.1235E-15, 1.1216E-15,  
     &      1.1351E-15, 1.1490E-15, 1.1633E-15, 1.1780E-15, 1.1931E-15,  
     &      1.1924E-15, 1.2083E-15, 1.2246E-15, 1.2415E-15, 1.2588E-15,  
     &      1.2593E-15, 1.2776E-15, 1.2964E-15, 1.3159E-15, 1.3360E-15,  
     &      1.3380E-15, 1.3592E-15, 1.3811E-15, 1.4038E-15, 1.4071E-15,  
     &      1.4311E-15, 1.4559E-15, 1.4816E-15, 1.4867E-15, 1.5139E-15,  
     &      1.5422E-15, 1.5715E-15, 1.5786E-15, 1.6098E-15, 1.6422E-15,  
     &      1.6509E-15, 1.6855E-15, 1.7215E-15, 1.7590E-15, 1.7706E-15,  
     &      1.8107E-15, 1.8526E-15, 1.8666E-15, 1.9115E-15, 1.9585E-15,  
     &      1.9754E-15, 2.0260E-15, 2.0447E-15, 2.0993E-15, 2.1565E-15,  
     &      2.1790E-15, 2.2411E-15, 2.2660E-15, 2.3335E-15, 2.4045E-15, 
     &      2.4346E-15, 2.5121E-15, 2.5456E-15, 2.6304E-15, 2.6678E-15, 
     &      2.7610E-15, 2.8028E-15, 2.9056E-15, 2.9525E-15, 3.0662E-15, 
     &      3.1191E-15, 3.1734E-15, 3.3053E-15, 3.3669E-15, 3.5144E-15, 
     &      3.5847E-15, 3.6572E-15, 3.8310E-15, 3.9146E-15, 4.0012E-15, 
     &      4.2082E-15 /                                                  

C Polynomial coefficients from Zhang et al. (2005)         
      DATA C0_COEFF / 28.376D0, -205.44D0, 653.37D0, -1031.7D0, 803.18D0, -247.08D0 /
      DATA X_COEFF  / 3.1657D0, -19.079D0,  55.72D0, -83.998D0, 63.436D0, -19.248D0 /

C Domain decomposition info from emission and meteorology files
      INTEGER, SAVE :: GXOFF, GYOFF             ! origin offset
      INTEGER, SAVE :: STRTCOL_O1, ENDCOL_O1, STRTROW_O1, ENDROW_O1
      INTEGER, SAVE :: STRTCOLMC3, ENDCOLMC3, STRTROWMC3, ENDROWMC3

      INTEGER, SAVE :: LOGDEV

      CONTAINS

C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
         FUNCTION SSEMIS_INIT( JDATE, JTIME, TSTEP ) RESULT ( SUCCESS )

C Revision History:
C   26 Feb 10 J.Young:  Eliminate BARRIER, etc. to prevent MPI race condition
C   16 Feb 11 S.Roselle: replaced I/O API include files with UTILIO_DEFN
C   04 Apr 11 S.Howard: dynamic allocation and generalization of indices for
C                       SSOUTD, WRSS_SPC, and SSSPC_MAP
C   11 May 11 D.Wong: incorporated twoway model implementation
C   31 Oct 14 D. Neumann: Look whether sea salt should be read in externally. 
C                           Sea salt emis input file is opened and mappings
C                           between netCDF variable IDs and SSSCP and AEROSPC
C                           are created.

            USE UTILIO_DEFN
            USE AERO_DATA

            IMPLICIT NONE

C Arguments:
            INTEGER, INTENT( IN ) :: JDATE, JTIME, TSTEP
            LOGICAL :: SUCCESS

C Includes:
            INCLUDE SUBST_FILES_ID  ! file name parameters

C External Functions:
            INTEGER, EXTERNAL :: SETUP_LOGDEV

C Local Variables:
            CHARACTER( 16 ), SAVE :: CTM_SSEMDIAG = 'CTM_SSEMDIAG' ! env var for
                                                     ! diagnostic file

            CHARACTER( 16 ), SAVE :: PNAME = 'SSEMIS_INIT'
            CHARACTER( 80 ) :: VARDESC
            CHARACTER( 96 ) :: XMSG = ' '
            INTEGER STATUS
            INTEGER S, N, L        

C NEW VARIABLES BY NEUMANND
            ! TEST WHETHER VARIABLES SHOULD BE READ IN FROM AN EXTERNAL FILE
            CHARACTER( 16 ), SAVE  :: CTM_READSSEMIS = 'CTM_READSSEMIS  '
            ! VARIABLE FOR FILENAME OF EXISTING SSEMIS FILE
            CHARACTER( 16 ), SAVE  :: SSEMIS_1 = 'SSEMIS_1'
            
            INTEGER, DIMENSION(:), ALLOCATABLE :: AEROSPC2SSSPC_TMP_MAP 
                        ! temporary mapping from aerosol species id (index in 
                        ! AEROSPC) to sea salt species id (index in SS_MAP_M 
                        ! bzw. AERO_MAP_M) (inversion of TMP_SSSPC_MAP and
                        ! SSSPC_MAP).
            LOGICAL, DIMENSION(:), ALLOCATABLE :: AEROSPC_PRESENT 
                        ! If AEROSPC_PRESENT at S is TRUE then the aerosol
                        ! species with id S is present in the netCDF input file
                        ! in the I, J, and/or K mode.
            INTEGER, DIMENSION(:), ALLOCATABLE :: TMP_SSSPC_MAP 
                        ! Temporary mapping from sea salt species id to aerosol
                        ! species id (inversion of AEROSPC2SSSPC_TMP_MAP). 
                        ! which SSSPC_MAP is created later. This array will
                        ! have the length N_AEROSPC while SSSPC_MAP will only
                        ! have the length NSSSPC. Therefore, we introduced this
                        ! seconda array
            INTEGER :: V
            INTEGER :: NMODES ! number of modes for which sea salt emissions
                              ! are read in or calculated:
                              !  == 2 => J and K modes
                              !  == 3 => I, J and K modes
            
C-----------------------------------------------------------------------

            LOGDEV = SETUP_LOGDEV()
            SUCCESS = .TRUE.

C Map aerospc data
            CALL MAP_AERO()

            
          ! START NEUMANND
          READSSEMIS = ENVYN( CTM_READSSEMIS,'Logical which indicates ', .FALSE., STATUS )
          ! Look whther environment variable CTM_READSSEMIS ist set. If it is 
          ! set and TRUE we try to read in sea salt emissins from an external
          ! file.
          
          IF ( READSSEMIS ) THEN
          ! OPEN FILE
            IF ( .NOT. OPEN3( SSEMIS_1, FSREAD3, PNAME ) ) THEN
               XMSG = 'Open failure for ' // SSEMIS_1
               CALL M3WARN( PNAME, JDATE, JTIME, XMSG )
               SUCCESS = .FALSE.; RETURN
            END IF
            
            IF ( .NOT. DESC3( SSEMIS_1 ) ) THEN
                XMSG = 'Description of file ' // trim(SSEMIS_1)
&                        // ' could not be retrieved.'
                CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, 1 )
            END IF
            
            ! ALLOCATE ARRAYS AND INITIALIZE VARIABLES
            ! these value is incremented below
            NSSSPC = 0 ! number of sea salt species
                       ! (e.g. ASO4J and ASO4K are ONE species)

            ALLOCATE(TMP_SSSPC_MAP(N_AEROSPC))
            TMP_SSSPC_MAP = -1
            ALLOCATE(AEROSPC2SSSPC_TMP_MAP(N_AEROSPC))
            AEROSPC2SSSPC_TMP_MAP = -1
            ALLOCATE(AEROSPC_PRESENT(N_AEROSPC))
            AEROSPC_PRESENT = .FALSE.
            ALLOCATE(AERO_MAP_M(N_AEROSPC, N_MODE))
            AERO_MAP_M = -1
            ALLOCATE(SS_MAP_N(N_MODE), SS_MAP_S(N_MODE))
            SS_MAP_N = -1
            SS_MAP_S = -1
            
            ! WE CREATE MAPPINGS FROM THE netCDF VARIABLES TO SEA SALT SPECIES
            ! AND AEROSOL SPECIES AND THE OTHER WAY ARROUND. THE MAPPINGS ARE 
            ! USER IN GET_SSEMIS LATER ON.
            ! iterate all input variables
            DO V = 1, NVARS3D
              WRITE(XMSG, '(A10,A16)') 'FIND VAR: ', VNAME3D( V )
              CALL M3MESG(XMSG)
              
              DO N = 1, N_MODE
                IF ( VNAME3D( V ) == AEROMODE( N )%NUM_NAME ) THEN
                  ! netCDF variable represents number emissions
                  CALL M3MESG('NUM')
                  SS_MAP_N( N ) = V
                ELSEIF ( VNAME3D( V ) == AEROMODE( N )%SRF_NAME ) THEN
                  ! netCDF variable represents surface emissions
                  CALL M3MESG('SRF')
                  SS_MAP_S( N ) = V
                ELSE
                  ! netCDF variable represents neither number nor surface
                  ! emissions.
                  CALL M3MESG('ELSE')
                  DO S = 1, N_AEROSPC
                    ! Iterate all aerosol species and have a look whether the
                    ! netCDF variable is an aerosol species (if not, we do 
                    ! nothing).
                    IF ( VNAME3D( V ) == AEROSPC( S )%NAME( N ) ) THEN
                      CALL M3MESG(AEROSPC( S )%NAME( N ))
                      ! INCREASE SPECIES COUNTER IF SPECIES NEW
                      IF ( .NOT. AEROSPC_PRESENT( S ) ) THEN
                        NSSSPC = NSSSPC + 1
                        AEROSPC_PRESENT( S ) = .TRUE.
                        AERO_MAP_M( NSSSPC, N ) = V
                        TMP_SSSPC_MAP( NSSSPC ) = S
                        AEROSPC2SSSPC_TMP_MAP( S ) = NSSSPC
                        IF (AEROSPC( S )%NAME(1) == 'AH2OI' ) KH2O = NSSSPC
                      ELSE
                        AERO_MAP_M( AEROSPC2SSSPC_TMP_MAP( S ), N ) = V
                        TMP_SSSPC_MAP( AEROSPC2SSSPC_TMP_MAP( S ) ) = S
                      END IF
                    END IF
                  END DO
                ENDIF
              END DO
            END DO
            
            ! Create and copy further mappings.
            ALLOCATE(SSSPC_MAP(NSSSPC))
            ALLOCATE(SS_MAP_M(NSSSPC, N_MODE)) ! LATER
            SS_MAP_M = -1
            
            DO S = 1, NSSSPC
              DO N = 1, N_MODE
                SS_MAP_M( S,N ) = AERO_MAP_M( S,N )
              END DO
              ! see variable declaration of TMP_SSSPC_MAP on why we do it this
              ! way and why we not directly write into SSSPC_MAP.
              SSSPC_MAP( S ) = TMP_SSSPC_MAP( S )
            END DO
            
            ! Allocate SS arrays using NSSSPC value
            ALLOCATE ( SSOUTM( NSSSPC,N_MODE,NCOLS,NROWS ),
     &                 SSOUTN( N_MODE,NCOLS,NROWS ),
     &                 SSOUTS( N_MODE,NCOLS,NROWS ), STAT = STATUS )
            IF ( STATUS .NE. 0 ) THEN
               XMSG = '*** SSOUTM, SSOUTN or SSOUTS memory allocation failed'
               CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
               SUCCESS = .FALSE.; RETURN
            END IF
            
            SSOUTM = 0.0
            SSOUTN = 0.0
            SSOUTS = 0.0

          ELSE
          ! END NEUMANND


C Determine number of sea-salt species
            NSSSPC = 0
            DO S = 1, N_AEROSPC
               IF ( AEROSPC_SSF( 2,S ) + AEROSPC_SSF( 3,S ) .GT. 0.0 ) NSSSPC = NSSSPC + 1
            END DO
C Include H2O
            NSSSPC = NSSSPC + 1

C Allocate mapping array
            ALLOCATE ( SSSPC_MAP( NSSSPC ), STAT = STATUS )
            IF ( STATUS .NE. 0 ) THEN
               XMSG = '*** SSSPC_MAP memory allocation failed'
               CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
               SUCCESS = .FALSE.; RETURN
            END IF

C build mapping array SSSPC_MAP
            NSSSPC = 0
            DO S = 1, N_AEROSPC
               IF ( AEROSPC_SSF( 2,S ) + AEROSPC_SSF( 3,S ) .GT. 0.0 ) THEN
                  NSSSPC = NSSSPC + 1
                  SSSPC_MAP( NSSSPC ) = S
               END IF
            END DO

C Include H2O to mapping
            NSSSPC = NSSSPC + 1      
            KH2O = NSSSPC
            SSSPC_MAP( NSSSPC ) = AH2O_IDX

C Allocate SS arrays using NSSSPC value
            ALLOCATE ( SSOUTM( NSSSPC,N_MODE,NCOLS,NROWS ),
     &                 SSOUTN( N_MODE,NCOLS,NROWS ),
     &                 SSOUTS( N_MODE,NCOLS,NROWS ), STAT = STATUS )
            IF ( STATUS .NE. 0 ) THEN
               XMSG = '*** SSOUTM, SSOUTN or SSOUTS memory allocation failed'
               CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
               SUCCESS = .FALSE.; RETURN
            END IF

            ALLOCATE ( OCEAN( NCOLS,NROWS ),
     &                 SZONE( NCOLS,NROWS ),
     &                 PRES ( NCOLS,NROWS ), 
     &                 TA   ( NCOLS,NROWS ), 
     &                 QV   ( NCOLS,NROWS ), 
     &                 WSPD10( NCOLS,NROWS ), STAT = STATUS )
            IF ( STATUS .NE. 0 ) THEN
               XMSG = '*** OCEAN, SZONE, PRES, TA, QV, or WSPD10 memory allocation failed'
               CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
               SUCCESS = .FALSE.; RETURN
            END IF

          ! START NEUMANND
          END IF ! IF READSSEMIS
          ! END NEUMANND
          
C Open the ocean file, which contains the ocean and surf-zone fractions

          IF ( .NOT. OPEN3( OCEAN_1, FSREAD3, PNAME ) ) THEN
             XMSG = 'Open failure for ' // OCEAN_1
             CALL M3WARN( PNAME, JDATE, JTIME, XMSG )
             SUCCESS = .FALSE.; RETURN
          END IF

          CALL SUBHFILE ( OCEAN_1, GXOFF, GYOFF,
     &                      STRTCOL_O1, ENDCOL_O1, STRTROW_O1, ENDROW_O1 )

C Open the 3D cross-point met file

          IF ( .NOT. OPEN3( MET_CRO_3D, FSREAD3, PNAME ) ) THEN
             XMSG = 'Open failure for ' // MET_CRO_3D
             CALL M3WARN( PNAME, JDATE, JTIME, XMSG )
             SUCCESS = .FALSE.; RETURN
          END IF

          CALL SUBHFILE ( MET_CRO_3D, GXOFF, GYOFF,
     &                      STRTCOLMC3, ENDCOLMC3, STRTROWMC3, ENDROWMC3 )


C Get sea-salt-emission diagnostic output file flag.

            SSEMDIAG = .FALSE.         ! default
            VARDESC = 'Flag for writing the sea-salt-emission diagnostic file'
            SSEMDIAG = ENVYN( CTM_SSEMDIAG, VARDESC, SSEMDIAG, STATUS )
            IF ( STATUS .NE. 0 ) WRITE( LOGDEV, '(5X, A)' ) VARDESC
            IF ( STATUS .EQ. 1 ) THEN
               XMSG = 'Environment variable improperly formatted'
               CALL M3WARN( PNAME, JDATE, JTIME, XMSG )
               SUCCESS = .FALSE.; RETURN
            ELSE IF ( STATUS .EQ. -1 ) THEN
               XMSG =
     &             'Environment variable set, but empty ... Using default:'
               WRITE( LOGDEV, '(5X, A, I9)' ) XMSG, JTIME
            ELSE IF ( STATUS .EQ. -2 ) THEN
               XMSG = 'Environment variable not set ... Using default:'
               WRITE( LOGDEV, '(5X, A, I9)' ) XMSG, JTIME
            END IF

            IF ( SSEMDIAG ) THEN    ! Open the sea-salt emission diagnostic file
C Determine number of sea-salt species For diagnostic file
            NSSDIAG = 0
               DO S = 1, N_AEROSPC
! NEW NEUMANND: I mode
                  IF ( AEROSPC_SSF( 1,S ) .GT. 0.0 ) NSSDIAG = NSSDIAG + 1 ! NEW
                  IF ( AEROSPC_SSF( 2,S ) .GT. 0.0 ) NSSDIAG = NSSDIAG + 1
                  IF ( AEROSPC_SSF( 3,S ) .GT. 0.0 ) NSSDIAG = NSSDIAG + 1
               END DO
               
               IF ( READSSEMIS ) THEN
                  NMODES = 3  ! Assumption: When ssemis are read in from a file
                              !  then they might be defined for all three modes.
                              !  If they are not, it is not further problematic
                              !  because non-existing netCDF variables will be
                              !  set to zero in the code.
               ELSE
                  NMODES = 2
               ENDIF
               
C Include H2O
! MODIFIED NEUMANND: generalized for variable number of modes
!              NSSDIAG = NSSDIAG + 2   !! count both J and K modes
!              NSSDIAG = NSSDIAG + 3   !! count both J and K modes
               NSSDIAG = NSSDIAG + NMODES
               
C Include NUMBER and SURFACE_AREA to diagnostic (J & K)
! MODIFIED NEUMANND: generalized for variable number of modes
!              NSSDIAG = NSSDIAG + 4
!              NSSDIAG = NSSDIAG + 6
               NSSDIAG = NSSDIAG + 2 * NMODES

C Allocate diagnostic arrays
               ALLOCATE ( SSOUTD( NSSDIAG ),
     &                    WRSS_SPC( NSSDIAG ), STAT = STATUS )
               IF ( STATUS .NE. 0 ) THEN
                  XMSG = '*** SSOUTD or WRSS_SPC memory allocation failed'
                  CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
                  SUCCESS = .FALSE.; RETURN
               END IF

C Build diagnostic file write buffer (WRSS_SPC) array
               NSSDIAG = 0
! MODIFIED NEUMANND: generalized for variable number of modes
!              DO N = 2, N_MODE
!              DO N = 1, N_MODE
               DO N = (3 + 1 - NMODES), N_MODE
                  DO S = 1, N_AEROSPC
                     IF ( AEROSPC_SSF( N,S ) .GT. 0.0 ) THEN
                        NSSDIAG = NSSDIAG + 1
                        WRSS_SPC( NSSDIAG ) = AEROSPC( S )%NAME( N )
                     END IF
                  END DO
C Include H2O
                  NSSDIAG = NSSDIAG + 1
                  WRSS_SPC( NSSDIAG ) = AEROSPC( AH2O_IDX )%NAME( N )
               END DO

C Include Mode NUMBER to diagnostic
! MODIFIED NEUMANND: generalized for variable number of modes
!              DO N = 2, N_MODE
!              DO N = 1, N_MODE
               DO N = (3 + 1 - NMODES), N_MODE
                  NSSDIAG = NSSDIAG + 1
                  WRSS_SPC( NSSDIAG ) = AEROMODE( N )%NUM_NAME
               END DO
 
C Include Mode SURFACE_AREA to diagnostic
! MODIFIED NEUMANND: generalized for variable number of modes
!              DO N = 2, N_MODE
!              DO N = 1, N_MODE
               DO N = (3 + 1 - NMODES), N_MODE
                  NSSDIAG = NSSDIAG + 1
                  WRSS_SPC( NSSDIAG ) = AEROMODE( N )%SRF_NAME
               END DO

C Open the sea-salt emission dignostic file
               SDATE = ENVINT( 'CTM_STDATE', 'Scenario Start (YYYYJJJ)', 0, STATUS )
               STIME = ENVINT( 'CTM_STTIME', 'Scenario Start (HHMMSS)', 0, STATUS )

               IF ( MYPE .EQ. 0 ) CALL OPSSEMIS ( SDATE, STIME, TSTEP, NSSDIAG, WRSS_SPC )

               ALLOCATE ( SSBF( NSSDIAG,MY_NCOLS,MY_NROWS ),
     &                    WRSS( NCOLS,NROWS ), STAT = STATUS )
               IF ( STATUS .NE. 0 ) THEN
                  XMSG = '*** SSBF or WRSS memory allocation failed'
                  CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
                  SUCCESS = .FALSE.; RETURN
               END IF

            END IF

         END FUNCTION SSEMIS_INIT

C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
         SUBROUTINE OPSSEMIS ( JDATE, JTIME, TSTEP, NSSDIAG, WRSS_SPC )

C   16 May 05 P.Bhave: original version, using OPDDEP as a template
C    4 Mar 10 J.Young: accomodate Steve Howard's reengineering
C   16 Feb 11 S.Roselle: replaced I/O API include files with UTILIO_DEFN

         USE GRID_CONF           ! horizontal & vertical domain specifications
         USE UTILIO_DEFN

         IMPLICIT NONE

         INCLUDE SUBST_FILES_ID  ! file name parameters

C Arguments:
         INTEGER      JDATE      ! current model date, coded YYYYDDD
         INTEGER      JTIME      ! current model time, coded HHMMSS
         INTEGER      TSTEP      ! output time step
         INTEGER      NSSDIAG
         CHARACTER( 16 ) :: WRSS_SPC( NSSDIAG )

C Local variables:
         CHARACTER( 16 ) :: PNAME = 'OPSSEMIS'
         CHARACTER( 96 ) :: XMSG = ' '

         INTEGER      V, L       ! loop induction variables

C-----------------------------------------------------------------------

C Try to open existing file for update
         IF ( .NOT. OPEN3( CTM_SSEMIS_1, FSRDWR3, PNAME ) ) THEN
            XMSG = 'Could not open CTM_SSEMIS_1 for update - '
     &           // 'try to open new'
            CALL M3MESG( XMSG )

C Set output file characteristics based on COORD.EXT and open diagnostic file
            FTYPE3D = GRDDED3
            SDATE3D = JDATE
            STIME3D = JTIME
            TSTEP3D = TSTEP
            CALL NEXTIME( SDATE3D, STIME3D, TSTEP3D ) !  start the next hour

            NVARS3D = NSSDIAG
            NCOLS3D = GL_NCOLS
            NROWS3D = GL_NROWS
            NLAYS3D = 1
            NTHIK3D = 1
            GDTYP3D = GDTYP_GD
            P_ALP3D = P_ALP_GD
            P_BET3D = P_BET_GD
            P_GAM3D = P_GAM_GD
            XORIG3D = XORIG_GD
            YORIG3D = YORIG_GD
            XCENT3D = XCENT_GD
            YCENT3D = YCENT_GD
            XCELL3D = XCELL_GD
            YCELL3D = YCELL_GD
            VGTYP3D = VGTYP_GD
            VGTOP3D = VGTOP_GD
!           VGTPUN3D = VGTPUN_GD ! currently, not defined
            DO L = 1, NLAYS3D + 1
               VGLVS3D( L ) = VGLVS_GD( L )
            END DO
!           GDNAM3D = GDNAME_GD
            GDNAM3D = GRID_NAME  ! from HGRD_DEFN

            DO V = 1, NSSDIAG
               VTYPE3D( V ) = M3REAL
               VNAME3D( V ) = WRSS_SPC( V )
               IF ( VNAME3D( V )(1:3) .EQ. 'NUM' ) THEN
                  UNITS3D( V ) = 'number/s'
               ELSE IF ( VNAME3D( V )(1:3) .EQ. 'SRF' ) THEN
                  UNITS3D( V ) = 'm2/s'
               ELSE
                  UNITS3D( V ) = 'g/s'
               END IF
               VDESC3D( V ) = 'hourly ' // TRIM( VNAME3D( V ) )
     &                     // ' sea-salt emission rate'
            END DO

            FDESC3D( 1 ) = 'hourly layer-1 sea-salt emission rates'
            DO L = 2, MXDESC3
               FDESC3D( L ) = ' '
            END DO

C Open sea-salt emissions diagnostic file
            IF ( .NOT. OPEN3( CTM_SSEMIS_1, FSNEW3, PNAME ) ) THEN
               XMSG = 'Could not create the CTM_SSEMIS_1 file'
               CALL M3EXIT( PNAME, SDATE3D, STIME3D, XMSG, XSTAT1 )
            END IF

         END IF

         RETURN

         END SUBROUTINE OPSSEMIS

C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
         SUBROUTINE GET_SSEMIS ( JDATE, JTIME, TSTEP, CELLVOL, RJACM, CELLHGT )

C GET_SSEMIS calculates the sea-salt emission rates in a grid cell
C given the fractional grid-cell area covered by open ocean and surf zone
 
C Key Subroutines/Functions Called:  NONE
 
C Revision History:
 
C   13 Jun 05 P.Bhave:  first version for public release
C   11 Apr 08 J.Kelly:  added code for (1) emission of coarse surface area,
C                       (2) emission of coarse water, (3) enhanced surf zone
C                       emissions, and (4) variable coarse std. deviation
C   20 Feb 10 J.Young:  move out of AERO_EMIS proper into this module
C   16 Feb 11 S.Roselle: replaced I/O API include files with UTILIO_DEFN
C   04 Apr 11 S.Howard: moved speciation to AERO_DATA and generalized species
C                       and modal references and density calculation
C   11 May 11 D.Wong: incorporated twoway model implementation
C   31 Oct 14 D. Neumann: 
 
C References:
C    Alduchov and Eskridge, "Improved Magnus Form Approximations of
C          Saturation Vapor Pressure," Jour. of Applied Meteorology, vol. 35,
C          pp 601-609, April, 1996.
C    Binkowski F.S., S.J. Roselle. Models-3 Community Multiscale Air Quality 
C          (CMAQ) model aerosol component 1: Model Description.  J. Geophys. 
C          Res., Vol 108, No D6, 4183, doi:10.1029/2001JD001409 (2003).
C    Seinfeld, J.H., S.N. Pandis.  "Atmospheric Chemistry and Physics - from
C          air pollution to climate change"  John Wiley & Sons, Inc. (1998).
C    Zhang, K.M., E.M. Knipping, A.S. Wexler, P.V. Bhave, G.S. Tonnesen
C          "Size distribution of sea-salt emissions as a function of relative
C          humidity"  Atmospheric Environment, 39, 3373-3379 (2005).
C    Clarke, A.D., S.R. Owens, and J. Zhou "An ultrafine sea-salt flux from
C          breaking waves: Implications for cloud condensation nuclei in the
C          remote marine atmosphere" J. Geophys. Res. (2006)
C-----------------------------------------------------------------------

            USE UTILIO_DEFN
            USE AERO_DATA
            USE AEROMET_DATA   ! Includes CONST.EXT

            IMPLICIT NONE

C Arguments
            INTEGER, INTENT( IN ) :: JDATE, JTIME, TSTEP( 3 )
            REAL,    INTENT( IN ) :: CELLVOL               ! grid-cell volume [m2*sigma]
            REAL,    INTENT( IN ) :: RJACM( NCOLS,NROWS )  ! reciprocal Jacobian [1/m]
            REAL,    INTENT( IN ) :: CELLHGT               ! grid-cell height [sigma]

C Includes:
            INCLUDE SUBST_FILES_ID  ! file name parameters

C Local Variables:
            REAL    :: OFRAC              ! fractional seawater cover
            REAL    :: SFRAC              ! fractional surf-zone cover
            REAL    :: BLKPRS             ! atmospheric pressure [Pa]
            REAL    :: BLKTA              ! air temperature [K]
            REAL    :: BLKQV              ! H2O mass mixing ratio [ kg/kg ]
            REAL    :: U10                ! wind speed at 10m [m/s]
            REAL    :: RLAY1HT            ! reciprocal of layer-1 height [1/m]
            REAL    :: AIRVOL             ! grid-cell volume [m3]

            CHARACTER( 16 ), SAVE :: PNAME = 'GET_SSEMIS'
            CHARACTER( 96 ) :: XMSG = ' '
            INTEGER C, R, L, N, V, S         ! loop indices

C Variables for calculating ambient RH as done in the aerosol module
            REAL, PARAMETER :: AA = 610.94      ! coefficients from Alduchov
            REAL, PARAMETER :: BB = 17.625      !  and Eskridge (1996)
            REAL, PARAMETER :: CC = 243.04
            REAL, PARAMETER :: EPSWATER = MWWAT / MWAIR
            REAL      :: ESAT             ! saturation vapor pressure
            REAL      :: BLKRH            ! fractional relative humidity

C Variables for calculating solute mass fraction & size-correction factors
            REAL      :: RH               ! BLKRH constrained to 45-99% range
            REAL( 8 ) :: DRH              ! double precision RH
            REAL( 8 ) :: C0               ! (Diam @ formation) / (Diam @ ambient RH)
            REAL( 8 ) :: C80              ! (Diam @ 80% RH) / (Diam @ ambient RH)    
            REAL( 8 ) :: XSW              ! fractional solute weight fraction

C Density of dry sea salt [kg/m3] using volume additivity
            REAL( 8 ) :: RHOSSDRY 
            REAL( 8 ) :: RHOSW            ! sea-salt particle density at ambient RH [g/m3]
            REAL( 8 ) :: RHSUM       

C Maximum and minimum diameters for emission-flux integrations [m]
            REAL, PARAMETER :: DPMAX = 20.0E-6 ! upper limit used for the
                                               ! numerical integrations
            REAL      :: DPMINO           ! minimum diameter for open-ocean flux
            REAL      :: DPMAXO           ! maximum diameter for open-ocean flux
   
C Wind-speed-dependent scaling factors for emissions flux functions
            REAL      :: WDSCALO          ! Gong-Monahan open-ocean flux

C Fraction of whitecap coverage
            REAL      :: WCAP

C Lognormal parameters and numerically-integrated fluxes at ambient RH
            INTEGER   :: IPOS             ! position of ambient RH in data arrays
            REAL      :: DGNJRH           ! geometric mean of accum-mode flux [m]
            REAL      :: DGNKRH           ! geometric mean of coarse-mode flux [m]
            REAL      :: SIGJRH           ! geom std deviation of accum-mode flux
            REAL      :: SIGKRH           ! geom std deviation of coarse-mode flux
            REAL      :: FNJRH            ! magnitude of accum-mode flux [1/m2/s]
            REAL      :: FNKRH            ! magnitude of coarse-mode flux [1/m2/s]
            REAL      :: VFLXRH           ! total particle volume [m3/m2/s]

C Alternate forms of geometric std devs stored for efficiency
            REAL      :: LNSGJ, LNSGK     ! natural log of std dev
            REAL      :: RDIVJ, RDIVK     ! reciprocal of sqrt(2)*lnsg

C Geometric mean diameters by volume/mass [m]
            REAL      :: DGVJRH           ! accumulation mode
            REAL      :: DGVKRH           ! coarse mode

C Variables for converting 3rd moment emission rates to M0 and M2 rates
            REAL      :: FACTNUMJ         ! conversion for accumulation mode M0
            REAL      :: FACTNUMK         ! conversion for coarse mode M0
            REAL      :: FACTM2J          ! conversion for accumulation mode M2
            REAL      :: FACTM2K          ! conversion for coarse mode M2 

C Variables for calculating fraction of mass emissions assigned to each mode
            REAL      :: FFINJ            ! finite integral for accumulation mode
            REAL      :: FFINK            ! finite integral for coarse mode
            REAL      :: FRACMJ           ! mass fraction assigned to accumulation mode
            REAL      :: FRACMK           ! mass fraction assigned to coarse mode

C Mode-specific sea-salt mass emission rates [g/m3/s]
            REAL( 8 ) :: ESEASJ           ! accumulation mode
            REAL( 8 ) :: ESEASK           ! coarse mode

C 3rd moment emission rates [m3/m3/s]
            REAL      :: M3J              ! accumulation mode, 3rd moment
            REAL      :: M3K              ! coarse mode, 3rd moment
         
C Variables for calculating dry surface-area emission rate for accum mode
            REAL      :: WET_M3J, DRY_M3J ! wet & dry 3rd moment emis rates [m3/m3/s] 
            REAL      :: WET_M2J, DRY_M2J ! wet & dry 2nd moment emis rates [m3/m3/s] 
            REAL      :: WET_M3K, DRY_M3K ! wet & dry 3rd moment emis rates [m3/m3/s] 
            REAL      :: WET_M2K, DRY_M2K ! wet & dry 2nd moment emis rates [m3/m3/s] 

            INTEGER, SAVE     :: WSTEP  = 0    ! local write counter
            INTEGER           :: MDATE, MTIME  ! SSEMIS write date&time

C Mathematical constants; statement function for ERF approximation
            REAL, PARAMETER :: ONE3 = 1.0 / 3.0
            REAL, PARAMETER :: TWO3 = 2.0 / 3.0
            REAL            :: ERF        ! ERF approx. statement function
            REAL            :: X          ! dummy argument for ERF
            ERF( X ) = SIGN( 1.0, X ) * SQRT( 1.0 - EXP( -4.0 * X * X / PI ) )
            
C NEW VARIABLES BY NEUMANND
            CHARACTER( 16 ), SAVE  :: SSEMIS_1 = 'SSEMIS_1'
            INTEGER                :: STATUS
            
            REAL, DIMENSION(:,:), ALLOCATABLE :: SS_TMP ! content of individual
                                         ! variables from the input netCDF file.
            
                    
C----------------------------- Begin calc ------------------------------

          ! START NEUMANND
          IF ( READSSEMIS ) THEN
            ! If sea salt emissions should be read in (READSSEMIS == TRUE) then
            ! we read them in from a netCDF file till "ELSE ! IF READSSEMIS"
            ! and skip the calculation part till "END IF ! IF READSSEMIS". The
            ! code after "ELSE ! IF READSSEMIS" is unchanged and normally
            ! processed.
            
            ! Allocate temporary array for reading in netCDF variables.
            ALLOCATE(SS_TMP( NCOLS,NROWS ))
            SS_TMP = 0.0
            
            ! GO THROUGH EMITTED MASSES BY MODE AND SPECIES
            DO S = 1, NSSSPC
              DO N = 1, N_MODE
                IF ( SS_MAP_M( S, N ) == -1 ) THEN
                  ! species does not exist in this mode
                  DO R = 1, MY_NROWS
                    DO C = 1, MY_NCOLS
                      SSOUTM( S,N,C,R ) = 0.0
                    END DO
                  END DO
                ELSE
                  ! species exists in this mode
                  ! Name of the current substance: AEROSPC(SPC_MODE_POS(V))%NAME(N)
                  V = SS_MAP_M( S, N )
                  IF ( .NOT. INTERPX( SSEMIS_1, AEROSPC(SSSPC_MAP(S))%NAME(N), PNAME,
&                                     STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3,
&                                     1,1, JDATE, JTIME, SS_TMP ) ) THEN
                     XMSG = 'Could not interpolate  from ' // SSEMIS_1
                     CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
                  END IF
                  
                  ! write input data into SSOUTM
                  DO R = 1, MY_NROWS
                    DO C = 1, MY_NCOLS
                      AIRVOL = CELLVOL / RJACM( C,R )
                      SSOUTM( S,N,C,R ) = SS_TMP( C,R ) / AIRVOL
                    END DO
                  END DO
                END IF
              END DO
            END DO
            
            ! GO THROUGH SURFACE AND NUMBER BY MODE
            DO N = 1, N_MODE
              ! Start with number emis
              IF ( SS_MAP_N( N ) == -1 ) THEN
                ! Number emis of mode N do not exist.
                DO R = 1, MY_NROWS
                  DO C = 1, MY_NCOLS
                    SSOUTN( N,C,R ) = 0.0
                  END DO
                END DO
              ELSE
                ! Number emis of mode N exist.
                V = SS_MAP_N( N )
                IF ( .NOT. INTERPX( SSEMIS_1, AEROMODE( N )%NUM_NAME, PNAME,
&                            STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3,
&                            1,1, JDATE, JTIME, SS_TMP ) ) THEN
                     XMSG = 'Could not interpolate  from ' // MET_CRO_3D
                     CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
                END IF
                
                ! write input data into SSOUTN
                DO R = 1, MY_NROWS
                  DO C = 1, MY_NCOLS
                    AIRVOL = CELLVOL / RJACM( C,R )
                    SSOUTN( N,C,R ) = SS_TMP( C,R ) / AIRVOL 
                  END DO
                END DO
              END IF
              
              ! Surface emis
              IF ( SS_MAP_S( N ) == -1 ) THEN
                ! Surface emis of mode N do not exist.
                DO R = 1, MY_NROWS
                  DO C = 1, MY_NCOLS
                    SSOUTS( N,C,R ) = 0.0
                  END DO
                END DO
              ELSE
                ! Surface emis of mode N exist.
                V = SS_MAP_S( N )
                IF ( .NOT. INTERPX( SSEMIS_1, AEROMODE( N )%SRF_NAME, PNAME,
&                            STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3,
&                            1,1, JDATE, JTIME, SS_TMP ) ) THEN
                     XMSG = 'Could not interpolate  from ' // MET_CRO_3D
                     CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
                END IF
                
                ! write input data into SSOUTS
                DO R = 1, MY_NROWS
                  DO C = 1, MY_NCOLS
                    AIRVOL = CELLVOL / RJACM( C,R )
                    SSOUTS( N,C,R ) = SS_TMP( C,R ) / AIRVOL 
                  END DO
                END DO
              END IF
            END DO
            
            
            ! Initialize sea-salt output buffer
            IF ( SSEMDIAG .AND. WSTEP .EQ. 0 ) SSBF = 0.0
            
            ! Update the SSBF array, for writing the diagnostic sea-salt-emission file.
            ! This code is just copied from the Readssemis-Else case. It needed
            ! to be copied because in the else case it is integrated in two
            ! large ROW and COL loops. In order to not duplicate these lines,
            ! I would have been forced to modify the original code parts which
            ! I did not want to do.
            DO R = 1, MY_NROWS
              DO C = 1, MY_NCOLS
                IF ( SSEMDIAG ) THEN
                   V = 0
                   DO N = 1, 3
                      DO S = 1, NSSSPC
                         IF ( S .EQ. KH2O .OR. AEROSPC_SSF( N,SSSPC_MAP( S ) ) .GT. 0 ) THEN
                            V = V + 1
                            SSOUTD( V ) = SSOUTM( S,N,C,R )
                         END IF
                      END DO
                   END DO

                   DO N = 1, 3
                      V = V + 1
                      SSOUTD( V ) = SSOUTN( N,C,R )
                   END DO

                   DO N = 1, 3
                      V = V + 1
                      SSOUTD( V ) = SSOUTS( N,C,R )
                   END DO

                   AIRVOL = CELLVOL / RJACM( C,R )

                   DO S = 1, NSSDIAG
                      SSBF( S,C,R ) = SSBF( S,C,R ) + SSOUTD( S ) * AIRVOL
&                                   * FLOAT( TIME2SEC ( TSTEP( 2 ) ) )
                   END DO
                END IF  ! SSEMDIAG
              END DO
            END DO
          ELSE ! IF READSSEMIS
          ! END NEUMANND
            WRITE(XMSG, '(A6,I4,A2,I4,A2,I4,A2,I4)') 'Test: ', STRTCOLMC3, ", ", ENDCOLMC3, ", ", STRTROWMC3, ", ", ENDROWMC3
            CALL M3MESG(XMSG)
          
            IF ( .NOT. INTERPX( MET_CRO_3D, 'PRES', PNAME,
&                               STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3,
&                               1,1, JDATE, JTIME, PRES ) ) THEN
               XMSG = 'Could not interpolate PRES from ' // MET_CRO_3D
               CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
            END IF

            IF ( .NOT. INTERPX( MET_CRO_3D, 'TA', PNAME,
     &                          STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3,
     &                          1,1, JDATE, JTIME, TA ) ) THEN
               XMSG = 'Could not interpolate TA from ' // MET_CRO_3D
               CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
            END IF

            IF ( .NOT. INTERPX( MET_CRO_3D, 'QV', PNAME,
     &                          STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3,
     &                          1,1, JDATE, JTIME, QV ) ) THEN
               XMSG = 'Could not interpolate QV from ' // MET_CRO_3D
               CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
            END IF

            IF ( .NOT. INTERPX( MET_CRO_2D, 'WSPD10', PNAME,
     &                          STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3,
     &                          1,1, JDATE, JTIME, WSPD10 ) ) THEN
               XMSG = 'Could not find WSPD10 in ' // MET_CRO_2D //
     &                ' search for WIND10.'
               CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
               IF ( .NOT. INTERPX( MET_CRO_2D, 'WIND10', PNAME,
     &                             STRTCOLMC3,ENDCOLMC3, STRTROWMC3,ENDROWMC3,
     &                             1,1, JDATE, JTIME, WSPD10 ) ) THEN
                  XMSG = 'Could not find WIND10 or WSPD10 in ' // MET_CRO_2D
                  CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
               END IF
            END IF
            
C Read fractional seawater and surf-zone coverage from the OCEAN file.
C Store results in the OCEAN and SZONE arrays.
            IF ( .NOT. INTERPX( OCEAN_1, 'OPEN', PNAME,
     &                          STRTCOL_O1,ENDCOL_O1, STRTROW_O1,ENDROW_O1,
     &                          1,1, JDATE, JTIME, OCEAN ) ) THEN
               XMSG = 'Could not interpolate OPEN from ' // OCEAN_1
               CALL M3EXIT ( PNAME, JDATE, JTIME, XMSG, XSTAT1 )
            END IF

            IF ( .NOT. INTERPX( OCEAN_1, 'SURF', PNAME,
     &                          STRTCOL_O1,ENDCOL_O1, STRTROW_O1,ENDROW_O1,
     &                          1,1, JDATE, JTIME, SZONE ) ) THEN
               XMSG = 'Could not interpolate SURF from ' // OCEAN_1 //
     &                '.  Surf-zone emissions set to zero.'
               CALL M3WARN ( PNAME, JDATE, JTIME, XMSG )
            END IF
          
C Initialize sea-salt output buffer
            IF ( SSEMDIAG .AND. WSTEP .EQ. 0 ) SSBF = 0.0

            DO R = 1, MY_NROWS
               DO C = 1, MY_NCOLS
                  IF ( OCEAN( C,R ) + SZONE( C,R ) .GT. 0.0 ) THEN
                     OFRAC   = OCEAN ( C,R )
                     SFRAC   = SZONE ( C,R )
                     BLKPRS  = PRES  ( C,R )
                     BLKTA   = TA    ( C,R )
                     BLKQV   = QV    ( C,R )
                     U10     = WSPD10( C,R )
                     RLAY1HT = RJACM ( C,R ) / CELLHGT
                  ELSE
                     DO N = 1, N_MODE
                        DO S = 1, NSSSPC
                           SSOUTM( S,N,C,R ) = 0.0
                        END DO
                        SSOUTN( N,C,R ) = 0.0
                        SSOUTS( N,C,R ) = 0.0
                     END DO
                     CYCLE
                  END IF

C Calculate fractional relative humidity in the grid cell
                  ESAT  = AA * EXP( BB * ( BLKTA - 273.15 ) / ( BLKTA - 273.15 + CC ) )
                  BLKRH = BLKPRS * BLKQV / ( ( EPSWATER + BLKQV ) *  ESAT )       
                  BLKRH = MIN( 0.99, MAX( 0.005, BLKRH ) )

C Calculate RH-dependent size-correction factors and solute mass fraction
C using equations from Zhang et al. (2005), which are valid over the
C 45-99% RH range.  Store powers of RH for polynomial calculations.
                  RH = BLKRH
                  RH = MAX( 0.45, MIN( 0.99, RH ) )
                  DRH = REAL( RH, 8 )

C Calculate size-correction factors from Equations 8 & 9 of Zhang et al.
                  C0 =         C0_COEFF( 1 )
     &               + DRH * ( C0_COEFF( 2 )
     &               + DRH * ( C0_COEFF( 3 )
     &               + DRH * ( C0_COEFF( 4 )
     &               + DRH * ( C0_COEFF( 5 )
     &               + DRH * ( C0_COEFF( 6 ) ) ) ) ) )
                  C80 = C0 / 1.97D0

C Calculate solute mass fraction using Equation 1 of Zhang et al.
                  XSW =         X_COEFF( 1 )
     &                + DRH * ( X_COEFF( 2 )
     &                + DRH * ( X_COEFF( 3 )
     &                + DRH * ( X_COEFF( 4 )
     &                + DRH * ( X_COEFF( 5 )
     &                + DRH * ( X_COEFF( 6 ) ) ) ) ) )

C Calculate dry sea-salt [kg/m3]
                  RHSUM = 0.0D0
                  DO S = 1, NSSSPC
                     IF ( S .EQ. KH2O ) CYCLE     !! do not include Water
                     RHSUM = RHSUM
     &                     + REAL( ( AEROSPC_SSF( 2,SSSPC_MAP( S ) ) 
     &                             / AEROSPC( SSSPC_MAP( S ) )%DENSITY ), 8 )
                  END DO
                  RHOSSDRY = 1.0D0 / RHSUM    

C Calculate sea-salt-particle density [g/m3] at ambient RH, assuming
C volume additivity of dry salt plus water
                  RHOSW = 1.0D+3 / ( XSW / RHOSSDRY
     &                  + ( 1.0D0 - XSW )
     &                  / REAL( AEROSPC( AH2O_IDX )%DENSITY, 8 ) )

C Set minimum and maximum diameters for integration using the size-
C correction factors of Zhang et al. (2005)

C Gong-Monahan flux function is valid from 0.03-5.0 um dry radius;
C multiply by 4 to get wet diameter at 80% RH
                  DPMINO = 1.2E-7 / REAL( C80, 4 )
                  DPMAXO = MIN( DPMAX, DPMAX / REAL( C80, 4 ) ) 

C deLeeuw flux function is valid from 1.6-20 um diameter at formation
!                 DPMINS = 1.6E-6 / C0 
!                 DPMAXS = MIN( DPMAX, DPMAX / C0 )
            
C Calculate scaling factors to correct the tabulated fluxes for the 10m
C wind speed in this cell.  Note: tabulated values are based on 1 m/s.
                  WDSCALO = MIN( U10, 20.0 ) ** 3.41        ! Gong-Monahan flux function
!                 WDSCALS = EXP( 0.23 * MIN( U10, 9.0 ) )   ! deLeeuw flux function

C Calculate whitecap coverage fraction 
                  WCAP = 3.84E-6 * WDSCALO ! Eq. 5, Clarke et al. (2006) JGR

C Find position in the data arrays that corresponds to ambient RH in this
C grid cell and set the appropriate data values
                  IF ( BLKRH .LE. 0.45 ) THEN
                     IPOS = 1
                  ELSE IF ( BLKRH .LE. 0.90 ) THEN
                     IPOS = NINT( 100.0 * BLKRH - 44.0 )
                  ELSE IF ( BLKRH .LE. 0.99 ) THEN
                     IPOS = NINT( 1000.0 * BLKRH - 854.0 )
                  ELSE
                     IPOS = 136
                  END IF

                  DGNJRH = DGNJ ( IPOS ) * 1.0E-6   ! convert to [m]
                  DGNKRH = DGNK ( IPOS ) * 1.0E-6   ! convert to [m]
                  SIGJRH = SIGJ ( IPOS )
                  SIGKRH = SIGK ( IPOS ) 
                  FNJRH  = FNJ  ( IPOS )
                  FNKRH  = FNK  ( IPOS )
                  VFLXRH = VFLUX( IPOS ) * ( OFRAC + SFRAC / WCAP )

C Save certain functions of the geometric standard deviations for later use
                  LNSGJ  = LOG( SIGJRH )
                  LNSGK  = LOG( SIGKRH )
                  RDIVJ  = 1.0 / ( SQRT( 2.0 ) * LNSGJ )
                  RDIVK  = 1.0 / ( SQRT( 2.0 ) * LNSGK )

C Calculate geometric-mean diameters by volume using Equation 7.52 of
C Seinfeld & Pandis (1998).
                  DGVJRH = DGNJRH * EXP( 3.0 * LNSGJ * LNSGJ )
                  DGVKRH = DGNKRH * EXP( 3.0 * LNSGK * LNSGK )

C Calculate modal volume fluxes [m3/m2/s] by evaluating finite integrals 
C from DPMIN to DPMAX over each lognormal distribution.  Use resulting
C values to calculate the fraction of the total number emissions to 
C assign to each mode.  See Equations 19 and 20 of Uma Shankar`s
C "Integration of Sea-Salt Fluxes" memo.
                  FFINJ  = 0.5 * FNJRH * DGVJRH ** 3 * EXP( -4.5 * LNSGJ * LNSGJ )
     &                         * ( ERF( LOG( DPMAXO / DGVJRH ) * RDIVJ )  
     &                           - ERF( LOG( DPMINO / DGVJRH ) * RDIVJ ) )
                  FFINK  = 0.5 * FNKRH * DGVKRH ** 3 * EXP( -4.5 * LNSGK * LNSGK )
     &                         * ( ERF( LOG( DPMAXO / DGVKRH ) * RDIVK )  
     &                           - ERF( LOG( DPMINO / DGVKRH ) * RDIVK ) )
                  FRACMJ = FFINJ / ( FFINJ + FFINK )
                  FRACMK = 1.0 - FRACMJ

C Calculate mode-specific mass emission rates [g/m3/s], by multiplying
C numerically-integrated volume fluxes by the modal volume fractions,
C scaling for wind speed, dividing by grid-cell height, and multiplying
C by particle density.  Multiply by chemical speciation factors to
C obtain speciated mass emissions.

                  ESEASJ = VFLXRH * FRACMJ * WDSCALO * RLAY1HT * REAL( RHOSW, 4 )
                  ESEASK = VFLXRH * FRACMK * WDSCALO * RLAY1HT * REAL( RHOSW, 4 )

                  DO S = 1, NSSSPC
                     SSOUTM( S,1,C,R ) = 0.0
                     IF ( S .EQ. KH2O ) THEN    ! Water
                        SSOUTM( S,2,C,R ) = ESEASJ * REAL( ( 1.0D0 - XSW ), 4 )
                        SSOUTM( S,3,C,R ) = ESEASK * REAL( ( 1.0D0 - XSW ), 4 )
                     ELSE
                        SSOUTM( S,2,C,R ) = ESEASJ * REAL( XSW, 4 )
     &                                    * AEROSPC_SSF( 2,SSSPC_MAP( S ) )
                        SSOUTM( S,3,C,R ) = ESEASK * REAL( XSW, 4 )
     &                                    * AEROSPC_SSF( 3,SSSPC_MAP( S ) )
                     END IF
                  END DO

C Calculate mode-specific 3rd moment emission rates [m3/m3/s]
                  M3J = ESEASJ * F6PI / REAL( RHOSW, 4 )
                  M3K = ESEASK * F6PI / REAL( RHOSW, 4 )

C Calculate factors for converting 3rd moment emission rates into number
C and 2nd moment emission rates.  See Equations 7b and 7c of Binkowski &
C Roselle (2003)
                  FACTNUMJ = EXP( 4.5 * LNSGJ * LNSGJ ) / DGVJRH ** 3
                  FACTNUMK = EXP( 4.5 * LNSGK * LNSGK ) / DGVKRH ** 3
                  FACTM2J  = EXP( 0.5 * LNSGJ * LNSGJ ) / DGVJRH
                  FACTM2K  = EXP( 0.5 * LNSGK * LNSGK ) / DGVKRH 

C Calculate mode-specific emission rates of particle number [1/s]
                  SSOUTN( 1,C,R ) = 0.0
                  SSOUTN( 2,C,R ) = M3J * FACTNUMJ
                  SSOUTN( 3,C,R ) = M3K * FACTNUMK

C Calculate mode-specific dry surface area emission rates [m2/m3/s].  
C Subtract water from 3rd moment to obtain dry 3rd moment emission rate.  
C Calculate dry 2nd moment while holding the standard deviation constant.
C Multiply dry 2nd moment by PI to obtain dry surface area emission rate.
                  WET_M3J = M3J 
                  WET_M2J = M3J * FACTM2J
                  DRY_M3J = WET_M3J - F6PI * SSOUTM( KH2O,2,C,R ) * 1.0E-06
                  DRY_M2J = WET_M2J * ( DRY_M3J / WET_M3J ) ** TWO3

                  WET_M3K = M3K
                  WET_M2K = M3K * FACTM2K
                  DRY_M3K = WET_M3K - F6PI * SSOUTM( KH2O,3,C,R ) * 1.0E-06
                  DRY_M2K = WET_M2K * ( DRY_M3K / WET_M3K ) ** TWO3

                  SSOUTS( 1,C,R ) = 0.0
                  SSOUTS( 2,C,R ) = PI * DRY_M2J
                  SSOUTS( 3,C,R ) = PI * DRY_M2K

C Update the SSBF array, for writing the diagnostic sea-salt-emission file.
                  IF ( SSEMDIAG ) THEN
                     V = 0
                     DO N = 2, 3
                        DO S = 1, NSSSPC
                           IF ( S .EQ. KH2O .OR. AEROSPC_SSF( N,SSSPC_MAP( S ) ) .GT. 0 ) THEN
                              V = V + 1
                              SSOUTD( V ) = SSOUTM( S,N,C,R )
                           END IF
                        END DO
                     END DO

                     DO N = 2, 3
                        V = V + 1
                        SSOUTD( V ) = SSOUTN( N,C,R )
                     END DO

                     DO N = 2, 3
                        V = V + 1
                        SSOUTD( V ) = SSOUTS( N,C,R )
                     END DO

                     AIRVOL = CELLVOL / RJACM( C,R )

                     DO S = 1, NSSDIAG
                        SSBF( S,C,R ) = SSBF( S,C,R ) + SSOUTD( S ) * AIRVOL
     &                                * FLOAT( TIME2SEC ( TSTEP( 2 ) ) )
                     END DO
                  END IF  ! SSEMDIAG

               END DO   ! C
            END DO   ! R

!           write( logdev,* ) 'sum ssoutm, ssoutn, ssouts: ',
!    &                         sum( ssoutm ), sum( ssoutn ), sum( ssouts )

            
          ! START NEUMANND
          END IF ! IF READSSEMIS
          ! END NEUMANND
          
C If last call this hour, write out the total sea-salt emissions [g/s].
C Then reset the sea-salt emissions array and local write counter.
            IF ( SSEMDIAG ) THEN

               WSTEP = WSTEP + TIME2SEC( TSTEP( 2 ) )

               IF ( WSTEP .GE. TIME2SEC( TSTEP( 1 ) ) ) THEN
                  IF ( .NOT. CURRSTEP( JDATE, JTIME, SDATE, STIME, TSTEP( 1 ),
     &                                 MDATE, MTIME ) ) THEN
                     XMSG = 'Cannot get step date and time'
                     CALL M3EXIT( PNAME, JDATE, JTIME, XMSG, XSTAT3 )
                  END IF
                  CALL NEXTIME( MDATE, MTIME, TSTEP( 1 ) )
                  DO V = 1, NSSDIAG
                     DO R = 1, MY_NROWS
                        DO C = 1, MY_NCOLS
                           WRSS( C,R ) = SSBF( V,C,R ) / FLOAT( WSTEP )
                        END DO
                     END DO
                     IF ( .NOT. WRITE3( CTM_SSEMIS_1, WRSS_SPC( V ),
     &                          MDATE, MTIME, WRSS ) ) THEN
                        XMSG = 'Could not write ' // CTM_SSEMIS_1 // ' file'
                        CALL M3EXIT( PNAME, MDATE, MTIME, XMSG, XSTAT1 )
                     END IF
                  END DO
                  WRITE( LOGDEV, '( /5X, 3( A, :, 1X ), I8, ":", I6.6 )' )
     &                  'Timestep written to', CTM_SSEMIS_1,
     &                  'for date and time', MDATE, MTIME
                  WSTEP = 0
                  SSBF = 0.0
               END IF

            END IF  ! SSEMDIAG

            RETURN
            
         END SUBROUTINE GET_SSEMIS

      END MODULE SSEMIS
