      program seasalt

!  This program generates gridded oceanic emissions for CAMx:
!      Sea salt PM component species ("NA", "PCL", "PSO4" g/hr)
!      Sea salt Cl component ("SS" for nitryl chloride gas-phase
!                             mechanism, mol/hr)
!      Methyl iodide ("CH3I" for iodine gas-phase mechanism, mol/hr)
!
!  NOTE: THIS PROGRAM READS CAMx V6 MET FILE FORMATS
!        USE METCONVERT FIRST IF YOU HAVE EARLIER MET FORMATS
!  
!  1. Sea salt emissions
!     The contribution from the open ocean is considered as well as the aerosol
!     flux from breaking waves in the surf zone.  Fluxes are parameterized as a
!     function of wet aerosol radius at RH=80%, which is roughly the same as dry
!     aerosol diameter.
!
!     For particles with dry radii > 4 microns, sea salt particle flux from the
!     open ocean is calculated using the Smith and Harrison (1998) 
!     parameterization.  For particles with dry radii < 4 microns, either the 
!     Gong (2003) or Monahan et al. (1986) parameterization may be used. 
!
!     The surf zone contribution is calculated using the parameterization of 
!     DeLeeuw et al. (2000) or the Gong (2003) open ocean parameterization
!     assuming 100% whitecap coverage.  
!
!     User must supply an input file that specifies cells containing ocean
!     and coastline.  More information on this is given in the documentation.
! 
!     The flux integral is performed over a default distribution of dry diameter
!     bins from Grini et al (2002) and between limits specified by the user in
!     the job script.  The user must provide the total size range over
!     which to integrate (e.g., 0-2.5 micros dry diameter).
!
!     The flux represents the rate of production of marine aerosol droplets
!     per unit area of the sea surface per increment of droplet radius.
!
!     Indirect production F0 is associated with the bursting of whitecap
!     bubbles.  Direct production F1 comes from mechanical disruption of wave
!     crests by the wind. Add F1 and F0 to get the total droplet production 
!     for each size bin.
!
!     From the particle flux, the emitted particle mass is calculated assuming
!     spherical particle geometry, and a dry density of 2250 kg m-3 (Grini et
!     al., 2002).
!
!     NOTE: ASSUMES OCEAN SURFACES ARE ICE-COVERED (0 flux) when Ts < -2C
!
!  2. Halo-methane emissions
!     10 halo-methane emissions are calculated, including CH3I, using
!     a satellite image of chlorophyll-a.  Fluxes of CH3I are tuned to yeild
!     global emissions consistent with Bell et al. (2002).  See documentation
!     about where to get chlorophyll-a data.
!
!  References:
!
!  Lewis, E. and Schwartz, S., 2004:  "Sea Salt Aerosol Production:
!  Mechanisms, Methods, and Models-A Critical Review".  American
!  Geophyscial Union, p.222.
!
!  Gong, S.L., 2003: "A parameterization for sea-salt aerosol 
!  source function for sub- and super-micron particle sizes".  
!  Biogeochemical Cycles, (17), 1097-1104.
!
!  Grini, A. et al. 2002: "Modeling the annual cycle of sea salt
!  in the Global 3D Model Oslo CTM-2: Concentrations, fluxes and
!  radiative impacts."  J. Climate, (15).
!
!  Smith, M.H. and Harrison, N.M., 1998: "The Sea Spray Generation
!  Function".  J. Aerosol. Sci. (29), S189-S190.
!
!  Monahan, E.C., et al. 1986:  "A model of marine aerosol 
!  generation via whitecaps and wave disruption", 'Oceanic
!  Whitecaps', p. 167-174.
!
!  DeLeeuw, G., Neele, F.P. Hill, M., Smith, M., and Vignati, E. (2000)
!  Production of sea spray aerosol in the surf zone.  J. Geophys. Res.
!  105, p. 29397.
!
!  Bell, N., Hsu L., Jacob., D. J., Schultz, M. G., Blake, D. R.,
!  Butler, J. H., King, D. B., Lobert, J. M., and Maier-Reimer, E.,
!  2002, "Methyl iodide: Atmospheric budget and use as a tracer of
!  marine convection in global models". J. Geophys. Res. 107, D17,
!  4340
!
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!c Copyright (C) 2012  ENVIRON International Corporation
!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 write to the Free Software Foundation, Inc.,
!c 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!
 
