      program oceanic

!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!c Copyright (C) 2012-2018  Ramboll
!c
!c This program is free software; you can redistribute it and/or
!c modify it under the terms of the GNU General Public License
!c as published by the Free Software Foundation; either version 2
!c of the License, or (at your option) any later version.
!c
!c This program is distributed in the hope that it will be useful,
!c but WITHOUT ANY WARRANTY; without even the implied warranty of
!c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!c GNU General Public License for more details.
!c
!c To obtain a copy of the GNU General Public License
!c go to the Free Software Foundation at http://www.fsf.org.
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

!  This program generates gridded oceanic emissions for CAMx:
!      Sea salt PM component species ("NA", "PCL", "PSO4" g/hr)
!      Sea salt gas component species ("SSCL", "SSBR", mol/hr)
!      Halo-methane compounds (for halogen gas-phase mechanism, mol/hr)
!      Dimethylsulfide (DMS, mol/hr)
!
!  NOTES: THIS PROGRAM READS CAMx V6 MET FILE FORMATS
!         USE METCONVERT FIRST IF YOU HAVE EARLIER MET FORMATS
!
!         SEA ICE IN THE CAMx 2-D SURFACE FILE IS USED TO LIMIT OR ZERO
!         EMISSIONS
!
!         CHLOROPHYLL-A FROM SATELLITE DATA IS USED FOR HALO-METHANE EMISSIONS
!
!         MONTHLY CLIMATOLOGICAL DMS OCEAN CONCENTRATIONS ARE USED FOR DMS
!         EMISSIONS
!
!         LAND/WATER (OCEAN) MASK CAN BE TEXT (O3MAP FORMAT) OR BINARY
!         (CAMx) FORMAT

      implicit none

!-----File units

      integer, parameter :: iout   = 10, &
                            icst   = 11, &
                            i3dmet = 12, &
                            i2dmet = 13, &
                            ilumet = 15, &
                            ifswm  = 17, &
                            ifchl  = 18, &
                            ifdms  = 19

!-----Constants

      real, parameter :: PI_6 = 0.52360  ! Pi / 6
      real, parameter :: LN10 = 2.30259  ! ln(10)
      real, parameter :: conv = 57.29578 ! degree/radian

!-----Size distribution parameters

      real, parameter :: D_MIN = 0.0     ! min integration diam [um], current RH
      real, parameter :: D_MAX = 2.5     ! max integration diam [um], current RH
      integer, parameter :: NBINS  = 10  ! number of integration sections

!-----Oceanic emission species composition

      integer, parameter :: MXSPEC = 19  ! number of supported output species

!-----Sea-salt speciation factors based on dry mass derived from
!     composition of seawater (Lewis and Schwartz, 2004, p.48):
!               mass ratio  MW
!       Cl(-)   19.35       35.45
!       SO4(2-)  2.71       96.06
!       Na(+)   10.78       22.99
!       Mg(2+)   1.28       24.31
!       Ca(2+)   0.41       40.08
!       K(+)     0.40       39.10
!     Non-sodium sea salt cations are represented by equivalent
!     concentrations of sodium to achieve electroneutrality with
!     the chloride and sulfate anions.

      real, parameter ::  SPCFAC_NA   = 0.387
      real, parameter ::  SPCFAC_CL   = 0.538
      real, parameter ::  SPCFAC_SO4  = 0.075
      real, parameter ::  SPCFAC_SSCL = 0.0282    ! g Cl/hr to mol Cl/hr
      real, parameter ::  SPCFAC_SSBR = 0.00154   ! mol Cl/hr to mol Br/hr

!-----Density of dry NaCl particle [kg/m3] (Grini et al., 2002)

      real, parameter :: dens_NaCl = 2200.

!-----Grid cell max concentration of chlorophyll-A

      real, parameter :: cap_chla = 1. ! [ug/m3]

!-----Halo-methane species index

      integer, parameter :: iMB3  = 1, &
                            iMB2  = 2, &
                            iMBC  = 3, &
                            iMB2C = 4, &
                            iMBC2 = 5, &
                            iCH3I = 6, &
                            iMIC  = 7, &
                            iMIB  = 8, &
                            iMI2  = 9

!-----Halo-methane speciation factors

      real, parameter, dimension(9) :: fhalo = (/ 2.00, & ! MB3
                                                  0.50, & ! MB2
                                                  0.08, & ! MBC
                                                  0.12, & ! MB2C
                                                  0.10, & ! MBC2
                                                  1.60, & ! CH3I
                                                  1.41, & ! MIC
                                                  0.42, & ! MIB
                                                  0.46 /) ! MI2

!-----Halo-methane diurnal profile

      real, parameter, dimension(24) :: fhrly = (/ &
         0.032, 0.032, 0.032, 0.033, 0.034, 0.036, & !  1 -  6
         0.039, 0.044, 0.051, 0.057, 0.062, 0.064, & !  7 - 12
         0.062, 0.057, 0.051, 0.044, 0.039, 0.036, & ! 13 - 18
         0.034, 0.033, 0.032, 0.032, 0.032, 0.032 /) ! 19 - 24

!-----DMS emissions related constants

      real, parameter :: mw_DMS = 62.13        ! DMS molecular weight [g mol^-1]
      real, parameter :: sqDMS = 0.53847613014 ! sqrt(mw_H2O/mw_DMS) 
      real, parameter :: R_gas = 0.082057338   ! R gas constant [L atm K^-1 mol^-1]

!-----Parameters for wind speed at 10 m

      real, parameter :: minspd = 2.0 ! min wind speed for flux estimates [m/s]

!-----Minimum salt water when lbin_lwmask = T

      real, parameter :: minsw = 0.0001

!-----Output file type identifier

      character(4), parameter, dimension(10) :: ifile = &
                          (/'E','M','I','S','S','I','O','N','S',' '/)

!-----Open ocean sea salt source function

      character(7), parameter :: OO_SSSF = &
                                'OO_OV14' ! Ovadnevaite et al., 2014
                              ! 'OO_GO03' ! Gong, 2003
                              ! 'OO_SO11' ! Sofiev et al., 2011

!-----Surf zone sea salt source function

      character(7), parameter :: SZ_SSSF = &
                                'SZ_GO03' ! Gong, 2003 with 100% whitecap
                                          !  (Kelly et al., 2010)
                              ! 'SZ_DL00' ! de Leeuw et al., 2000

