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

module 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, 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
   
   
   
   implicit none
   private

! !PUBLIC MEMBER FUNCTIONS:

   public :: path_tracer_cnt,        &
             path_init,              &
             path_set_interior,      &
             path_tavg_forcing

!EOP
!BOC

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

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

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

   integer(int_kind), parameter :: &
       pa_ind = 1,      &     ! Pa index
       th_ind = 2          ! Th index

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

   type(ind_name_pair), dimension(path_tracer_cnt) :: &
       ind_name_table = (/ ind_name_pair(pa_ind, 'Pa'),ind_name_pair(th_ind, 'Th') /)

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

   integer (int_kind) :: &
      
       tavg_PA_D,                &
       tavg_PA_P,                &
       tavg_TH_D,                &
       tavg_TH_P,                &
       tavg_PA_SOURCE,           &
       tavg_PA_SINK,             &
       tavg_TH_SOURCE,           &
       tavg_TH_SINK,             &
       tavg_PA_RESET_TEND,        &
       tavg_TH_RESET_TEND,       &
       tavg_PA_TOP,			&
       tavg_PA_BOT,			&
       tavg_TH_TOP,			&
       tavg_TH_BOT
      
     
   real(r8),dimension(:,:,:,:),allocatable :: &
       r_path_calcite,                             & ! calcite mass concentration ratio
       r_path_poc,                                 &  ! poc mass concentration ratio	
       r_path_opal,                                &  ! opal mass concentration ratio
       r_path_dust,                                & ! dust mass concentration ratio
       prodk_cndk_pa,                            &  ! denominator in equation 14(rempfer2011)
       prodk_cndk_th,							 &
       pa_d,                               &
       pa_p,                               &
       th_d,                               &
       th_p,                               &
       pa_source,                          &
       pa_sink,                            &
       th_source,                          &
       th_sink,							  &
	   pa_top,							  &
	   pa_bot,					   		  &
	   th_top,							  &
	   th_bot

   real(r8)          :: &
       path_w,			&!sinking velocity for particles
       beta_pa,			&
       beta_th,			&
       lambda_pa,		&
       lambda_th 


!EOC
!*****************************************************************************

contains

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