!----Declare variables

      implicit none  

      integer,parameter :: MXSPEC = 16
      integer,parameter :: MAXBINS = 16
      real,parameter    :: PI = 3.14159

      integer i, j, k, r, l, m, n
      integer ihr, idt, num_bins
      integer jdate, nx, ny, nz
      integer i2dmet, i3dmet, icrmet, ilumet, iout, icst, ifswm
      integer n3dmet, n2dmet, ncrmet, nlumet
      integer itimez, itzon, ibd, ied
      integer upperbin, lowerbin
      integer nspecs
      integer iutm, iproj, istag, idum
      integer inx, iny
      integer izz, ipp, itt, iuu, ivv, iww, its, icd
      integer ifchl
      integer nxchl, nychl                                     
      integer mm1, mm2, dd1, dd2, year1, year2                 
      integer yy1, yy2, idate1, idate2                         
      integer ios                                              
      integer col, row                                         
      integer maxi,mini,maxj,minj                              
      integer ilochr                                           
      integer iMB3,iMB2,iMBC,iMB2C,iMBC2,iCH3I,iMIC,iMIB,iMI2  
      integer*4 today(3), now(3) 

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

      integer,allocatable,dimension(:,:) :: col_loc,row_loc
      real hr, xorg, yorg
      real plon, plat, tlat1, tlat2
      real bt, et
      real,allocatable,dimension(:,:,:) :: pp,zz,tt,uu,vv,rmet3d
      real,allocatable,dimension(:,:) :: water,ts,spd10m,sea_ice, &
                 fcloud,cldtrns,coastline_len,srfzone_width,rlat,rlon,rmet2d
      real term1, term2, term3, term4, B, minspd
      real density
      real bin_ctr(MAXBINS),bin_width(MAXBINS)
      real binctr_diam(MAXBINS), binwdt_diam(MAXBINS)
      real radius(MAXBINS)
	  
      real tempc, tempk
      real henry, rgas
      real schmidt, kw, ch3i_g, ch3i_aq, doc, alpha
      real coszen, tau, solflux
      real prod, kcl, depth, Mcl
      real Mw
      real aersl_flux, dry_mass, area

      real,allocatable,dimension(:,:,:) :: flux,dflux0,dflux1,dflux_srf, &
                                           mass_crs_fn
      real,allocatable,dimension(:,:) ::  mass,rchla_conc,area2d

      real bin_min,bin_max
      real dry_srf_mass
      real srf_upper_limit, srf_lower_limit
      real theta, theta0, press0, dtheta
      real thetabar, wind, ri, deltaz, uwnd, vwnd
      real zscale, cm, fm, ustar2, ustar, wr
      real eval_r, dk
      real term_0, A_g, theta_g, dens_NaCl, upper_limit, lower_limit
      real hhmm, zenith
      real R01, R02, f1, f2, A1, A2 
      real xloc, yloc, dx, dy, clen
      real vk, grav, gamma, z0

      real dxchl, dychl, x1, y1, lat1, lon1, dmin
      real x2, y2, dmin1, dmin2, lat2, lon2
      real f(9), fhrly(24), dxdy, centlon
      real jlat, conv
      real cap_chla
      real, allocatable :: chla_conc(:,:)
      real, allocatable :: lon_chl(:,:), lat_chl(:,:)
     
      logical ldark, lss, lnhalo

      character*4 ifile(10), fname(10), mspec(10,MXSPEC), namvar(10)
      character*4 fileiden(60), metnote(60)
      character*10 c10, metspc
      character*60 c60
      character*80 inline, method, method_s
      character*120 ipath, icstpath, ichla

!----Sea-salt speciation factors based on dry mass derived from
!    composition of seawater, and requiring charge balance
 
      REAL, PARAMETER ::  SPCFAC_NA   = 0.528
      REAL, PARAMETER ::  SPCFAC_CL   = 0.414
      REAL, PARAMETER ::  SPCFAC_SO4  = 0.0569
      REAL, PARAMETER ::  SPCFAC_SS   = 0.028169  !Converts g Cl/hr to mol Cl/hr
 
!----Data Statements
      
      data num_bins /MAXBINS/
      data binctr_diam /0.0385, 0.0585, 0.09, 0.135, 0.205, 0.31, 0.47, 0.72, &
                        1.085, 1.65, 2.55, 3.9, 5.9, 8.95, 13.6, 20.7/
      data binwdt_diam /0.017, 0.023,  0.04, 0.05,  0.09,  0.12, 0.20, 0.3, &
                        0.43,  0.7,  1.1,  1.6, 2.4, 3.7,   5.6,  8.6/

      data dens_NaCl /2250./     !density of dry NaCl kg/m3
      data density   /1150./     !density of aerosol at RH=80% kg/m3

      data minspd  /2.0/         !min wind speed for flux estimates (m/s)

      data vk/0.4/, grav/9.8/, gamma/0.286/
      data rgas/8.314/, alpha/1.71e-11/, depth/5.e3/
      data Mcl/0.54/, Mw/141.939/
      data conv/57.29578/ ! degree/radian 
      data cap_chla /1./  ! 1 ug/m3 cap on chlorophyll-A
      data iout   /10/
      data icst   /11/
      data i3dmet /12/
      data i2dmet /13/
      data icrmet /14/
      data ilumet /15/
      data ifswm  /17/
      data ifchl  /18/

      data iMB3   /1/
      data iMB2   /2/                                           
      data iMBC   /3/                                           
      data iMB2C  /4/                                           
      data iMBC2  /5/                                           
      data iCH3I  /6/                                           
      data iMIC   /7/                                           
      data iMIB   /8/                                           
      data iMI2   /9/

!             MB3  MB2  MBC  MB2C MBC2 CH3I MIC  MIB  MI2      
      data f /2.00,0.50,0.08,0.12,0.10,1.60,1.41,0.42,0.46/    

!                     1      2      3      4      5      6     
!                     7      8      9     10     11     12     
!                    13     14     15     16     17     18     
!                    19     20     21     22     23     24     
      data fhrly /0.032, 0.032, 0.032, 0.033, 0.034, 0.036, &  
                  0.039, 0.044, 0.051, 0.057, 0.062, 0.064, &  
                  0.062, 0.057, 0.051, 0.044, 0.039, 0.036, &  
                  0.034, 0.033, 0.032, 0.032, 0.032, 0.032/

      data ifile /'E','M','I','S','S','I','O','N','S',' '/

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

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

      read(*,'(20x,i)') jdate
      write(*,*) 'Julian date    :',jdate
      read(*,'(20x,i)') itimez
      write(*,*) 'Time zone      :',itimez
      centlon = - (itimez * 15.) ! timezone central longitude 
      read(*,'(20x,a)') method
      write(*,*) 'Open ocean     :',trim(method)
      read(*,'(20x,a)') method_s
      write(*,*) 'Surf zone      :',trim(method_s)
      read(*,'(20x,a)') icstpath     
      write(*,*) 'Coastline length filename:',trim(icstpath)
      read(*,'(20x,f)') lower_limit
      write(*,*) 'Min SS size    :',lower_limit
      read(*,'(20x,f)') upper_limit
      write(*,*) 'Max SS size    :',upper_limit

