
>>>>>>>>>>>>>>>>>>> FILE: ../traj/earth.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE EARTH( T, TREF, XE)
      DOUBLE PRECISION  T, TREF, XE(6)
C
C     SUBROUTINE MOON ( XM)
      DOUBLE PRECISION  XM(6)
C
C-----------------------------------------------------------------------
C  EARTH (EARTH state) computes the inertial-frame heliocentric state
C  (position and velocity) of the Earth at a specified time.
C
C  Inputs:
C   T         Time (ephemeris days w.r.t. TREF).
C   TREF      Reference time (JED).
C
C  Output:
C   XE(1..6)  Inertial-frame position and velocity components of the
C             Earth w.r.t. the Sun at time T (AU, AU/d).
C
C----------------------------------------------------------------------
C  MOON (MOON state) returns the inertial-frame geocentric state
C  (position and velocity) of the Moon at the time specified in the
C  last call to EARTH.  EARTH must be called before MOON.
C
C  Outputs:
C   XM(1..6)  Inertial-frame position and velocity components of the
C             Moon w.r.t. the Earth at the time specified in last call
C             to EARTH (AU, AU/d).
C
C$---------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../traj/intgrt.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE INTGRT( T, SGLSTP, TI, X, NCA)
C
      USE DERIVS_MOD, ONLY: DERIVS_F, DERIVS_CASRCH, DERIVS_INI
      USE INTPAR_MOD, ONLY: DVX, DVY, DVZ, DVT
      USE ESTPAR_MOD, ONLY: DVX_EST, DVX_COL, DVY_EST, DVY_COL, 
     .     DVZ_EST, DVZ_COL
C
      DOUBLE PRECISION   T, TI, X(6,*)
      LOGICAL            SGLSTP
      INTEGER            NCA
C
C     SUBROUTINE INTINI( T, TEND, TREF, X, NC, WRTF10, PRTSUM, CAPMON)
      DOUBLE PRECISION   TEND, TREF
      INTEGER            NC
      LOGICAL            WRTF10, PRTSUM, CAPMON
C
C     SUBROUTINE INTTOL( TOL)
      DOUBLE PRECISION   TOL(3)
C
C     SUBROUTINE INTF10( CSTEP, TSPAN)
      DOUBLE PRECISION   CSTEP, TSPAN(2)
C
C     SUBROUTINE INTDIN( NSTEP)
      INTEGER            NSTEP
C
C     SUBROUTINE INTDHD( DHDR)
      CHARACTER*(*)      DHDR
C
C     SUBROUTINE INTINT( T, N, X)
      INTEGER            N
C
C     SUBROUTINE INTSSZ( SSZ)
      DOUBLE PRECISION   SSZ
C
C     SUBROUTINE INTFIN
C
C     SUBROUTINE INTSTA( NSTP, NREC)
      INTEGER            NSTP, NREC