! !DESCRIPTION:
!  Initialize nd 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(path_tracer_cnt), intent(inout) :: &
       tracer_d_module   ! descriptors for each tracer

   real(r8), dimension(nx_block,ny_block,km,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 = 'path_mod:path_init'

   character(char_len) :: &
       init_path_option,           & ! option for initialization of nd
       init_path_init_file,        & ! filename for option 'file'
       init_path_init_file_fmt,    & ! file format for option 'file'
  
       prescribed_filename_path,       & ! fielname for dust deposition
       r_path_calcite_varname,        & ! variable name for calcite R
       r_path_poc_varname,            & ! variable name for poc R
       r_path_opal_varname,           & ! variable name for opal R
       r_path_dust_varname              ! variable name for dust R

   real(r8)          :: &
       k_poc_pa,                    &
       k_calcite_pa,                &
       k_opal_pa,                   &
       k_dust_pa,					&
       k_poc_th,                    &
       k_calcite_th,                &
       k_opal_th,                   &
       k_dust_th

  
   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(path_tracer_cnt) :: &
       path_init_ext        ! namelist variable for initializing tracers

   namelist /path_nml/ &
       init_path_option, init_path_init_file, path_init_ext, &
       init_path_init_file_fmt, prescribed_filename_path, &
       r_path_calcite_varname,r_path_poc_varname,       &
       r_path_opal_varname, r_path_dust_varname,        &
       k_poc_pa,k_calcite_pa,k_opal_pa,k_dust_pa,         &
       k_poc_th,k_calcite_th,k_opal_th,k_dust_th,path_w,beta_pa,beta_th,lambda_pa,lambda_th


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

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


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

   init_path_option = 'unknown'
   init_path_init_file = 'unknown'
   init_path_init_file_fmt = 'bin'

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

   r_path_calcite_varname = 'unknown'
   r_path_poc_varname =     'unknown'
   r_path_opal_varname  =   'unknown'
   r_path_dust_varname  =   'unknown'
  
   k_poc_pa = c0
   k_calcite_pa = c0
   k_opal_pa = c0
   k_dust_pa = c0
   
   k_poc_th = c0
   k_calcite_th = c0
   k_opal_th = c0
   k_dust_th = c0
   
   path_w = c0
   beta_pa = c0
   beta_th = c0
   lambda_pa = c0
   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
       else
           nml_error =  1      
       endif
       do while (nml_error > 0)
           read(nml_in, nml=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, 'path_nml not found')
       call exit_POP(sigAbort, 'stopping in ' /&
                           &/ subname)
   endif

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

   call broadcast_scalar(init_path_option , master_task)
   call broadcast_scalar(init_path_init_file, master_task)
   call broadcast_scalar(init_path_init_file_fmt, master_task)

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


   call broadcast_scalar(prescribed_filename_path, master_task)
   
  
   call broadcast_scalar(r_path_calcite_varname, master_task)
   call broadcast_scalar(r_path_poc_varname, master_task)
   call broadcast_scalar(r_path_opal_varname, master_task)
   call broadcast_scalar(r_path_dust_varname, master_task)
   
  
   call broadcast_scalar(k_poc_pa, master_task)
   call broadcast_scalar(k_calcite_pa, master_task)
   call broadcast_scalar(k_opal_pa, master_task)
   call broadcast_scalar(k_dust_pa, master_task)
   call broadcast_scalar(k_poc_th, master_task)
   call broadcast_scalar(k_calcite_th, master_task)
   call broadcast_scalar(k_opal_th, master_task)
   call broadcast_scalar(k_dust_th, master_task)
   
   
   call broadcast_scalar(path_w, master_task)
   call broadcast_scalar(beta_pa,master_task)
   call broadcast_scalar(beta_th,master_task)
   call broadcast_scalar(lambda_pa,master_task)
   call broadcast_scalar(lambda_th,master_task)
   
!-----------------------------------------------------------------------
!  initialize tracers
!-----------------------------------------------------------------------

   select case (init_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' )
      path_restart_filename = char_blank
      if (init_path_init_file == 'same_as_TS') then
        if (read_restart_filename == 'undefined') then
            call document(subname, 'no restart file to read PaTh from')
            call exit_POP(sigAbort, 'stopping in ' /&
                                 &/ subname)
        endif
        path_restart_filename = read_restart_filename
        init_path_init_file_fmt = init_ts_file_fmt

      else  ! do not read from TS restart file
        path_restart_filename = trim(init_path_init_file)

      endif
      call rest_read_tracer_block(init_path_init_file_fmt, &
                                  path_restart_filename,   &
                                  tracer_d_module,         &
                                  TRACER_MODULE)

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

      call file_read_tracer_block(init_path_init_file_fmt, &
                                  init_path_init_file,     &
                                  tracer_d_module,         &
                                  ind_name_table,          &
                                  path_init_ext,         &
                                  TRACER_MODULE)
 
      if (n_topo_smooth > 0) then
        do n = 1,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, &
                  'path_init: error in fill_points for PaTh')
               return
            endif
          end do
         end do
         
      endif
    case default
      call document(subname, 'init_path_option = ', init_path_option)
      call exit_POP(sigAbort, 'stopping in ' /&
                           &/ subname)

   end select

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

   do iblock=1,nblocks_clinic
      do n = 1,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_PA_RESET_TEND, 'PA_RESET_TEND',2,  &
                          long_name='surface reset tendency of PA', &
                          units='dpm/m3/s', grid_loc='2110',           &
                          coordinates='TLONG TLAT time')
   call define_tavg_field(tavg_TH_RESET_TEND, 'TH_RESET_TEND',2,  &
                          long_name='surface reset tendency of TH', &
                          units='dpm/m3/s', grid_loc='2110',           &
                          coordinates='TLONG TLAT time')
   
    
   call define_tavg_field(tavg_PA_D, 'PA_D',3,&
                          long_name='PA disolved',&
                          units='dpm/m^3', grid_loc='3111')
     
     
   call define_tavg_field(tavg_PA_P, 'PA_P',3,&
                          long_name='PA particle related',&
                          units='dpm/m^3', grid_loc='3111')
     
     
   call define_tavg_field(tavg_TH_D, 'TH_D',3,&
                          long_name='TH disolved',&
                          units='dpm/m^3', grid_loc='3111')
     
   call define_tavg_field(tavg_TH_P, 'TH_P',3,&
                          long_name='TH particle related',&
                           units='dpm/m^3', grid_loc='3111')
     
   call define_tavg_field(tavg_PA_SOURCE, 'PA_SOURCE',3,&
                          long_name='231Pa source term',&
                          grid_loc='3111',coordinates='TLONG TLAT z_t time')
  
   call define_tavg_field(tavg_PA_SINK, 'PA_SINK',3,&
                          long_name='231Pa sink term',&
                          grid_loc='3111',coordinates='TLONG TLAT z_t time')
             

   call define_tavg_field(tavg_TH_SOURCE, 'TH_SOURCE',3,&
                          long_name='230 Th source term',&
                          grid_loc='3111',coordinates='TLONG TLAT z_t time')
  
   call define_tavg_field(tavg_TH_SINK, 'TH_SINK',3,&
                          long_name='230 Th sink term',&
                          grid_loc='3111',coordinates='TLONG TLAT z_t time')
    
   						  
   
   call define_tavg_field(tavg_PA_TOP, 'PA_TOP',3,				&
   						  long_name='PA_TOP',		&
   						  grid_loc='3111',coordinates='TLONG TLAT z_t time')
   call define_tavg_field(tavg_PA_BOT, 'PA_BOT',3,				&
   						  long_name='PA_BOT',		&
   						  grid_loc='3111',coordinates='TLONG TLAT z_t time')
   call define_tavg_field(tavg_TH_TOP, 'TH_TOP',3,				&
   						  long_name='TH_TOP',		&
   						  grid_loc='3111',coordinates='TLONG TLAT z_t time')
   call define_tavg_field(tavg_TH_BOT, 'TH_BOT',3,				&
   						  long_name='TH_BOT',		&
   						  grid_loc='3111',coordinates='TLONG TLAT z_t time')			
   						  
   						  						  
   write(stdout,delim_fmt)
   write(stdout,*) ' successful define tavg_PATH' 
   write(stdout,delim_fmt)