!-----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 identifier   :',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.
      lnhalo = .false.
      do i = 1,nspecs
        read(*,'(20x,a10)') c10
        write(*,*) c10
        if (trim(c10).ne.'PCL'  .and. trim(c10).ne.'NA'   .and. &
            trim(c10).ne.'PSO4' .and. trim(c10).ne.'SS'   .and. &
            trim(c10).ne.'WIND' .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') then
!           trim(c10).ne.'WIND' .and. trim(c10).ne.'SOLFLUX' .and. &
          write(*,'(2a)')'Unrecognized output species name: ',c10
          write(*,*)' Allowed names are:'
          write(*,*)'PCL'
          write(*,*)'NA'
          write(*,*)'PSO4'
          write(*,*)'SS'
          write(*,*)'MB3'
          write(*,*)'MB2'
          write(*,*)'MBC'
          write(*,*)'MB2C'
          write(*,*)'MBC2'
          write(*,*)'CH3I'
          write(*,*)'MIC'
          write(*,*)'MIB'
          write(*,*)'MI2'
          write(*,*)'I2'
          write(*,*)'WIND (diagnostic only)'
!         write(*,*)'SOLFLUX (diagnostic only)'
          write(*,*)'CHLA (diagnostic only)'
          stop
        endif
        if (trim(c10).eq.'PCL'  .or. trim(c10).eq.'NA' .or. &
            trim(c10).eq.'PSO4' .or. trim(c10).eq.'SS' .or. &
            trim(c10).eq.'WIND') 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'   .or. &
            trim(c10).eq.'CHLA') then
          lnhalo = .true.
        endif
        do j = 1,10
          mspec(j,i)(1:1) = c10(j:j)
        enddo
      enddo

!-----CAMx meteorological 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(icrmet,file=ipath,form='unformatted',status='old')
!     write(*,*) 'Input CAMx cloud/rain 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
      open(ifswm,file=ipath,form='formatted', status='old')
      write(*,*) 'Salt water mask file       :',trim(ipath)

!-----Halo-methane algorithm

      if (lnhalo) then
        read(*,'(20x,A)') ichla
        open(ifchl,file=ichla,form='formatted',status='old')
        write(*,*) 'Chlorophyll-A file       :',trim(ichla)
!     else
!       read(*,'(20x,a)') inline
!       read(inline,*) doc,ch3i_g
!       write(*,*) 'DOC and CH3I(g):',doc,ch3i_g
      endif

!-----Read input met file headers

      read(i3dmet) fname,metnote,itzon,n3dmet,ibd,bt,ied,et
      read(i2dmet) fname,metnote,itzon,n2dmet,ibd,bt,ied,et
!     read(icrmet) fname,metnote,itzon,ncrmet,ibd,bt,ied,et
      read(ilumet) fname,metnote,idum, nlumet,ibd,bt,ied,et
      read(i3dmet) plon,plat,iutm,xorg,yorg,dx,dy,nx,ny,nz,iproj,istag, &
                   tlat1,tlat2
      read(i2dmet)
!     read(icrmet)
      read(ilumet)
      read(i3dmet)
      read(i2dmet)
!     read(icrmet)
      read(ilumet)
      read(i3dmet) 
      read(i2dmet) 
!     read(icrmet) 
      read(ilumet) 

      if (itimez.ne.itzon) then
        write(*,*) 'Time zone mismatch!'
        write(*,*) 'User Time zone : ',itimez
        write(*,*) 'Met Time zone  : ',itzon
        stop
      endif
      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 new output emissions file header

      write(*,1005) jdate,0.,jdate,24.
 1005 format(5x,2(i10,f10.2))

      write(iout) ifile,fileiden,1,nspecs,jdate,0.,jdate,24.
      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(rmet2d(nx,ny))
      allocate(rmet3d(nx,ny,nz))
      allocate(salt_water(nx,ny))
      allocate(water(nx,ny))
      allocate(ts(nx,ny))
      allocate(pp(nx,ny,nz))
      allocate(zz(nx,ny,nz))
      allocate(tt(nx,ny,nz))
      allocate(uu(nx,ny,nz))
      allocate(vv(nx,ny,nz))
      allocate(spd10m(nx,ny))
      allocate(sea_ice(nx,ny))
!     allocate(fcloud(nx,ny))
!     allocate(cldtrns(nx,ny))
      allocate(coastline_len(nx,ny))
      allocate(srfzone_width(nx,ny))
      allocate(rlat(nx,ny))
      allocate(rlon(nx,ny))
      allocate(flux(nx,ny,MXSPEC))
      allocate(dflux0(nx,ny,MAXBINS))
      allocate(dflux1(nx,ny,MAXBINS))
      allocate(dflux_srf(nx,ny,MAXBINS))
      allocate(mass(nx,ny))
      allocate(mass_crs_fn(nx,ny,MXSPEC))
      allocate(rchla_conc(nx,ny))
      allocate(col_loc(nx,ny))
      allocate(row_loc(nx,ny))
      allocate(area2d(nx,ny))