C
C=======================================================================
C This package of routines interfaces with the integrator DIVA:
C
C-----------------------------------------------------------------------
C  INTGRT (INTeGRaTe) performs the integration, returning to the user
C  when the integration time TI has reached or just passed a specified
C  target time T, or after a single step has been taken, if the user has
C  set the flag SGLSTP, or when a close approach has occurred and the
C  user has requested monitoring of close approaches in the initializa-
C  tion call to INTINI.  INTGRT returns the current integration time TI,
C  the state (position and velocity) of the object at that time, and the
C  number of close approaches NCA which occurred in the latest step.
C  If NCA > 0, the user should call DERCAP to obtain the close approach
C  information.  If the ephemeris file (file 10) was requested in the
C  call to INTINI, INTGRT writes to this file as the integration
C  proceeds.  File 10 contains the integrator state at uniform time
C  intervals.  INTINI must always be called before INTGRT, and if file
C  10 is to be written, INTF10 must be called before INTINI.
C
C  Inputs:
C   T        Target time (d wrt reference time TREF).
C   SGLSTP   Set TRUE to return after only a single step.
C
C  Outputs:
C   TI          Integrator time (d wrt TREF).
C   X(1..6, 1)  State (position and velocity) at time TI.  Elements are
C               ordered X, Y, Z, VX, VY, VZ, where V is velocity.
C   NCA         Number of close approaches in latest integration step.
C
C-----------------------------------------------------------------------
C  INTINI (INTegrator INItialize) initializes the integrator data areas.
C  The integrator state includes the variational partials if requested
C  via the input argument NC.  If the writing of file 10 is requested,
C  it must already have been opened as an unformatted file, and INTF10
C  must already have been called.
C
C  Inputs:
C   T              Start time for the integration (d wrt ref time TREF).
C   TEND           End time for the integration, (d wrt ref time TREF).
C   TREF           Reference time (JED).
C   X(1..6, 1..NC) State and partials at start time T.  Column 1
C                  contains the position (AU) and velocity (AU/d) in the
C                  order X, Y, Z, VX, VY, VZ, where V is velocity.
C                  If NC > 1, columns 2 through NC contain the partials.
C   NC             Number of columns to integrate.  If NC = 1, only the
C                  state is integrated; if NC > 1, the variational
C                  partials are integrated along with the state.
C   WRTF10         Set TRUE to have file 10 written.
C   PRTSUM         Set TRUE to have an integration summary printed.
C   CAPMON         Set TRUE to have close approaches monitored.
C
C-----------------------------------------------------------------------
C  INTTOL (INTegrator TOLerance) sets the integrator tolerances on the
C  state, state partials, and non-grav partials.  The default values are
C  1.D-13, 1.D-8, and 1.D-6.  INTTOL must be called before INTINI.
C
C  Input:
C   TOL(1..3)  Integrator tolerances on state, state partials, and
C              non-grav partials.
C
C-----------------------------------------------------------------------
C  INTF10 (INTegrator File 10) sets the step size and time span to be
C  used when writing file 10.  INTF10 must be called before INTINI if
C  file 10 is requested.
C
C  Inputs:
C   CSTEP        File 10 step size (d).
C   TSPAN(1..2)  File 10 time span: start and stop times (JED).
C
C-----------------------------------------------------------------------
C  INTDIN (INTegration Display INitialize) sets the number of steps
C  between updates of the integration display.  A value of 0 disables
C  the integration display; the default is 100.  INTDIN must be called
C  before INTINI.
C
C  Input:
C   NSTEP    Number of steps between integration display updates.
C
C-----------------------------------------------------------------------
C  INTDHD (INTegration Display HeaDer) sets the header to be used on the
C  integration display.  The default is 'Current integration step:'.
C  The special header value '^' suppresses the writing of the header but
C  does not disable the integration display.  INTDHD must be called
C  before INTINI.
C
C  Input:
C   DHDR     Header for integration display.
C
C-----------------------------------------------------------------------
C  INTINT (INTegrator INTerpolate) interpolates for the state and
C  variational partials at a specified time within the latest
C  integration step.  If N = 1, only the state is interpolated.
C  INTINT can only be called after INTGRT is called.
C
C  Inputs:
C   T        Time at which the state and partials are desired (d wrt
C            ref time).
C   N        Number of columns of state and partials to be interpolated.
C            If N = 1, only the state is computed.
C
C  Outputs:
C   X(1..6, 1..N)  State and partials at time T (AU, AU/d).  Elements
C                  are ordered X, Y, Z, VX, VY, VZ, where V is velocity.
C
C-----------------------------------------------------------------------
C  INTFIN (INTegrator FINish) is used to finish off the integration.
C  This involves writing a software end of file on file 10 and rewinding
C  it, if this file has been requested, and printing the integration
C  summary, if it has been requested.  INTFIN can only be called after
C  INTINI.
C
C-----------------------------------------------------------------------
C  INTSTA (INTegrator STAtus) is used to retrieve the current inte-
C  gration status (number of steps taken and number of records written
C  to file 10).
C
C  Outputs:
C   NSTP  Number of steps.
C   NREC  Number of records written to file 10.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../traj/opngst.f <<<<<<<<<<<<<<<<<<<

      SUBROUTINE OPNGST( U )
      INTEGER            U
C
C-----------------------------------------------------------------------
C  OPNGST (OPeN Galilean Satellite Theory file) finds the filename for
C  the GST file from the environment variable and/or varlist input,
C  ensures that the file exists, and then opens it.
C
C  Input:
C   U    Unit number to use for GST file.
C
C$----------------------------------------------------------------------

>>>>>>>>>>>>>>>>>>> FILE: ../traj/intpar_mod.f90 <<<<<<<<<<<<<<<<<<<

module intpar_mod

use constants_mod, only: wp, gm_conv, rad_per_deg, deg_per_rad, au_per_day
use derivs_mod, only: sb_list, sb_mass, pl_mass, gs_mass, &
   tng_start, tng_stop, ng_a1, ng_a2, ng_a3, ng_dt, &
   ng_r0, ng_aln, ng_nm, ng_nn, ng_nk, &
   ng_aj1, ng_aj2, ng_alf, ng_del, ng_et1, ng_et2, ng_dth, &
   yarko, yrk_lgk, yrk_rho, &
   jg_sep, two_body, relativity, em_sep, jup_j2, jup_j4, area_mass_ratio, &
   derivs_calim
use estpar_mod, only: estpar_name, n_est_sb, sb_navio_idx, sb_idx, &
   dvx_idx, dvy_idx, dvz_idx, io_conversion
use yarko_mod, only: yarko_ini, yarko_msg, yrk_tag