!-----------------------------------------------------------------------
!  read in particle concentration ratio 
!-----------------------------------------------------------------------

   allocate(r_path_calcite(nx_block,ny_block,km,max_blocks_clinic))
   r_path_calcite = c0
   call read_field_3D('nc',prescribed_filename_path,r_path_calcite_varname,r_path_calcite)

   allocate(r_path_poc(nx_block,ny_block,km,max_blocks_clinic))
   r_path_poc = c0
   call read_field_3D('nc',prescribed_filename_path,r_path_poc_varname,r_path_poc)

   allocate(r_path_opal(nx_block,ny_block,km,max_blocks_clinic))
   r_path_opal = c0
   call read_field_3D('nc',prescribed_filename_path,r_path_opal_varname,r_path_opal)

   allocate(r_path_dust(nx_block,ny_block,km,max_blocks_clinic))
   r_path_dust = c0
   call read_field_3D('nc',prescribed_filename_path,r_path_dust_varname,r_path_dust)

!-----------------------------------------------------------------------
!  calculate prodk_cndk
!-----------------------------------------------------------------------
   allocate(prodk_cndk_pa(nx_block,ny_block,km,max_blocks_clinic))
   allocate(prodk_cndk_th(nx_block,ny_block,km,max_blocks_clinic))

   prodk_cndk_pa = c0
   prodk_cndk_th = c0	



   do iblock = 1,nblocks_clinic
     do nx = 1,nx_block
       do ny = 1,ny_block 
            do k = 1,km
                
                prodk_cndk_pa(nx,ny,k,iblock) = c1 + k_poc_pa*r_path_poc(nx,ny,k,iblock) +k_calcite_pa*r_path_calcite(nx,ny,k,iblock) + &
                                             k_opal_pa*r_path_opal(nx,ny,k,iblock)+ k_dust_pa*r_path_dust(nx,ny,k,iblock)
            
                prodk_cndk_th(nx,ny,k,iblock) = c1 + k_poc_th*r_path_poc(nx,ny,k,iblock) +k_calcite_th*r_path_calcite(nx,ny,k,iblock) + &
                                             k_opal_th*r_path_opal(nx,ny,k,iblock)+ k_dust_th*r_path_dust(nx,ny,k,iblock)
            
            end do
      end do
     end do
   end do

  