!-----Determine number of bins to use for specified size range

      if (lss) then
        do r = 1,num_bins
          bin_min = binctr_diam(r) - 0.5*binwdt_diam(r)
          bin_max = binctr_diam(r) + 0.5*binwdt_diam(r)
          if (lower_limit.ge.bin_min .and. lower_limit.le.bin_max) then
            lowerbin = r
            if (lowerbin .ne. 1) then
              binwdt_diam(lowerbin) = bin_max - lower_limit
              binctr_diam(lowerbin) = lower_limit + (bin_max - lower_limit)/2.
            endif 
          endif
          if (r .eq. 1 .and. lower_limit .lt. bin_min) then
            lowerbin = 1
            binwdt_diam(lowerbin) = bin_max - lower_limit
            binctr_diam(lowerbin) = lower_limit + (bin_max - lower_limit)/2.
          endif
          if (upper_limit.ge.bin_min .and. upper_limit.le.bin_max) then
            upperbin = r
            if (upperbin .ne. num_bins) then
              binwdt_diam(upperbin) = upper_limit - bin_min
              binctr_diam(upperbin) = bin_min + (upper_limit - bin_min)/2.
            endif
          endif
          if (r .eq. num_bins .and. upper_limit .gt. bin_max) then
            upperbin = num_bins
            binwdt_diam(upperbin) = upper_limit - bin_min
            binctr_diam(upperbin) = bin_min + (upper_limit - bin_min)/2.
          endif
        enddo

