module atmos_tracer_driver_mod
! <CONTACT EMAIL="GFDL.Climate.Model.Info@noaa.gov">
!   William Cooke
! </CONTACT>

! <REVIEWER EMAIL="GFDL.Climate.Model.Info@noaa.gov">
!   Matt Harrison
! </REVIEWER>

! <REVIEWER EMAIL="GFDL.Climate.Model.Info@noaa.gov">
!   Bruce Wyman
! </REVIEWER>


! <OVERVIEW>
!     This code allows the user to easily add tracers to the FMS framework.
! </OVERVIEW>

! <DESCRIPTION>
!
!     This code allows a user to easily implement tracer code in the FMS
!     framework.  The tracer and tracer tendency arrays are supplied along
!     with longtitude,  latitude, wind, temperature, and pressure
!     information which allows a  user to implement sources and sinks of the
!     tracer which depend on these parameters.
!
!     In the following example, radon being implemented in the atmosphere
!     will be used as an example of how to implement a tracer in the FMS
!     framework.
!
!     Within the global scope of tracer_driver_mod a use
!     statement should be inserted for each tracer to be added.
!<PRE>      use radon_mod, only : radon_sourcesink, radon_init, radon_end </PRE>
!     
!     An integer parameter, which  will be used as an identifier for the
!     tracer, should be assigned.
!<PRE>
!      integer :: nradon 
!</PRE>
!     Within tracer_driver_init a call to the tracer manager is needed in
!     order  to identify which tracer it has set the tracer as.
!<PRE>
!      nradon = get_tracer_index(MODEL_ATMOS,'radon')
!</PRE>
!     Here MODEL_ATMOS is a parameter defined in field_manager. 
!          'radon' is the name of the tracer within the field_table.
!     
!     If the tracer exists then the integer returned will be positive and it
!     can be used to call the initialization routines for the individual
!     tracers.
!<PRE>
!      if (nradon > 0) then
!           call radon_init(Argument list)
!      endif
!</PRE>     
!
!     Within tracer_driver the user can also use the identifier to surround
!     calls to the source-sink routines for the tracer of interest.
!
!<PRE>
!      if (nradon > 0 .and. nradon <= nt) then
!          call radon_sourcesink (Argument list)
!          rdt(:,:,:,nradon)=rdt(:,:,:,nradon)+rtnd(:,:,:)
!      endif
!</PRE>
!
!     It is the users responsibility to add the tendency generated by the
!     sourcesink routine.
!      
!     Within tracer_driver_end the user can add calls to the
!     terminators for the appropriate source sink routines.
!
!<PRE>      call radon_end</PRE>
!
!     This may simply be a deallocation statement or a routine to send
!     output to the logfile stating that the termination routine has been
!     called.
!
! <NOTE>
! This code has been modified by Paul.Ginoux in 2005
! to include the following aerosols:
!  - SO4 with a simplified SOx chemistry
!  - sea salt with a size distributon splitted into 5 bins
!  - mineral dust with a size distributon splitted into 5 bins
!  - carbonaceous aerosols with new emissions by Shekar.Reddy
! </NOTE>
!
! </DESCRIPTION>


!-----------------------------------------------------------------------

use mpp_mod,               only : input_nml_file 
use fms_mod,               only : file_exist, close_file,&
                                  open_namelist_file, check_nml_error, &
                                  write_version_number, &
                                  error_mesg, &
                                  FATAL, &
                                  mpp_pe, &
                                  mpp_root_pe, &
                                  stdlog, &
                                  mpp_clock_id, &
                                  mpp_clock_begin, &
                                  mpp_clock_end, &
                                  CLOCK_MODULE
use time_manager_mod,      only : time_type, &
                                  get_date, get_date_julian, &
                                  real_to_time_type
use diag_manager_mod,      only : register_diag_field, send_data
use astronomy_mod,         only : astronomy_init, diurnal_solar
use tracer_manager_mod,    only : get_tracer_index,   &
                                  get_number_tracers, &
                                  get_tracer_names,   &
                                  get_tracer_indices, &
                                  adjust_positive_def
use field_manager_mod,     only : MODEL_ATMOS
use atmos_tracer_utilities_mod, only :                     &
                                  dry_deposition,     &
                                  dry_deposition_time_vary, &
                                  dry_deposition_endts,     &
                                  atmos_tracer_utilities_init, &
                                  get_rh, get_w10m, get_cldf, &
                                  sjl_fillz
use constants_mod,         only : grav, WTMAIR
use atmos_radon_mod,       only : atmos_radon_sourcesink,   &
                                  atmos_radon_init,         &
                                  atmos_radon_end
use atmos_carbon_aerosol_mod, only : &
                                  atmos_carbon_aerosol_time_vary,  &
                                  atmos_carbon_aerosol_endts,      &
                                  atmos_carbon_aerosol_driver,  &
                                  atmos_carbon_aerosol_init,&
                                  atmos_carbon_aerosol_end
use atmos_convection_tracer_mod,only: &
                                  atmos_convection_tracer_init, &
                                  atmos_cnvct_tracer_sourcesink, &
                                  atmos_convection_tracer_end
use atmos_sulfur_hex_mod,  only : atmos_sf6_sourcesink,     &
                                  atmos_sulfur_hex_init,    &
                                  atmos_sulfur_hex_end
use atmos_ch3i_mod,        only : atmos_ch3i, &
                                  atmos_ch3i_init, &
                                  atmos_ch3i_time_vary, &
                                  atmos_ch3i_endts, &
                                  atmos_ch3i_end
use atmos_sea_salt_mod,    only : atmos_sea_salt_sourcesink,     &
                                  atmos_sea_salt_init,    &
                                  atmos_sea_salt_end
use atmos_dust_mod,        only : atmos_dust_sourcesink,     &
                                  atmos_dust_init,    &
                                  atmos_dust_time_vary,    &
                                  atmos_dust_endts,        &
                                  atmos_dust_end
use atmos_sulfate_mod,     only : atmos_sulfate_init, &
                                  atmos_sulfate_time_vary, &
                                  atmos_sulfate_endts,     &
                                  atmos_sulfate_end, &
                                  atmos_DMS_emission, &
                                  atmos_SOx_emission, &
                                  atmos_SOx_chem
use atmos_soa_mod,         only : atmos_SOA_init, &
                                  atmos_SOA_time_vary, &
                                  atmos_SOA_endts, &
                                  atmos_SOA_end, &
                                  atmos_SOA_chem
use tropchem_driver_mod,   only : tropchem_driver, &
                                  tropchem_driver_time_vary, &
                                  tropchem_driver_endts, &
                                  tropchem_driver_init
use strat_chem_driver_mod, only : strat_chem, strat_chem_driver_init
use atmos_age_tracer_mod,  only : atmos_age_tracer_init, atmos_age_tracer, &
                                  atmos_age_tracer_end
use atmos_co2_mod,         only : atmos_co2_sourcesink,   &
                                  atmos_co2_emissions,          &
                                  atmos_co2_gather_data,        &
                                  atmos_co2_flux_init,          &
                                  atmos_co2_init,               &
                                  atmos_co2_end

use interpolator_mod,      only : interpolate_type

implicit none
private
!-----------------------------------------------------------------------
!----- interfaces -------

public  atmos_tracer_driver,            &
        atmos_tracer_driver_init,       &
        atmos_tracer_driver_time_vary,       &
        atmos_tracer_driver_endts,       &
        atmos_tracer_driver_end,        &
        atmos_tracer_flux_init,         &
        atmos_tracer_driver_gather_data

!-----------------------------------------------------------------------
!----------- namelist -------------------
logical :: prevent_flux_through_ice = .false.  
                               ! when true, tracers will only be fluxed 
                               ! through the non-ice-covered portions of
                               ! ocean grid boxes
                                             

namelist /atmos_tracer_driver_nml / prevent_flux_through_ice
!-----------------------------------------------------------------------
!
!  When initializing additional tracers, the user needs to make the
!  following changes.
!
!  Add an integer variable below for each additional tracer. 
!  This should be initialized to zero. 
!
!-----------------------------------------------------------------------

! Indices for timing the various chemistry routines
integer :: radon_clock = 0
integer :: convect_clock = 0
integer :: age_tracer_clock = 0
integer :: stratozone_clock = 0
integer :: tropchem_clock = 0
integer :: carbon_clock = 0
integer :: dust_clock = 0
integer :: seasalt_clock = 0
integer :: sulfur_clock = 0
integer :: SOA_clock = 0
integer :: sf6_clock = 0
integer :: ch3i_clock = 0
integer :: co2_clock = 0

logical :: do_tropchem = .false.  ! Do tropospheric chemistry?
logical :: do_coupled_stratozone = .FALSE. !Do stratospheric chemistry?

integer, dimension(6) :: itime   ! JA's time (simpler than model time) 
integer :: nsphum  ! Specific humidity parameter

integer :: no3 = 0
integer :: no3ch = 0
integer :: nextinct = 0
integer :: naerosol = 0
integer :: nbcphobic =0
integer :: nbcphilic =0
integer :: nomphobic =0
integer :: nomphilic =0
integer :: nclay     =0
integer :: nsilt     =0
integer :: nseasalt1 =0
integer :: nseasalt2 =0
integer :: nseasalt3 =0
integer :: nseasalt4 =0
integer :: nseasalt5 =0
integer :: ndust1    =0
integer :: ndust2    =0
integer :: ndust3    =0
integer :: ndust4    =0
integer :: ndust5    =0
integer :: nsf6      =0
integer :: nDMS      =0
integer :: nSO2      =0
integer :: nSO4      =0
integer :: nMSA      =0
integer :: nSOA      =0
integer :: nH2O2     =0
integer :: nch3i     =0
integer :: nage      =0
integer :: nco2      =0
integer :: nNH4NO3   =0
integer :: nNH4      =0
integer :: nDMS_cmip =0
integer :: nSO2_cmip =0
integer :: noh       =0