implicit none
save
private

! Public procedures
public :: intpar_ini, intpar_flags, intpar_dump, intpar_def

! Public variables
public :: dvx, dvy, dvz, dvt

! Settable model parameters
real (kind=wp) :: dvt = 0._wp, dvx = 0._wp, dvy = 0._wp, dvz = 0._wp

!===============================================================================
! The intpar_mod module for initializing propagation options.
!
! Public routines:
!    intpar_ini   - Load propagator VARLIST options into integrator and 
!                   force model.
!    intpar_flags - Load command line switches.
!    intpar_dump  - Output integration parameters to orb file.
!    intpar_def   - VLDEF the relevant VARLIST options and specify valid 
!                   command line switches.
!
! The normal calling order should be intpar_def, intpar_flags, intpar_ini.
! The estpar_mod must be initialized and the small body ephemeris must be
! opened before calling intpar_ini.
!
!-------------------------------------------------------------------------------
! Written       Mar 06, 2002       Steve Chesley
! Modified      May 30, 2002       Steve Chesley
!   - Fixed bug in intpar_dump relating to sbsize call when sbeph not opened.
! Modified      Mar 06, 2003       Steve Chesley
!   - Added Yarkovsky stuff
! Modified      Nov 04, 2003       Steve Chesley
!   - Updated documentation.
!   - Restored '-q' option for backwards compatibility.
!   - Turned off dumping of nonstandard small body masses in intpar_dump.
! Modified      Nov 12, 2003       Steve Chesley
!   - Fixed bug in VLDEF call for TSPAN.
! Modified      Nov 23, 2004       Steve Chesley
!   - Let lone argument to INTPAR_DEF become optional.
! Modified      Aug 06, 2007       Steve Chesley
!   - Add parameters for impulsive delta-V
! Modified      Sep 04, 2007       Steve Chesley
!   - More tweaking of DV stuff. Added On-Off switch for comet nongravs.
! Modified      Sep 19, 2008       Steve Chesley
!   - Add output of g(r) parameters
!
! $Id: intpar_mod.f90,v 1.8 2008/09/19 23:30:38 chesley Exp $
!$==============================================================================

!^########################## INTPAR_INI ########################################
subroutine intpar_ini(tref, small_bodies, num_small_bodies, tspan_found)

real (kind=wp), intent(in) :: tref
character (len=*), intent(inout) :: small_bodies(:)
integer, intent(out) :: num_small_bodies
logical, intent(out) :: tspan_found

!===============================================================================
! intpar_ini loads propagation options. Some options pertain to the integrator 
! and others to the force model.
!
! Input:
!    tref - Reference time for delta-V epoch
!
! Input/Output:
!    small_bodies - array of names for requested small body perturbers. 
!       Updated on output with array of "official" names in the standard 
!       order. If first element is 'T' then all bodies in SBNAVIO file will 
!       be used. If first element is 'F' then no bodies will be used.
!
! Outputs:
!    num_small_bodies - The total number of perturbing small bodies to be
!        used in the force model.
!    tspan_found - Indicates that FILE10 is requested.
!$==============================================================================

!^########################## INTPAR_FLAGS ######################################
subroutine intpar_flags(flags)

character(len=*), intent(in) :: flags

!===============================================================================
! intpar_flags parses the command line switches and sets the relevant flags 
!    in the propagator.
!
! Input:
!    flags - command line flags, starting with a '-', e.g., '-q -nr'
!$==============================================================================

!^########################## INTPAR_DEF ########################################
subroutine intpar_def(valid_flags)

character (len=*), optional, intent(out) :: valid_flags

!===============================================================================
! intpar_def VLDEFs the VARLIST variables controlling the propagation model 
!    and specifies the valid command line switches to control the propagator.
!
! Output:
!    flags - command line flags, starting with a '-', e.g., '-q -nr'
!
! NOTE: The previous -q and -nr are replaced by the VARLIST variables QUIET and 
!    RELATIVITY, respectively. The -q option is preserved for backwards 
!    compatibility.
!$==============================================================================

!^########################## INTPAR_DUMP #######################################
subroutine intpar_dump(lin, k)

character (len=*), intent(inout) :: lin
integer, intent(inout) :: k

!===============================================================================
! intpar_dump dumps the user-definable integration parameters in ASCII
!    form to a file, ordinarily the orb file, using the LININI package. It
!    is C assumed that this package has already been initialized.
!
! Input/Output:
!    lin - Line buffer.
!    k   - Current line length.
!
! NOTE: Dumping of nonstandard perturber masses is disabled until VARLIST
!    issues can be resolved.
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../traj/estpar_mod.f90 <<<<<<<<<<<<<<<<<<<

module estpar_mod

use constants_mod, only: wp, gm_conv, deg_per_rad, sec_per_day, au_per_day

