!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

module ecosys_path_mod

!BOP
! !MODULE: path_mod
!
! !DESCRIPTION:Module for Pa/Th 
!  The activity of 231Pa and 230Th are treated as tracers.
!  Unit for Nd143 and Nd144 is dpm/m^3
!
! !REVISION HISTORY:
!  SVN:$Id: path_mod.F90 26603 2011-01-28 23:09:02Z njn01 $

! !USES:

   use POP_KindsMod
   use POP_IOUnitsMod
   use POP_ErrorMod

   use blocks, only: nx_block, ny_block, block
   use domain_size, only: max_blocks_clinic, km, nx_global, ny_global
   use domain, only: nblocks_clinic
   use exit_mod, only: sigAbort, exit_POP
   use communicate, only: my_task, master_task
   use prognostic, only: tracer_field
   use kinds_mod
   use constants, only: c0, c1000,c10000,c1, p5, char_blank, delim_fmt, field_type_scalar
   use io, only: data_set
   use io_types, only: stdout, nml_in, nml_filename, datafile, io_dim,  &
       io_field_desc, rec_type_dbl, construct_file, construct_io_dim,   &
       construct_io_field, destroy_file, destroy_io_field
   use io_tools, only: document
   use tavg, only: define_tavg_field, accumulate_tavg_field, accumulate_tavg_now
   use passive_tracer_tools, only: ind_name_pair, tracer_read, &
       rest_read_tracer_block, file_read_tracer_block, read_field
   
   use ecosys_share
   use ecosys_parms   
   
   implicit none
   private

! !PUBLIC MEMBER FUNCTIONS:

   public :: ecosys_path_tracer_cnt,        &
             ecosys_path_init,              &
             ecosys_path_set_interior,      &
             ecosys_path_tavg_forcing

!EOP
!BOC

!-----------------------------------------------------------------------
!  module variables required by passive_tracer
!-----------------------------------------------------------------------

   integer(int_kind), parameter :: &
       ecosys_path_tracer_cnt = 2

!-----------------------------------------------------------------------
!  relative tracer indices
!-----------------------------------------------------------------------

   integer(int_kind), parameter :: &
       ecosys_pa_ind = 1,      &     ! Pa index
       ecosys_th_ind = 2          ! Th index

!-----------------------------------------------------------------------
!  derived type & parameter for tracer index lookup
!-----------------------------------------------------------------------

   type(ind_name_pair), dimension(ecosys_path_tracer_cnt) :: &
       ind_name_table = (/ ind_name_pair(ecosys_pa_ind, 'eco_Pa'),ind_name_pair(ecosys_th_ind, 'eco_Th') /)

!-----------------------------------------------------------------------
!  tavg ids for non-standard tavg variables
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      
       tavg_eco_PA_D,                &
       tavg_eco_PA_P,                &
       tavg_eco_TH_D,                &
       tavg_eco_TH_P,                &
       tavg_eco_PA_RESET_TEND,        &
       tavg_eco_TH_RESET_TEND,		&
       tavg_eco_POC_flux_path,            &
       tavg_eco_OPAL_flux_path,           &
       tavg_eco_calcite_flux_path,        &
       tavg_eco_dust_flux_path,			  &
       tavg_betapa
	 
           
   real(r8),dimension(:,:,:),allocatable :: &
   	   betapa
   	   
   real(r8),dimension(:,:,:,:),allocatable :: &
    
       prodk_eco_pa,                            &  ! denominator in equation 14(rempfer2011)
       prodk_eco_th,							 &
       prodk_eco_pa_in,                            &  ! denominator in equation 14(rempfer2011)
       prodk_eco_th_in,							 &
       prodk_eco_pa_out,                            &  ! denominator in equation 14(rempfer2011)
       prodk_eco_th_out,							 &
       eco_pa_d,                               &
       eco_pa_p,                               &
       eco_th_d,                               &
       eco_th_p,                               &
       c_POC_path,								  & !POC concentration
   	   c_calcite_path,							  & !calcite concentration
   	   c_SiO2_path,								  & !SiO2 concentration
   	   c_dust_path,								  & !dust concentration
   	   eco_poc_flux_path,   &
	   eco_opal_flux_path,   &
	   eco_calcite_flux_path,   &
	   eco_dust_flux_path,     &
	   eco_poc_flux_in_path,		&
	   eco_poc_flux_out_path,	&
	   eco_opal_flux_in_path,	&
	   eco_opal_flux_out_path,   &
	   eco_calcite_flux_in_path,	&
	   eco_calcite_flux_out_path,&
	   eco_dust_flux_in_path,	&
	   eco_dust_flux_out_path


   real(r8)          :: &
       eco_path_w,			&!sinking velocity for particles
       eco_beta_pa,			&
       eco_beta_th,			&
       eco_lambda_pa,		&
       eco_lambda_th,		&
       k_eco_poc_pa,                    &
       k_eco_calcite_pa,                &
       k_eco_opal_pa,                   &
       k_eco_dust_pa,					&
       k_eco_poc_th,                    &
       k_eco_calcite_th,                &
       k_eco_opal_th,                   &
       k_eco_dust_th 


!EOC
!*****************************************************************************
contains

!*****************************************************************************
!BOP
! !IROUTINE: ecosys_path_init
! !INTERFACE:
subroutine ecosys_path_init(init_ts_file_fmt, read_restart_filename, &
                      tracer_d_module, TRACER_MODULE, errorCode)

! !DESCRIPTION:
!  Initialize pa_th tracer module. This involves setting metadata, reading
!  the module namelist and setting initial conditions.
!
! !REVISION HISTORY:
!  same as module