real    :: ozon(11,48),cosp(14),cosphc(48),photo(132,14,11,48),   &
           solardata(1801),chlb(90,15),ozb(144,90,12),tropc(151,9),  &
           dfdage(90,48,8),anoy(90,48)


integer, dimension(:), pointer :: nradon
integer, dimension(:), pointer :: nconvect

integer :: nt     ! number of activated tracers
integer :: ntp    ! number of activated prognostic tracers

logical :: use_tau=.false.

character(len=6), parameter :: module_name = 'tracer'
character(len=7), parameter :: mod_name = 'tracers'

logical :: module_is_initialized = .FALSE.

type(interpolate_type), allocatable :: drydep_data(:)

integer, allocatable :: local_indices(:) 
! This is the array of indices for the local model. 
! local_indices(1) = 5 implies that the first local tracer is the fifth
! tracer in the tracer_manager.
  
integer :: id_landfr, id_seaicefr, id_snowfr, id_vegnfr, id_vegnlai
integer :: id_om_ddep, id_bc_ddep, id_ssalt_ddep, id_dust_ddep, &
           id_nh4_ddep_cmip
integer :: id_ssalt_emis, id_dust_emis
integer :: id_nh4no3_col, id_nh4_col
integer :: id_nh4no3_cmip, id_nh4_cmip
integer :: id_nh4no3_cmipv2, id_nh4_cmipv2
integer :: id_so2_cmip, id_dms_cmip
integer :: id_so2_cmipv2, id_dms_cmipv2

!-----------------------------------------------------------------------
type(time_type) :: Time

!---- version number -----
character(len=128) :: version = '$Id: atmos_tracer_driver.F90,v 20.0 2013/12/13 23:24:10 fms Exp $'
character(len=128) :: tagname = '$Name: tikal $'
!-----------------------------------------------------------------------

contains

!#######################################################################

! <SUBROUTINE NAME="atmos_tracer_driver">
!   <OVERVIEW>
!     A routine which allows tracer code to be called.
!   </OVERVIEW>
!   <DESCRIPTION>
!     This subroutine calls the source sink routines for atmospheric
!     tracers. This is the interface between the dynamical core of the 
!     model and the tracer code. It should supply all the necessary 
!     information to a user that they need in order to calculate the 
!     tendency of that tracer with respect to emissions or chemical losses.
!
!   </DESCRIPTION>
!   <TEMPLATE>
!     call atmos_tracer_driver (is, ie, js, je, Time, lon, lat, land, phalf, pfull, r,  &
!                           u, v, t, q, u_star, rdt, rm, rdiag, kbot)
!   </TEMPLATE>
!   <IN NAME="is, ie, js, je" TYPE="integer">
!     Local domain boundaries.
!   </IN>
!   <IN NAME="Time" TYPE="type(time_type)">
!     Model time.
!   </IN>
!   <IN NAME="lon" TYPE="real" DIM="(:,:)">
!     Longitude of the centre of the model gridcells
!   </IN>
!   <IN NAME="lat" TYPE="real" DIM="(:,:)">
!     Latitude of the centre of the model gridcells
!   </IN>
!   <IN NAME="land" TYPE="logical" DIM="(:,:)">
!     Land/sea mask.
!   </IN>
!   <IN NAME="phalf" TYPE="real" DIM="(:,:,:)">
!     Pressures on the model half levels.
!   </IN>
!   <IN NAME="pfull" TYPE="real" DIM="(:,:,:)">
!     Pressures on the model full levels.
!   </IN>
!   <IN NAME="r" TYPE="real" DIM="(:,:,:,:)">
!     The tracer array in the component model.
!   </IN>
!   <IN NAME="u" TYPE="real" DIM="(:,:,:)">
!     Zonal wind speed.
!   </IN>
!   <IN NAME="v" TYPE="real" DIM="(:,:,:)">
!     Meridonal wind speed.
!   </IN>
!   <IN NAME="t" TYPE="real" DIM="(:,:,:)">
!     Temperature.
!   </IN>
!   <IN NAME="q" TYPE="real" DIM="(:,:,:)">
!     Specific humidity. This may also be accessible as a
!                        portion of the tracer array.
!   </IN>
!   <IN NAME="u_star" TYPE="real" DIM="(:,:)">
!     Friction velocity :: 
!     The magnitude of the wind stress is density*(ustar**2)
!     The drag coefficient for momentum is u_star**2/(u**2+v**2)
!   </IN>
!   <INOUT NAME="rdt" TYPE="real" DIM="(:,:,:,:)">
!     The tendency of the tracer array in the compenent
!     model. The tendency due to sources and sinks computed
!     in the individual tracer routines should be added to
!     this array before exiting tracer_driver.
!   </INOUT>
!   <IN NAME="rm" TYPE="real" DIM="(:,:,:,:)">
!     The tracer array in the component model for the previous timestep.
!   </IN>
!++amf
!   <IN NAME="flux_sw_down_vis_dir" TYPE="real" DIM="(:,:)">
!     Visible direct radiation at the surface in W / m2
!   </IN>
!   <IN NAME="flux_sw_down_vis_dif" TYPE="real" DIM="(:,:)">
!     Visible diffuse radiation at the surface in W / m2
!   </IN>
!--amf
!   <INOUT NAME="rdiag" TYPE="real" DIM="(:,:,:,:)">
!     The array of diagnostic tracers. As these may be changed within the
!     tracer routines for diagnostic purposes, they need to be writable.
!   </INOUT>
!   <IN NAME="kbot" TYPE="integer, optional" DIM="(:,:)">
!     Integer array describing which model layer intercepts the surface.
!   </IN>
 subroutine atmos_tracer_driver (is, ie, js, je, Time, lon, lat,  &
                           area, z_pbl, rough_mom, &
                           frac_open_sea, &
                           land, phalf, pfull,     &
                           u, v, t, q, r,          &
                           rm, rdt, dt,            &
                           u_star, b_star, q_star, &
                           z_half, z_full,         &
                           t_surf_rad, albedo,     &
                           Time_next,              &
                           flux_sw_down_vis_dir,   &
                           flux_sw_down_vis_dif,   &
                           mask,                   &
                           kbot)

!-----------------------------------------------------------------------
integer, intent(in)                           :: is, ie, js, je
type(time_type), intent(in)                   :: Time
real, intent(in),    dimension(:,:)           :: lon, lat
real, intent(in),    dimension(:,:)           :: u_star, b_star, q_star
real, intent(in),    dimension(:,:)           :: land
real, intent(in),    dimension(:,:)           :: area, z_pbl, rough_mom
real, intent(in),    dimension(:,:)           :: frac_open_sea
real, intent(in),    dimension(:,:,:)         :: phalf, pfull
real, intent(in),    dimension(:,:,:)         :: u, v, t, q
real, intent(inout), dimension(:,:,:,:)       :: r
real, intent(inout), dimension(:,:,:,:)       :: rm
real, intent(inout), dimension(:,:,:,:)       :: rdt
real, intent(in)                              :: dt !timestep(used in chem_interface)
real, intent(in),    dimension(:,:,:)         :: z_half !height in meters at half levels
real, intent(in),    dimension(:,:,:)         :: z_full !height in meters at full levels
real, intent(in),    dimension(:,:)           :: t_surf_rad !surface temperature
real, intent(in),    dimension(:,:)           :: albedo
real, intent(in), dimension(:,:)              :: flux_sw_down_vis_dir
real, intent(in), dimension(:,:)              :: flux_sw_down_vis_dif
type(time_type), intent(in)                   :: Time_next
integer, intent(in), dimension(:,:), optional :: kbot
real, intent(in), dimension(:,:,:),  optional :: mask

!-----------------------------------------------------------------------
! Local variables
!-----------------------------------------------------------------------
real, dimension(size(r,1),size(r,2),size(r,3)) :: rtnd, pwt, ozone, o3_prod, &
                                                  aerosol, rho
real, dimension(size(r,1),size(r,2),size(r,3)) :: rtndso2, rtndso4
real, dimension(size(r,1),size(r,2),size(r,3)) :: rtnddms, rtndmsa, rtndh2o2
real, dimension(size(r,1),size(r,2),size(r,3)) :: rtndbcphob, rtndbcphil
real, dimension(size(r,1),size(r,2),size(r,3)) :: rtndomphob, rtndomphil
real, dimension(size(r,1),size(r,2),size(r,3)) :: rtndco2, rtndco2_emis
real, dimension(size(r,1),size(r,2),size(rdt,4)) :: dsinku
real, dimension(size(r,1),size(r,2)) ::  w10m_ocean, w10m_land
integer :: year,month,day,hour,minute,second
integer :: jday
real, dimension(size(rdt,1),size(rdt,2),size(rdt,3),size(rdt,4)) :: chem_tend
real, dimension(size(r,1),size(r,2))           :: coszen, fracday, half_day
real :: rrsun
real, dimension(size(r,1),size(r,2),size(r,3)) :: cldf ! cloud fraction
real, dimension(size(r,1),size(r,2),size(r,3)) :: rh  ! relative humidity
real, dimension(size(r,1),size(r,2),size(r,3)) :: lwc ! liq water content
real, dimension(size(r,1),size(r,2),size(r,3),size(r,4)) :: tracer
real, dimension(size(r,1),size(r,3)) :: dp, temp
real, dimension(size(r,1),size(r,2),5) ::  ssalt_settl, dust_settl              
real, dimension(size(r,1),size(r,2),5) ::  ssalt_emis, dust_emis                
real, dimension(size(r,1),size(r,2)) ::  all_salt_settl, all_dust_settl         
real, dimension(size(r,1),size(r,2)) ::  suma, ocn_flx_fraction
real, dimension(size(r,1),size(r,2)) ::  frland, frsnow, frsea, frice

