      subroutine _name_(z,wavenm,ncomp, acomp,bcomp,atotal,btotal)
c
c     user template on file h41\template\ab.txt
c
c     called by INISHAMP, BEAMc1, BEAMcz, and RHOTAU via ABSCAT
c
c     calls pureh2o, astar, bstar, and abmodel
c
c     This template shows the form required for subroutines that 
c     return the absorption (a) and scattering (b) coefficients given the
c     depth and wavelength (the "ab" routines).
c
c     The subroutines on files h41\maincode\abconst.f, abcase1.f,
c     abac9.f, etc are specific examples of how to implement this 
c     template.
c
c     The ab routine is called by Hydrolight core routines to obtain
c     the absorption and scattering coefficients.  The ab routine is 
c     incorporated into the Hydrolight core routines via include statements 
c     of the form
c        INCLUDE 'ABSCAT.INC'
c     The include file ABSCAT.INC is created by the front-end program
c     when the user selects the "ab" subroutine to be used in the run.
c
c     The required format of the subroutine is
c
c     subroutine _name_(z,wavenm,ncomp, acomp,bcomp,atotal,btotal)
c
c     where
c
c     _name_ is replaced by the desired name of the routine, e.g. myab
c
c     Input:
c
c     z = the depth in meters
c     wavenm = the wavelength in nanometers
c     ncomp = the number of IOP components expected by Hydrolight
c
c     Output:
c
c     acomp(i) = the absorption coefficient in 1/m of the ith component, 
c                i = 1 to ncomp
c     bcomp(i) = the scattering coefficient in 1/m of the ith component, 
c                i = 1 to ncomp
c     atotal = the total absorption coefficient in 1/m (= sum of the acomp(i)
c              values)
c     btotal = the total scattering coefficient in 1/m (= sum of the bcomp(i)
c              values)
c
c     A component's specific absorption and scattering can be calculated via
c     a call to the astar and bstar routines.  Both astar and bstar are
c     functions which can be called in the form:
c     	astarchl = astar(icomp, wavenm)
c     The GUI specified options/params for astar/bstar are passed
c     via a common block from initial (where they are read) directly to
c     the function. 
c
c     The a* specifications are:
c         iastropt = -1:  Subroutine called for aCDOM (a* routine not called)
c                    -2:  aCDOM solved as a percentage of aw+achl in "ab" routine (a* not called)
c		        0:  user-supplied data file read to get a* values
c                     1:  Pope & Fry absorption model used (pure water only)
c                     2:  Smith and Baker absorption model used (pure water only)
c                     3:  Prieur-Sathyenranath-Morel model used (Chlorophyll only)
c                     4:  Exponential model used (CDOM only)
c                     5:  Prieur-Sathyendranath-Morel model used (CDOM only)
c     The b* specifications are:
c         ibstropt =  0:  user-supplied data file read to get b* values
c                     1:  Power Law used
c                     2:  Linear model
c                     3:  Gordon-Morel power-law model used (Chlorophyll only)
c                     4:  Gould-Arnone-Morel linear model used (Chlorophyll only)
cc
c**** WARNING:  In this routine, the user can define the components to be 
c     whatever is desired, e.g. component 1 is pure water, component 2 is CDOM, 
c     component 3 is particles, etc.  IT IS THE USER'S RESPONSIBILITY TO MAKE 
c     SURE THAT AN APPROPRIATE PHASE FUNCTION IS USED FOR EACH COMPONENT
c     WHEN ENTERING THE NAMES OF THE FILES OF DISCRETIZED PHASE FUNCTIONS
c     IN THE FRONT-END ROUTINE.
c
c**** Note:  after creating a new routine from this template,
c        1.  save the file with a .f extension, e.g. myab.f
c        2.  place the new file in the maincode directory
c        3.  add the new filename to the makefile in the maincode directory
c
      INCLUDE 'DIMENS.INC'
c
      real compconc(mxcomp), astarcomp(mxcomp), bstarcomp(mxcomp)
      real astar1(mxcomp), bstar1(mxcomp)

      dimension acomp(mxcomp),bcomp(mxcomp)
      common/cmisc/ imisc(30),fmisc(30)
c
c     Specific ab variables: only iastropt and ibstropt used here!
      Common /Castar/ iastropt, astarRef, astar0, asgamma
      Common /Cbstar/ ibstropt, bstarRef, bstar0, PLn, GAMm, GAMi
      Integer iastropt(mxcomp+2), ibstropt(mxcomp+2)
      Real astarRef(mxcomp+2), astar0(mxcomp+2)
      Real asgamma(mxcomp+2)
      Real bstarRef(mxcomp+2), bstar0(mxcomp+2)
      Real PLn(mxcomp+2), GAMm(mxcomp+2),GAMi(mxcomp+2)

c	itype indicates the method used to specify the component concentrations
c     = 0 if constant; =1 if subroutine, =2 if datafile
	common /Cconc/  itype(mxcomp)  
        integer itype
c
      data kall/0/, waveold/0.0/
c
c     The number of components in the model may be specified either here
c     or below in the code
c      data numcomp/4/
c
c     Note on dynamic storage:  MUST save data values between subroutine 
c     calls:
      save

c
      if(kall.eq.0) then