! !USES:

   use broadcast, only: broadcast_scalar
   use prognostic, only: curtime, oldtime
   use grid, only: KMT, n_topo_smooth, fill_points, dz
   use time_management, only: seconds_in_year

! !INPUT PARAMETERS:

   character (*), intent(in) ::  &
       init_ts_file_fmt,       &   ! format (bin or nc) for input file
       read_restart_filename       ! file name for restart file

    
! !INPUT/OUTPUT PARAMETERS:

   type (tracer_field), dimension(ecosys_path_tracer_cnt), intent(inout) :: &
       tracer_d_module   ! descriptors for each tracer

   real(r8), dimension(nx_block,ny_block,km,ecosys_path_tracer_cnt,3,max_blocks_clinic), &
       intent(inout) :: TRACER_MODULE

! !OUTPUT PARAMETERS:

   integer (POP_i4), intent(out) :: &
       errorCode            ! returned error code

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   character(*), parameter :: subname = 'ecosys_path_mod:ecosys_path_init'

   character(char_len) :: &
       init_ecosys_path_option,           & ! option for initialization of nd
       init_ecosys_path_init_file,        & ! filename for option 'file'
       init_ecosys_path_init_file_fmt       ! file format for option 'file'
  




  
   logical(log_kind) :: &
       lnml_found             ! Was path_nml found ?

   integer(int_kind) :: &
       n,                   & ! index for looping over trlsacers
       k,                   & ! index for looping over depth levels
       nx,                  & ! index for looping over x direction
       ny,                  & ! index for looping over y direction
       iblock,              & ! index for looping over blocks
       nml_error              ! namelist i/o error flag

!     l,                   & ! index for looping over time levels
   type(tracer_read), dimension(ecosys_path_tracer_cnt) :: &
       ecosys_path_init_ext        ! namelist variable for initializing tracers

   namelist /ecosys_path_nml/ &
       init_ecosys_path_option, init_ecosys_path_init_file, ecosys_path_init_ext, &
       init_ecosys_path_init_file_fmt, &
       k_eco_poc_pa,k_eco_calcite_pa,k_eco_opal_pa,k_eco_dust_pa,         &
       k_eco_poc_th,k_eco_calcite_th,k_eco_opal_th,k_eco_dust_th,eco_path_w,eco_beta_pa,eco_beta_th,eco_lambda_pa,eco_lambda_th


   character (char_len) ::  &
       eco_path_restart_filename  ! modified file name for restart file
      
  
!-----------------------------------------------------------------------
!  initialize tracer_d values
!-----------------------------------------------------------------------

   errorCode = POP_Success
  
   tracer_d_module(ecosys_pa_ind)%short_name = 'eco_Pa'
   tracer_d_module(ecosys_pa_ind)%long_name  = 'eco_Pa'
   tracer_d_module(ecosys_pa_ind)%units      = 'dpm/(m^3)'
   tracer_d_module(ecosys_pa_ind)%tend_units = 'dpm/(m^3s)'
!   tracer_d_module(pa_ind)%flux_units = 'cm years/s'


   tracer_d_module(ecosys_th_ind)%short_name = 'eco_Th'
   tracer_d_module(ecosys_th_ind)%long_name  = 'eco_Th'
   tracer_d_module(ecosys_th_ind)%units      = 'dpm/(m^3)'
   tracer_d_module(ecosys_th_ind)%tend_units = 'dpm/(m^3s)'
!  tracer_d_module(th_ind)%flux_units = 'cm years/s'
!-----------------------------------------------------------------------
!  default namelist settings
!-----------------------------------------------------------------------

   init_ecosys_path_option = 'unknown'
   init_ecosys_path_init_file = 'unknown'
   init_ecosys_path_init_file_fmt = 'bin'

   do n = 1,ecosys_path_tracer_cnt
       ecosys_path_init_ext(n)%mod_varname  = 'unknown'
       ecosys_path_init_ext(n)%filename     = 'unknown'
       ecosys_path_init_ext(n)%file_varname = 'unknown'
       ecosys_path_init_ext(n)%scale_factor = c1
       ecosys_path_init_ext(n)%default_val  = c0
       ecosys_path_init_ext(n)%file_fmt     = 'nc'
   end do
   

   k_eco_poc_pa = c0
   k_eco_calcite_pa = c0
   k_eco_opal_pa = c0
   k_eco_dust_pa = c0
   
   k_eco_poc_th = c0
   k_eco_calcite_th = c0
   k_eco_opal_th = c0
   k_eco_dust_th = c0
   
   eco_path_w = c0
   eco_beta_pa = c0
   eco_beta_th = c0
   eco_lambda_pa = c0
   eco_lambda_th = c0
  
   if (my_task == master_task) then
       open (nml_in, file=nml_filename, status='old',iostat=nml_error)
       if (nml_error /= 0) then  
           nml_error = -1
           call document(subname,'nml_error = -1')
       else
           nml_error =  1  
           call document(subname,'nml_error = 1')    
       endif
       do while (nml_error > 0)
           read(nml_in, nml=ecosys_path_nml,iostat=nml_error)
       end do
       if (nml_error == 0) close(nml_in)
   endif

   call broadcast_scalar(nml_error, master_task)
   if (nml_error /= 0) then
       call document(subname, 'ecosys_path_nml not found')
       call exit_POP(sigAbort, 'stopping in path_nml' /&
                           &/ subname)
   endif