integer :: j, k, id, jd, kd, nt
integer :: nqq  ! index of specific humidity
integer :: nql  ! index of cloud liquid specific humidity
integer :: nqi  ! index of cloud ice water specific humidity
integer :: nqa  ! index of cloud amount
integer :: n, nnn
logical :: used


!-----------------------------------------------------------------------

!   <ERROR MSG="tracer_driver_init must be called first." STATUS="FATAL">
!     Tracer_driver_init needs to be called before tracer_driver.
!   </ERROR>
      if (.not. module_is_initialized)  &
      call error_mesg ('Tracer_driver','tracer_driver_init must be called first.', FATAL)

!-----------------------------------------------------------------------
     id=size(r,1);jd=size(r,2);kd=size(r,3); nt=size(r,4)

      nqa = get_tracer_index(MODEL_ATMOS,'cld_amt')
      nqi = get_tracer_index(MODEL_ATMOS,'ice_wat')
      nql = get_tracer_index(MODEL_ATMOS,'liq_wat')
      nqq = get_tracer_index(MODEL_ATMOS,'sphum')

      ssalt_settl = 0.
      ssalt_emis = 0.
      dust_settl = 0.
      dust_emis = 0.
!------------------------------------------------------------------------
! Make local copies of all the tracers
!------------------------------------------------------------------------
!++lwh
      if (use_tau) then
        do n = 1,nt
          tracer(:,:,:,n) = r(:,:,:,n)
!------------------------------------------------------------------------
! For tracers other than specific humdity, cloud amount, ice water and &
! liquid water fill eventual negative values
!------------------------------------------------------------------------
!        does tracer need to be adjusted to remain positive definite?
!         if (n /= nqq .and. n/=nqa .and. n/=nqi .and. n/=nql) then
          if ( adjust_positive_def(MODEL_ATMOS,n) ) then
!
            do j=1,jd
              do k=1,kd
                temp(:,k) = tracer(:,j,k,n)
                dp(:,k) = phalf(:,j,k+1)-phalf(:,j,k)
              enddo
              call sjl_fillz(id,kd,1,temp,dp)
              do k=1,kd
                tracer(:,j,k,n) = temp(:,k)
              enddo
            enddo
            if (n <= ntp) then
               rdt(:,:,:,n) = rdt(:,:,:,n) + (tracer(:,:,:,n)-r(:,:,:,n))/dt
            end if
          endif
        end do
      else
        do n = 1,nt
          if (n <= ntp) then
             tracer(:,:,:,n)=rm(:,:,:,n)+rdt(:,:,:,n)*dt
          else
             tracer(:,:,:,n)=r(:,:,:,n)
          end if
!        does tracer need to be adjusted to remain positive definite?
!         if (n /= nqq .and. n/=nqa .and. n/=nqi .and. n/=nql) then
          if ( adjust_positive_def(MODEL_ATMOS,n) ) then
            do j=1,jd
              do k=1,kd
                temp(:,k) = tracer(:,j,k,n)
                dp(:,k) = phalf(:,j,k+1)-phalf(:,j,k)
              enddo
              call sjl_fillz(id,kd,1,temp,dp)
              do k=1,kd
                tracer(:,j,k,n) = temp(:,k)
              enddo
            enddo
            if (n <= ntp) then
               rdt(:,:,:,n) = (tracer(:,:,:,n) - rm(:,:,:,n)) /dt
            end if
          end if
        end do
      end if
!--lwh

!------------------------------------------------------------------------
! Rediagnose meteoroligical variables. Note these parameterizations
! are not consistent with those used elsewhere in the GCM
!------------------------------------------------------------------------

!-----------------------------------------------------------------------
!-----------Calculate relative humidity
!-----------------------------------------------------------------------
      call get_rh(t,q,pfull,rh,mask)
!-----------------------------------------------------------------------
!--------- Calculate wind speed at 10 meters
!-----------------------------------------------------------------------
      call get_w10m(z_full(:,:,kd) - z_half(:,:,kd+1), &
                    u(:,:,kd), v(:,:,kd), &
                    rough_mom, u_star, b_star, q_star, &
                    w10m_ocean, w10m_land, Time_next, is, js)
!-----------------------------------------------------------------------
!------Cloud liquid water content
!-----------------------------------------------------------------------
      if (nqi > 0) then
        lwc(:,:,:)=max(tracer(:,:,:,nqi),0.) 
      else
        lwc(:,:,:) = 0.0
      endif
      if (nql > 0) lwc(:,:,:) = lwc(:,:,:) + max(tracer(:,:,:,nql),0.) 
!-----------------------------------------------------------------------
!--------- Cloud fraction -----------------------------
!-----------------------------------------------------------------------
      if (nqa > 0 ) then
!-----  cloud fraction is a prognostic variable ------------
        cldf(:,:,:)= max(0.,min(1.,tracer(:,:,:,nqa) ))
      else
!-----   cloud fraction estimated from RH-------------
        call get_cldf(phalf(:,:,kd+1),pfull,rh,cldf)
      endif
!-----------------------------------------------------------------------
!--------- Get Julian data
!-----------------------------------------------------------------------
      call get_date_julian(Time, year, month, jday, hour, minute, second)
!-----------------------------------------------------------------------
!--------- Get current date
!-----------------------------------------------------------------------
      call get_date(Time, year, month, day, hour, minute, second)

! Calculate cosine of solar zenith angle
!
       call diurnal_solar( lat, lon, Time, coszen, fracday, &
                           rrsun, dt_time=real_to_time_type(dt), &
                           half_day_out=half_day )

!------------------------------------------------------------------------
! Get air mass in layer (in kg/m2), equal to dP/g
!------------------------------------------------------------------------
      do k=1,kd
         pwt(:,:,k)=(phalf(:,:,k+1)-phalf(:,:,k))/grav
         rho(:,:,k) = pwt(:,:,k)/(z_half(:,:,k) - z_half(:,:,k+1))
      enddo

!------------------------------------------------------------------------
!    define various land fractions needed for dry deposition calculation.
!------------------------------------------------------------------------
      frland(:,:) = min(1., max(0.,     land(:,:) ) )
      frsea(:,:)  = min(1., max(0., 1. - frland(:,:) ) )
!      frsnow(:,:) = min(frland(:,:), max(0., snow_area(:,:) ) )
      frice(:,:)  = min(frsea(:,:),   &
                               max(0., frsea(:,:) - frac_open_sea(:,:) ) ) 

!------------------------------------------------------------------------
!    output land fraction information, if desired.
!------------------------------------------------------------------------
   used = send_data ( id_landfr, frland, Time_next, is_in =is,js_in=js)
   used = send_data ( id_seaicefr, frice, Time_next, is_in =is,js_in=js)
!   used = send_data ( id_snowfr, frsnow, Time_next, is_in =is,js_in=js)
!   used = send_data ( id_vegnfr, vegn_cover, Time_next, is_in =is,js_in=js)
!   used = send_data ( id_vegnlai, vegn_lai, Time_next, is_in =is,js_in=js)

!------------------------------------------------------------------------
! For tracers other than specific humdity, cloud amount, ice water and &
! liquid water calculate flux at surface due to dry deposition
!------------------------------------------------------------------------
!++lwh
      do n=1,ntp
         if (n /= nqq .and. n/=nqa .and. n/=nqi .and. n/=nql) then
            call dry_deposition( n, is, js, u(:,:,kd), v(:,:,kd), t(:,:,kd), &
                                 pwt(:,:,kd), pfull(:,:,kd), &
                                 z_half(:,:,kd)-z_half(:,:,kd+1), u_star, &
                                 (land > 0.5), dsinku(:,:,n), &
                                 tracer(:,:,kd,n), Time, Time_next, &
                                 lon, half_day, &
                                 drydep_data(n))!, frland, frice, frsnow, &
!                                 vegn_cover, vegn_lai, &
!                                 b_star, z_pbl, rough_mom)
            rdt(:,:,kd,n) = rdt(:,:,kd,n) - dsinku(:,:,n)
         end if
      enddo

      if (id_om_ddep > 0) then
        used  = send_data (id_om_ddep,  &
         pwt(:,:,kd)*(dsinku(:,:,nomphilic) + dsinku(:,:,nomphobic)),  &
                                              Time_next, is_in=is, js_in=js)
      endif
      if (id_bc_ddep > 0) then
        used  = send_data (id_bc_ddep,  &
         pwt(:,:,kd)*(dsinku(:,:,nbcphilic) + dsinku(:,:,nbcphobic)),  &
                                               Time_next, is_in=is, js_in=js)
      endif
      if (id_nh4_ddep_cmip > 0) then
        used  = send_data (id_nh4_ddep_cmip,  &
        0.018*1.0e03*pwt(:,:,kd)*(dsinku(:,:,nNH4NO3) + dsinku(:,:,nNH4))/WTMAIR,  &
                                              Time_next, is_in=is, js_in=js)
      endif