!-----DMS water-side velocity scheme

      character(4), parameter :: DMS_KW = &
                                'LM86'    ! Liss & Merlivat (1986)
                              ! 'NI00'    ! Nightingale et al. (2000)

      real, parameter :: theta_GO03 = 30. ! theta for GO03 parameterization
      real, parameter :: sclfac_Hs = 0.75 ! factor to scale down fully developed
                                          ! wind wave Hs
      real, parameter, dimension(5) :: Dg_OV14 = &
                  (/ 0.018, 0.041, 0.09, 0.23, 0.83 /)
      real, parameter, dimension(5) :: sig_OV14 = &
                  (/ 1.37, 1.5, 1.42, 1.53, 1.85 /)
      real, parameter, dimension(5) :: A_OV14 = &
                  (/ 104.5, 0.0442, 149.6, 2.96, 0.51 /)
      real, parameter, dimension(5) :: B_OV14 = &
                  (/ 1.e5, 1.e5, 1.e5, 1.e5, 2.e5 /)
      real, parameter, dimension(5) :: C_OV14 = &
                  (/ 0.556, 1.08, 0.545, 0.79, 0.87 /)

      integer :: i, j, k, l, m, n, ii, jj, iii, jjj
      integer :: ibin, imode, idt
      integer :: itzon, jdate, ied, idum
      integer :: idt1,idt2
      integer :: nspecs, n3dmet, n2dmet, nlumet, nlwmet
      integer :: iutm, nx, ny, nz, iproj, istag
      integer :: inx, iny
      integer :: ipp, itt, iwv, iuu, ivv, iww, isi, its, iow, isf
      integer :: nxchl, nychl
      integer :: nxdms, nydms, idms, nsum
      integer :: mm1, mm2, dd1, dd2, year1, year2
      integer :: yy1, yy2, idate1, idate2
      integer :: ios
      integer :: ilochr
      integer :: addday, subday

      real :: factor
      real :: efactor
      real :: sal

      real :: plon, plat, xorg, yorg, dx, dy, tlat1, tlat2
      real :: bt, et, hr1, hr2
      real :: xloc, yloc
      real :: area, areas, jlat, clen
      real :: centlon
      real :: dxchl, dychl

      real :: qwatr, ev, es

      real :: r80(0:NBINS), wt(0:NBINS)
      real :: logwidth
      real :: r80_min, r80_min_oo, r80_min_sz
      real :: r80_max, r80_max_oo, r80_max_sz
      real :: r80_lower_limit_oo, r80_lower_limit_sz
      real :: r80_upper_limit_oo, r80_upper_limit_sz
      real :: sum_vol, sum_vol_sz, dF_dDp, dF_dlogDp, dV_dlogDp
      real :: term1, term2, term3
      real :: A_GO03, B_GO03
      real :: whtcap, aa, bb
      real :: ftw, ftw298, ftw288, ftw278, ftw271
      real :: fsw, fsw33, fsw09
      real :: Cd, Hs, Vis, ReHw, xerf, xx1, xx2
      real :: spd10
      real :: alpha,ka,kw,k600,Tc,SC,gamma_a,sumdms

      integer, allocatable, dimension(:,:) :: salt_water
      integer, allocatable, dimension(:,:) :: iloc, jloc

      real, allocatable, dimension(:,:) :: rlat, rlon
      real, allocatable, dimension(:,:) :: rmet2d, water, sea_ice, surf
      real, allocatable, dimension(:,:) :: ts, spd10m, relhum
      real, allocatable, dimension(:,:) :: u10, v10
      real, allocatable, dimension(:,:) :: coastline_len, srfzone_width
      real, allocatable, dimension(:,:) :: area2d, areasf, rchla_conc
      real, allocatable, dimension(:,:) :: mass,massdms
      real, allocatable, dimension(:,:) :: chla_conc
      real, allocatable, dimension(:,:) :: rdms_conc,dms_conc

      real, allocatable, dimension(:,:,:) :: rmet3d, pp, tt, wv
      real, allocatable, dimension(:,:,:) :: flux, mass_crs_fn
     
      logical :: lss, lhalo, lbin_lwmask, ldms, ldmssulf

      character(4) :: fname(10), mspec(10,MXSPEC), namvar(10)
      character(4) :: fileiden(60), metnote(60)
      character(10) :: c10, metspc
      character(60) :: c60
      character(80) :: inline
      character(256) :: ipath, icstpath, ichla, idmspath

!===============================================================================
!-----Begin main

      lbin_lwmask = .false.

!-----Get environment variable toggle for binary land/water mask

      call getenv('BIN_LWMASK', inline)
      if (inline.eq.'') then
        write(*,*) 'Environment variable BIN_LWMASK not defined:'
        write(*,*) '  Set to default - ', lbin_lwmask
      else
        inline = ADJUSTL(inline)
        if (inline(1:1).eq.'Y' .or. inline(1:1).eq.'y' .or. &
            inline(1:1).eq.'T' .or. inline(1:1).eq.'t') then
          lbin_lwmask = .true.
        else if (inline(1:1).eq.'N' .or. inline(1:1).eq.'n' .or. &
                 inline(1:1).eq.'F' .or. inline(1:1).eq.'f') then
          lbin_lwmask = .false.
        else
          stop'ERROR: Invalid value for BIN_LWMASK'
        endif
        write(*,*) 'Environment variable BIN_LWMASK set to', lbin_lwmask
      endif

!-----Read input data from the script

      read(*,'(20x,f10.0)') efactor
      write(*,*) 'Emis scaling factor:',efactor
      read(*,'(20x,f10.0)') sal
      write(*,*) 'Ocean salinity [g/kg]:',sal
      read(*,'(20x,a)') icstpath     
      write(*,*) 'Coastline length filename:',trim(icstpath)

!-----Information for the output emissions file

      read(*,'(20x,a)') ipath
      open(iout,file=ipath,form='unformatted',status='unknown')
      write(*,*) 'Output emissions filename:',trim(ipath)
      read(*,'(20x,a)') c60               
      write(*,*) 'Output file note         :',trim(c60)
      do j = 1,60
        fileiden(j)(1:1) = c60(j:j)
      enddo
      read(*,'(20x,i5)') nspecs
      write(*,*) 'Ouput species  : ',nspecs
      if (nspecs .gt. MXSPEC) then
        write(*,'(2a)') 'You are asking for more species than the ', &
                        'program can generate'
        write(*,*) 'Max is ', MXSPEC
        stop
      endif

      lss = .false.
      lhalo = .false.
      ldms = .false.
      do i = 1,nspecs
        read(*,'(20x,a10)') c10
        write(*,*) c10
        if (trim(c10).ne.'NA'   .and. trim(c10).ne.'PCL'  .and. &
            trim(c10).ne.'PSO4' .and. trim(c10).ne.'SSCL' .and. &
            trim(c10).ne.'SSBR' .and. trim(c10).ne.'MB3'  .and. &
            trim(c10).ne.'MB2'  .and. trim(c10).ne.'MBC'  .and. &
            trim(c10).ne.'MB2C' .and. trim(c10).ne.'MBC2' .and. &
            trim(c10).ne.'CH3I' .and. trim(c10).ne.'MIC'  .and. &
            trim(c10).ne.'MIB'  .and. trim(c10).ne.'MI2'  .and. &
            trim(c10).ne.'I2'   .and. trim(c10).ne.'CHLA' .and. &
            trim(c10).ne.'DMS'  .and. trim(c10).ne.'DMSWC'.and. &
            trim(c10).ne.'WIND' .and. trim(c10).ne.'RELH' .and. &
            trim(c10).ne.'SRFT') then
          write(*,'(2a)')'Unrecognized output species name: ',c10
          write(*,*)' Allowed names are:'
          write(*,*)'NA'
          write(*,*)'PCL'
          write(*,*)'PSO4'
          write(*,*)'SSCL'
          write(*,*)'SSBR'
          write(*,*)'MB3'
          write(*,*)'MB2'
          write(*,*)'MBC'
          write(*,*)'MB2C'
          write(*,*)'MBC2'
          write(*,*)'CH3I'
          write(*,*)'MIC'
          write(*,*)'MIB'
          write(*,*)'MI2'
          write(*,*)'I2'
          write(*,*)'DMS'
          write(*,*)'DMSWC (diagnostic only)'
          write(*,*)'CHLA (diagnostic only)'
          write(*,*)'WIND (diagnostic only)'
          write(*,*)'RELH (diagnostic only)'
          write(*,*)'SRFT (diagnostic only)'
          stop
        endif
        if (trim(c10).eq.'NA'   .or. trim(c10).eq.'PCL'  .or. &
            trim(c10).eq.'PSO4' .or. trim(c10).eq.'SSCL' .or. &
            trim(c10).eq.'SSBR') then
          lss = .true.
        endif

        if (trim(c10).eq.'MB3'  .or. trim(c10).eq.'MB2'  .or. &
            trim(c10).eq.'MBC'  .or. trim(c10).eq.'MB2C' .or. &
            trim(c10).eq.'MBC2' .or. trim(c10).eq.'CH3I' .or. &
            trim(c10).eq.'MIC'  .or. trim(c10).eq.'MIB'  .or. &
            trim(c10).eq.'MI2'  .or. trim(c10).eq.'I2') then
          lhalo = .true.
        endif

        if (trim(c10).eq.'DMS') then
          ldms = .true.
          idms = i
        endif

        read(c10,'(10a1)') (mspec(j,i),j=1,10)
      enddo