!-----------------------------------------------------------------------
!  broadcast all namelist variables
!-----------------------------------------------------------------------

   call broadcast_scalar(init_ecosys_path_option , master_task)
   call broadcast_scalar(init_ecosys_path_init_file, master_task)
   call broadcast_scalar(init_ecosys_path_init_file_fmt, master_task)

   do n = 1,ecosys_path_tracer_cnt
      call broadcast_scalar(ecosys_path_init_ext(n)%mod_varname, master_task)
      call broadcast_scalar(ecosys_path_init_ext(n)%filename, master_task)
      call broadcast_scalar(ecosys_path_init_ext(n)%file_varname, master_task)
      call broadcast_scalar(ecosys_path_init_ext(n)%scale_factor, master_task)
      call broadcast_scalar(ecosys_path_init_ext(n)%default_val, master_task)
      call broadcast_scalar(ecosys_path_init_ext(n)%file_fmt, master_task)
   end do

  
   call broadcast_scalar(k_eco_poc_pa, master_task)
   call broadcast_scalar(k_eco_calcite_pa, master_task)
   call broadcast_scalar(k_eco_opal_pa, master_task)
   call broadcast_scalar(k_eco_dust_pa, master_task)
   call broadcast_scalar(k_eco_poc_th, master_task)
   call broadcast_scalar(k_eco_calcite_th, master_task)
   call broadcast_scalar(k_eco_opal_th, master_task)
   call broadcast_scalar(k_eco_dust_th, master_task)
   
   
   call broadcast_scalar(eco_path_w, master_task)
   call broadcast_scalar(eco_beta_pa,master_task)
   call broadcast_scalar(eco_beta_th,master_task)
   call broadcast_scalar(eco_lambda_pa,master_task)
   call broadcast_scalar(eco_lambda_th,master_task)
   
!-----------------------------------------------------------------------
!  initialize tracers
!-----------------------------------------------------------------------

   select case (init_ecosys_path_option)

   case ('ccsm_startup', 'zero', 'ccsm_startup_spunup')
      TRACER_MODULE = c0
      if (my_task == master_task) then
          write(stdout,delim_fmt)
          write(stdout,*) ' Initial 3-d  isotope ratio set to all zeros' 
          write(stdout,delim_fmt)
          call POP_IOUnitsFlush(POP_stdout) ; call POP_IOUnitsFlush(stdout)
      endif
       
   case ('restart', 'ccsm_continue', 'ccsm_branch', 'ccsm_hybrid' )
      eco_path_restart_filename = char_blank
      if (init_ecosys_path_init_file == 'same_as_TS') then
        if (read_restart_filename == 'undefined') then
            call document(subname, 'no restart file to read eco PaTh from')
            call exit_POP(sigAbort, 'stopping in case 11111' /&
                                 &/ subname)
        endif
        eco_path_restart_filename = read_restart_filename
        init_ecosys_path_init_file_fmt = init_ts_file_fmt

      else  ! do not read from TS restart file
        eco_path_restart_filename = trim(init_ecosys_path_init_file)

      endif
      call rest_read_tracer_block(init_ecosys_path_init_file_fmt, &
                                  eco_path_restart_filename,   &
                                  tracer_d_module,         &
                                  TRACER_MODULE)

   case ('file')
      call document(subname, 'eco PaTh being read from separate file')

      call file_read_tracer_block(init_ecosys_path_init_file_fmt, &
                                  init_ecosys_path_init_file,     &
                                  tracer_d_module,         &
                                  ind_name_table,          &
                                  ecosys_path_init_ext,         &
                                  TRACER_MODULE)
 
      if (n_topo_smooth > 0) then
        do n = 1,ecosys_path_tracer_cnt
         do k=1,km
            call fill_points(k,TRACER_MODULE(:,:,k,n,curtime,:), errorCode)

            if (errorCode /= POP_Success) then
               call POP_ErrorSet(errorCode, &
                  'ecosys_path_init: error in fill_points for ecoPaTh (curtime)')
               return
            endif
            call fill_points(k,TRACER_MODULE(:,:,k,n,oldtime,:), errorCode)

            if (errorCode /= POP_Success) then
               call POP_ErrorSet(errorCode, &
                  'ecosys_path_init: error in fill_points for eco PaTh (oldime)')
               return
            endif
            
          end do
         end do
         
      endif
    case default
      call document(subname, 'init_ecosys_path_option = ', init_ecosys_path_option)
      call exit_POP(sigAbort, 'stopping in  case default' /&
                           &/ subname)

   end select