!----------------------------------------------------------------------
!   output the nh4no3 and nh4 loads.
!----------------------------------------------------------------------
      if(id_nh4_col > 0) then
        suma = 0.
        do k=1,kd
          suma(:,:) = suma(:,:) + pwt(:,:,k)*(tracer(:,:,k,nNH4) + &
                           tracer(:,:,k,nNH4NO3))
        end do
        used  = send_data (id_nh4_col,  &
               0.018*1.0e03*suma(:,:)/WTMAIR,  &
                                              Time_next, is_in=is, js_in=js)
      endif
      if(id_nh4no3_col > 0) then
        suma = 0.
        do k=1,kd
          suma(:,:) = suma(:,:) + pwt(:,:,k)*tracer(:,:,k,nNH4NO3)
        end do
        used  = send_data (id_nh4no3_col,  &
                  0.062*1.0e03*suma(:,:)/WTMAIR,  &
                                              Time_next, is_in=is, js_in=js)
      endif

!----------------------------------------------------------------------
!   output the tracer fields needed for CMIP.
!----------------------------------------------------------------------
      if (id_nh4_cmip > 0) then
        used  = send_data (id_nh4_cmip,  &
               0.018*1.0e03* (tracer(:,:,:,nNH4NO3) + &
                              tracer(:,:,:,nNH4)) /WTMAIR,  &
                                          Time_next, is_in=is, js_in=js, ks_in=1)
      endif
      if (id_nh4_cmipv2 > 0) then
        used  = send_data (id_nh4_cmipv2,  &
               0.018*1.0e03*rho(:,:,:)* (tracer(:,:,:,nNH4NO3) + &
                              tracer(:,:,:,nNH4)) /WTMAIR,  &
                                          Time_next, is_in=is, js_in=js, ks_in=1)
      endif
      if(id_nh4no3_cmip > 0) then
        used  = send_data (id_nh4no3_cmip,  &
                0.062*1.0e03*tracer(:,:,:,nNH4NO3)/WTMAIR,  &
                                         Time_next, is_in=is, js_in=js, ks_in=1)
     endif
      if(id_nh4no3_cmipv2 > 0) then
        used  = send_data (id_nh4no3_cmipv2,  &
                0.062*1.0e03*rho(:,:,:)*tracer(:,:,:,nNH4NO3)/WTMAIR,  &
                                         Time_next, is_in=is, js_in=js, ks_in=1)
     endif
     if(id_so2_cmip > 0) then
       used  = send_data (id_so2_cmip,  &
                 0.064*1.0e03*tracer(:,:,:,nSO2_cmip)/WTMAIR,  &
                                         Time_next, is_in=is, js_in=js, ks_in=1)
     endif
     if(id_so2_cmipv2 > 0) then
       used  = send_data (id_so2_cmipv2,  &
               0.064*1.0e03*rho(:,:,:)*tracer(:,:,:,nSO2_cmip)/WTMAIR,  &
                                         Time_next, is_in=is, js_in=js, ks_in=1)
     endif
     if(id_dms_cmip > 0) then
       used  = send_data (id_dms_cmip,  &
                0.062*1.0e03*tracer(:,:,:,nDMS_cmip)/WTMAIR,  &
                                         Time_next, is_in=is, js_in=js, ks_in=1)
     endif
     if(id_dms_cmipv2 > 0) then
       used  = send_data (id_dms_cmipv2,  &
                0.062*1.0e03*rho(:,:,:)*tracer(:,:,:,nDMS_cmip)/WTMAIR,  &
                                         Time_next, is_in=is, js_in=js, ks_in=1)
     endif

!------------------------------------------------------------------------
! Compute radon source-sink tendency
!------------------------------------------------------------------------
    call mpp_clock_begin (radon_clock)
    do nnn = 1, size(nradon(:))
     if (nradon(nnn) > 0) then
       if (nradon(nnn) > nt) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for radon', FATAL)
         call atmos_radon_sourcesink (lon,lat,land,pwt,tracer(:,:,:,nradon(nnn)),  &
                                 rtnd, Time, kbot)
       rdt(:,:,:,nradon(nnn))=rdt(:,:,:,nradon(nnn))+rtnd(:,:,:)
    endif
 
   end do
   call mpp_clock_end (radon_clock)

!------------------------------------------------------------------------
! Compute convection tracer source-sink tendency
!------------------------------------------------------------------------
   call mpp_clock_begin (convect_clock)
   do nnn = 1, size(nconvect(:))
     if (nconvect(nnn) > 0) then
       if (nconvect(nnn) > nt) call error_mesg ('Tracer_driver', &
           'Number of tracers .lt. number for convection tracer', FATAL)
       call atmos_cnvct_tracer_sourcesink (lon,lat,land,pwt,  &
                                       tracer(:,:,:,nconvect(nnn)),  &
                                       rtnd, Time, is, ie, js, je,kbot)
       rdt(:,:,:,nconvect(nnn))=rdt(:,:,:,nconvect(nnn))+rtnd(:,:,:)
     endif
   end do
   call mpp_clock_end (convect_clock)

!------------------------------------------------------------------------
! Stratospheric chemistry
!------------------------------------------------------------------------
  if(do_coupled_stratozone) then
    call mpp_clock_begin (stratozone_clock)
    itime(:) = 0
    call get_date(time,itime(1),itime(2),itime(3),itime(4),itime(5),itime(6))

    call strat_chem(lon,lat,r,rdt,pfull,t,itime,is,ie,js,je,dt,coszen,  &
!    ozon,cosp,cosphc,photo,solardata,chlb,ozb,dfdage,tropc,anoy,        &
    nsphum,chem_tend,ozone,o3_prod,aerosol,mpp_pe())
!
! The water vapour tendency is included here, partially coupling the 
! chemistry to the radiation
! 
      rdt(:,:,:,:) = rdt(:,:,:,:) + chem_tend(:,:,:,1:ntp) 
      if(nt.gt.(ntp+1))  then
! Modify the diagnostic tracers.
        tracer(:,:,:,no3)      = ozone(:,:,:) 
        tracer(:,:,:,no3ch)    = o3_prod(:,:,:) 
        if (naerosol > 0 ) tracer(:,:,:,naerosol) = aerosol(:,:,:) 
      endif
    call mpp_clock_end (stratozone_clock)
  endif

!-----------------------------------------------------------------------
!   define sfc area through which tracer flux from ocean is allowed; it
!   is area fraction of sea unless prevent_flux_through_ice is .true.,
!   in which case it is the fractional area of open sea.
!-----------------------------------------------------------------------
   if (prevent_flux_through_ice) then
     ocn_flx_fraction = frac_open_sea
   else
     ocn_flx_fraction = 1. - land       
   endif

!------------------------------------------------------------------------
! Tropospheric chemistry
!------------------------------------------------------------------------
   if ( do_tropchem ) then
!------------------------------------------------------------------------
! Compute age tracer source-sink tendency
!------------------------------------------------------------------------
      if (nage > 0) then
        call mpp_clock_begin (age_tracer_clock)
        if (nage > nt) call error_mesg ('Tracer_driver', &
           'Number of tracers .lt. number for age tracer', FATAL)
        call atmos_age_tracer( lon, lat, pwt,  &
                               tracer(:,:,:,nage),  &
                               rtnd, Time, kbot)
          rdt(:,:,:,nage)=rdt(:,:,:,nage)+rtnd(:,:,:)
        call mpp_clock_end (age_tracer_clock)
      endif

      call mpp_clock_begin (tropchem_clock)
      call tropchem_driver( lon, lat, land, ocn_flx_fraction, pwt, &
                            tracer(:,:,:,1:ntp),chem_tend, &
                            Time, phalf, pfull, t, is, ie, js, je, dt, &
                            z_half, z_full, q, t_surf_rad, albedo, coszen, rrsun, &
                            area, w10m_ocean, &
                            flux_sw_down_vis_dir, flux_sw_down_vis_dif, & 
                            half_day, &
                            Time_next, tracer(:,:,:,MIN(ntp+1,nt):nt), kbot)
      rdt(:,:,:,:) = rdt(:,:,:,:) + chem_tend(:,:,:,:)
      call mpp_clock_end (tropchem_clock)
   endif

!! RSH 4/8/04
!! note that if there are no diagnostic tracers, that argument in the
!! call to sourcesink should be made optional and omitted in the calls
!! below. note the switch in argument order to make this argument
!! optional.
    if (nt == ntp) then  ! implies no diagnostic tracers
!   if(do_tropchem) then
!      if(present(kbot)) then
!        call sourcesink(lon,lat,land,pwt,r+rdt*dt,chem_tend,Time,phalf,pfull,t,is,js,je,dt,&
!                          z_half, z_full,q,t_surf_rad,albedo,coszen, Time_next,&
!                          u,v,u_star,&
!                          kbot)
!      else
!       call sourcesink(lon,lat,land,pwt,r+rdt*dt,chem_tend,Time,phalf,pfull,t,is,js,je,dt, &
!                          z_half, z_full,q, t_surf_rad, albedo, coszen, Time_next, &
!                          u,v,u_star,&
!                          )
!      endif
!      rdt(:,:,:,:) = rdt(:,:,:,:) + chem_tend(:,:,:,:)
!   endif        
      else   ! case of diagnostic tracers being present