!----Convert dry diameters to wet radii at RH=80% (Rwet = 2*Rdry = Ddry)

        do r = lowerbin,upperbin
          bin_ctr(r)   = binctr_diam(r) 
          bin_width(r) = binwdt_diam(r)
        enddo
      endif

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

      iww = 0
      read(ilumet,end=900)
      do 10 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
          do j = 1,ny
            do i = 1,nx
              water(i,j) = rmet2d(i,j)
            enddo
          enddo
        endif
 10   continue
      if (iww.eq.0) goto 900
      close(ilumet)

      do j = ny,1,-1
        read(ifswm,'(500(:,i1))') (salt_water(i,j),i=1,nx)
      enddo
      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. 
          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).ne.0) then
            if (iproj.ne.0) then
              area = water(i,j)*dx*dy
              clen = 0.5*(dx + dy)
            else
              jlat = yorg + (dy*j)-(dy/2.)
              area = abs(dx*111134*cos(jlat/conv)*dy*111134)
              area = water(i,j)*area
              clen = 0.5*(abs(dx*111134*cos(jlat/conv)) + dy*111134)
            endif
            area2d(i,j) = area
          endif
          if (salt_water(i,j).gt.1) coastline_len(i,j) = clen !meters
        enddo
      enddo
      if (trim(icstpath) .ne. 'NONE') then
        open (unit=icst,file=icstpath,form='formatted')
        do i = 1,nx*ny 
          read(icst,*) inx, iny, coastline_len(inx,iny)
        enddo
        close(icst)
      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,plat,plon,tlat1,tlat2,xloc,yloc,rlon(i,j),rlat(i,j))
            elseif (iproj.eq.3) then
              call pspgeo(1,plon,plat,xloc,yloc,rlon(i,j),rlat(i,j))
            endif
          enddo
        enddo
      endif

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

      if (lnhalo) 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), lon_chl(nxchl,nychl), &
                 lat_chl(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

        ! Define LAT/LON of global grid cells

        dxchl = 360./nxchl
        dychl = 180./nychl
        do j = 1, nychl
          lon_chl(1,j) = -180.0 + dxchl/2.
        enddo
        do i = 1, nxchl
          lat_chl(i,1) = -90.0 + dychl/2.
        enddo
        do j = 1, nychl-1
          do i = 1, nxchl-1
            lat_chl(i,j+1) = lat_chl(i,j) + dychl
            lon_chl(i+1,j) = lon_chl(i,j) + dxchl
          enddo
        enddo
        do i = 1, nxchl-1
          lon_chl(i+1,nychl) = lon_chl(i,nychl) + dxchl
        enddo
        do j = 1, nychl-1
          lat_chl(nxchl,j+1) = lon_chl(nxchl,j) + dxchl
        enddo

1000    format('Date ', i2.2, '/', i2.2, '/', i2.2 '; time ', &
               i2.2, ':', i2.2, ':', i2.2)

        ! Mapping

!       write(*,*) 'Scanning/Mapping south/north border lines'
!       call idate(today) ! today(1)=month, (2)=day, (3)=year
!       call itime(now)   ! now(1)=hour,(2)= minute, (3)=second
!       write(*,1000) today(1),today(2),today(3),now
        do col = 1, nx ! Scan south and north border lines simultaneously
          lat1 = rlat(col,1)  !lat of south border line cells
          lon1 = rlon(col,1)  !lon of south border line cells
          lat2 = rlat(col,ny) !lat of north border line cells
          lon2 = rlon(col,ny) !lon of north border line cells
          dmin1 = 1.0e+30
          dmin2 = 1.0e+30
          maxj = -999
          minj = 99999
          do i = 1, nxchl
            do j = 1, nychl
              x1 = (lat1-lat_chl(i,j))**2
              y1 = (lon1-lon_chl(i,j))**2
              x2 = (lat2-lat_chl(i,j))**2
              y2 = (lon2-lon_chl(i,j))**2
              if ((x1+y1).lt.dmin1) then
                dmin1 = x1 + y1
                col_loc(col,1) = i
                row_loc(col,1) = j
              endif
              if ((x2+y2).lt.dmin2) then
                dmin2 = x2 + y2
                col_loc(col,ny) = i
                row_loc(col,ny) = j
              endif
            enddo
          enddo
          if (row_loc(col,1).lt.minj) minj = row_loc(col,1)
          if (row_loc(col,ny).gt.maxj) maxj = row_loc(col,ny)
        enddo

!       write(*,*) 'Scanning/Mapping west/east border lines'
!       call idate(today) ! today(1)=month, (2)=day, (3)=year
!       call itime(now)   ! now(1)=hour,(2)= minute, (3)=second
!       write(*,1000) today(1),today(2),today(3),now
        do row = 1, ny ! Scan west and east border lines simultaneously
          lat1 = rlat(1,row)  !lat of west border line cells
          lon1 = rlon(1,row)  !lon of west border line cells
          lat2 = rlat(nx,row) !lat of east border line cells
          lon2 = rlon(nx,row) !lon of east border line cells
          dmin1 = 1.0e+30
          dmin2 = 1.0e+30
          maxi = -999
          mini = 99999
          do i = 1, nxchl
            do j = minj, maxj
              x1 = (lat1-lat_chl(i,j))**2
              y1 = (lon1-lon_chl(i,j))**2
              x2 = (lat2-lat_chl(i,j))**2
              y2 = (lon2-lon_chl(i,j))**2
              if ((x1+y1).lt.dmin1) then
                dmin1 = x1 + y1
                col_loc(1,row) = i
                row_loc(1,row) = j
              endif
              if ((x2+y2).lt.dmin2) then
                dmin2 = x2 + y2
                col_loc(nx,row) = i
                row_loc(nx,row) = j
              endif
            enddo
          enddo
          if (col_loc(1,row).lt.mini) mini = col_loc(1,row)
          if (col_loc(nx,row).gt.maxi) maxi = col_loc(nx,row)
        enddo
!       write(*,*)'mini,maxi,minj,maxj'
!       write(*,*)mini,maxi,minj,maxj

!       write(*,*) 'Mapping is processing'
!       call idate(today) ! today(1)=month, (2)=day, (3)=year
!       call itime(now)   ! now(1)=hour,(2)= minute, (3)=second
!       write(*,1000) today(1),today(2),today(3),now
        do col = 2, nx-1
          do row = 2, ny-1
            lat1 = rlat(col,row)
            lon1 = rlon(col,row)
            dmin = 1.0e+30
            if (mini <= maxi) then
              do j = minj, maxj
                do i = mini, maxi
                    x1 = (lat1-lat_chl(i,j))**2
                    y1 = (lon1-lon_chl(i,j))**2
                    if ((x1+y1).lt.dmin) then
                      dmin = x1 + y1
                      col_loc(col,row) = i
                      row_loc(col,row) = j
                    endif
                enddo
              enddo
            else ! The regional domain crosses the International Date Line
              do j = minj, maxj
                do i = mini, nxchl
                    x1 = (lat1-lat_chl(i,j))**2
                    y1 = (lon1-lon_chl(i,j))**2
                    if ((x1+y1).lt.dmin) then
                      dmin = x1 + y1
                      col_loc(col,row) = i
                      row_loc(col,row) = j
                    endif
                enddo
                do i = 1, maxi
                    x1 = (lat1-lat_chl(i,j))**2
                    y1 = (lon1-lon_chl(i,j))**2
                    if ((x1+y1).lt.dmin) then
                      dmin = x1 + y1
                      col_loc(col,row) = i
                      row_loc(col,row) = j
                    endif
                enddo
              enddo
            endif
          enddo
        enddo
!       call idate(today) ! today(1)=month, (2)=day, (3)=year
!       call itime(now)   ! now(1)=hour,(2)= minute, (3)=second
!       write(*,1000) today(1),today(2),today(3),now

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

        !Initialize

        do j = 1,ny
          do i = 1,nx
            rchla_conc(i,j) = 0.
          enddo
        enddo

        !Mapping Global chlorophyll-A conc to regional grid

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

!-----Loop over hours in a day

      do ihr = 1,24

!-----Initialize...

        do j = 1,ny
          do i = 1,nx
            sea_ice(i,j) = 0.
            mass(i,j) = 0.
            do r = 1,MAXBINS
              dflux0(i,j,r) = 0.            
              dflux1(i,j,r) = 0.
              dflux_srf(i,j,r) = 0.
            enddo
            do l = 1,MXSPEC
              mass_crs_fn(i,j,l) = 0.
              flux(i,j,l) = 0.
            enddo
          enddo
        enddo

!-----Read the met

        its = 0
        read(i2dmet,end=901) idt,hr
        do 20 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
            do j = 1,ny
              do i = 1,nx
                ts(i,j) = rmet2d(i,j)
              enddo
            enddo
          endif
 20     continue
        if (its.eq.0) goto 901

        izz = 0
        ipp = 0
        itt = 0
        iuu = 0
        ivv = 0
        read(i3dmet,end=902) idt,hr
        do 30 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.'ZGRID_M') then
            izz = 1
            do k = 1,nz
              do j = 1,ny
                do i = 1,nx
                  zz(i,j,k) = rmet3d(i,j,k)
                enddo
              enddo
            enddo
          endif
          if (metspc.eq.'PRESS_MB') then
            ipp = 1
            do k = 1,nz
              do j = 1,ny
                do i = 1,nx
                  pp(i,j,k) = rmet3d(i,j,k)
                enddo
              enddo
            enddo
          endif
          if (metspc.eq.'TEMP_K') then
            itt = 1
            do k = 1,nz
              do j = 1,ny
                do i = 1,nx
                  tt(i,j,k) = rmet3d(i,j,k)
                enddo
              enddo
            enddo
          endif
          if (metspc.eq.'UWIND_MpS') then
            iuu = 1
            do k = 1,nz
              do j = 1,ny
                do i = 1,nx
                  uu(i,j,k) = rmet3d(i,j,k)
                enddo
              enddo
            enddo
          endif
          if (metspc.eq.'VWIND_MpS') then
            ivv = 1
            do k = 1,nz
              do j = 1,ny
                do i = 1,nx
                  vv(i,j,k) = rmet3d(i,j,k)
                enddo
              enddo
            enddo
          endif
 30     continue
        if (izz.eq.0) goto 902
        if (ipp.eq.0) goto 903
        if (itt.eq.0) goto 904
        if (iuu.eq.0) goto 905
        if (ivv.eq.0) goto 906