!-----------------------------------------------------------------------
!  apply land mask to tracers
!-----------------------------------------------------------------------

   do iblock=1,nblocks_clinic
      do n = 1,ecosys_path_tracer_cnt
         do k = 1,km
            where (k > KMT(:,:,iblock))
                 TRACER_MODULE(:,:,k,n,curtime,iblock) = c0
                 TRACER_MODULE(:,:,k,n,oldtime,iblock) = c0
            end where
         end do
      end do
   enddo

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

   call define_tavg_field(tavg_eco_PA_RESET_TEND, 'eco_PA_RESET_TEND',2,  &
                          long_name='surface reset tendency of eco_PA', &
                          units='dpm/m3/s', grid_loc='2110',           &
                          coordinates='TLONG TLAT time')
   call define_tavg_field(tavg_eco_TH_RESET_TEND, 'eco_TH_RESET_TEND',2,  &
                          long_name='surface reset tendency of eco_TH', &
                          units='dpm/m3/s', grid_loc='2110',           &
                          coordinates='TLONG TLAT time')
   
    
   call define_tavg_field(tavg_eco_PA_D, 'eco_PA_D',3,&
                          long_name='eco_PA disolved',&
                          units='dpm/m^3', grid_loc='3111')
     
     
   call define_tavg_field(tavg_eco_PA_P, 'eco_PA_P',3,&
                          long_name='eco_PA particle related',&
                          units='dpm/m^3', grid_loc='3111')
     
     
   call define_tavg_field(tavg_eco_TH_D, 'eco_TH_D',3,&
                          long_name='eco_TH disolved',&
                          units='dpm/m^3', grid_loc='3111')
     
   call define_tavg_field(tavg_eco_TH_P, 'eco_TH_P',3,&
                          long_name='eco_TH particle related',&
                           units='dpm/m^3', grid_loc='3111')
     
     
     
   call define_tavg_field(tavg_eco_POC_flux_path, 'eco_POC_flux_path',3,				&
   						  long_name=' ',		&
   						  grid_loc='3111',coordinates='TLONG TLAT z_t time', units = 'base unit/cm^2/s')	
   						  						  
   call define_tavg_field(tavg_eco_OPAL_flux_path, 'eco_OPAL_flux_path',3,				&
   						  long_name='eco_OPAL_flux',		&
   						  grid_loc='3111',coordinates='TLONG TLAT z_t time', units = 'base unit/cm^2/s')	

   call define_tavg_field(tavg_eco_calcite_flux_path, 'eco_calcite_flux_path',3,				&
   						  long_name='eco_calcite_flux',		&
   						  grid_loc='3111',coordinates='TLONG TLAT z_t time', units = 'base unit/cm^2/s')	

   call define_tavg_field(tavg_eco_dust_flux_path, 'eco_dust_flux_path',3,				&
   						  long_name='eco_dust_flux',		&
   						  grid_loc='3111',coordinates='TLONG TLAT z_t time', units = 'base unit/cm^2/s')	
  
   						  						  

   	call define_tavg_field(tavg_betapa, 'beta_pa',2,                   &
   						  long_name='beta_pa',&
   						  grid_loc='2110')



  

!-----------------------------------------------------------------------
!  allocate space for pa_d,pa_p,th_d,th_p
!-----------------------------------------------------------------------

 allocate(c_POC_path(nx_block,ny_block,km,max_blocks_clinic))
   c_POC_path = c0

   allocate(c_calcite_path(nx_block,ny_block,km,max_blocks_clinic))
   c_calcite_path = c0

   allocate(c_SiO2_path(nx_block,ny_block,km,max_blocks_clinic))
   c_SiO2_path = c0

   allocate(c_dust_path(nx_block,ny_block,km,max_blocks_clinic))
   c_dust_path = c0
   
   
allocate(eco_pa_d(nx_block,ny_block,km,max_blocks_clinic))
allocate(eco_pa_p(nx_block,ny_block,km,max_blocks_clinic))
allocate(eco_th_d(nx_block,ny_block,km,max_blocks_clinic))
allocate(eco_th_p(nx_block,ny_block,km,max_blocks_clinic))

eco_pa_d = c0
eco_pa_p = c0
eco_th_d = c0
eco_th_p = c0


allocate(prodk_eco_pa(nx_block,ny_block,km,max_blocks_clinic))
prodk_eco_pa = c0
allocate(prodk_eco_th(nx_block,ny_block,km,max_blocks_clinic))
prodk_eco_th = c0

allocate(prodk_eco_pa_in(nx_block,ny_block,km,max_blocks_clinic))
prodk_eco_pa_in = c0
allocate(prodk_eco_th_in(nx_block,ny_block,km,max_blocks_clinic))
prodk_eco_th_in = c0

allocate(prodk_eco_pa_out(nx_block,ny_block,km,max_blocks_clinic))
prodk_eco_pa_out = c0
allocate(prodk_eco_th_out(nx_block,ny_block,km,max_blocks_clinic))
prodk_eco_th_out = c0


allocate(eco_poc_flux_path(nx_block,ny_block,km,max_blocks_clinic))
allocate(eco_opal_flux_path(nx_block,ny_block,km,max_blocks_clinic))
allocate(eco_calcite_flux_path(nx_block,ny_block,km,max_blocks_clinic))
allocate(eco_dust_flux_path(nx_block,ny_block,km,max_blocks_clinic))


allocate(eco_poc_flux_in_path(nx_block,ny_block,km,max_blocks_clinic))
allocate(eco_opal_flux_in_path(nx_block,ny_block,km,max_blocks_clinic))
allocate(eco_calcite_flux_in_path(nx_block,ny_block,km,max_blocks_clinic))
allocate(eco_dust_flux_in_path(nx_block,ny_block,km,max_blocks_clinic))


allocate(eco_poc_flux_out_path(nx_block,ny_block,km,max_blocks_clinic))
allocate(eco_opal_flux_out_path(nx_block,ny_block,km,max_blocks_clinic))
allocate(eco_calcite_flux_out_path(nx_block,ny_block,km,max_blocks_clinic))
allocate(eco_dust_flux_out_path(nx_block,ny_block,km,max_blocks_clinic))


eco_poc_flux_path = c0
eco_opal_flux_path = c0
eco_calcite_flux_path = c0
eco_dust_flux_path = c0


eco_poc_flux_in_path = c0
eco_opal_flux_in_path = c0
eco_calcite_flux_in_path = c0
eco_dust_flux_in_path = c0


eco_poc_flux_out_path = c0
eco_opal_flux_out_path = c0
eco_calcite_flux_out_path = c0
eco_dust_flux_out_path = c0


allocate(betapa(nx_block,ny_block,max_blocks_clinic))
betapa = c0
!EOC

 end subroutine ecosys_path_init
 
 !***********************************************************************
!BOP
! !IROUTINE: ecosys_path_set_interior
! !INTERFACE:

 
subroutine ecosys_path_set_interior(k, TRACER_MODULE_OLD, TRACER_MODULE_CUR, &
                           DTRACER_MODULE,this_block)