!   if(do_tropchem) then
!      if(present(kbot)) then
!        call sourcesink(lon,lat,land,pwt,r+rdt*dt,chem_tend,Time,phalf,pfull,t,is,js,je,dt,&
!                          z_half, z_full,q,t_surf_rad,albedo,coszen, Time_next,&
!                          u,v,u_star,&
!       rdiag=rm(:,:,:,nt+1:ntp),  &  ! (the diagnostic tracers)
!                   kbot=kbot)
!      else
!       call sourcesink(lon,lat,land,pwt,r+rdt*dt,chem_tend,Time,phalf,pfull,t,is,js,je,dt, &
!                          z_half, z_full,q, t_surf_rad, albedo, coszen, Time_next, &
!                          u,v,u_star,&
!        rdiag=rm(:,:,:,nt+1:ntp),  & ! (the diagnostic tracers)
!                          )
!      endif
!      rdt(:,:,:,:) = rdt(:,:,:,:) + chem_tend(:,:,:,:)
!   endif        
   endif  ! (no diagnostic tracers)

!------------------------------------------------------------------------
!   carbonaceous aerosols
!------------------------------------------------------------------------
   if (nbcphobic > 0 .and. nbcphilic > 0 .and. &
       nomphobic > 0 .and. nomphilic > 0) then
         if (nbcphobic > ntp .or. nbcphilic > ntp) &
            call error_mesg ('Tracer_driver', &
            'Number of tracers .lt. number for black carbon', FATAL)
     call mpp_clock_begin (carbon_clock)
     call atmos_carbon_aerosol_driver(lon,lat,ocn_flx_fraction, pfull,phalf,z_half,z_pbl, &
                                      t_surf_rad, w10m_ocean, &
                                      T, pwt, &
                                      tracer(:,:,:,nbcphobic), rtndbcphob, &
                                      tracer(:,:,:,nbcphilic), rtndbcphil, &
                                      tracer(:,:,:,nomphobic), rtndomphob, &
                                      tracer(:,:,:,nomphilic), rtndomphil, &
                                      tracer(:,:,:,noh),    &
                                      Time_next,is,ie,js,je)
      rdt(:,:,:,nbcphobic)=rdt(:,:,:,nbcphobic)+rtndbcphob(:,:,:)
      rdt(:,:,:,nbcphilic)=rdt(:,:,:,nbcphilic)+rtndbcphil(:,:,:)
      rdt(:,:,:,nomphobic)=rdt(:,:,:,nomphobic)+rtndomphob(:,:,:)
      rdt(:,:,:,nomphilic)=rdt(:,:,:,nomphilic)+rtndomphil(:,:,:)
      call mpp_clock_end (carbon_clock)
   endif