!-----Decide if DMS or SO2 is written out

      read(*,'(20x,l10)') ldmssulf          
      if (ldms) then
        write(*,*) 'DMS output as SO2:',ldmssulf
        if (ldmssulf) then
          c10 = 'SO2       '
          read(c10,'(10a1)') (mspec(j,idms),j=1,10)
        endif
      endif

!-----CAMx meteorological and land/water mask data

      read(*,'(20x,A)') ipath
      open(i3dmet,file=ipath,form='unformatted',status='old')
      write(*,*) 'Input CAMx 3D file         :',trim(ipath)
      read(*,'(20x,A)') ipath
      open(i2dmet,file=ipath,form='unformatted',status='old')
      write(*,*) 'Input CAMx 2D file       :',trim(ipath)
      read(*,'(20x,A)') ipath
      open(ilumet,file=ipath,form='unformatted',status='old')
      write(*,*) 'Input CAMx landuse file    :',trim(ipath)
      read(*,'(20x,A)') ipath
      if (lbin_lwmask) then
        open(ifswm,file=ipath,form='unformatted',status='old')
      else
        open(ifswm,file=ipath,form='formatted',status='old')
      endif
      write(*,*) 'Salt water mask file       :',trim(ipath)

!-----Halo-methane algorithm

      if (lhalo) then
        read(*,'(20x,A)') ichla
        open(ifchl,file=ichla,form='formatted',status='old')
        write(*,*) 'Chlorophyll-A file       :',trim(ichla)
      endif

!-----DMS algorithm

      if (ldms) then
        read(*,'(20x,A)') idmspath
        open(ifdms,file=idmspath,form='formatted',status='old')
        write(*,*) 'DMS Concentration file   :',trim(idmspath)
      endif

!-----Read input met file headers (and binary lwmask file if used)

      read(i3dmet) fname,metnote,itzon,n3dmet,jdate,bt,ied,et
      read(i2dmet) fname,metnote,idum, n2dmet
      read(ilumet) fname,metnote,idum, nlumet
      if (lbin_lwmask) &
        read(ifswm) fname,metnote,idum,nlwmet

      read(i3dmet) plon,plat,iutm,xorg,yorg,dx,dy,nx,ny,nz,iproj,istag, &
                   tlat1,tlat2
      read(i2dmet)
      read(ilumet)
      if (lbin_lwmask) read(ifswm)
      read(i3dmet)
      read(i2dmet)
      read(ilumet)
      if (lbin_lwmask) read(ifswm)
      read(i3dmet) 
      read(i2dmet) 
      read(ilumet) 
      if (lbin_lwmask) read(ifswm)

      if (et.eq.0.) then
        et = 24.
        ied = subday(ied)
      endif

      write(*,1005) 'Time/date range: ',jdate,bt,ied,et
 1005 format(1x,a,2(i10,f10.2))
      write(*,*) 'Time zone      : ',itzon
      write(*,*) 'Map projection : ',iproj
      write(*,*) 'UTM zone       : ',iutm
      write(*,*) 'Pole lon/lat   : ',plon,plat
      write(*,*) 'True lats      : ',tlat1,tlat2
      write(*,*) 'Grid origin    : ',xorg,yorg
      write(*,*) 'Grid cell size : ',dx,dy
      write(*,*) 'Grid dimensions: ',nx,ny,nz
      write(*,*)
      centlon = -itzon*15.       ! timezone central longitude 

!-----Write new output emissions file header

      write(iout) ifile,fileiden,1,nspecs,jdate,bt,ied,et
      write(iout) plon,plat,iutm,xorg,yorg,dx,dy,nx,ny,1,iproj,istag, &
                  tlat1,tlat2,0.
      write(iout) 1,1,nx,ny
      write(iout) ((mspec(i,j),i=1,10),j=1,nspecs)

!-----Allocate local variables

      allocate(salt_water(nx,ny), &
               iloc(nx,ny), &
               jloc(nx,ny), &
               rlat(nx,ny), &
               rlon(nx,ny), &
               rmet2d(nx,ny), &
               water(nx,ny), &
               sea_ice(nx,ny), &
               surf(nx,ny), &
               coastline_len(nx,ny), &
               srfzone_width(nx,ny), &
               area2d(nx,ny), &
               areasf(nx,ny), &
               rchla_conc(nx,ny), &
               rdms_conc(nx,ny), &
               ts(nx,ny), &
               u10(nx,ny), &
               v10(nx,ny), &
               rmet3d(nx,ny,nz), &
               pp(nx,ny,nz), &
               tt(nx,ny,nz), &
               wv(nx,ny,nz), &
               spd10m(nx,ny), &
               relhum(nx,ny), &
               flux(nx,ny,MXSPEC), &
               mass(nx,ny), &
               massdms(nx,ny), &
               mass_crs_fn(nx,ny,MXSPEC), &
               stat=ios)
      if (ios.ne.0) stop'ERROR: memory allocation failed'

      salt_water = 0
      surf = 0.
      if (lss) then

!-----Set weighting factors for uniform grid trapezoidal rule

        wt(0)  = 0.5
        do ibin = 1, NBINS-1
          wt(ibin)  = 1.0
        enddo
        wt(NBINS)  = 0.5

!-----Set valid particle size upper limit for the selected parameterization

        if (OO_SSSF .eq. 'OO_GO03') then ! see de Leeuw et al., 2011;
                                         !  0.1 < r80 < 8
          r80_lower_limit_oo = 0.1
          r80_upper_limit_oo = 8.0
        else if (OO_SSSF .eq. 'OO_SO11') then ! see Sofiev et al., 2011;
                                              !  0.01 < r80 < 10
          r80_lower_limit_oo = 0.01
          r80_upper_limit_oo = 10.0
        else if (OO_SSSF .eq. 'OO_OV14') then ! see Ovadnevaite et al., 2014;
                                              !  0.015 < r80 < 6
          r80_lower_limit_oo = 0.015
          r80_upper_limit_oo = 6.0
        else
          stop'ERROR: Unrecognized open ocean parameteriziation'
        endif

        if (SZ_SSSF .eq. 'SZ_DL00') then ! see de Leeuw et al., 2000;
                                         !  1.6 < D98 (= 4 x r80) < 20
          r80_lower_limit_sz = 0.4
          r80_upper_limit_sz = 5.0
        else if (SZ_SSSF .eq. 'SZ_GO03') then ! same as OO_GO03
          r80_lower_limit_sz = 0.1
          r80_upper_limit_sz = 8.0
        else
          stop'ERROR: Unrecognized surf zone parameteriziation'
        endif

      endif