c
c******* Note:  Be sure to define the number of components used in the model
c        via a statment or in the data statement above
c           numcomp = ...
c
c
c        Perform a check to see if the number of components used in the
c        ab model corresponds with what is expected by Hydrolight:
c
         if(numcomp.ne.ncomp) then
            write(6,100) ncomp, numcomp
  100       format(//' Error in ab subroutine:  ncomp = ',i2,
     1      ' must equal ',i2,' for this ab routine')
            stop
         endif
c
c        Initialization on the first call:
c
c        In this section, perform any initialization calculations that
c        need to be done only once, such as computing the parameters of
c        a model or reading in a data file of a and b values.
c
c------- Add user initialization code here:  ---------------------------
c
c        call ABMODEL first time to trigger messages
         depth = 0.0
	  call abmodel(depth, compconc)
c
c        *** You may need to normalize astars and bstars
c            here depending on what models you are using!
c        *** You may need to normalize astars and bstars
c            here depending on what models you are using!
         Do ii=2,ncomp
           	conc = compconc(ii)
	   If(iastropt(ii).eq.3) then	!use Prieur-Sathyendranath-Morel model
         	astar1(ii) = astar(ii, astarRef(ii))
              If(abs(astar1(ii)).lt.1.e-7) then
                write(10,'(1x,a,i2,a,f4.0,a)') 'Astar for component ',
     1                      ii,'is zero at reference wavelength ',
     2                      astarRef(ii),' nm.'
                write(10,'(5x,a)') '(used to normalize absorption)'
     1                      'for component ',ii,' is zero.'
                write(10,*) 'Absorption for component ',ii,' will be ',
     1                      'set to zero for all depths and wavelengths!'
                write(10,*)
              Endif
          Endif
c         ** no model currently uses a normalized bstar
c        	bstar1(ii) = bstar(ii, bstarRef(ii))
         Enddoc
c------- End user initialization code ----------------------------------
c
         kall = 1
      endif
c
c     End of initialization ----------------------------------------------
c
c     Subsequent calls to the ab routine execute only the following code.
c     Because the ab routine is called many times when solving the RTE
c     (in particular, when integrating the Riccati equations), use care
c     to make this code as efficient as possible.  
c
c---- Add user code here: ------------------------------------------------
c
c       Save the routine input into dummy arrays
      depth = z
      wavelen = wavenm
c
c****** Note:  For multiwavelength ab models, it can be efficient to check
c       to see if the call is at the same wavelength as the previous call
c       or at a new one, and then perform any calculations that are needed
c       only once for a given wavelength.  This can be done via statements
c       of the form
c
      if(wavenm.gt.waveold) then
c     	 new-wavelength calculations here....
c
c	 call routine pureH2O to get the pure water values
c       for use with some a and b models
        call pureH2O(wavelen, awater,bwater,cwater)
c
c****** Note:  if pure water is a component of your model, you
c	 may need to treat it separately (see abcase2 for example)
c
c       call astar and bstar routines for each component
c       using the form:  Xstar(icomp, wavenm)
        Do ii=1,ncomp
           astarcomp(ii) = astar(ii, wavelen)
           bstarcomp(ii) = bstar(ii, wavelen)
        Enddo
        waveold = wavenm
      endif
c
c******* Then perform the remainder of the a and b calculations for the current
c        depth and wavelength....
c
c        A call to abmodel will provide the concentration info for each
c        component:  compconc(1:ncomp)
	  call abmodel(depth, compconc)
c
c        If needed, these concentration values can be overwritten
c        here...

c        Call the astar and bstar routines to define the specific
c        scattering and absorption coefficients
c        Then use the concentrations to get the total a and b for
c	  each component
c
      Do ii=1,ncomp
c
c     **** if concentration is zero (or less than), set a and b to zero
       If(compconc(ii).le.0) then
              acomp(ii) = 0.0
              bcomp(ii) = 0.0
       Else
c     ----- calculate absorption ----
	   If(iastropt(ii).eq.3) then	!use Prieur-Sathyendranath-Morel model
            If(abs(astar1(ii)).ge.1e-7) then   !make sure astar1 is nonzero
c             Get the Chl absorption, by the Morel model (L&W eq 3.27).
        	acomp(ii) = 0.06 * astarcomp(ii)/astar1(ii)
     1                         * compconc(ii)**(0.65)
            Else
              acomp(ii) = 0.0
            Endif
      	   Else	!for all other options, use astarchl
        	acomp(ii) = astarcomp(ii) * compconc(ii)
	   Endif
c
c     ----- calculate scattering ----
          If(ibstropt(ii).le.0) then    !data file or nonscattering
c             Use the general form, b = conc * bstar
            	bcomp(ii) = compconc(ii) * bstarcomp(ii)
          Else
c             Use the general form, b = conc**N * bstar
            	bcomp(ii) = compconc(ii)**compN(ii) * bstarcomp(ii)
          Endif

       Endif

      Enddo
cc
c
c******* WARNING:  This section of the code MUST define the a and b value
c        of each component.  Thus the code must include statements of the
c        form (or something equivalent)
c
c        acomp(1) = ...
c          ...
c        acomp(ncomp) = ...
c        bcomp(1) = ...
c          ...
c        bcomp(ncomp) = ...
c
c---- End user code section ---------------------------------------------
c
c     Now compute the total a and b values from the component values:
      atotal = 0.0
      btotal = 0.0
      do i=1,ncomp
         atotal = atotal + acomp(i)
         btotal = btotal + bcomp(i)
      enddo
c
      return
      end