! !DESCRIPTION:
!  set interior source/sink term for pa_th 
!
! !REVISION HISTORY:
!  same as module

! !USES:

   use time_management, only: seconds_in_year
   use grid, only : dz,KMT,REGION_MASK,dzwr
   

! !INPUT PARAMETERS:

   integer(int_kind), intent(in) :: &
      k                   ! vertical level index
      
   real (r8), dimension(:,:,:,:), intent(in) :: &
      TRACER_MODULE_OLD, &! old tracer values
      TRACER_MODULE_CUR   ! current tracer values
      
   type (block), intent(in) :: &
      this_block          ! block info for the current block

! !OUTPUT PARAMETERS:

   real(r8), dimension(nx_block,ny_block,ecosys_path_tracer_cnt), intent(out) :: &
      DTRACER_MODULE      ! computed source/sink term
    
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------
   
      
   integer (int_kind)       :: &
      bid  ,        &               ! local_block id
      nx,           &
      ny,           &
      ntracer
   
   logical(log_kind),dimension(nx_block,ny_block) :: &
      mask
  
   real(r8), dimension(nx_block,ny_block) :: &
      p_pa_top,            &
      p_pa_bot,            &
      p_th_top,            &
      p_th_bot,            &
      pa_cur,           &
      th_cur,           &
      pa_upper,         &
      th_upper,         &
      pa_lower,         &
      th_lower
   
!EOP
!BOC

    bid = this_block%local_id
    

    betapa = eco_beta_pa
!-----------------------------------------------------------------------
    DTRACER_MODULE = c0


!-----------------------------------------------------------------------
!  apply production of new Pa or Th from the decay of 235U and 234U
!-----------------------------------------------------------------------
    DTRACER_MODULE(:,:,ecosys_pa_ind) = DTRACER_MODULE(:,:,ecosys_pa_ind)+eco_beta_pa/seconds_in_year
    DTRACER_MODULE(:,:,ecosys_th_ind) = DTRACER_MODULE(:,:,ecosys_th_ind)+eco_beta_th/seconds_in_year
    
    
!-----------------------------------------------------------------------
!  apply radiodecay for Pa and Th
!-----------------------------------------------------------------------    
    pa_cur = p5*(TRACER_MODULE_OLD(:,:,k,ecosys_pa_ind)+TRACER_MODULE_CUR(:,:,k,ecosys_pa_ind))
    th_cur = p5*(TRACER_MODULE_OLD(:,:,k,ecosys_th_ind)+TRACER_MODULE_CUR(:,:,k,ecosys_th_ind))
    DTRACER_MODULE(:,:,ecosys_pa_ind) = DTRACER_MODULE(:,:,ecosys_pa_ind) - eco_lambda_pa*pa_cur/seconds_in_year
    DTRACER_MODULE(:,:,ecosys_th_ind) = DTRACER_MODULE(:,:,ecosys_th_ind) - eco_lambda_th*th_cur/seconds_in_year
    
    
!-----------------------------------------------------------------------
!  sink term
!-----------------------------------------------------------------------  