implicit none
save
private

! Public procedures
public :: estpar, estpar_dump, estpar_def

! Public variables
public :: estpar_msg
public :: estpar_name, io_conversion, apr_sigs, apr_vals
public :: n_est, n_est_dp, n_est_ndp, n_con, n_con_dp, n_con_ndp, &
   n_param, n_param_p1, n_dyn_param, n_nondyn_param, n_col, n_utm, n_utm_rhs, &
   last_est_col, n_elems, &
   a1_est, a1_col, a1_idx, &
   a2_est, a2_col, a2_idx, &
   a3_est, a3_col, a3_idx, &
   dt_est, dt_col, dt_idx, &
   aj1_est, aj1_col, aj1_idx, &
   aj2_est, aj2_col, aj2_idx, &
   alf_est, alf_col, alf_idx, &
   del_est, del_col, del_idx, &
   et1_est, et1_col, et1_idx, &
   et2_est, et2_col, et2_idx, &
   dth_est, dth_col, dth_idx, &
   lgk_est, lgk_col, lgk_idx, &
   rho_est, rho_col, rho_idx, &
   s0_est, s0_idx, &
   n_est_sb, n_est_pl, n_est_gs, &
   sb_idx, sb_col, sb_navio_idx, &
   pl_idx, pl_col, pl_navio_idx, &
   gs_idx, gs_col, gsthry_idx, &
   dvx_est, dvx_col, dvx_idx, &
   dvy_est, dvy_col, dvy_idx, &
   dvz_est, dvz_col, dvz_idx

! Module parameters (private)
integer, parameter :: max_est_sb = 50
integer, parameter :: n_pl = 10
integer, parameter :: n_gs = 4
integer, parameter :: max_est = max_est_sb + n_pl + n_gs + 4
integer, parameter :: max_con = max_est
integer, parameter :: max_param = 6 + max_est

! Public module variables
character (len=256) :: estpar_msg ! Diagnostic message

! Estimation parameter counts. Initialize with generic values
integer :: n_est=0, n_est_dp=0, n_est_ndp=0
integer :: n_con=0, n_con_dp=0, n_con_ndp=0
integer :: n_elems=6, n_param=6, n_param_p1=7, n_utm=21, n_utm_rhs=28
integer :: n_dyn_param=6, n_nondyn_param=0, n_col=7, last_est_col=7
integer :: n_est_sb=0, n_est_pl=0, n_est_gs=0

! Estimation flags, column numbers (_col) and input indices (_idx)
logical :: s0_est
integer :: s0_idx
logical :: a1_est, a2_est, a3_est, dt_est
integer :: a1_idx, a2_idx, a3_idx, dt_idx
integer :: a1_col, a2_col, a3_col, dt_col
logical :: aj1_est, aj2_est, alf_est, del_est, et1_est, et2_est, dth_est
integer :: aj1_idx, aj2_idx, alf_idx, del_idx, et1_idx, et2_idx, dth_idx
integer :: aj1_col, aj2_col, alf_col, del_col, et1_col, et2_col, dth_col
logical :: lgk_est, rho_est
integer :: lgk_idx, rho_idx
integer :: lgk_col, rho_col
integer :: sb_idx(max_est_sb), sb_col(max_est_sb), sb_navio_idx(max_est_sb)
integer :: pl_idx(n_pl), pl_col(n_pl), pl_navio_idx(n_pl)
integer :: gs_idx(n_gs), gs_col(n_gs), gsthry_idx(n_gs)
logical :: dvx_est, dvy_est, dvz_est
integer :: dvx_idx, dvy_idx, dvz_idx
integer :: dvx_col, dvy_col, dvz_col

! Allocatable public arrays: a priori sigmas and parameter names
real (kind=wp), allocatable :: apr_sigs(:)
real (kind=wp), allocatable :: apr_vals(:)
real (kind=wp), allocatable :: io_conversion(:)
character (len=3), allocatable :: estpar_name(:)