!-----Calculate lat/lon at cell centroids
 
      if (iproj.eq.0) then
        do j = 1,ny
          yloc = yorg + dy*(float(j) - 0.5)
          do i = 1,nx
            xloc = xorg + dx*(float(i) - 0.5)
            rlat(i,j) = yloc
            rlon(i,j) = xloc
          enddo
        enddo
      else
        do j = 1,ny
          yloc = (yorg + dy*(float(j) - 0.5))/1000.
          do i = 1,nx
            xloc = (xorg + dx*(float(i) - 0.5))/1000.
            if (iproj.eq.1) then
              call utmgeo(1,iutm,xloc,yloc,rlon(i,j),rlat(i,j))
            elseif (iproj.eq.2) then
              call lcpgeo(1,plon,plat,tlat1,tlat2,xloc,yloc,rlon(i,j),rlat(i,j))
            elseif (iproj.eq.3) then
              call rpsgeo(1,plon,plat,xloc,yloc,rlon(i,j),rlat(i,j))
            elseif (iproj.eq.4) then
              call pspgeo(1,plon,plat,tlat1,xloc,yloc,rlon(i,j),rlat(i,j))
            elseif (iproj.eq.5) then
              call mrcgeo(1,plon,plat,tlat1,xloc,yloc,rlon(i,j),rlat(i,j))
            endif
          enddo
        enddo
      endif

!-----Read land use file and salt water mask to determine distribution of
!     ocean/coastline grid cells.

      iww = 0
      isi = 0
      read(ilumet,end=900)
      do n = 1,nlumet
        read(ilumet) idum,(namvar(i),i=1,10),((rmet2d(i,j),i=1,nx),j=1,ny)
        write(metspc,'(10a1)') (namvar(i),i=1,10)
        if (metspc.eq.'WATER') then
          iww = 1
          water(:,:) = rmet2d(:,:)
        endif
        if (metspc.eq.'ICE') then
          isi = 1
          sea_ice(:,:) = rmet2d(:,:)
        endif
      enddo
      if (iww.eq.0) goto 901
      if (isi.eq.0) goto 902
      close(ilumet)

!-----Read binary or text land/water mask file

      if (lbin_lwmask) then
        iow = 0
        isf = 0
        read(ifswm,end=903)
        do n = 1,nlwmet
          read(ifswm) idum,(namvar(i),i=1,10),((rmet2d(i,j),i=1,nx),j=1,ny)
          write(metspc,'(10a1)') (namvar(i),i=1,10)
          if (metspc.eq.'OPEN') then
            iow = 1
            water(:,:) = max(0.,rmet2d(:,:)) ! override water fraction
          endif
          if (metspc.eq.'SURF') then
            isf = 1
            surf(:,:) = max(0.,rmet2d(:,:))
            surf(:,:) = 0.5 * surf(:,:) ! surf zone width: 50m->25m
                                        ! (Gantt et al., 2015)
          endif
        enddo
        if (iow.eq.0) goto 904
        if (isf.eq.0) goto 905
      else
        do j = ny,1,-1
          read(ifswm,'(500(:,i1))') (salt_water(i,j),i=1,nx)
        enddo
      endif
      close(ifswm)

!-----Set coastline length, surfzone width, and ocean area

      do j = 1,ny
        do i = 1,nx
          coastline_len(i,j) = 0.
          srfzone_width(i,j) = 0. 
          if (iproj.ne.0) then
            area = water(i,j)*dx*dy
            areas = surf(i,j)*dx*dy
          else
            jlat = yorg + dy*(j - 0.5)
            area = abs(dx*111134.*cos(jlat/conv)*dy*111134.)
            area = water(i,j)*area
            areas = surf(i,j)*area
          endif
          if (lbin_lwmask) then
            area2d(i,j) = (1. - sea_ice(i,j))*area
            areasf(i,j) = (1. - sea_ice(i,j))*areas
          else ! areasf will be calculated in the next big do loop.
            area2d(i,j) = area
          endif

          if (.not.lbin_lwmask) then
            if (salt_water(i,j).eq.0) area2d(i,j) = 0.
            if (salt_water(i,j).eq.2) srfzone_width(i,j) = 10.
            if (salt_water(i,j).eq.3) srfzone_width(i,j) = 20.
            if (salt_water(i,j).eq.4) srfzone_width(i,j) = 50.
            if (salt_water(i,j).eq.5) srfzone_width(i,j) = 100.
            if (salt_water(i,j).gt.1) then
              if (iproj.ne.0) then
                clen = 0.5*(dx + dy)
              else
                clen = 0.5*(abs(dx*111134.*cos(jlat/conv)) + dy*111134.)
              endif
              coastline_len(i,j) = clen ! meters
            endif
          endif
        enddo
      enddo

      if (.not.lbin_lwmask ) then
        if (trim(icstpath).ne.'NONE') then
          open (unit=icst,file=icstpath,form='formatted')
          do i = 1,nx*ny 
            read(icst,*,end=99) inx,iny,coastline_len(inx,iny)
          enddo
 99       close(icst)
        endif
        areasf(:,:) = (1. - sea_ice(:,:)) * &
                      srfzone_width(:,:) * coastline_len(:,:)
      endif

!-----Read a halo-methane file, hard-coded for a specific format

      if (lhalo) then
        do i = 1,2 ! skip 2 lines
          read(ifchl,*)
        enddo
        read(ifchl,'(12x,i4)') nychl
        read(ifchl,'(12x,i4)') nxchl
        write(*,*) 'Chlorophyll-A file nx, ny: ',nxchl, nychl
        allocate(chla_conc(nxchl,nychl), stat=ios)
        if (ios.ne.0) then
          write(*,'(2a)') 'Allocation of 2D variables related with ', &
                          'Chlorophyll-A file failed'
          stop
        endif
        do i = 1,23 ! skip 23 lines
          read(ifchl,*)
        enddo
        read(ifchl,'(24x,i2,1x,i2,1x,i4)') mm1,dd1,year1
        read(ifchl,'(23x,i2,1x,i2,1x,i4)') mm2,dd2,year2
        write(*,'(a,i2,a,i2,a,i4)') 'Start date of Chl-A file:', &
             mm1,'/',dd1,'/',year1
        write(*,'(a,i2,a,i2,a,i4)') 'End date of Chl-A file:', &
             mm2,'/',dd2,'/',year2
        yy1 = mod(year1,100)
        yy2 = mod(year2,100)
        idate1 = yy1*10000 + mm1*100 + dd1
        idate2 = yy2*10000 + mm2*100 + dd2
        call juldate(idate1)
        call juldate(idate2)
        if ((jdate.lt.idate1).or.(jdate.gt.idate2)) then
          write(*,'(2a)')'WARNING - Chlorophyll-A file is outside', &
                         ' of the processing date'
        endif
        do i = 1,4 ! skip 4 lines
          read(ifchl,*)
        enddo
        read(ifchl,*) ((chla_conc(i,j),i=1,nxchl),j=nychl,1,-1)
        close(ifchl)