!-----------------------------------------------------------------------
!  allocate space for pa_d,pa_p,th_d,th_p
!-----------------------------------------------------------------------
allocate(pa_d(nx_block,ny_block,km,max_blocks_clinic))
allocate(pa_p(nx_block,ny_block,km,max_blocks_clinic))
allocate(th_d(nx_block,ny_block,km,max_blocks_clinic))
allocate(th_p(nx_block,ny_block,km,max_blocks_clinic))

allocate(pa_source(nx_block,ny_block,km,max_blocks_clinic))
allocate(pa_sink(nx_block,ny_block,km,max_blocks_clinic))
allocate(th_source(nx_block,ny_block,km,max_blocks_clinic))
allocate(th_sink(nx_block,ny_block,km,max_blocks_clinic))


allocate(pa_top(nx_block,ny_block,km,max_blocks_clinic))
allocate(pa_bot(nx_block,ny_block,km,max_blocks_clinic))
allocate(th_top(nx_block,ny_block,km,max_blocks_clinic))
allocate(th_bot(nx_block,ny_block,km,max_blocks_clinic))

pa_top = c0
th_top = c0
pa_bot = c0
th_bot = c0


pa_d = c0
pa_p = c0
th_d = c0
th_p = c0

pa_source = c0
pa_sink = c0
th_source = c0
th_sink = c0

!EOC

 end subroutine path_init

!***********************************************************************
!BOP
! !IROUTINE: path_set_interior
! !INTERFACE:

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

! !DESCRIPTION:
!  set interior source/sink term for nd isotope ratio tracer
!
! !REVISION HISTORY:
!  same as module

! !USES:

   use time_management, only: seconds_in_year
   use grid, only : dz,KMT,REGION_MASK,dzwr
   
   use forcing_fields, only : ROFF_F
! !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,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_k_upper,        &
      p_pa_k,              &
      p_pa_k_lower,        &
      p_th_k_upper,        &
      p_th_k,              &
      p_th_k_lower,        &
      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
    

    
!-----------------------------------------------------------------------
    DTRACER_MODULE = c0


!-----------------------------------------------------------------------
!  apply production of new Pa or Th from the decay of 235U and 234U
!-----------------------------------------------------------------------
    DTRACER_MODULE(:,:,pa_ind) = DTRACER_MODULE(:,:,pa_ind)+beta_pa/seconds_in_year
    DTRACER_MODULE(:,:,th_ind) = DTRACER_MODULE(:,:,th_ind)+beta_th/seconds_in_year
    
    