!===============================================================================
! The estpar_mod module is the area where information is computed and stored 
! concerning the number and types of estimated and considered parameters. The 
! following VARLIST variables are used:
!     EST(*) - List of estimated parameters. Possible values are A1, A2, A3,
!              S0, DT, AJ, and M01-M99. Example: The parameter M34 indicates
!              the mass of the 34th object in the small body perturber file.
!     CON(*) - List of consider parameters. Possible values include those 
!              listed for EST(*), plus P01-P10 and GS1-GS4, which refer to 
!              planetary and Galilean satellite masses, respectively.
!     APR_EST_SIGS(*) - Apriori uncertainty values for the parameters listed in 
!              EST(*). The ordering is assumed to be consistent with the
!              EST(*) ordering.
!     APR_CON_SIGS(*) - Apriori uncertainty values for the parameters listed in 
!              CON(*). The ordering is assumed to be consistent with the
!              CON(*) ordering.
!     APR_EST_VALS(*) - Apriori parameter values for the parameters listed in 
!              EST(*). The ordering is assumed to be consistent with the
!              EST(*) ordering. If a value is not present then the nominal 
!              value will be used.
!     APR_ELE_SIGS(6) - Apriori uncertainty values for the six orbital elements.
!     APR_ELE_VALS(6) - Apriori parameter values for the six orbital elements.
!
! The module keeps track of the _index_ for each parameter, which is the 
! location in the estpar_name array, as well as within the elements vector,
! covariance, etc. It also maintains the _column_ for each (dynamic) parameter.
! This is the column within the combined (6 x n_col) state and partials matrix
! that is to be propagated.
!
! Public routines: 
!    estpar      - Loads the various VARLIST options and sets the flags and 
!        values for use by other routines.
!    estpar_def  - VLDEFs the various VARLIST options used by estpar.
!    estpar_dump - Dump EST, CON, and APR_CON_SIGS to orbfile.
!
! Public variables:
!    estpar_msg - Diagnostic message in case of error
!    estpar_name - Array of parameter names
!    io_conversion - Array of conversion constants. 
!          io_conversion * internal_array converts the internal_array to 
!          units used for output, e.g., rad -> deg or AU^3/d^2 -> km_3/s^2.
!    apr_sigs - Array of a priori uncertainty values associated with the
!        parameters in estpar_name.
!    apr_vals - Array of a priori parameter values associated with the
!        parameters in estpar_name.
!
!------------------ Parameter counts and associated constants ------------------
!
!    n_est - Number estimated parameters (from EST list)
!    n_est_dp - Number of estimated dynamical parameters 
!    n_est_ndp - Number of estimated nondynamical parameters
!    n_con - Number consider parameters (from CON list)
!    n_con_dp - Number of dynamical consider parameters
!    n_con_ndp - Number of nondynamical consider parameters
!    n_param - Total number of parameters, including elements, estimated and 
!        considered (6+n_est+n_con). This is the dimension of estpar_name.
!    n_param_p1 - n_param + 1
!    n_dyn_param - Total number of dynamical parameters, including elements, 
!        estimated and considered (6+n_est_dp+n_con_dp).
!    n_nondyn_param - Totla number of nondynamical parameters 
!        (n_est_ndp+n_con_ndp).
!    n_elems - Number of elements in the parameter list that are
!        actually estimated, as opposed to considered.
!    n_col - Number of columns in X, the combined state and partials matrix
!        used in the propagation.
!    n_utm - Number of entries in a vector-stored covariance or information
!        matrix (n_param*(n_param+1)/2)
!    n_utm_rhs - Number of entries in a vector-stored matrix information
!        matrix with RHS (n_utm + n_param_p1).
!    last_est_col - Column location of the last estimated (as opposed to
!        considered) parameter in the X matrix.
!
!------------------ Flags, indices and column numbers --------------------------
!
!    The numerous estimable quantities all have the associated variables:
!         *_est - flag to indicate the variable is to be estimated/considered.
!         *_idx - integer to indicate the index within the parameter list.
!         *_col - integer to indicate the column location in X. (This quantity 
!                 is not present for nondynamical parameters, e.g., S0.)
!
!    The available prefixes for these paramters are
!         Center-of-light offset
!           s0 - offset in km at 1 AU from sun
!         Classical comet non-gravs
!           a1, a2, a3 - RTN acceleration paramters
!           dt - Time delay, in days AFTER perihelion of the peak accel.
!         Comet rotating jet nongravs (1,2 refer to jet number(
!           aj1, aj2 - Jet strengths
!           et1, et2 - Latitude of comet jet
!           alf, del - RA, DEC of comet pole
!         Yarkovsky acceleration:
!           lgk - Log10 of K, the thermal conductivity
!           rho - The bulk density
!
!    Similar parameters exist in array form for the various perturber masses:
!           n_est_sb - Number of estimated or considered small body masses
!           n_est_pl - Number of considered planetary masses
!           n_est_gs - Number of considered Galilean satallite masses
!           sb_idx, pl_idx, gs_idx - Array of param. list locations for est/con
!               small body masses, planetary masses and galilean sat. masses.
!           sb_col, pl_col, gs_col - Array of column locs. within the X matrix
!               for est/con small body, planetary, and galilean sat. masses.
!           sb_navio_idx, pl_navio_idx, gsthry_idx - Array of locs. within the
!               NAVIO/GSTRHY file for estimated or considered small body,
!               planetary, and galilean sat. masses.
!
!-------------------------------------------------------------------------------
! References: 
!-------------------------------------------------------------------------------
! Written       Mar 06, 2002       Steve Chesley
! Modified      Mar 18, 2002       Steve Chesley
!   - Added io_conversion array.
!   - Check for missing APR_CON_SIGS.
! Modified      Apr 25, 2002       Steve Chesley
!   - Added comet jetting stuff.
! Modified      Mar 06, 2003       Steve Chesley
!   - Added Yarkovsky stuff.
! Modified      Sep 17, 2003       Steve Chesley
!   - Added documentation.
! Modified      Dec 01, 2004       Steve Chesley
!   - Initialized n_param, etc. for generic OD case.
! Modified      Jan 06, 2005       Steve Chesley
!   - Corrected initial values for n_est and n_est_dp.
!   - Check and warn for neglected EST and CON entries.
!
! $Id: estpar_mod.f90,v 1.6 2007/08/14 22:09:16 chesley Exp $
!===============================================================================

contains

!########################### ESTPAR ############################################
subroutine estpar(errcod)

integer, intent(out) :: errcod

!===============================================================================
! estpar loads the various VARLIST options and sets the flags and values for
! use by other routines.
!
! Outputs:
!    errcod - Zero indicates success.
!$==============================================================================

!^########################## ESTPAR_DUMP #######################################
subroutine estpar_dump(lin, k)

character (len=*), intent(inout) :: lin
integer, intent(inout) :: k

!===============================================================================
! estpar_dump dumps the user-definable integration
!  parameters in ASCII form to a file using the LININI package.  It is
!  assumed that this package has already been initialized.
!
! Output the following items, which are needed for mapping:
!        EST, CON, APR_CON_SIGS
! The following items are not recorded since they are not needed for mapping.
!        APR_EST_SIGS, APR_EST_VALS, APR_ELE_SIGS, APR_ELE_VALS
! (We have the a posteriori values in the SRC.)
!
!  Input/Outputs:
!   lin - Line buffer.
!   k   - Current line length.
!
!$==============================================================================

!^########################## ESTPAR_DEF ########################################
subroutine estpar_def()

!===============================================================================
! estpar_def VLDEFs the various VARLIST options used by estpar.
!$==============================================================================

>>>>>>>>>>>>>>>>>>> FILE: ../traj/derivs_mod.f90 <<<<<<<<<<<<<<<<<<<

module derivs_mod

use constants_mod, only: wp, au, au_per_day, v_light, &
   gm_sun, gm_conv, earth_moon_ratio, earth_moon_ratio_p1
use estpar_mod, only: n_col, a1_col, a1_est, a2_col, a2_est, a3_col, a3_est, &
   dt_col, dt_est, aj1_col, aj1_est, aj2_col, aj2_est, alf_col, alf_est, &
   del_col, del_est, et1_col, et1_est, et2_col, et2_est, dth_col, dth_est, &
   lgk_col, lgk_est, rho_col, rho_est, &
   n_est_sb, sb_col, sb_navio_idx, &
   n_est_pl, pl_col, pl_navio_idx, &
   n_est_gs, gs_col, gsthry_idx
use jet_mod, only: jet, jet_part_r, jet_part_p
use yarko_mod, only: yarko_accel, yarko_part, yarko_msg

implicit none
save
private

! Public procedures
public :: derivs_f, derivs_ini, derivs_calim, derivs_casrch, derivs_ca

! Public variables
public :: derivs_msg
public :: two_body, em_sep, jg_sep, relativity, jup_j2, jup_j4
public :: area_mass_ratio
public :: tng_start, tng_stop, ng_a1, ng_a2, ng_a3, ng_dt
public :: ng_r0, ng_aln, ng_nm, ng_nn, ng_nk
public :: ng_aj1, ng_aj2, ng_alf, ng_del, ng_et1, ng_et2, ng_dth
public :: yarko, yrk_lgk, yrk_rho
public :: sb_list, pl_mass, gs_mass, sb_mass

! Private module parameters
integer, parameter :: max_n_col = 100
integer, parameter :: idx0_sb = 14
integer, parameter :: max_n_sb = 128
integer, parameter :: max_n_pert = idx0_sb + max_n_sb

! Public module variables
character (len=256) :: derivs_msg ! Diagnostic message

! Model options
logical :: two_body    = .false.  ! Keplerian model
logical :: relativity  = .true.   ! Include relativistic acceleration terms
logical :: em_sep      = .true.   ! Separate Earth and Moon
logical :: jg_sep      = .false.  ! Separate Galileamn sats. from Jupiter
logical :: jup_j2      = .false.  ! Include J_2 of jupiter
logical :: jup_j4      = .false.  ! Include J_4 of Jupiter
logical :: yarko       = .false.  ! Include Yarkovsky acceleration

! Settable model parameters
real (kind=wp) :: area_mass_ratio = 0._wp
real (kind=wp) :: tng_start = -9.e99_wp
real (kind=wp) :: tng_stop = 9.e99_wp
real (kind=wp) :: ng_a1  = 0._wp
real (kind=wp) :: ng_a2  = 0._wp
real (kind=wp) :: ng_a3  = 0._wp
real (kind=wp) :: ng_dt  = 0._wp
real (kind=wp) :: ng_r0  = 2.808_wp
real (kind=wp) :: ng_aln = 0.1112620426_wp
real (kind=wp) :: ng_nm  = 2.15_wp
real (kind=wp) :: ng_nn  = 5.093_wp
real (kind=wp) :: ng_nk  = 4.6142_wp
real (kind=wp) :: ng_aj1  = 0._wp
real (kind=wp) :: ng_aj2  = 0._wp
real (kind=wp) :: ng_alf = 0._wp
real (kind=wp) :: ng_del = 0._wp
real (kind=wp) :: ng_et1 = 0._wp
real (kind=wp) :: ng_et2 = 0._wp
real (kind=wp) :: ng_dth = 0._wp
real (kind=wp) :: yrk_lgk = -2._wp   ! log10(W/m/K)
real (kind=wp) :: yrk_rho = 2500._wp ! kg/m^3

! Perturber request lists, masses and close approach distances
real (kind=wp) :: pl_mass(10) = -9.e99_wp
real (kind=wp) :: gs_mass(4) = -9.e99_wp
integer :: sb_list(max_n_sb) = 0
real (kind=wp) :: sb_mass(max_n_sb) = -9.e99_wp

!===============================================================================
! The derivs_mod module. It does a lot.
!
! Public routines: 
!    derivs_f      - Computes the accelerations and derivatives of the 
!                    variational equations.
!    derivs_ini    - Initializes force model. Must be called before derivs_f.
!    derivs_calim  - Sets the modeule limits for close approach monitoring.
!    derivs_casrch - Searches for all close approaches in the current time step.
!    derivs_ca     - Returns the bodies, times and ranges of all close 
!                    approaches in the current timestep.
!
! Public variables:
!     derivs_msg - Diagnostic message in case of error.
!   Force modeling:
!     two_body   - Turn off all forces except the solar acceleration.
!     relativity - Use relativity in force model.
!     area_mass_ratio - Solar radiation pressure parameter.
!   Comet Nongravs:
!     ng_a1      - A1, radial acceleration parameter
!     ng_a2      - A2, transverse acceleration parameter
!     ng_a3      - A3, normal acceleration parameter
!     ng_r0      - g(r): r_0
!     ng_aln     - g(r): alpha, normalizing factor
!     ng_nm      - g(r): m exponent
!     ng_nn      - g(r): n exponent
!     ng_nk      - g(r): k exponent
!     ng_dt      - g(r): DT, peak acceleration delay
!     ng_aj1     - Jet model: AJ1, intrinsic strength of jet 1
!     ng_aj2     - Jet model: AJ2, intrinsic strength of jet 2
!     ng_alf     - Jet model: alpha, RA of pole
!     ng_del     - Jet model: delta, DEC of pole
!     ng_et1     - Jet model: eta1, colatitude of jet 1
!     ng_et2     - Jet model: eta2, colatitude of jet 2
!     ng_dth     - Jet model: diurnal lag angle
!   Yarkovsky:
!     yarko      - Apply Yarkovsky accels.
!     yrk_lgk    - log10 K, where K is surface thermal conductivity.
!     yrk_rho    - bulk density
!   Perturbers:
!     em_sep     - Separate the gravity fields of the Earth and Moon.
!     jg_sep     - Separate the gravityt fields of Jupter and the Galileans.
!     jup_j2     - Apply J2 of Jovian gravity field.
!     jup_j4     - Apply J4 of Jovian gravity field.
!     sb_list    - List of small body perturbers used
!     pl_mass    - Planetary masses
!     gs_mass    - Galilean masses
!     sb_mass    - Small body masses
!
! References: 
! R.A. Jacobson, "The Gravity Field of the Jovian System and the Orbits of 
!    the Regular Jovian Satellites" BAAS 33(3) 1039, 2001.
!-------------------------------------------------------------------------------
! Original F77 version written by  Paul Chodas, July, 1990      
!     - Modeled after earlier Unisys routine.
! Translated to F90 by Steve Chesley, March, 2002
!     - Added more flexible parameter estimation capability and small body 
!       mass partials.
! Modified      Nov 04, 2003       Steve Chesley
!     - Many changes since March, 2003!
!     - A3 and DT partials
!     - Jet model
!     - Yarkovsky acceleration
! Modified      Jan 06, 2005       Steve Chesley
!     - Formalized initialization of rho_lim(:10).
! Modified      Sep 04, 2008       Steve Chesley
!     - Tweaked to properly handle zero mass planets
!     - And output perturber masses when nonstd masses are used
!
! $Id: derivs_mod.f90,v 1.8 2008/09/19 23:29:01 chesley Exp $
!$==============================================================================

!^########################## DERIVS_F ##########################################
subroutine derivs_f(t, y, f, kord_1)

real (kind=wp), intent(in)  :: t
real (kind=wp), intent(in)  :: y(2,3,n_col)
integer, intent(in) :: kord_1
real (kind=wp), intent(out) :: f(3,n_col)

!===============================================================================
! derivs_f computes the inertial-frame components of the acceleration
!    with respect to the Sun of an orbiting body at a given time, given the
!    heliocentric position and velocity of the orbiting body. The routine
!    also computes the derivative of the variational partials if requested
!    in the initialization call to derivs_ini. The equations of motion
!    include point-mass perturbations due to all the planets, perturbations
!    due to cometary outgassing, and relativistic effects. Subroutine
!    derivs_ini must always be called after changing force model parameters
!    and before calling derivs_f.
!
! Inputs:
!    t      - Current integration time, relative to tref.
!    y      - 2 x 3 x n_col matrix containing the current state and partials. 
!             The indices are coded as follows:
!                Index 1: 1, refers to position and 2 refers to velocity
!                Index 2: X,Y,Z components, respectively.
!                Index 3: 1 contains the actual state, 2-7 contain the state 
!                  partials, and higher values contain dynamic param. partials.
!    kord_1 - Flag indicating reason for call:
!                1 - Predict state,
!                3 - Correct state, predict variational partials if present
!                2 - If n_int_col > 1 then correct variational partials only
!                        (with state unchanged)
!                  - If n_int_col = 1 then correct state (partials not present).
!
! Outputs:
!    f      - Column 1 contains acceleration at time T (AU/d^2), other columns
!                contain the time derivatives of the variational partials.
!$==============================================================================

!^########################## DERIVS_INI ########################################
subroutine derivs_ini(tref, n_intgrt_col, forward)

real (kind=wp), intent(in) :: tref
integer, intent(in) :: n_intgrt_col
logical, intent(in) :: forward

!===============================================================================
! derivs_ini initializes the force model by reading constants from the
!    planetary ephemeris and computing various other constants used by DERIVS.
!
! Inputs:
!    tref         - Time reference point (JED).
!    n_intgrt_col - Number of columns in the state and partials matrix. Acts 
!       as a flag for computing varaitional equations: 1 for state only, higher
!       numbers if propagation of partials is requested.
!    forward      - Flag to indicate the time direction of the integration. 
!       Used during close approach search.
!$==============================================================================

!^########################## DERIVS_CALIM ######################################
subroutine derivs_calim(pl_calim, gs_calim, sb_calim)

real (kind=wp), intent(in) :: pl_calim(:)
real (kind=wp), intent(in) :: gs_calim(:)
real (kind=wp), intent(in) :: sb_calim(:)

!===============================================================================
! derivs_calim is used to specify the range limits to be used when
!    checking for close approaches to the 9 planets, the Moon, Galilean
!    satelites, and specified small bodies.  The routine derivs_casrch,
!    below, checks for close approaches only when the range from a body is
!    less than the range limit for that body.  The default range limits for
!    the 9 planets, in AU, are 0.1, 3*0.5, 4*2, and 1; the defaults for the
!    small bodies are 0.1 AU.
!
! Inputs:
!    pl_calim - Range limits for the 9 planets and Earth's moon.
!    gs_calim - Range limits for the Galilean moons J1-J4.
!    sb_calim - Range limits for selected small bodies (AU).
!$==============================================================================

!^########################## DERIVS_CASRCH #####################################
subroutine derivs_casrch(n_cap)

integer, intent(out) :: n_cap

!===============================================================================
! derivs_casrch is used to perform the search for close approaches.
!    It should be called after derivs_ini and after every integrator step.
!    The number of close approaches found in the current step is returned
!    to the calling routine.  The user should call derivs_ca that number of
!    times to get the exact close approach times and ranges before
!    proceeding with the integration.
!
! Output:
!    n_cap - Number of close approaches in current integration step.
!$==============================================================================

!^########################## DERIVS_CA #########################################
subroutine derivs_ca(i_body, t_ca, rho_ca)

integer, intent(out) :: i_body
real (kind=wp), intent(out) :: t_ca
real (kind=wp), intent(out) :: rho_ca

!===============================================================================
! derivs_cap is used to obtain the body index of a close approach, and
!    the precise time and close approach range.
!
! Outputs:
!    i_body - Body index of close approach.
!    t_ca   - Time of close approach (d wrt TREF).
!    rho_ca - Range of close approach (AU).
!$==============================================================================