!------------------------------------------------------------------------
! Mineral Dust 
!------------------------------------------------------------------------
   call mpp_clock_begin (dust_clock)
   if (ndust1 > 0) then
         if (ndust1 > ntp ) call error_mesg ('Tracer_driver', &
                   'Number of tracers .lt. number for dust', FATAL)
         call atmos_dust_sourcesink (&
              1, 0.1e-6, 1.0e-6, 0.75e-6, 2500., &
              lon,lat,land,pwt, &
              z_half, pfull, w10m_land, t, rh, &
              tracer(:,:,:,ndust1), rtnd, dust_emis(:,:,1), &
              dust_settl(:,:,1), Time, Time_next, &
              is,ie,js,je, kbot)
      rdt(:,:,:,ndust1)=rdt(:,:,:,ndust1)+rtnd(:,:,:)
   endif

   if (ndust2 > 0) then
         if (ndust2 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for dust', FATAL)
         call atmos_dust_sourcesink (&
              2, 1.e-6, 2.e-6, 1.5e-6, 2650., &
              lon,lat,land,pwt, &
              z_half, pfull, w10m_land, t, rh, &
              tracer(:,:,:,ndust2), rtnd, dust_emis(:,:,2), &
              dust_settl(:,:,2), Time, Time_next, &
              is,ie,js,je, kbot)
      rdt(:,:,:,ndust2)=rdt(:,:,:,ndust2)+rtnd(:,:,:)
   endif

   if (ndust3 > 0) then
         if (ndust3 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for dust', FATAL)
         call atmos_dust_sourcesink (&
              3, 2.e-6, 3.e-6, 2.5e-6, 2650., &
              lon,lat,land,pwt, &
              z_half, pfull, w10m_land, t, rh, &
              tracer(:,:,:,ndust3), rtnd, dust_emis(:,:,3), &
              dust_settl(:,:,3), Time, Time_next, &
              is,ie,js,je, kbot)
      rdt(:,:,:,ndust3)=rdt(:,:,:,ndust3)+rtnd(:,:,:)
   endif

   if (ndust4 > 0) then
         if (ndust4 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for dust', FATAL)
         call atmos_dust_sourcesink (&
              4, 3.e-6, 6.e-6, 4.5e-6, 2650., &
              lon,lat,land,pwt, &
              z_half, pfull, w10m_land, t, rh, &
              tracer(:,:,:,ndust4), rtnd, dust_emis(:,:,4), &
              dust_settl(:,:,4), Time, Time_next, &
              is,ie,js,je, kbot)
      rdt(:,:,:,ndust4)=rdt(:,:,:,ndust4)+rtnd(:,:,:)
   endif

   if (ndust5 > 0) then
         if (ndust5 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for dust', FATAL)
         call atmos_dust_sourcesink (&
              5, 6.e-6, 10.e-6, 8.e-6, 2650., &
              lon,lat,land,pwt, &
              z_half, pfull, w10m_land, t, rh, &
              tracer(:,:,:,ndust5), rtnd, dust_emis(:,:,5), &
              dust_settl(:,:,5), Time, Time_next, &
              is,ie,js,je, kbot)
      rdt(:,:,:,ndust5)=rdt(:,:,:,ndust5)+rtnd(:,:,:)
   endif
   if (id_dust_ddep > 0) then
     all_dust_settl(:,:) = dust_settl(:,:,1) +  &
                dust_settl(:,:,2) +  dust_settl(:,:,3) +  &
                dust_settl(:,:,4) +  dust_settl(:,:,5)
     used  = send_data (id_dust_ddep,  all_dust_settl(:,:) + &
        pwt(:,:,kd)*(dsinku(:,:,ndust1) + dsinku(:,:,ndust2) + &
                     dsinku(:,:,ndust3) + dsinku(:,:,ndust4) + &
                      dsinku(:,:,ndust5)), Time_next, is_in=is, js_in=js)
   endif
   if (id_dust_emis > 0) then
     used  = send_data (id_dust_emis,  &
                 dust_emis(:,:,1) + dust_emis(:,:,2) + &
                 dust_emis(:,:,3) + dust_emis(:,:,4) + &
                 dust_emis(:,:,5), Time_next, is_in=is, js_in=js)
   endif
   call mpp_clock_end (dust_clock)

!------------------------------------------------------------------------
!sea salt
!------------------------------------------------------------------------
   call mpp_clock_begin (seasalt_clock)
   if (nseasalt1 > 0) then
         if (nseasalt1 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for sea salt', FATAL)
         rtnd(:,:,:) = 0.
         call atmos_sea_salt_sourcesink ( &
              1, 0.1e-6,0.5e-6,0.3e-6, 2200., &
              lon,lat,ocn_flx_fraction,pwt, &
              z_half, pfull, w10m_ocean, t, rh, &
              tracer(:,:,:,nseasalt1), rtnd, dt, &
               ssalt_settl(:,:,1), ssalt_emis(:,:,1), &
              Time,Time_next,is,ie,js,je, kbot)
      rdt(:,:,:,nseasalt1)=rdt(:,:,:,nseasalt1)+rtnd(:,:,:)
   endif
   if (nseasalt2 > 0) then
         if (nseasalt2 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for sea salt', FATAL)
         rtnd(:,:,:) = 0.
         call atmos_sea_salt_sourcesink ( &
              2, 0.5e-6,1.0e-6,0.75e-6, 2200., &
              lon,lat,ocn_flx_fraction,pwt, &
              z_half, pfull, w10m_ocean, t, rh, &
              tracer(:,:,:,nseasalt2), rtnd, dt, &
               ssalt_settl(:,:,2), ssalt_emis(:,:,2), &
              Time,Time_next,is,ie,js,je, kbot)
      rdt(:,:,:,nseasalt2)=rdt(:,:,:,nseasalt2)+rtnd(:,:,:)
   endif
   if (nseasalt3 > 0) then
         if (nseasalt3 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for sea salt', FATAL)
         rtnd(:,:,:) = 0.
         call atmos_sea_salt_sourcesink ( &
              3, 1.e-6,2.5e-6,1.75e-6, 2200., &
              lon,lat,ocn_flx_fraction,pwt, &
              z_half, pfull, w10m_ocean, t, rh, &
              tracer(:,:,:,nseasalt3), rtnd, dt, &
               ssalt_settl(:,:,3), ssalt_emis(:,:,3), &
              Time,Time_next,is,ie,js,je, kbot)
      rdt(:,:,:,nseasalt3)=rdt(:,:,:,nseasalt3)+rtnd(:,:,:)
   endif
   if (nseasalt4 > 0) then
         if (nseasalt4 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for sea salt', FATAL)
         rtnd(:,:,:) = 0.
         call atmos_sea_salt_sourcesink ( &
              4, 2.5e-6,5.0e-6,3.75e-6, 2200., &
              lon,lat,ocn_flx_fraction,pwt, &
              z_half, pfull, w10m_ocean, t, rh, &
              tracer(:,:,:,nseasalt4), rtnd, dt, &
               ssalt_settl(:,:,4), ssalt_emis(:,:,4), &
              Time,Time_next,is,ie,js,je, kbot)
      rdt(:,:,:,nseasalt4)=rdt(:,:,:,nseasalt4)+rtnd(:,:,:)
   endif
   if (nseasalt5 > 0) then
         if (nseasalt5 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for sea salt', FATAL)
         rtnd(:,:,:) = 0.
         call atmos_sea_salt_sourcesink ( &
              5, 5.e-6,10.0e-6,7.5e-6, 2200., &
              lon,lat,ocn_flx_fraction,pwt, &
              z_half, pfull, w10m_ocean, t, rh, &
              tracer(:,:,:,nseasalt5), rtnd, dt, &
               ssalt_settl(:,:,5), ssalt_emis(:,:,5), &
              Time,Time_next,is,ie,js,je, kbot)
      rdt(:,:,:,nseasalt5)=rdt(:,:,:,nseasalt5)+rtnd(:,:,:)
   endif
   if (id_ssalt_ddep > 0) then
     all_salt_settl(:,:) = ssalt_settl(:,:,1) + ssalt_settl(:,:,2) +  &
                         ssalt_settl(:,:,3) + ssalt_settl(:,:,4) + &
                         ssalt_settl(:,:,5)
     used  = send_data (id_ssalt_ddep, all_salt_settl(:,:) + &
         pwt(:,:,kd)*(dsinku(:,:,nseasalt1) + dsinku(:,:,nseasalt2) + &
                      dsinku(:,:,nseasalt3) + dsinku(:,:,nseasalt4) + &
                      dsinku(:,:,nseasalt5)), Time_next, is_in=is, js_in=js)
   endif
   if (id_ssalt_emis > 0) then
     used  = send_data (id_ssalt_emis,  &
                    ssalt_emis(:,:,1) + ssalt_emis(:,:,2) + &
                    ssalt_emis(:,:,3) + ssalt_emis(:,:,4) + &
                    ssalt_emis(:,:,5), Time_next, is_in=is, js_in=js)
   endif
   call mpp_clock_end (seasalt_clock)

!------------------------------------------------------------------------
! Sulfur chemistry
!------------------------------------------------------------------------
   if (nDMS > 0 .and. nSO2 > 0 .and. nSO4 > 0 .and. nMSA > 0 ) then
      if (nDMS > ntp ) call error_mesg ('Tracer_driver', &
                     'Number of tracers .lt. number for DMS', FATAL)
      if (nSO2 > ntp ) call error_mesg ('Tracer_driver', &
                     'Number of tracers .lt. number for SO2', FATAL)
      if (nSO4 > ntp ) call error_mesg ('Tracer_driver', &
                     'Number of tracers .lt. number for SO4', FATAL)

      call mpp_clock_begin (sulfur_clock)
      call atmos_DMS_emission(lon, lat, area, ocn_flx_fraction, t_surf_rad, &
             w10m_ocean, pwt, rtnddms, Time, Time_next, is,ie,js,je,kbot)
      rdt(:,:,kd,nDMS) = rdt(:,:,kd,nDMS) + rtnddms(:,:,kd)
      call atmos_SOx_emission(lon, lat, area, land, &
             z_pbl, z_half, phalf, pwt, rtndso2, rtndso4, &
             Time, Time_next, is,ie,js,je,kbot)
      rdt(:,:,:,nSO2) = rdt(:,:,:,nSO2) + rtndso2(:,:,:)
      rdt(:,:,:,nSO4) = rdt(:,:,:,nSO4) + rtndso4(:,:,:)
      call atmos_SOx_chem( pwt, t, pfull, phalf, dt, lwc, &
                jday,hour,minute,second,lat,lon,    &
                tracer(:,:,:,nSO2), tracer(:,:,:,nSO4), tracer(:,:,:,nDMS), &
                tracer(:,:,:,nMSA), tracer(:,:,:,nH2O2), &
                tracer(:,:,:,noh), &
                rtndso2, rtndso4, rtnddms, rtndmsa, rtndh2o2, &
                Time,Time_next, is,ie,js,je,kbot)
      rdt(:,:,:,nSO2) = rdt(:,:,:,nSO2) + rtndso2(:,:,:)
      rdt(:,:,:,nSO4) = rdt(:,:,:,nSO4) + rtndso4(:,:,:)
      rdt(:,:,:,nDMS) = rdt(:,:,:,nDMS) + rtnddms(:,:,:)
      rdt(:,:,:,nMSA) = rdt(:,:,:,nMSA) + rtndmsa(:,:,:)
      rdt(:,:,:,nH2O2) = rdt(:,:,:,nH2O2) + rtndh2o2(:,:,:)
      call mpp_clock_end (sulfur_clock)
   endif

!------------------------------------------------------------------------
! Secondary organic aerosols
!------------------------------------------------------------------------
   if (nSOA > 0 ) then
      if (nSOA > ntp ) call error_mesg ('Tracer_driver', &
                     'Number of tracers .lt. number for SOA', FATAL)

      call mpp_clock_begin (SOA_clock)
      call atmos_SOA_chem(pwt,t,pfull,phalf,dt, &
                jday,hour,minute,second,lat,lon,    &
                tracer(:,:,:,nSOA),rtnd, Time,Time_next,is,ie,js,je,kbot )

      rdt(:,:,:,nSOA)=rdt(:,:,:,nSOA)+rtnd(:,:,:)
      call mpp_clock_end (SOA_clock)

   endif

!------------------------------------------------------------------------
! Sulfur hexafluoride (SF6)
!------------------------------------------------------------------------
   if (nsf6 > 0) then
         if (nsf6 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers .lt. number for sulfur hexafluoride', FATAL)
         call mpp_clock_begin (sf6_clock)
         call atmos_sf6_sourcesink (lon,lat,land,pwt,r(:,:,:,nsf6),  &
                                 rtnd, Time,is,ie,js,je, kbot)
      rdt(:,:,:,nsf6)=rdt(:,:,:,nsf6)+rtnd(:,:,:)
         call mpp_clock_end (sf6_clock)
   endif

   if (nch3i > 0) then
         if (nch3i > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers < number for ch3i', FATAL)
         call mpp_clock_begin (ch3i_clock)
         call atmos_ch3i( lon, lat, land, pwt, &
                          r(:,:,:,nch3i)+dt*rdt(:,:,:,nch3i), rtnd,  &
                          Time, phalf, pfull, t, is, js, dt, &
                          z_half, z_full,q, t_surf_rad, albedo, coszen, &
                          Time_next, kbot)
         rdt(:,:,:,nch3i) = rdt(:,:,:,nch3i) + rtnd(:,:,:)
         call mpp_clock_end (ch3i_clock)
   endif

!------------------------------------------------------------------------
! Compute CO2 source-sink tendency:  tracer(:,:,:,nco2) is in moist mmr
!  atmos_co2_sourcesink will convert to dry vmr and return rtnd as moist mmr
!------------------------------------------------------------------------
   if (nco2 > 0 ) then
         if (nco2 > ntp ) call error_mesg ('Tracer_driver', &
                            'Number of tracers < number for co2', FATAL)
         call mpp_clock_begin (co2_clock)
         call atmos_co2_emissions (is, ie, js, je, Time, Time_next, dt, pwt, tracer(:,:,:,nco2),     &
                                   tracer(:,:,:,nsphum), rtndco2_emis, kbot)
         rdt(:,:,:,nco2)=rdt(:,:,:,nco2)+rtndco2_emis(:,:,:)

         call atmos_co2_sourcesink (is, ie, js, je, Time, Time_next, dt, pwt, tracer(:,:,:,nco2),     &
                                    tracer(:,:,:,nsphum), rtndco2)
         rdt(:,:,:,nco2)=rdt(:,:,:,nco2)+rtndco2(:,:,:)
         call mpp_clock_end (co2_clock)
   endif

!------------------------------------------------------------------------
! Save diagnostic tracer concentrations back to tracer array
!------------------------------------------------------------------------
   if (nt > ntp) then
      r(:,:,:,ntp+1:nt) = tracer(:,:,:,ntp+1:nt)
   end if

 end subroutine atmos_tracer_driver
! </SUBROUTINE>

!#######################################################################

! <SUBROUTINE NAME="atmos_tracer_driver_init">
!   <OVERVIEW>
!     Subroutine to initialize the tracer driver module.
!   </OVERVIEW>
!   <DESCRIPTION>
!   The purpose of the arguments here are for passing on to the individual
!   tracer code. The user may wish to provide initial values which can be
!   implemented in the initialization part of the tracer code. Remember that
!   the tracer manager will provide a simple fixed or exponential profile if
!   the user provides data for this within the field table. However if a more
!   complicated profile is required then it should be set up in the
!   initialization section of the user tracer code.
!   </DESCRIPTION>
!   <TEMPLATE>
!     call atmos_tracer_driver_init (lonb,latb, r, mask, axes, Time)
!   </TEMPLATE>
!   <IN NAME="lonb" TYPE="real" DIM="(:,:)">
!     The longitude corners for the local domain.
!   </IN>
!   <IN NAME="latb" TYPE="real" DIM="(:,:)">
!     The latitude corners for the local domain.
!   </IN>
!   <IN NAME="mask" TYPE="real, optional" DIM="(:,:,:)">
!      optional mask (0. or 1.) that designates which grid points
!           are above (=1.) or below (=0.) the ground dimensioned as
!           (nlon,nlat,nlev).
!   </IN>
!   <IN NAME="Time" TYPE="type(time_type)">
!     Model time.
!   </IN>
!   <IN NAME="axes" TYPE="integer" DIM="(4)">
!     The axes relating to the tracer array dimensioned as
!      (nlon, nlat, nlev, ntime)
!   </IN>
!   <INOUT NAME="r" TYPE="real" DIM="(:,:,:,:)">
!     Tracer fields dimensioned as (nlon,nlat,nlev,ntrace). 
!   </INOUT>
 subroutine atmos_tracer_driver_init (lonb, latb, r, axes, Time, phalf, mask)

!-----------------------------------------------------------------------
           real, intent(in),    dimension(:,:)             :: lonb, latb
           real, intent(inout), dimension(:,:,:,:)         :: r
type(time_type), intent(in)                                :: Time
        integer, intent(in)                                :: axes(4)
           real, intent(in),    dimension(:,:,:)           :: phalf
           real, intent(in),    dimension(:,:,:), optional :: mask

!-----------------------------------------------------------------------
! Local variables
!-----------------------------------------------------------------------
      integer :: nbr_layers
      integer :: unit, ierr, io, logunit
  
!-----------------------------------------------------------------------
!
!  When initializing additional tracers, the user needs to make changes 
!
!-----------------------------------------------------------------------

      if (module_is_initialized) return

!------------------------------------------------------------------------
!   read namelist.
!------------------------------------------------------------------------
      if ( file_exist('input.nml')) then
#ifdef INTERNAL_FILE_NML
        read (input_nml_file, nml=atmos_tracer_driver_nml, iostat=io)
        ierr = check_nml_error(io,'atmos_tracer_driver_nml')
#else
        unit =  open_namelist_file ( )
        ierr=1; do while (ierr /= 0)
        read  (unit, nml=atmos_tracer_driver_nml, iostat=io, end=10)
        ierr = check_nml_error(io, 'atmos_tracer_driver_nml')
        end do
10      call close_file (unit)
#endif
      endif
!--------- write version and namelist to standard log ------------
      call write_version_number(version, tagname)
      logunit=stdlog()
      if ( mpp_pe() == mpp_root_pe() ) &
                         write ( logunit, nml=atmos_tracer_driver_nml )

!---------------------------------------------------------------------
!  make sure that astronomy_mod has been initialized (if radiation
!  not being called in this run, it will not have previously been 
!  initialized).
!---------------------------------------------------------------------
      call astronomy_init

!If we wish to automatically register diagnostics for wet and dry 
! deposition, do it now.
      call atmos_tracer_utilities_init(lonb, latb, axes, Time)

!----- set initial value of radon ------------

      call atmos_radon_init(r, axes, Time, nradon, mask)
      radon_clock = mpp_clock_id( 'Tracer: Radon', &
           grain=CLOCK_MODULE )

!----- initialize the convection tracer ------------

      call atmos_convection_tracer_init(r, phalf, axes, Time, &
                                        nconvect,  mask)
      convect_clock = mpp_clock_id( 'Tracer: Convection tracer', &
           grain=CLOCK_MODULE )

!----- initialize the age tracer ------------

      call atmos_age_tracer_init( r, axes, Time, nage, lonb, latb, phalf, mask)
      age_tracer_clock = mpp_clock_id( 'Tracer: Age tracer', grain=CLOCK_MODULE )

      call get_number_tracers (MODEL_ATMOS, num_tracers=nt, &
                               num_prog=ntp)


!------------------------------------------------------------------------
! Initialize stratospheric chemistry
!------------------------------------------------------------------------
      nsphum   = get_tracer_index(MODEL_ATMOS,'sphum')
      no3      = get_tracer_index(MODEL_ATMOS,'O3')
      no3ch    = get_tracer_index(MODEL_ATMOS,'O3_chem')
      nextinct = get_tracer_index(MODEL_ATMOS,'Extinction')
      naerosol = get_tracer_index(MODEL_ATMOS,'Aerosol')
!
!  Set up photolysis rates etc.
!
      do_coupled_stratozone = strat_chem_driver_init()
      stratozone_clock = mpp_clock_id( 'Tracer: Stratospheric Ozone', &
           grain=CLOCK_MODULE )
         

!------------------------------------------------------------------------
! Initialize tropospheric chemistry
!------------------------------------------------------------------------
      allocate( drydep_data(nt) )
      do_tropchem = tropchem_driver_init(r,mask,axes,Time,lonb,latb,phalf,drydep_data)
      tropchem_clock = mpp_clock_id( 'Tracer: Tropospheric chemistry', &
           grain=CLOCK_MODULE )

! ----------Interactive traceres--------------------
! If any of the interactive tracers are activated, get the 
! tracer number and initialize it.
!      ntraceer= get_tracer_index(MODEL_ATMOS,'tracer_name')
!      if (ntracer > 0) then
!        call {tracer}_init(lonb, latb, r, axes, Time, mask)
!      endif


!    get tracer indices
      nbcphobic = get_tracer_index(MODEL_ATMOS,'bcphob')
      nbcphilic = get_tracer_index(MODEL_ATMOS,'bcphil')
      nomphobic = get_tracer_index(MODEL_ATMOS,'omphob')
      nomphilic = get_tracer_index(MODEL_ATMOS,'omphil')
      ndust1    = get_tracer_index(MODEL_ATMOS,'dust1')
      ndust2    = get_tracer_index(MODEL_ATMOS,'dust2')
      ndust3    = get_tracer_index(MODEL_ATMOS,'dust3')
      ndust4    = get_tracer_index(MODEL_ATMOS,'dust4')
      ndust5    = get_tracer_index(MODEL_ATMOS,'dust5')
      nseasalt1 = get_tracer_index(MODEL_ATMOS,'ssalt1')
      nseasalt2 = get_tracer_index(MODEL_ATMOS,'ssalt2')
      nseasalt3 = get_tracer_index(MODEL_ATMOS,'ssalt3')
      nseasalt4 = get_tracer_index(MODEL_ATMOS,'ssalt4')
      nseasalt5 = get_tracer_index(MODEL_ATMOS,'ssalt5')
      nDMS      = get_tracer_index(MODEL_ATMOS,'simpleDMS')
      nSO2      = get_tracer_index(MODEL_ATMOS,'simpleSO2')
      nSO4      = get_tracer_index(MODEL_ATMOS,'simpleSO4')
      nMSA      = get_tracer_index(MODEL_ATMOS,'simpleMSA')
      nH2O2     = get_tracer_index(MODEL_ATMOS,'simpleH2O2')
      nnH4NO3   = get_tracer_index(MODEL_ATMOS,'nh4no3')
      nnH4      = get_tracer_index(MODEL_ATMOS,'nh4')
      nSOA      = get_tracer_index(MODEL_ATMOS,'SOA')
      nsf6      = get_tracer_index(MODEL_ATMOS,'sf6')
      nch3i     = get_tracer_index(MODEL_ATMOS,'ch3i')
      nco2      = get_tracer_index(MODEL_ATMOS,'co2')
      nDMS_cmip = get_tracer_index(MODEL_ATMOS,'DMS')
      nSO2_cmip = get_tracer_index(MODEL_ATMOS,'so2')
      noh       = get_tracer_index(MODEL_ATMOS,'oh')

! Number of vertical layers
      nbr_layers=size(r,3)

! initialize the tracers
!carbonaceous aerosols
      if (nbcphobic > 0 .or. nbcphilic >0 .or.  &
          nomphobic > 0 .or. nomphilic >0 ) then
        call atmos_carbon_aerosol_init(lonb, latb, axes, Time, mask)
        carbon_clock = mpp_clock_id( 'Tracer: Carbonaceous aerosol', &
                       grain=CLOCK_MODULE )

      endif
!dust aerosols
      if (ndust1 > 0.or.ndust2 > 0.or.ndust3 > 0 &
        .or.ndust4 > 0.or.ndust5 > 0 ) then
        call atmos_dust_init (lonb, latb, axes, Time, mask)
        dust_clock = mpp_clock_id( 'Tracer: Dust aerosol', &
                     grain=CLOCK_MODULE )
      endif
!sea salt
      if (nseasalt1 > 0.or.nseasalt2 > 0.or.nseasalt3 > 0 &
        .or.nseasalt4 > 0.or.nseasalt5 > 0 ) then
        call atmos_sea_salt_init (lonb, latb, axes, Time, mask)
        seasalt_clock = mpp_clock_id( 'Tracer: Seasalt aerosol', &
                        grain=CLOCK_MODULE )
      endif
!sulfur cycle
      if (nDMS > 0 .or. nSO2 > 0 .or. nSO4 > 0 &
                   .or. nMSA > 0 .or. nH2O2 > 0 ) then
        call atmos_sulfate_init ( lonb, latb, nbr_layers, axes, Time, mask)
        sulfur_clock = mpp_clock_id( 'Tracer: Sulfur', &
                       grain=CLOCK_MODULE )
      endif
!SOA
      if ( nSOA > 0 ) then
        call atmos_SOA_init ( lonb, latb, nbr_layers, axes, Time, mask)
        SOA_clock = mpp_clock_id( 'Tracer: SOA', &
                    grain=CLOCK_MODULE )
      endif
!sf6
      if (nsf6 > 0) then
        call atmos_sulfur_hex_init (lonb, latb, r, axes, Time, mask)
        sf6_clock = mpp_clock_id( 'Tracer: SF6', &
                    grain=CLOCK_MODULE )
      endif
!ch3i
      if (nch3i > 0) then
        call atmos_ch3i_init (lonb, latb, axes, Time, mask)
        ch3i_clock = mpp_clock_id( 'Tracer: CH3I', &
                     grain=CLOCK_MODULE )
      endif

!co2
      if (nco2 > 0) then
      call atmos_co2_init ( Time, axes(1:3))
        co2_clock = mpp_clock_id( 'Tracer: CO2', &
                    grain=CLOCK_MODULE )
      endif

      call get_number_tracers (MODEL_ATMOS, num_tracers=nt, &
                               num_prog=ntp)

     id_landfr = register_diag_field ( mod_name,                    &
            'landfr_atm', axes(1:2), Time,               &
            'land fraction',                                 &
            'fraction', missing_value=-999.     )
     id_seaicefr = register_diag_field ( mod_name,                    &
            'seaicefr_atm', axes(1:2), Time,               &
            'seaice fraction',                                 &
            'fraction', missing_value=-999.     )
     id_snowfr = register_diag_field ( mod_name,                    &
            'snowfr_atm', axes(1:2), Time,               &
            'snow cover fraction',                                 &
            'fraction', missing_value=-999.     )
     id_vegnfr = register_diag_field ( mod_name,                    &
            'vegnfr_atm', axes(1:2), Time,               &
            'vegetation cover fraction',                                 &
            'fraction', missing_value=-999.     )
     id_vegnlai = register_diag_field ( mod_name,                    &
            'vegnlai_atm', axes(1:2), Time,               &
            'vegetation leaf area index',                                &
            'fraction', missing_value=-999.     )

      id_om_ddep = register_diag_field (mod_name, &
          'om_ddep', axes(1:2), Time, &
          'total dry deposition of om', 'kg/m2/s')

      id_bc_ddep = register_diag_field (mod_name, &
          'bc_ddep', axes(1:2), Time, &
          'total dry deposition of bc', 'kg/m2/s')

      id_ssalt_ddep = register_diag_field (mod_name, &
          'ssalt_ddep', axes(1:2), Time, &
          'dry deposition and settling of seasalt', 'kg/m2/s')

      id_ssalt_emis = register_diag_field (mod_name, &
          'ssalt_emis', axes(1:2), Time, &
          'total emission of seasalt', 'kg/m2/s')

      id_dust_ddep = register_diag_field (mod_name, &
          'dust_ddep', axes(1:2), Time, &
          'dry deposition and settling of dust', 'kg/m2/s')

      id_dust_emis = register_diag_field (mod_name, &
          'dust_emis', axes(1:2), Time, &
          'total emission of dust', 'kg/m2/s')

      id_nh4_ddep_cmip = register_diag_field (mod_name, &
          'tot_nh4_ddep_cmip', axes(1:2), Time, &
          'total dry deposition of ammonium', 'kg/m2/s')

      id_nh4no3_col = register_diag_field (mod_name, &
          'tot_no3_col', axes(1:2), Time, &
          'total aerosol load of nitrate', 'kg/m2')

      id_nh4_col  = register_diag_field (mod_name, &
          'tot_nh4_col', axes(1:2), Time, &
          'total aerosol load of ammonium', 'kg/m2')

      id_nh4no3_cmip = register_diag_field (mod_name, &
          'tot_no3', axes(1:3), Time, &
          'total nitrate', 'kg/m3')

      id_nh4_cmip  = register_diag_field (mod_name, &
          'tot_nh4', axes(1:3), Time, &
          'total ammonium', 'kg/m3')

      id_nh4no3_cmipv2 = register_diag_field (mod_name, &
          'tot_no3v2', axes(1:3), Time, &
          'total nitrate', 'kg/m3')

      id_nh4_cmipv2  = register_diag_field (mod_name, &
          'tot_nh4v2', axes(1:3), Time, &
          'total ammonium', 'kg/m3')

      id_so2_cmip  = register_diag_field (mod_name, &
          'so2_cmip', axes(1:3), Time, &
          'SO2', 'kg/m3')

      id_dms_cmip  = register_diag_field (mod_name, &
          'dms_cmip', axes(1:3), Time, &
          'DMS', 'kg/m3')

      id_so2_cmipv2  = register_diag_field (mod_name, &
          'so2_cmipv2', axes(1:3), Time, &
          'SO2', 'kg/m3')

      id_dms_cmipv2  = register_diag_field (mod_name, &
          'dms_cmipv2', axes(1:3), Time, &
          'DMS', 'kg/m3')

      module_is_initialized = .TRUE.

 end subroutine atmos_tracer_driver_init



!#####################################################################

subroutine atmos_tracer_driver_time_vary (Time)

type(time_type), intent(in) :: Time

      if (nbcphobic > 0 .and. nbcphilic > 0 .and. &
          nomphobic > 0 .and. nomphilic > 0) then
        call atmos_carbon_aerosol_time_vary (Time)
      endif

      if (nch3i > 0) then
        call atmos_ch3i_time_vary (Time)
      endif

      if (ndust1 > 0.or.ndust2 > 0.or.ndust3 > 0 &
          .or.ndust4 > 0.or.ndust5 > 0 ) then
        call atmos_dust_time_vary (Time)
      endif

      if (nSOA > 0 ) then
        call atmos_SOA_time_vary (Time)
      endif

      if (nDMS > 0 .or. nSO2 > 0 .or. nSO4 > 0 &
                   .or. nMSA > 0 .or. nH2O2 > 0 ) then
        call atmos_sulfate_time_vary (Time)
      endif

      call dry_deposition_time_vary (drydep_data, Time)

      if (do_tropchem) then
        call tropchem_driver_time_vary (Time)
      endif



end subroutine atmos_tracer_driver_time_vary




!#####################################################################

subroutine atmos_tracer_driver_endts

      if (do_tropchem) then
        call tropchem_driver_endts
      endif
      if (nbcphobic > 0 .and. nbcphilic > 0 .and. &
          nomphobic > 0 .and. nomphilic > 0) then
        call atmos_carbon_aerosol_endts
      endif
      if (nch3i > 0) then
        call atmos_ch3i_endts
      endif
      if (ndust1 > 0.or.ndust2 > 0.or.ndust3 > 0 &
          .or.ndust4 > 0.or.ndust5 > 0 ) then
        call atmos_dust_endts
      endif
      if (nSOA > 0 ) then
        call atmos_soa_endts
      endif
      if (nDMS > 0 .or. nSO2 > 0 .or. nSO4 > 0 &
                   .or. nMSA > 0 .or. nH2O2 > 0 ) then
        call atmos_sulfate_endts
      endif
      call dry_deposition_endts (drydep_data)

end subroutine atmos_tracer_driver_endts




!#####################################################################


! </SUBROUTINE>

!#######################################################################

! <SUBROUTINE NAME="atmos_tracer_driver_end">
!   <OVERVIEW>
!     Subroutine to terminate the tracer driver module.
!   </OVERVIEW>
!   <DESCRIPTION>
!     Termination routine for tracer_driver. It should also call
!     the destructors for the individual tracer routines.
!   </DESCRIPTION>
!   <TEMPLATE>
!     call atmos_tracer_driver_end
!   </TEMPLATE>
 subroutine atmos_tracer_driver_end

!-----------------------------------------------------------------------
integer :: logunit
      if (mpp_pe() /= mpp_root_pe()) return

      logunit=stdlog()
      write (logunit,'(/,(a))') 'Exiting tracer_driver, have a nice day ...'

      call atmos_radon_end
      call atmos_sulfur_hex_end
      call atmos_convection_tracer_end
      call atmos_dust_end     
      call atmos_sea_salt_end 
      call atmos_sulfate_end
      call atmos_SOA_end      
      call atmos_carbon_aerosol_end
      if (nch3i > 0) then
        call atmos_ch3i_end
      endif
      call atmos_age_tracer_end      
      call atmos_co2_end

      module_is_initialized = .FALSE.

!-----------------------------------------------------------------------

 end subroutine atmos_tracer_driver_end
! </SUBROUTINE>

!#######################################################################
! <SUBROUTINE NAME="atmos_tracer_flux_init">
!   <OVERVIEW>
!     Subroutine to initialize the ocean-atmosphere gas flux modules
!   </OVERVIEW>
!   <DESCRIPTION>
!     Subroutine to initialize the ocean-atmosphere gas flux modules
!   </DESCRIPTION>

subroutine atmos_tracer_flux_init

call atmos_co2_flux_init

return

end subroutine atmos_tracer_flux_init
! </SUBROUTINE>

!######################################################################
! <SUBROUTINE NAME="atmos_tracer_driver_gather_data">
!   <OVERVIEW>
!     Subroutine to terminate the tracer driver module.
!   </OVERVIEW>
!   <DESCRIPTION>
!     Termination routine for tracer_driver. It should also call
!     the destructors for the individual tracer routines.
!   </DESCRIPTION>
!   <TEMPLATE>
!     call atmos_tracer_driver_gather_data
!   </TEMPLATE>
 subroutine atmos_tracer_driver_gather_data(gas_fields, tr_bot)

use coupler_types_mod, only: coupler_2d_bc_type

type(coupler_2d_bc_type), intent(inout) :: gas_fields
real, dimension(:,:,:), intent(in)      :: tr_bot

!-----------------------------------------------------------------------

  call atmos_co2_gather_data(gas_fields, tr_bot)

!-----------------------------------------------------------------------

 end subroutine atmos_tracer_driver_gather_data
! </SUBROUTINE>

!######################################################################


end module atmos_tracer_driver_mod