!-----------------------------------------------------------------------
!  apply radiodecay for Pa and Th
!-----------------------------------------------------------------------    
    pa_cur = p5*(TRACER_MODULE_OLD(:,:,k,pa_ind)+TRACER_MODULE_CUR(:,:,k,pa_ind))
    th_cur = p5*(TRACER_MODULE_OLD(:,:,k,th_ind)+TRACER_MODULE_CUR(:,:,k,th_ind))
    DTRACER_MODULE(:,:,pa_ind) = DTRACER_MODULE(:,:,pa_ind) - lambda_pa*pa_cur/seconds_in_year
    DTRACER_MODULE(:,:,th_ind) = DTRACER_MODULE(:,:,th_ind) - 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)
        pa_d(:,:,k,bid) = pa_cur/prodk_cndk_pa(:,:,k,bid)
        pa_p(:,:,k,bid) = pa_cur - pa_d(:,:,k,bid)
    
        th_d(:,:,k,bid) = th_cur/prodk_cndk_th(:,:,k,bid)
        th_p(:,:,k,bid) = th_cur - th_d(:,:,k,bid)
    end where  
    
    
    
    where(pa_d(:,:,k,bid).lt.c0)
        pa_d(:,:,k,bid) = c0
    end where
    
    where(pa_p(:,:,k,bid).lt.c0)
        pa_p(:,:,k,bid) = c0
    end where
    
    where(th_d(:,:,k,bid).lt.c0)
        th_d(:,:,k,bid) = c0
    end where
    
    where(th_p(:,:,k,bid).lt.c0)
        th_p(:,:,k,bid) = c0
    end where
                  
    if(k==1) then
        p_pa_k = pa_p(:,:,k,bid)
        pa_lower =  p5*(TRACER_MODULE_OLD(:,:,k+1,pa_ind)+ TRACER_MODULE_CUR(:,:,k+1,pa_ind))
        th_lower =  p5*(TRACER_MODULE_OLD(:,:,k+1,th_ind)+ TRACER_MODULE_CUR(:,:,k+1,th_ind))              
        p_pa_k_lower = pa_lower*(1-c1/prodk_cndk_pa(:,:,k+1,bid))                  
        
        p_pa_bot = p_pa_k + (p_pa_k_lower-p_pa_k)*dzwr(k)*p5*dz(k)
        
        
        
        p_th_k       = th_p(:,:,k,bid)        
        p_th_k_lower = th_lower*(1-c1/prodk_cndk_th(:,:,k+1,bid))
    
        p_th_bot    = p_th_k + (p_th_k_lower-p_th_k)*dzwr(k)*p5*dz(k)
                        
    	
    
    else if(k >1) then
       where((KMT(:,:,bid))>k)
    
            pa_upper  = p5*(TRACER_MODULE_OLD(:,:,k-1,pa_ind)+ TRACER_MODULE_CUR(:,:,k-1,pa_ind))
            pa_lower =  p5*(TRACER_MODULE_OLD(:,:,k+1,pa_ind)+ TRACER_MODULE_CUR(:,:,k+1,pa_ind))
     
            p_pa_k_upper = pa_upper*(1-c1/prodk_cndk_pa(:,:,k-1,bid))                         
   
            p_pa_k       = pa_p(:,:,k,bid)
    
            p_pa_k_lower = pa_lower*(1-c1/prodk_cndk_pa(:,:,k+1,bid))
    
    		p_pa_top = p_pa_k_upper + (p_pa_k - p_pa_k_upper)*dzwr(k-1)*p5*dz(k-1)
    		p_pa_bot = p_pa_k + (p_pa_k_lower - p_pa_k)*dzwr(k)*p5*dz(k)
    		
    
            th_upper  = p5*(TRACER_MODULE_OLD(:,:,k-1,th_ind)+ TRACER_MODULE_CUR(:,:,k-1,th_ind))
            th_lower =  p5*(TRACER_MODULE_OLD(:,:,k+1,th_ind)+ TRACER_MODULE_CUR(:,:,k+1,th_ind))
    
            p_th_k_upper = th_upper*(1-c1/prodk_cndk_th(:,:,k-1,bid))                       
        
            p_th_k       = th_p(:,:,k,bid)
    
            p_th_k_lower = th_lower*(1-c1/prodk_cndk_th(:,:,k+1,bid))
                
    		p_th_top = p_th_k_upper + (p_th_k - p_th_k_upper)*dzwr(k-1)*p5*dz(k-1)
    		p_th_bot = p_th_k + (p_th_k_lower - p_th_k)*dzwr(k)*p5*dz(k)

       end where
       
       
       where((KMT(:,:,bid)).eq.k)
            pa_upper  = p5*(TRACER_MODULE_OLD(:,:,k-1,pa_ind)+ TRACER_MODULE_CUR(:,:,k-1,pa_ind))
            p_pa_k_upper = pa_upper*(1-c1/prodk_cndk_pa(:,:,k-1,bid))
            p_pa_k       = pa_p(:,:,k,bid)
            p_pa_k_lower = c0
            
            p_pa_top = p_pa_k_upper + (p_pa_k - p_pa_k_upper)*dzwr(k-1)*p5*dz(k-1)
			p_pa_bot = p_pa_k
    
            th_upper  = p5*(TRACER_MODULE_OLD(:,:,k-1,th_ind)+ TRACER_MODULE_CUR(:,:,k-1,th_ind))
            p_th_k_upper = th_upper*(1-c1/prodk_cndk_th(:,:,k-1,bid)) 
            p_th_k       = th_p(:,:,k,bid)
            p_th_k_lower = c0
    
            p_th_top = p_th_k_upper + (p_th_k - p_th_k_upper)*dzwr(k-1)*p5*dz(k-1)
			p_th_bot = p_th_k
    
       end where
    

    end if
    
    pa_source(:,:,k,bid) =  DTRACER_MODULE(:,:,pa_ind)
    th_source(:,:,k,bid) =  DTRACER_MODULE(:,:,pa_ind)
    
    pa_sink(:,:,k,bid) = path_w*(p_pa_bot - p_pa_top)/dz(k)
    th_sink(:,:,k,bid) = path_w*(p_th_bot - p_th_top)/dz(k)
    

    
    DTRACER_MODULE(:,:,pa_ind) = DTRACER_MODULE(:,:,pa_ind) - pa_sink(:,:,k,bid)
    DTRACER_MODULE(:,:,th_ind) = DTRACER_MODULE(:,:,th_ind) - th_sink(:,:,k,bid)
 
    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(:,:,pa_ind) = c0
              DTRACER_MODULE(:,:,th_ind) = c0
    end where
      
    
!-----------------------------------------------------------------------
!EOC

 end subroutine path_set_interior

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

!***********************************************************************
!BOP
! !IROUTINE: path_tavg_forcing
! !INTERFACE:

 subroutine 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
   

       do k = 1,km
          call accumulate_tavg_field(pa_d(:,:,k,iblock),tavg_PA_D,iblock,k)
          call accumulate_tavg_field(pa_p(:,:,k,iblock),tavg_PA_P,iblock,k)
          call accumulate_tavg_field(th_d(:,:,k,iblock),tavg_TH_D,iblock,k)
          call accumulate_tavg_field(th_p(:,:,k,iblock),tavg_TH_P,iblock,k)
          call accumulate_tavg_field(pa_source(:,:,k,iblock),tavg_pa_source,iblock,k)
          call accumulate_tavg_field(pa_sink(:,:,k,iblock),tavg_pa_sink,iblock,k)
          call accumulate_tavg_field(th_source(:,:,k,iblock),tavg_th_source,iblock,k)
          call accumulate_tavg_field(th_sink(:,:,k,iblock),tavg_th_sink,iblock,k)
    	  		
    	 

      end do  
   
   end do

   !$OMP END PARALLEL DO

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

 end subroutine 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 path_mod

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