!-----Map regional domain grid cells to the global chlorophyll-A grid

        dxchl = 360./nxchl
        dychl = 180./nychl
        rchla_conc = 0.
        iloc(:,:) = int((rlon(:,:) - (-180.))/dxchl) + 1
        jloc(:,:) = int((rlat(:,:) - (-90.))/dychl) + 1

        write(*,*)
        write(*,'(a)') '=== GLOBAL CHL-A GRID CELL RANGE ==='
        write(*,'(5x,4(a,i4),a)') &
                       '('    ,iloc(1,ny),',',jloc(1,ny), &
                       ') - (',iloc(nx,ny),',',jloc(nx,ny),')'
        write(*,'(5x,a)') '    |           |'
        write(*,'(5x,4(a,i4),a)') &
                       '('    ,iloc(1,1),',',jloc(1,1), &
                       ') - (',iloc(nx,1),',',jloc(nx,1),')'
        write(*,'(a)') '===================================='
        write(*,*)

!-----Load global chlorophyll-A conc to regional grid

        do i = 1,nx
          do j = 1,ny
            if (chla_conc(iloc(i,j),jloc(i,j)).ge.0.0) then
              rchla_conc(i,j) = chla_conc(iloc(i,j),jloc(i,j))
              rchla_conc(i,j) = min(cap_chla,rchla_conc(i,j))
            endif
          enddo
        enddo
      endif

!-----Read DMS concentration file, hard-coded for a specific format

      if (ldms) then
        read(ifdms,'(1x,i3)') nxdms
        read(ifdms,'(1x,i3)') nydms
        write(*,*) 'DMS file nx, ny: ',nxdms, nydms
        allocate(dms_conc(nxdms,nydms), stat=ios)
        if (ios.ne.0) then
          write(*,'(2a)') 'Allocation of 2D variables related with ', &
                          'DMS file failed'
          stop
        endif
        do j = 1,nydms
         read(ifdms,*) dms_conc(:,j)
        end do
        close(ifdms)

!-----Map regional domain grid cells to the global DMS grid

        dxchl = 360./nxdms
        dychl = 180./nydms
        rdms_conc = 0.
        iloc(:,:) = int((rlon(:,:) - (-180.))/dxchl) + 1
        jloc(:,:) = int((90. - rlat(:,:))/dychl) + 1

        write(*,*)
        write(*,'(a)') '=== GLOBAL DMS GRID CELL RANGE ==='
        write(*,'(5x,4(a,i4),a)') &
                       '('    ,iloc(1,ny),',',jloc(1,ny), &
                       ') - (',iloc(nx,ny),',',jloc(nx,ny),')'
        write(*,'(5x,a)') '    |           |'
        write(*,'(5x,4(a,i4),a)') &
                       '('    ,iloc(1,1),',',jloc(1,1), &
                       ') - (',iloc(nx,1),',',jloc(nx,1),')'
        write(*,'(a)') '===================================='
        write(*,*)

!-----Load global DMS conc to regional grid
!     Fill in poorly-resolved coastline DMS 1x1 deg grid with local average

        do i = 1,nx
          do j = 1,ny
            if (dms_conc(iloc(i,j),jloc(i,j)).ge.0.0) then
              rdms_conc(i,j) = dms_conc(iloc(i,j),jloc(i,j))
            else
              sumdms = 0.
              nsum = 0
              do ii = -1,1
                do jj = -1,1
                  iii = iloc(i,j) + ii
                  jjj = jloc(i,j) + jj
                  if (iii .gt. nxdms) iii = 1
                  if (iii .lt. 1) iii = nxdms
                  if (jjj .gt. nydms) jjj = nydms - 1
                  if (jjj .lt. 1) jjj = 2
                  if (dms_conc(iii,jjj).lt.0.) cycle
                  sumdms = sumdms + dms_conc(iii,jjj)
                  nsum = nsum + 1
                enddo
              enddo
              if (nsum.gt.0) rdms_conc(i,j) = sumdms/float(nsum)
            endif
          enddo
        enddo
      endif

!-----Loop over hours in the 2D/3D met files

 100  continue

!-----Read the met

        its = 0
        iuu = 0
        ivv = 0
        read(i2dmet,end=906) idt1,hr1,idt2,hr2
        do n = 1,n2dmet
          read(i2dmet) idum,(namvar(i),i=1,10),((rmet2d(i,j),i=1,nx),j=1,ny)
          write(metspc,'(10a1)') (namvar(i),i=1,10)
          if (metspc.eq.'TSURF_K') then
            its = 1
            ts(:,:) = rmet2d(:,:)
          endif
          if (metspc.eq.'U10_MpS') then
            iuu = 1
            u10(:,:) = rmet2d(:,:)
          endif
          if (metspc.eq.'V10_MpS') then
            ivv = 1
            v10(:,:) = rmet2d(:,:)
          endif
        enddo
        if (its.eq.0) goto 907
        if (iuu.eq.0) goto 908
        if (ivv.eq.0) goto 909

        ipp = 0
        itt = 0
        iwv = 0
        read(i3dmet,end=910) idt1,hr1,idt2,hr2
        do n = 1,n3dmet
          do k = 1,nz
            read(i3dmet) idum,(namvar(i),i=1,10),((rmet3d(i,j,k),i=1,nx),j=1,ny)
          enddo
          write(metspc,'(10a1)') (namvar(i),i=1,10)
          if (metspc.eq.'PRESS_MB') then
            ipp = 1
            pp(:,:,:) = rmet3d(:,:,:)
          endif
          if (metspc.eq.'TEMP_K') then
            itt = 1
            tt(:,:,:) = rmet3d(:,:,:)
          endif
          if (metspc.eq.'HUMID_PPM') then
            iwv = 1
            wv(:,:,:) = rmet3d(:,:,:)
          endif
        enddo
        if (ipp.eq.0) goto 911
        if (itt.eq.0) goto 912
        if (iwv.eq.0) goto 913

!-----Check if this met hour extends past end date/time

        if (idt1.gt.ied .or. (idt1.eq.ied .and. hr1.gt.et)) goto 1000
        write(*,*)'Read met data at: ',idt1,hr1

!-----Calculate the 10 meter wind speed (m/s) and relative humidity
!     (decimal fraction)

        do j = 1,ny
          do i = 1,nx
            spd10m(i,j) = max(0.1,sqrt(u10(i,j)**2 + v10(i,j)**2))
            qwatr = 1.e-6*wv(i,j,1)*18./28.8
            ev = qwatr*pp(i,j,1)/(qwatr + 0.622)
            es = 6.11*exp((2.5e6/461.)*(1./273. - 1./tt(i,j,1)))
            relhum(i,j) = MIN(0.99,ev/es)
          enddo
        enddo

!-----Calculate sea salt emissions

        if (lss) then