!       icd = 0
!       read(icrmet,end=907) idt,hr
!       do 40 n = 1,ncrmet
!         do k = 1,nz
!           read(icrmet,err=907,end=907) &
!             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.'CLOUDOD') then
!           icd = 1
!           do j = 1,ny
!             do i = 1,nx
!               tau = rmet3d(i,j,nz)
!               do k = nz-1,1,-1
!                 tau = tau + rmet3d(i,j,k)
!               enddo
!               if (tau.lt.5.) then
!                 cldtrns(i,j) = 1.
!                 fcloud(i,j) = 0.
!               else
!                 cldtrns(i,j) = (5. - exp(-tau))/(4. + 0.42*tau)
!                 fcloud(i,j) = 1.
!               endif
!             enddo
!           enddo
!         endif
!40     continue
!       if (icd.eq.0) goto 907

        write(*,*)'Read met data at: ',hr,idt

!-----Calculate the 10 meter wind speed (m/s) 

        if (lss) then
          do j = 1,ny
            do 100 i = 1,nx
              spd10m(i,j) = minspd
              if (salt_water(i,j).eq.0) goto 100
              deltaz = zz(i,j,1)/2.
              press0 = pp(i,j,1) - &
                       2.*deltaz*(pp(i,j,2) - pp(i,j,1))/zz(i,j,2)
              if (istag.eq.0 .or. i.eq.1 .or. j.eq.1) then
                wind = amax1(0.1,sqrt(uu(i,j,1)**2 + vv(i,j,1)**2))
              else
                uwnd = (uu(i,j,1) + uu(i-1,j,1))/2.
                vwnd = (vv(i,j,1) + uu(i,j-1,1))/2.
                wind = amax1(0.1,sqrt(uwnd**2 + vwnd**2))
              endif
              theta  = tt(i,j,1)*(1000./pp(i,j,1))**gamma
              theta0 = ts(i,j)*(1000./press0)**gamma
              dtheta = theta - theta0
              thetabar = (theta + theta0)/2.
              ri = (grav/thetabar)*deltaz*dtheta/wind**2
 
              z0 = amax1(1.e-4,2.e-6*wind**2.5)
              zscale = vk/alog(deltaz/z0)
              if (ri.lt.0.) then
                cm    = 69.56*sqrt(deltaz/z0)*zscale**2
                fm    = 1. - 9.4*ri/(1. + cm*sqrt(abs(ri)))
              else
                fm = 1./((1. + 4.7*ri)**2)
              endif

              ustar2 = fm*(wind*zscale)**2
              ustar2 = amax1(1.e-10,ustar2)
              ustar = sqrt(ustar2)
              wr = ustar/(vk*sqrt(fm))*alog(10./z0)
              spd10m(i,j) = max(minspd,wr)
 100        continue
          enddo

!-----Calculate d(flux)/d(radius) from bubbles bursting for
!     dry radii <= 4 microns (Rwet = Ddry <= 8 microns)

          do j = 1,ny
            do i = 1,nx
              do 200 r = lowerbin,upperbin
                dflux_srf(i,j,r) = 0.
                dflux1(i,j,r)    = 0.
                dflux0(i,j,r)    = 0.
                if (salt_water(i,j).eq.0) goto 200

                radius(r) = bin_ctr(r)
                if (radius(r) .le. 8.) then
                  eval_r = bin_ctr(r) 
                  dk = bin_width(r)

!----Original Monahan et al. (1986) parameterization

                  if (trim(method) .eq. 'MONAHAN') then
                    B = (0.380 - log10(eval_r))/0.650
                    term1 = 1.373*(spd10m(i,j)**3.41)*(eval_r**(-3.))
                    term2 = 1. + 0.057*eval_r**1.05
                    term3 = 1.19*exp(-(B**2))
                    term4 = 10.**(term3)

!----Gong et al. (2003) parameterization.  Reduces emissions at submicron
!    sizes relative to Monahan.

                  elseif (trim(method) .eq. 'GONG') then
                    theta_g = 30.
                    B       = (0.433 - log10(eval_r))/0.433
                    term_0  = -0.017*(eval_r**(-1.44))
                    A_g     = 4.7*(1.0 + theta_g*eval_r)**term_0
                    term1   = 1.373*(spd10m(i,j)**3.41)*(eval_r**(-A_g))
                    term2   = 1. + 0.057*(eval_r**3.45)
                    term3   = 1.607*exp(-(B**2))
                    term4   = 10.**(term3)
                  else
                    write(*,'(2a)')'Unrecognized open ocean ', &
                                   'parameterization--stopping.'
                    stop
                  endif
                  dflux0(i,j,r) =  term1*term2*term4*dk