!  calculate pa_d,pa_p,th_d,th_p and p_pa_top, p_pa_bot,p_th_top,p_th_bot
    p_pa_top = c0
    p_pa_bot = c0
    p_th_top = c0
    p_th_bot = c0
  
  
  where((KMT(:,:,bid))>k .and. k.eq.1)
    eco_poc_flux_path(:,:,k,bid) = (POC%sflux_out(:,:,bid) + POC%hflux_out(:,:,bid))
    eco_opal_flux_path(:,:,k,bid) = (P_SiO2%sflux_out(:,:,bid) + P_SiO2%hflux_out(:,:,bid))
    eco_calcite_flux_path(:,:,k,bid) = (P_CaCO3%sflux_out(:,:,bid) + P_CaCO3%hflux_out(:,:,bid))
    eco_dust_flux_path(:,:,k,bid) = p5*(dust%sflux_in(:,:,bid) + dust%hflux_in(:,:,bid) + dust%sflux_out(:,:,bid) + dust%hflux_out(:,:,bid))
  end where
 
  where((KMT(:,:,bid))>k .and. k.gt.1)
    eco_poc_flux_path(:,:,k,bid) = p5*(POC%sflux_in(:,:,bid) + POC%hflux_in(:,:,bid) + POC%sflux_out(:,:,bid) + POC%hflux_out(:,:,bid))
    eco_opal_flux_path(:,:,k,bid) = p5*(P_SiO2%sflux_in(:,:,bid) + P_SiO2%hflux_in(:,:,bid) + P_SiO2%sflux_out(:,:,bid) + P_SiO2%hflux_out(:,:,bid))
    eco_calcite_flux_path(:,:,k,bid) = p5*(P_CaCO3%sflux_in(:,:,bid) + P_CaCO3%hflux_in(:,:,bid) + P_CaCO3%sflux_out(:,:,bid) + P_CaCO3%hflux_out(:,:,bid))
    eco_dust_flux_path(:,:,k,bid) = p5*(dust%sflux_in(:,:,bid) + dust%hflux_in(:,:,bid) + dust%sflux_out(:,:,bid) + dust%hflux_out(:,:,bid))
  end where
   
   where((KMT(:,:,bid)).eq.k)
   eco_poc_flux_path(:,:,k,bid) = POC%sflux_in(:,:,bid) + POC%hflux_in(:,:,bid)
    eco_opal_flux_path(:,:,k,bid) = P_SiO2%sflux_in(:,:,bid) + P_SiO2%hflux_in(:,:,bid)
    eco_calcite_flux_path(:,:,k,bid) = P_CaCO3%sflux_in(:,:,bid) + P_CaCO3%hflux_in(:,:,bid) 
    eco_dust_flux_path(:,:,k,bid) = dust%sflux_in(:,:,bid) + dust%hflux_in(:,:,bid)
   end where
   
    eco_poc_flux_in_path(:,:,k,bid) = (POC%sflux_in(:,:,bid) + POC%hflux_in(:,:,bid))
    eco_opal_flux_in_path(:,:,k,bid) = (P_SiO2%sflux_in(:,:,bid) + P_SiO2%hflux_in(:,:,bid))
    eco_calcite_flux_in_path(:,:,k,bid) = (P_CaCO3%sflux_in(:,:,bid) + P_CaCO3%hflux_in(:,:,bid))
    eco_dust_flux_in_path(:,:,k,bid) = (dust%sflux_in(:,:,bid) + dust%hflux_in(:,:,bid))
    
    
    eco_poc_flux_out_path(:,:,k,bid) = (POC%sflux_out(:,:,bid) + POC%hflux_out(:,:,bid))
    eco_opal_flux_out_path(:,:,k,bid) = (P_SiO2%sflux_out(:,:,bid) + P_SiO2%hflux_out(:,:,bid))
    eco_calcite_flux_out_path(:,:,k,bid) = (P_CaCO3%sflux_out(:,:,bid) + P_CaCO3%hflux_out(:,:,bid))
    eco_dust_flux_out_path(:,:,k,bid) = (dust%sflux_out(:,:,bid) + dust%hflux_out(:,:,bid))
   
    c_POC_path(:,:,k,bid) =  eco_poc_flux_path(:,:,k,bid)*POC%mass*1.0e-3_r8/eco_path_w
    c_calcite_path(:,:,k,bid) = eco_calcite_flux_path(:,:,k,bid)*P_CaCO3%mass*1.0e-3_r8/eco_path_w
    c_SiO2_path(:,:,k,bid) =  eco_opal_flux_path(:,:,k,bid)*P_SiO2%mass*1.0e-3_r8/eco_path_w
    c_dust_path(:,:,k,bid) = eco_dust_flux_path(:,:,k,bid)*1.0e6_r8/eco_path_w
    
    
    prodk_eco_pa(:,:,k,bid) = c1 + k_eco_poc_pa*c_POC_path(:,:,k,bid)/(1024.5_r8*c1000)+k_eco_calcite_pa*c_calcite_path(:,:,k,bid)/(1024.5_r8*c1000) + &
						    k_eco_opal_pa*c_SiO2_path(:,:,k,bid)/(1024.5_r8*c1000) + k_eco_dust_pa*c_dust_path(:,:,k,bid)/(1024.5_r8*c1000)
	

    prodk_eco_th(:,:,k,bid) = c1 + k_eco_poc_th*c_POC_path(:,:,k,bid)/(1024.5_r8*c1000)+k_eco_calcite_th*c_calcite_path(:,:,k,bid)/(1024.5_r8*c1000) + &
						    k_eco_opal_th*c_SiO2_path(:,:,k,bid)/(1024.5_r8*c1000) + k_eco_dust_th*c_dust_path(:,:,k,bid)/(1024.5_r8*c1000)
  
  
    prodk_eco_pa_in(:,:,k,bid) = c1 +   k_eco_poc_pa*(eco_poc_flux_in_path(:,:,k,bid)*POC%mass*1.0e-3_r8/eco_path_w)/(1024.5_r8*c1000)+ &
    									k_eco_calcite_pa*(eco_calcite_flux_in_path(:,:,k,bid)*P_CaCO3%mass*1.0e-3_r8/eco_path_w)/(1024.5_r8*c1000) + &
						    			k_eco_opal_pa*(eco_opal_flux_in_path(:,:,k,bid)*P_SiO2%mass*1.0e-3_r8/eco_path_w)/(1024.5_r8*c1000) + &
						    			k_eco_dust_pa*(eco_dust_flux_in_path(:,:,k,bid)*1.0e6_r8/eco_path_w)/(1024.5_r8*c1000)
  
  
    prodk_eco_pa_out(:,:,k,bid) = c1 +   k_eco_poc_pa*(eco_poc_flux_out_path(:,:,k,bid)*POC%mass*1.0e-3_r8/eco_path_w)/(1024.5_r8*c1000)+ &
    									k_eco_calcite_pa*(eco_calcite_flux_out_path(:,:,k,bid)*P_CaCO3%mass*1.0e-3_r8/eco_path_w)/(1024.5_r8*c1000) + &
						    			k_eco_opal_pa*(eco_opal_flux_out_path(:,:,k,bid)*P_SiO2%mass*1.0e-3_r8/eco_path_w)/(1024.5_r8*c1000) + &
						    			k_eco_dust_pa*(eco_dust_flux_out_path(:,:,k,bid)*1.0e6_r8/eco_path_w)/(1024.5_r8*c1000)
  
  
  
    prodk_eco_th_in(:,:,k,bid) = c1 +   k_eco_poc_th*(eco_poc_flux_in_path(:,:,k,bid)*POC%mass*1.0e-3_r8/eco_path_w)/(1024.5_r8*c1000)+ &
    									k_eco_calcite_th*(eco_calcite_flux_in_path(:,:,k,bid)*P_CaCO3%mass*1.0e-3_r8/eco_path_w)/(1024.5_r8*c1000) + &
						    			k_eco_opal_th*(eco_opal_flux_in_path(:,:,k,bid)*P_SiO2%mass*1.0e-3_r8/eco_path_w)/(1024.5_r8*c1000) + &
						    			k_eco_dust_th*(eco_dust_flux_in_path(:,:,k,bid)*1.0e6_r8/eco_path_w)/(1024.5_r8*c1000)
  
  
    prodk_eco_th_out(:,:,k,bid) = c1 +   k_eco_poc_th*(eco_poc_flux_out_path(:,:,k,bid)*POC%mass*1.0e-3_r8/eco_path_w)/(1024.5_r8*c1000)+ &
    									k_eco_calcite_th*(eco_calcite_flux_out_path(:,:,k,bid)*P_CaCO3%mass*1.0e-3_r8/eco_path_w)/(1024.5_r8*c1000) + &
						    			k_eco_opal_th*(eco_opal_flux_out_path(:,:,k,bid)*P_SiO2%mass*1.0e-3_r8/eco_path_w)/(1024.5_r8*c1000) + &
						    			k_eco_dust_th*(eco_dust_flux_out_path(:,:,k,bid)*1.0e6_r8/eco_path_w)/(1024.5_r8*c1000)
  
  

    where(KMT(:,:,bid)>=k)
        eco_pa_d(:,:,k,bid) = pa_cur/prodk_eco_pa(:,:,k,bid)
        eco_pa_p(:,:,k,bid) = pa_cur - eco_pa_d(:,:,k,bid)
    
        eco_th_d(:,:,k,bid) = th_cur/prodk_eco_th(:,:,k,bid)
        eco_th_p(:,:,k,bid) = th_cur - eco_th_d(:,:,k,bid)
    end where  
    
    
    where(eco_pa_d(:,:,k,bid).lt.c0)
        eco_pa_d(:,:,k,bid) = c0
    end where
    
    where(eco_pa_p(:,:,k,bid).lt.c0)
        eco_pa_p(:,:,k,bid) = c0
    end where
    
    where(eco_th_d(:,:,k,bid).lt.c0)
        eco_th_d(:,:,k,bid) = c0
    end where
    
    where(eco_th_p(:,:,k,bid).lt.c0)
        eco_th_p(:,:,k,bid) = c0
    end where
                  
    if(k==1) then

        
        pa_lower =  p5*(TRACER_MODULE_OLD(:,:,k+1,ecosys_pa_ind)+ TRACER_MODULE_CUR(:,:,k+1,ecosys_pa_ind))
        th_lower =  p5*(TRACER_MODULE_OLD(:,:,k+1,ecosys_th_ind)+ TRACER_MODULE_CUR(:,:,k+1,ecosys_th_ind))              
        
	    p_pa_bot = (pa_cur + (pa_lower-pa_cur)*dzwr(k)*p5*dz(k))*(1-c1/prodk_eco_pa_out(:,:,k,bid))               
        p_th_bot = (th_cur + (th_lower-th_cur)*dzwr(k)*p5*dz(k))*(1-c1/prodk_eco_th_out(:,:,k,bid))
    
    else if(k >1) then
       where((KMT(:,:,bid))>k)
    
            pa_upper  = p5*(TRACER_MODULE_OLD(:,:,k-1,ecosys_pa_ind)+ TRACER_MODULE_CUR(:,:,k-1,ecosys_pa_ind))
            pa_lower =  p5*(TRACER_MODULE_OLD(:,:,k+1,ecosys_pa_ind)+ TRACER_MODULE_CUR(:,:,k+1,ecosys_pa_ind))
            p_pa_top = (pa_upper + (pa_cur-pa_upper)*dzwr(k-1)*p5*dz(k-1))*(1-c1/prodk_eco_pa_in(:,:,k,bid))
            p_pa_bot = (pa_cur + (pa_lower-pa_cur)*dzwr(k)*p5*dz(k))*(1-c1/prodk_eco_pa_out(:,:,k,bid))                       
   
            th_upper  = p5*(TRACER_MODULE_OLD(:,:,k-1,ecosys_th_ind)+ TRACER_MODULE_CUR(:,:,k-1,ecosys_th_ind))
            th_lower =  p5*(TRACER_MODULE_OLD(:,:,k+1,ecosys_th_ind)+ TRACER_MODULE_CUR(:,:,k+1,ecosys_th_ind))
     	    p_th_top = (th_upper + (th_cur-th_upper)*dzwr(k-1)*p5*dz(k-1))*(1-c1/prodk_eco_th_in(:,:,k,bid))
            p_th_bot = (th_cur + (th_lower-th_cur)*dzwr(k)*p5*dz(k))*(1-c1/prodk_eco_th_out(:,:,k,bid))          
            

       end where
       
       
       where((KMT(:,:,bid)).eq.k)
            pa_upper  = p5*(TRACER_MODULE_OLD(:,:,k-1,ecosys_pa_ind)+ TRACER_MODULE_CUR(:,:,k-1,ecosys_pa_ind))
            p_pa_top = (pa_upper + (pa_cur-pa_upper)*dzwr(k-1)*p5*dz(k-1))*(1-c1/prodk_eco_pa_in(:,:,k,bid))
			p_pa_bot = eco_pa_p(:,:,k,bid)
    
            th_upper  = p5*(TRACER_MODULE_OLD(:,:,k-1,ecosys_th_ind)+ TRACER_MODULE_CUR(:,:,k-1,ecosys_th_ind))
            p_th_top = (th_upper + (th_cur-th_upper)*dzwr(k-1)*p5*dz(k-1))*(1-c1/prodk_eco_th_in(:,:,k,bid))
			p_th_bot = eco_th_p(:,:,k,bid)
    
       end where
    

    end if
    
    
    DTRACER_MODULE(:,:,ecosys_pa_ind) = DTRACER_MODULE(:,:,ecosys_pa_ind) - eco_path_w*(p_pa_bot - p_pa_top)/dz(k)
    DTRACER_MODULE(:,:,ecosys_th_ind) = DTRACER_MODULE(:,:,ecosys_th_ind) - eco_path_w*(p_th_bot - p_th_top)/dz(k)
 
    where(REGION_MASK(:,:,bid).eq.-14 .or. REGION_MASK(:,:,bid).eq.-13 .or. REGION_MASK(:,:,bid).eq.-12 .or.    &
           REGION_MASK(:,:,bid).eq.-5 .or. REGION_MASK(:,:,bid).eq.4 .or. REGION_MASK(:,:,bid).eq.11) 
              DTRACER_MODULE(:,:,ecosys_pa_ind) = c0
              DTRACER_MODULE(:,:,ecosys_th_ind) = c0
    end where
      
    