!-----Calculate dry particle mass emissions [g/hr]

          mass = 0.
          do j = 1,ny
            do i = 1,nx
              if (area2d(i,j).lt.100. .and. areasf(i,j).eq.100.) CYCLE
              spd10 = max(minspd,spd10m(i,j))

!-----Convert D_MAX to r80 based on Eq.(2) of Zhang et al., 2006
!     Eq.(2) is valid over RH of 0.45-0.99, but marine atmospheric RH is almost
!     always greater than 0.45 (Lewis and Schwartz, 2004)

              r80_min = 1.82*((1. - relhum(i,j))/(2. - relhum(i,j)))**0.33333 &
                        *D_MIN/2.
              r80_max = 1.82*((1. - relhum(i,j))/(2. - relhum(i,j)))**0.33333 &
                        *D_MAX/2.
              r80_min_oo = MAX(r80_min, r80_lower_limit_oo)
              r80_max_oo = MIN(r80_max, r80_upper_limit_oo)
              if (r80_min_oo .gt. r80_max_oo) then
                write(*,*) 'ERROR: invalid integration range for OO SSSF' &
                          ,' at (',i,',',j,'):'
                write(*,*) '[',r80_min_oo,' - ',r80_max_oo,']'
                stop
              endif
              sum_vol = 0.0    ! total dry particle volume flux for open ocean
                               !  [um^3 m^-2 s^-1]

!-----Gong (2003) parameterization
!     Numerical integration using uniform grid trapezoidal rule

              if (OO_SSSF .eq. 'OO_GO03') then
                call get_bin_size(NBINS,r80_min_oo,r80_max_oo,logwidth,r80)

                term1 = 1.373*(spd10**3.41)
                do ibin = 0,NBINS
                  A_GO03 = 4.7*(1. + theta_GO03*r80(ibin)) &
                               **(-0.017*r80(ibin)**(-1.44))
                  B_GO03 = (0.433 - log10(r80(ibin)))/0.433
                  dF_dDp = (1. + 0.057*r80(ibin)**3.45)/r80(ibin)**A_GO03 &
                           *10.**(1.607*exp(-B_GO03*B_GO03))
                  dF_dlogDp = LN10*r80(ibin)*dF_dDp
                  dV_dlogDp = dF_dlogDp*PI_6*r80(ibin)**3.
                  sum_vol = sum_vol + wt(ibin)*dV_dlogDp
                enddo
                sum_vol = term1*logwidth*sum_vol

!-----Sofiev et al. (2011) parameterization, with Ts and salinity adjustment
!     Numerical integration using uniform grid trapezoidal rule

              else if (OO_SSSF .eq. 'OO_SO11') then
                call get_bin_size(NBINS,r80_min_oo,r80_max_oo,logwidth,r80)

                whtcap = 3.84e-6*(spd10**3.41)
                do ibin = 0,NBINS
                  aa     = -0.09/(r80(ibin) + 3.e-3)
                  term1  = exp(aa)/(2. + exp(-5./r80(ibin)))
                  term2  = (1. + 0.05*r80(ibin)**1.05)/r80(ibin)**3.
                  bb     = (0.27 - log10(r80(ibin)))/1.1
                  term3  = 10.**(1.05*exp(-bb*bb))

                  ftw298 = 1.
                  ftw288 = 0.480*r80(ibin)**-0.36
                  ftw278 = 0.150*r80(ibin)**-0.88
                  ftw271 = 0.092*r80(ibin)**-0.96
                  ftw = ftw298
                  if (ts(i,j).ge.288. .and. ts(i,j).lt.298.) then
                    ftw = (log10(ftw298) - log10(ftw288))/10.*(ts(i,j) - 288.)
                    ftw = ftw288*10.**ftw
                  elseif (ts(i,j).ge.278. .and. ts(i,j).lt.288.) then
                    ftw = (log10(ftw288) - log10(ftw278))/10.*(ts(i,j) - 278.)
                    ftw = ftw278*10.**ftw
                  elseif (ts(i,j).ge.271. .and. ts(i,j).lt.278.) then
                    ftw = (log10(ftw278) - log10(ftw271))/7.*(ts(i,j) - 271.)
                    ftw = ftw271*10.**ftw
                  elseif (ts(i,j).lt.271.) then
                    ftw = ftw271
                  endif

                  fsw33 = 1.
                  fsw09 = 0.12*r80(ibin)**-0.71
                  fsw = fsw33
                  if (sal.ge.9.2 .and. sal.lt.33.) then
                    fsw = (log10(fsw33) - log10(fsw09))/23.8*(sal - 9.2)
                    fsw = fsw09*10.**fsw
                  elseif (sal.lt.9.2) then
                    fsw = 0.
                  endif

                  dF_dDp = 1.e6*term1*term2*term3*ftw*fsw
                  dF_dlogDp = LN10*r80(ibin)*dF_dDp
                  dV_dlogDp = dF_dlogDp*PI_6*r80(ibin)**3.
                  sum_vol = sum_vol + wt(ibin)*dV_dlogDp
                enddo
                sum_vol = whtcap*logwidth*sum_vol

!-----Ovadnevaite et al. (2014) parameterization - Ts and salinity dependence
!     through Reynolds number.  Solve integration using err function

              else if (OO_SSSF .eq. 'OO_OV14') then
                Cd = (0.8 + 0.065*spd10)*1.e-3 ! drag coefficient
                                                     !  (Wu, 1980)
                Hs = sclfac_Hs*0.025*spd10**2. ! significant wave height
                                                     !  (Alves et al., 2003)
                ReHw = sqrt(Cd)*spd10*Hs/Vis(ts(i,j), sal)
                do imode = 1,5
                  term1 = 3.*LN10*log10(sig_OV14(imode))
                  term2 = log10(Dg_OV14(imode)) + term1*log10(sig_OV14(imode))
                  xx1 = (log10(r80_min_oo) - term2)/ &
                        (sqrt(2.)*log10(sig_OV14(imode)))
                  xx2 = (log10(r80_max_oo) - term2)/ &
                        (sqrt(2.)*log10(sig_OV14(imode)))
                  term3 = 3.*LN10*log10(Dg_OV14(imode)) + term1*term1/2.
                  sum_vol = sum_vol + (PI_6/2.)*exp(term3)*(A_OV14(imode)* &
                            MAX(ReHw - B_OV14(imode),0.)**C_OV14(imode))* &
                            (xerf(xx2) - xerf(xx1))
                enddo
              endif

!-----Surf zone sea salt generation

              sum_vol_sz = 0.0 ! total dry particle volume flux for surf zone
                               !  [um^3 m^-2 s^-1]
              if (salt_water(i,j).gt.1 .or. surf(i,j).gt.minsw) then
                r80_min_sz = MAX(r80_min, r80_lower_limit_sz)
                r80_max_sz = MIN(r80_max, r80_upper_limit_sz)
                if (r80_min_sz .gt. r80_max_sz) then
                  write(*,*) 'ERROR: invalid integration range for SZ SSSF' &
                            ,' at (',i,',',j,'):'
                  write(*,*) '[',r80_min_sz,' - ',r80_max_sz,']'
                  stop
                endif