!----Calculate d(flux)/d(radius) as in Smith and Harrison (1998) from spume 
!    and bubbles for radius > 4 microns.   

                elseif (radius(r) .gt. 8.) then  
                  dk = bin_width(r)
                  eval_r = bin_ctr(r)
                  R01 = 3.0
                  R02 = 30. 
                  f1  = 1.5
                  f2  = 1.
                  A1  = 0.2*(spd10m(i,j)**3.5)
                  A2  = 6.8*(10.**( -3))*(spd10m(i,j)**3) 
  
                  term1 = alog((eval_r/R01)**2)
                  term2 = exp(-f1*term1)
                  dflux0(i,j,r) = A1*term2*dk
  
                  term1 = alog((eval_r/R02)**2)
                  term2 = exp(-f2*term1)
                  dflux1(i,j,r) = A2*term2*dk
                endif

!----DeLeeuw Surf Zone parameterization
!    Evaluating results of integral of DeLeeuw (2000) equation (4)
!    See for example, Lewis and Schwartz, (2004), p. 222

                if (salt_water(i,j).gt.1) then
                  if (trim(method_s) .eq. 'DELEEUW') then
                    term1 = (1.1e7/(-0.65))*exp(0.23*min(9.0,spd10m(i,j)))
                    srf_lower_limit = (bin_ctr(r) - (bin_width(r)/2.))*2.
                    srf_upper_limit = (bin_ctr(r) + (bin_width(r)/2.))*2.
                    if (srf_lower_limit .ge. 1.6 .and. &
                        srf_upper_limit .le. 20.) then
                      term2 = (srf_upper_limit)**(-0.65)
                      term3 = (srf_lower_limit)**(-0.65)
                      dflux_srf(i,j,r) = term1*(term2 - term3)
                    else
                      dflux_srf(i,j,r) = 0.
                    endif
                  elseif (trim(method_s) .eq. 'GONG') then
                    dk = bin_width(r)
                    eval_r = bin_ctr(r)
                    theta_g = 30.
                    B       = (0.433 - log10(eval_r))/0.433
                    term_0  = -0.017*(eval_r**(-1.44))
                    A_g     = 4.7*(1.0 + theta_g*eval_r)**term_0
  
!-----Whitecap coverage fraction = 1 

                    term1   = 3.576e5*(eval_r**(-A_g)) 
                    term2   = 1. + 0.057*((eval_r)**3.45)
                    term3   = 1.607*exp(-(B**2))
                    term4   = 10.**(term3)
                    dflux_srf(i,j,r) = term1*term2*term4*dk
                  else 
                    write(*,'(2a)') 'Unrecognized surf zone ', &
                                    'parameteriziation--stopping.'
                    stop
                  endif
                endif
 200          continue           
            enddo
          enddo 

!-----Add contributions from direct and indirect terms to get total aerosol
!     production. Units of the fluxes are particles / m^2-s.  
!    
!     Need to convert from number of particles to mass/hr going  into the
!     gridcell.  The number of particles produced is converted to a mass of
!     particles using:  Mass = Number x (particle density) x (4/3)pi radius^3
!     This gives a mass flux (kg/m^2-s) so still need to multiply by the
!     area of the grid box to get the mass per hour.

          do j = 1,ny
            do 300 i = 1,nx
              if (salt_water(i,j).eq.0) goto 300
              mass(i,j) = 0.
              do r = lowerbin,upperbin
                aersl_flux = dflux0(i,j,r) + dflux1(i,j,r) 
                dry_mass = aersl_flux*dens_NaCl* &
                           (4./3.)*PI*((radius(r)/2.)*1.e-6)**3
                dry_mass = dry_mass*area2d(i,j)*1000.*3600.
   
                dry_srf_mass = dflux_srf(i,j,r)*dens_NaCl* &
                               (4./3.)*PI*((radius(r)/2.)*1.e-6)**3
                area = srfzone_width(i,j)*coastline_len(i,j)
                dry_srf_mass = dry_srf_mass*area*1000.*3600.
   
   
                mass(i,j) = mass(i,j) + &
                            (1. - sea_ice(i,j))*(dry_srf_mass + dry_mass)
              enddo 
 300        continue
          enddo 
        endif

!-----Calculate halo-methane emissions
!     (Only occurs for sea temperature > -2C)
 
        if (lnhalo) then
          do j = 1,ny
            do 400 i = 1,nx
              if (salt_water(i,j) .eq. 0) goto 400
              if (sea_ice(i,j).ne.0 .or. ts(i,j).lt.271.) goto 400

!-----Solar flux

!             hhmm = ihr*100.
!             call getznth(rlat(i,j),rlon(i,j),hhmm,jdate,itimez,zenith,ldark)
!             coszen = cos(DBLE(zenith)*DBLE(PI)/180.)
!             solflux = (990.*coszen - 30.)* &
!                       (1. - fcloud(i,j)*(1. - cldtrns(i,j)))
!             solflux = dmax1(DBLE(0.),DBLE(solflux))