!-----------------------------------------------------------------------
!EOC

 end subroutine ecosys_path_set_interior



!***********************************************************************

!***********************************************************************
!BOP
! !IROUTINE: ecosys_path_tavg_forcing
! !INTERFACE:

 subroutine ecosys_path_tavg_forcing
 


! !DESCRIPTION:
!  Make accumulation calls for forcing related tavg fields. This is
!  necessary because the forcing routines are called before tavg flags
!  are set.

! !REVISION HISTORY:
!  same as module

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   integer (int_kind) :: &
      iblock,   &      ! block loop index
      k

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

   !$OMP PARALLEL DO PRIVATE(iblock)

   do iblock = 1, nblocks_clinic
   
       call accumulate_tavg_field(betapa(:,:,iblock),tavg_betapa,iblock,1)
       
       do k = 1,km
          call accumulate_tavg_field(eco_pa_d(:,:,k,iblock),tavg_eco_PA_D,iblock,k)
          call accumulate_tavg_field(eco_pa_p(:,:,k,iblock),tavg_eco_PA_P,iblock,k)
          call accumulate_tavg_field(eco_th_d(:,:,k,iblock),tavg_eco_TH_D,iblock,k)
          call accumulate_tavg_field(eco_th_p(:,:,k,iblock),tavg_eco_TH_P,iblock,k)
          call accumulate_tavg_field(eco_poc_flux_path(:,:,k,iblock),tavg_eco_POC_flux_path,iblock,k)
    	  call accumulate_tavg_field(eco_opal_flux_path(:,:,k,iblock),tavg_eco_OPAL_flux_path,iblock,k)
    	  call accumulate_tavg_field(eco_calcite_flux_path(:,:,k,iblock),tavg_eco_calcite_flux_path,iblock,k)
    	  call accumulate_tavg_field(eco_dust_flux_path(:,:,k,iblock),tavg_eco_dust_flux_path,iblock,k)
         
    	 

      end do  
   
   end do

   !$OMP END PARALLEL DO