!-----de Leeuw et al. (2000) parameterization for surf zone:
!       dF/dDf = 1.1e7 * exp(0.23 * U10) / Df^1.65
!     Zhang et al. (2006):
!       dF/dr80 = 2 * 3.7 * [((1-0.8)/(2-0.8))(35/sal)]^(1/3) * dF/dDf
!     Dry particle volume flux:
!       dV/dDp = dF/dDp * (Pi/6) * Dp^3; dF/dDp = dF/dr80
!     Solve integration analytically

                if (SZ_SSSF .eq. 'SZ_DL00') then
                  term1 = 1.1e7*exp(0.23*min(spd10,9.0))
                  term2 = 4.072377*(35./sal)**0.33333
                  sum_vol_sz = term1*term2**(-0.65)*PI_6* &
                               (r80_max_sz**2.35 - r80_min_sz**2.35)/2.35

!-----Gong (2003) with 100% whitecap
!     Numerical integration using uniform grid trapezoidal rule

                elseif (SZ_SSSF .eq. 'SZ_GO03') then
                  call get_bin_size(NBINS,r80_min_sz,r80_max_sz,logwidth,r80)

                  term1 = 3.576e5
                  do ibin = 0,NBINS
                    A_GO03 = 4.7*(1. + theta_GO03*r80(ibin)) &
                             **(-0.017*r80(ibin)**(-1.44))
                    B_GO03 = (0.433 - log10(r80(ibin)))/0.433
                    dF_dDp = (1. + 0.057*r80(ibin)**3.45)/r80(ibin)**A_GO03 &
                             *10.**(1.607*exp(-B_GO03 * B_GO03))
                    dF_dlogDp = LN10*r80(ibin)*dF_dDp
                    dV_dlogDp = dF_dlogDp*PI_6*r80(ibin)**3.
                    sum_vol_sz = sum_vol_sz + wt(ibin)*dV_dlogDp
                  enddo
                  sum_vol_sz = term1*logwidth*sum_vol_sz
                endif
              endif

!-----Convert volume flux to mass emission: [um^3 m^-2 s^-1] -> [g/hr]

              mass(i,j) = (sum_vol*area2d(i,j) + sum_vol_sz*areasf(i,j)) &
                          *dens_NaCl*1.e-15*3600.
            enddo
          enddo
        endif

!-----Calculate halo-methane emissions
 
        if (lhalo) then
          flux = 0.
          do j = 1,ny
            do i = 1,nx
              if (area2d(i,j).lt.100.) CYCLE
              ilochr = 1 + int(hr1) - int((centlon + 7.5 - rlon(i,j))/15.)
              if (ilochr.lt.1) ilochr = 24 + ilochr
              if (ilochr.gt.24) ilochr = ilochr - 24
              factor = 4.31e-8*area2d(i,j)*fhrly(ilochr)*rchla_conc(i,j)
              flux(i,j,1) = factor*fhalo(iCH3I)   ! mol/hr-cell
              flux(i,j,2) = factor*fhalo(iMB3)    ! mol/hr-cell
              flux(i,j,3) = factor*fhalo(iMB2)    ! mol/hr-cell
              flux(i,j,4) = factor*fhalo(iMBC)    ! mol/hr-cell
              flux(i,j,5) = factor*fhalo(iMB2C)   ! mol/hr-cell
              flux(i,j,6) = factor*fhalo(iMBC2)   ! mol/hr-cell
              flux(i,j,7) = factor*fhalo(iMIC)    ! mol/hr-cell
              flux(i,j,8) = factor*fhalo(iMIB)    ! mol/hr-cell
              flux(i,j,9) = factor*fhalo(iMI2)    ! mol/hr-cell
              flux(i,j,10) = 2.39e-8*area2d(i,j)  ! I2 mol/hr-cell from
                                                  !  4.E8 molecules/cm2/s
              flux(i,j,11) = rchla_conc(i,j) ! mg/m3
            enddo
          enddo
        endif

!-----Calculate DMS emissions; oceanic DMS flux parameterized according to
!     Lana et al. (2011)

        if (ldms) then

          massdms=0.0
          do j = 1,ny
            do i = 1,nx
             if (area2d(i,j).lt.100.) CYCLE

             ka = 659.*spd10m(i,j)*sqDMS
             alpha = R_gas*ts(i,j)/exp(12.64 - 3547./ts(i,j))
             Tc = min(47.,ts(i,j) - 273.15)
             SC = 2674. + Tc*(-147.12 + Tc*(3.726 - Tc*0.038))
             if (DMS_KW .eq. 'NI00') then
               k600 = 0.222*spd10m(i,j)*spd10m(i,j) + 0.333*spd10m(i,j)
             else
               if (spd10m(i,j).le.3.6) then
                 k600 = 0.17*spd10m(i,j)
               elseif (spd10m(i,j).le.13.) then
                 k600 = 2.85*spd10m(i,j) - 9.65
               else
                 k600 = 5.9*spd10m(i,j) - 49.3
               endif
             endif
             kw = k600*sqrt(600./SC)
             if (kw .le. 1.e-4) then 
               gamma_a = 0.012
             else
               gamma_a = 1./(1. + ka/(alpha*kw))
             endif

!-----Convert DMS oceanic concentrations in nM to mass emission -> [mol/hr]
!     [1e-2 m h^-1] [1e-6 mol m^-3][m^2] -> [g/hr]

             massdms(i,j) = kw*(1. - gamma_a)*rdms_conc(i,j)*1.e-8*area2d(i,j)
            enddo
          enddo
        endif 

!-----Partition emissions into component species, if necessary

        do l = 1,nspecs
          write(c10,'(10a1)') (mspec(m,l),m=1,10)
          if (trim(c10) .eq. 'CH3I') then
            mass_crs_fn(:,:,l) = flux(:,:,1)
          endif
          if (trim(c10) .eq. 'MB3') then
            mass_crs_fn(:,:,l) = flux(:,:,2)
          endif
          if (trim(c10) .eq. 'MB2') then
            mass_crs_fn(:,:,l) = flux(:,:,3)
          endif
          if (trim(c10) .eq. 'MBC') then
            mass_crs_fn(:,:,l) = flux(:,:,4)
          endif
          if (trim(c10) .eq. 'MB2C') then
            mass_crs_fn(:,:,l) = flux(:,:,5)
          endif
          if (trim(c10) .eq. 'MBC2') then
            mass_crs_fn(:,:,l) = flux(:,:,6)
          endif
          if (trim(c10) .eq. 'MIC') then
            mass_crs_fn(:,:,l) = flux(:,:,7)
          endif
          if (trim(c10) .eq. 'MIB') then
            mass_crs_fn(:,:,l) = flux(:,:,8)
          endif
          if (trim(c10) .eq. 'MI2') then
            mass_crs_fn(:,:,l) = flux(:,:,9)
          endif
          if (trim(c10) .eq. 'I2') then
            mass_crs_fn(:,:,l) = flux(:,:,10)
          endif
          if (trim(c10) .eq. 'CHLA') then
            mass_crs_fn(:,:,l) = flux(:,:,11)
          endif
          if (trim(c10) .eq. 'NA') then
            mass_crs_fn(:,:,l) = mass(:,:)*efactor*SPCFAC_NA
          endif
          if (trim(c10) .eq. 'PCL') then
            mass_crs_fn(:,:,l) = mass(:,:)*efactor*SPCFAC_CL
          endif
          if (trim(c10) .eq. 'PSO4') then
            mass_crs_fn(:,:,l) = mass(:,:)*efactor*SPCFAC_SO4
          endif
          if (trim(c10) .eq. 'SSCL') then
            mass_crs_fn(:,:,l) = mass(:,:)*efactor*SPCFAC_CL*SPCFAC_SSCL
          endif
          if (trim(c10) .eq. 'SSBR') then
            mass_crs_fn(:,:,l) = mass(:,:)*efactor*SPCFAC_CL*SPCFAC_SSCL* &
                                 SPCFAC_SSBR
          endif
          if (trim(c10) .eq. 'WIND') then
            mass_crs_fn(:,:,l) = spd10m(:,:)
          endif
          if (trim(c10) .eq. 'RELH') then
            mass_crs_fn(:,:,l) = relhum(:,:)
          endif
          if (trim(c10) .eq. 'SRFT') then
            mass_crs_fn(:,:,l) = ts(:,:)
          endif
          if (trim(c10) .eq. 'DMS' .OR. trim(c10) .eq. 'SO2') then
            mass_crs_fn(:,:,l) = massdms(:,:)
          endif
          if (trim(c10) .eq. 'DMSWC') then
            mass_crs_fn(:,:,l) = rdms_conc(:,:)
          endif
        enddo 