!-----Halo-methane emission process

              ilochr = ihr - int((centlon+7.5-rlon(i,j))/15.)
              if (ilochr .le.0) ilochr = 24 + ilochr
              flux(i,j,1) = 4.31e-8 * f(iCH3I) * area2d(i,j) * fhrly(ilochr) * &
                            rchla_conc(i,j) ! mol/hr-cell
              flux(i,j,2) = 4.31e-8 * f(iMB3) * area2d(i,j) * fhrly(ilochr) * &
                            rchla_conc(i,j) ! mol/hr-cell
              flux(i,j,3) = 4.31e-8 * f(iMB2) * area2d(i,j) * fhrly(ilochr) * &
                            rchla_conc(i,j) ! mol/hr-cell
              flux(i,j,4) = 4.31e-8 * f(iMBC) * area2d(i,j) * fhrly(ilochr) * &
                            rchla_conc(i,j) ! mol/hr-cell
              flux(i,j,5) = 4.31e-8 * f(iMB2C) * area2d(i,j) * fhrly(ilochr) * &
                            rchla_conc(i,j) ! mol/hr-cell
              flux(i,j,6) = 4.31e-8 * f(iMBC2) * area2d(i,j) * fhrly(ilochr) * &
                            rchla_conc(i,j) ! mol/hr-cell
              flux(i,j,7) = 4.31e-8 * f(iMIC) * area2d(i,j) * fhrly(ilochr) * &
                            rchla_conc(i,j) ! mol/hr-cell
              flux(i,j,8) = 4.31e-8 * f(iMIB) * area2d(i,j) * fhrly(ilochr) * &
                            rchla_conc(i,j) ! mol/hr-cell
              flux(i,j,9) = 4.31e-8 * f(iMI2) * area2d(i,j) * fhrly(ilochr) * &
                            rchla_conc(i,j) ! mol/hr-cell
              flux(i,j,10) = 2.39e-8 * area2d(i,j) ! mol/hr-cell from
                                                   ! 4.E8 molecules/cm2/s
              flux(i,j,11) = rchla_conc(i,j) ! mg/m3

! old CH3I mechanism
!             tempk = ts(i,j)
!             tempc = tempk - 273.15
!             henry = 1.38e-3*rgas*tempk* &
!                     exp(-14.43**((tempk - 298.)/tempk))
!             schmidt = 1.095*(2004. - 93.5*tempc + 1.39*tempc**2)
!             kw = 0.31*spd10m(i,j)**2*sqrt(660./schmidt)
!             prod = alpha*(solflux*doc)
!             kcl = 2.80e17*exp(-13518./tempk)
!             ch3i_aq = (prod + (kw/depth)*henry*ch3i_g)/ &
!                       ((kw/depth) + kcl*Mcl)
!             flux(i,j,1) = 0.24*kw*ch3i_aq*water(i,j)  ! ug/m2-day
!             flux(i,j,1) = (4.1e-8*dx*dy/Mw)*flux(i,j,1) ! mol/hr-cell
!             flux(i,j,3) = solflux                       ! W m-2
 400        continue
          enddo
        endif

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

        do j = 1,ny
          do i = 1,nx
            do l = 1,nspecs
              do m = 1,10
                c10(m:m) = mspec(m,l)(1:1) 
              enddo
              if (trim(c10) .eq. 'CH3I') then
                mass_crs_fn(i,j,l) = flux(i,j,1)
              endif
              if (trim(c10) .eq. 'MB3') then
                mass_crs_fn(i,j,l) = flux(i,j,2)
              endif
              if (trim(c10) .eq. 'MB2') then
                mass_crs_fn(i,j,l) = flux(i,j,3)
              endif
              if (trim(c10) .eq. 'MBC') then
                mass_crs_fn(i,j,l) = flux(i,j,4)
              endif
              if (trim(c10) .eq. 'MB2C') then
                mass_crs_fn(i,j,l) = flux(i,j,5)
              endif
              if (trim(c10) .eq. 'MBC2') then
                mass_crs_fn(i,j,l) = flux(i,j,6)
              endif
              if (trim(c10) .eq. 'MIC') then
                mass_crs_fn(i,j,l) = flux(i,j,7)
              endif
              if (trim(c10) .eq. 'MIB') then
                mass_crs_fn(i,j,l) = flux(i,j,8)
              endif
              if (trim(c10) .eq. 'MI2') then
                mass_crs_fn(i,j,l) = flux(i,j,9)
              endif
              if (trim(c10) .eq. 'I2') then
                mass_crs_fn(i,j,l) = flux(i,j,10)
              endif
              if (trim(c10) .eq. 'CHLA') then
                mass_crs_fn(i,j,l) = flux(i,j,11)
              endif
              if (trim(c10) .eq. 'NA') then
                mass_crs_fn(i,j,l) = mass(i,j)*SPCFAC_NA
              endif
              if (trim(c10) .eq. 'PCL') then
                mass_crs_fn(i,j,l) = mass(i,j)*SPCFAC_CL
              endif
              if (trim(c10) .eq. 'PSO4') then
                mass_crs_fn(i,j,l) = mass(i,j)*SPCFAC_SO4
              endif
              if (trim(c10) .eq. 'SS') then
                mass_crs_fn(i,j,l) = mass(i,j)*SPCFAC_CL*SPCFAC_SS
              endif
              if (trim(c10) .eq. 'WIND') then
                mass_crs_fn(i,j,l) = spd10m(i,j)
              endif
!             if (trim(c10) .eq. 'SOLFLUX') then
!               mass_crs_fn(i,j,l) = flux(i,j,3)
!             endif
            enddo 
          enddo
        enddo

!-----Write out emission records

        write(*,1005) jdate,float(ihr-1),jdate,float(ihr)
        write(iout) jdate,float(ihr-1),jdate,float(ihr)
        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

      enddo     !end of the ihr loop for this day
      stop

!-----Met file error handling

 900  continue
      write(*,*)'Did not find WATER lanudse'
      stop

 901  continue
      write(*,*)'Did not find TSURF'
      stop

 902  continue
      write(*,*)'Did not find ZGRID'
      stop
 903  continue
      write(*,*)'Did not find PRESS'
      stop
 904  continue
      write(*,*)'Did not find TEMP'
      stop
 905  continue
      write(*,*)'Did not find UWIND'
      stop
 906  continue
      write(*,*)'Did not find VWIND'
      stop

 907  continue
      write(*,*)'Did not find CLOUDOD'
      stop

      end