!-----------------------------------------------------------------------
!EOC

 end subroutine ecosys_path_tavg_forcing


!***********************************************************************

!***********************************************************************
!BOP
! !IROUTINE: read_field_3D
! !INTERFACE:

 subroutine read_field_3D(fmt, filename, fieldname, FIELD, record_length)

! !DESCRIPTION:
!  read 3D field from a file
!  Assumes the field is (nx_global,ny_global), cell centered, and scalar.
!  The length of the 3rd dimension is determined by the dimension of FIELD.
!  For binary files, the default external precision is double precision.
!  This can be overridden by passing the desired precision into record_length.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

   character (*), intent(in) ::  &
      fmt,                 & ! format (bin or nc)
      filename,            & ! file to read from
      fieldname              ! field to be read

   integer(int_kind), intent(in), optional :: &
      record_length          ! record length type for binary files

! !INPUT/OUTPUT PARAMETERS:

   real(r8), dimension(:,:,:,:), intent(inout), target :: &
      FIELD                  ! field to be read in

!EOP
!BOC
!-----------------------------------------------------------------------
!  local variables
!-----------------------------------------------------------------------

   character(*), parameter :: &
      subname = 'passive_tracer_tools:read_field_3D'

   integer(int_kind) :: &
      record_length_loc    ! record length type for binary files

   type (io_field_desc) :: &
      FIELD_DESC           ! IO field descriptors for FIELD

   type (datafile) :: &
      restart_file         ! io file descriptor

   type (io_dim) :: &
      i_dim, j_dim, k_dim  ! dimension descriptors

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

   call document(subname, 'reading ' /&
                       &/ trim(fieldname) /&
                       &/ ' from ' /&
                       &/ trim(filename))

   if (present(record_length)) then
      record_length_loc = record_length
   else
      record_length_loc = rec_type_dbl
   endif

   restart_file =                                     &
      construct_file(fmt,                             &
                     full_name=trim(filename),        &
                     record_length=record_length_loc, &
                     recl_words=nx_global*ny_global)

   call data_set(restart_file, 'open_read')

   i_dim = construct_io_dim('i', nx_global)
   j_dim = construct_io_dim('j', ny_global)
   k_dim = construct_io_dim('k', size(FIELD,3))

   FIELD_DESC =                                       &
      construct_io_field(trim(fieldname),             &
                         dim1=i_dim,                  &
                         dim2=j_dim,                  &
                         dim3=k_dim,                  &
                         grid_loc ='3111',            &
                         d3d_array = FIELD)

   call data_set (restart_file, 'define', FIELD_DESC)

   call data_set (restart_file, 'read', FIELD_DESC)

   call destroy_io_field (FIELD_DESC)

   call data_set (restart_file, 'close')

   call destroy_file (restart_file)

!-----------------------------------------------------------------------
!EOC

 end subroutine read_field_3D


end module ecosys_path_mod

!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