!-----Write out emission records

        hr2 = hr2 + 1.
        if (hr2.gt.24.) then
          hr2 = hr2 - 24.
          idt2 = addday(idt2)
        endif
        write(*,1005) 'Writing emissions: ',idt1,hr1,idt2,hr2
        write(iout) idt1,hr1,idt2,hr2
        do l = 1,nspecs
          write(iout) 1,(mspec(m,l),m=1,10),((mass_crs_fn(i,j,l),i=1,nx),j=1,ny)
        enddo

      goto 100

!-----Met file error handling

 900  continue
      write(*,*)'Unexpected end of CAMx LU file'
      stop
 901  continue
      write(*,*)'Did not find WATER lanudse'
      stop
 902  continue
      write(*,*)'Did not find ICE lanudse'
      stop

 903  continue
      write(*,*)'Unexpected end of binary watermask file'
      stop
 904  continue
      write(*,*)'Did not find OPEN from lwmask'
      stop
 905  continue
      write(*,*)'Did not find SURF from lwmask'
      stop

 906  continue
      write(*,*)'Unexpected end of CAMx 2D file'
      stop
 907  continue
      write(*,*)'Did not find TSURF'
      stop
 908  continue
      write(*,*)'Did not find U10'
      stop
 909  continue
      write(*,*)'Did not find V10'
      stop

 910  continue
      write(*,*)'Unexpected end of CAMx 3D file'
      stop
 911  continue
      write(*,*)'Did not find PRESS'
      stop
 912  continue
      write(*,*)'Did not find TEMP'
      stop
 913  continue
      write(*,*)'Did not find HUMID'
      stop

 1000 write(*,*)'End of met data'
      write(*,*)'Normal Completion'
      end

!===============================================================================

      subroutine get_bin_size(ns,s_min,s_max,swidth,sa)
      implicit none
!
!     Set log-uniform distribution of size points
!
      integer :: ns, is
      real :: s_min, s_max, swidth
      real :: sa(0:ns)

      swidth = (log10(s_max) - log10(s_min)) / REAL(ns)
      sa(0) = s_min
      do is = 1, ns-1
        sa(is) = 10.**(log10(sa(is-1)) + swidth)
      enddo
      sa(ns) = s_max

      return
      end

!===============================================================================

      real function Vis(sst, sal)
      implicit none
!
!     Calculates kinematic viscosity of sea water
!     http://web.mit.edu/seawater/ (Version: 8 July 2016)
!
      real :: sst ! sea surface temperature [K]
      real :: sal ! ocean salinity [g/kg]

      real, parameter, dimension(10) :: a1 = (/ 1.5700386464E-01, &
                                                6.4992620050E+01, &
                                               -9.1296496657E+01, &
                                                4.2844324477E-05, &
                                                1.5409136040E+00, &
                                                1.9981117208E-02, &
                                               -9.5203865864E-05, &
                                                7.9739318223E+00, &
                                               -7.5614568881E-02, &
                                                4.7237011074E-04 /)

      real, parameter, dimension(5)  :: a2 = (/ 9.9992293295E+02, &
                                                2.0341179217E-02, &
                                               -6.1624591598E-03, &
                                                2.2614664708E-05, &
                                               -4.6570659168E-08 /)

      real, parameter, dimension(5)  :: b2 = (/ 8.0200240891E+02, &
                                               -2.0005183488E+00, &
                                                1.6771024982E-02, &
                                               -3.0600536746E-05, &
                                               -1.6132224742E-05 /)

      real :: tc ! SST in degree C
      real :: sal0 ! salinity in [g/g]
      real :: mu_w, AA, BB, mu
      real :: rho_w, D_rho, rho_sw_sharq

      tc = sst - 273.
      sal0 = sal / 1000.

      mu_w = a1(4) + 1./(a1(1) * (tc + a1(2)) * (tc + a1(2)) + a1(3))
      AA   = a1(5) + a1(6) * tc + a1(7) * tc * tc
      BB   = a1(8) + a1(9) * tc + a1(10) * tc * tc
      mu   = mu_w * (1. + AA * sal0 + BB * sal0 * sal0);

      rho_w = a2(1) + (a2(2) + (a2(3) + (a2(4) + a2(5) * tc) * tc) * tc) * tc
      D_rho = b2(1) * sal0 + &
             (b2(2) + (b2(3) + b2(4) * tc + b2(5) * sal0) * tc) * sal0 * tc
      rho_sw_sharq = rho_w + D_rho

      Vis = mu / rho_sw_sharq

      return
      end

!===============================================================================

      real function xerf(x)
      implicit none
!
!     Returns the error function erf(x)
!     Numerical Recipes in Fortran, 2nd Ed., 1992
!
      real erfcc, x
      real t, z
      z=abs(x)
      t=1./(1.+0.5*z)
      erfcc=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+ &
           t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+ &
           t*(1.48851587+t*(-.82215223+t*.17087277)))))))))
      if (x.lt.0.) erfcc=2.-erfcc

      xerf = 1.-erfcc ! return erf(x)

      return
      end

!===============================================================================

      integer function addday(jdate)
      implicit none
      integer jdate,iyr,idy
      iyr = jdate/1000
      idy = jdate - iyr*1000
      if ((mod(iyr,4).eq.0 .and. idy.eq.366) .or. &
          (mod(iyr,4).ne.0 .and. idy.eq.365)) then
        iyr = iyr + 1
        if (iyr.gt.99) iyr = 0
        addday = iyr*1000 + 1
      else
        addday = jdate + 1
      endif
      end

!===============================================================================

      integer function subday(jdate)
      implicit none
      integer jdate,iyr,idy
      iyr = jdate/1000
      idy = jdate - iyr*1000
      if (idy.eq.1) then
        iyr = iyr - 1
        if (iyr.lt.0) iyr = 99
        if (mod(iyr,4).eq.0) then
          idy = 366
        else
          idy = 365
        endif
        subday = iyr*1000 + idy
      else
        subday = jdate - 1
      endif
      end
