jacchia_roberts_module – jacchia-roberts-fortran

Jacchia-Roberts atmospheric density model.

The Jacchia-Roberts atmosphere model is valid for altitudes from 100 km to 2500 km.

Reference:

Roberts, C. E., Jr., "An Analytic Model for Upper Atmosphere Densities Based Upon Jacchia's 1970 Models", Celestial Mechanics, Vol. 4, pp. 368-377, 1971.
See GMAT: General Mission Analysis Tool files: JacchiaRobertsAtmosphere.cpp and SolarFluxReader.cpp.
See also: https://github.com/jacobwilliams/INPE-atmosphere-models

Uses

- jacchia_roberts_kinds
- space_weather_module
Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.

Variables

Type	Visibility	Attributes		Name		Initial
real(kind=dp),	private,	parameter	::	PI	=	acos(-1.0_dp)	Pi constant
real(kind=dp),	private,	parameter	::	RAD_PER_DEG	=	PI/180.0_dp	degree to radian conversion factor
real(kind=dp),	private,	parameter	::	RHO_ZERO	=	3.46e-9_dp	Low altitude density in g/cm^3
real(kind=dp),	private,	parameter	::	TZERO	=	183.0_dp	Temperature in degrees Kelvin at height of 90 km
real(kind=dp),	private,	parameter	::	G_ZERO	=	9.80665_dp	Earth gravitational constant m/s^2
real(kind=dp),	private,	parameter	::	GAS_CON	=	8.31432_dp	Gas constant (joules/(degK-mole))
real(kind=dp),	private,	parameter	::	AVOGADRO	=	6.022045e23_dp	Avogadro's number
real(kind=dp),	private,	parameter	::	nominalKp	=	3.0_dp	Nominal Kp index
real(kind=dp),	private,	parameter	::	nominalF107	=	150.0_dp	Nominal F10.7 value (solar flux units)
real(kind=dp),	private,	parameter	::	nominalF107a	=	150.0_dp	Nominal F10.7a value (solar flux units)
real(kind=dp),	private,	parameter	::	CON_C(5)	=	[-89284375.0_dp, 3542400.0_dp, -52687.5_dp, 340.5_dp, -0.8_dp]	Constants for series expansion
real(kind=dp),	private,	parameter	::	CON_L(5)	=	[0.1031445e5_dp, 0.2341230e1_dp, 0.1579202e-2_dp, -0.1252487e-5_dp, 0.2462708e-9_dp]	Constants for series expansion
real(kind=dp),	private,	parameter	::	MZERO	=	28.82678_dp	Constants required between 90 km and 100 km
real(kind=dp),	private,	parameter	::	M_CON(7)	=	[-435093.363387_dp, 28275.5646391_dp, -765.33466108_dp, 11.043387545_dp, -0.08958790995_dp, 0.00038737586_dp, -0.000000697444_dp]	Constants for series expansion of M(z) function
real(kind=dp),	private,	parameter	::	S_CON(6)	=	[3144902516.672729_dp, -123774885.4832917_dp, 1816141.096520398_dp, -11403.31079489267_dp, 24.36498612105595_dp, 0.008957502869707995_dp]	Constants for series expansion of S(z) function
real(kind=dp),	private,	parameter	::	S_BETA(6)	=	[-52864482.17910969_dp, -16632.50847336828_dp, -1.308252378125_dp, 0.0_dp, 0.0_dp, 0.0_dp]	Constants for series expansion of S(z) function - temperature part
real(kind=dp),	private,	parameter	::	OMEGA	=	-0.94585589_dp	Constants required for altitudes between 100 km and 125 km
real(kind=dp),	private,	parameter	::	ZETA_CON(7)	=	[0.1985549e-10_dp, -0.1833490e-14_dp, 0.1711735e-17_dp, -0.1021474e-20_dp, 0.3727894e-24_dp, -0.7734110e-28_dp, 0.7026942e-32_dp]	Constants for series expansion
real(kind=dp),	private,	parameter	::	MOL_MASS(6)	=	[28.0134_dp, 39.948_dp, 4.0026_dp, 31.9988_dp, 15.9994_dp, 1.00797_dp]	Molecular masses of atmospheric constituents in grams/mole
real(kind=dp),	private,	parameter	::	NUM_DENS(5)	=	[0.78110_dp, 0.93432e-2_dp, 0.61471e-5_dp, 0.161778_dp, 0.95544e-1_dp]	Number density divided by Avogadro's number of atmospheric constituents
real(kind=dp),	private,	parameter	::	CON_DEN(5,7)	=	reshape([0.1093155e2_dp, 0.1186783e-2_dp, -0.1677341e-5_dp, 0.1420228e-8_dp, -0.7139785e-12_dp, 0.1969715e-15_dp, -0.2296182e-19_dp, 0.8049405e1_dp, 0.2382822e-2_dp, -0.3391366e-5_dp, 0.2909714e-8_dp, -0.1481702e-11_dp, 0.4127600e-15_dp, -0.4837461e-19_dp, 0.7646886e1_dp, -0.4383486e-3_dp, 0.4694319e-6_dp, -0.2894886e-9_dp, 0.9451989e-13_dp, -0.1270838e-16_dp, 0.0_dp, 0.9924237e1_dp, 0.1600311e-2_dp, -0.2274761e-5_dp, 0.1938454e-8_dp, -0.9782183e-12_dp, 0.2698450e-15_dp, -0.3131808e-19_dp, 0.1097083e2_dp, 0.6118742e-4_dp, -0.1165003e-6_dp, 0.9239354e-10_dp, -0.3490739e-13_dp, 0.5116298e-17_dp, 0.0_dp], [5, 7], order=[2, 1])	Polynomial coefficients for constituent densities of each atmospheric gas

Derived Types

type, public :: geoparms_type

Geomagnetic parameters

Components

Type	Visibility	Attributes		Name		Initial
real(kind=dp),	public		::	tkp	=	0.0_dp	Geomagnetic index Kp
real(kind=dp),	public		::	xtemp	=	0.0_dp	Exospheric temperature (K)

type, public :: jacchia_roberts_type

Jacchia-Roberts atmosphere model type

Components

Type	Visibility		Name		Initial
real(kind=dp),	private	::	cb_polar_radius	=	0.0_dp	Central body polar radius (km)
real(kind=dp),	private	::	cb_polar_squared	=	0.0_dp	Polar radius squared (km^2)
real(kind=dp),	private	::	root1	=	0.0_dp	Auxiliary temperature root
real(kind=dp),	private	::	root2	=	0.0_dp	Auxiliary temperature root
real(kind=dp),	private	::	x_root	=	0.0_dp	Complex root real part
real(kind=dp),	private	::	y_root	=	0.0_dp	Complex root imaginary part (absolute value)
real(kind=dp),	private	::	t_infinity	=	0.0_dp	Exospheric temperature at infinity
real(kind=dp),	private	::	tx	=	0.0_dp	Temperature at boundary
real(kind=dp),	private	::	sum	=	0.0_dp	Intermediate temperature sum
logical,	private	::	use_sw_file	=	.true.	Flag to control whether to use space weather file data or specified values.
real(kind=dp),	private	::	f107_override	=	nominalF107	Override value for F10.7 if not using space weather file
real(kind=dp),	private	::	f107a_override	=	nominalF107a	Override value for F10.7a if not using space weather file
real(kind=dp),	private	::	kp_override	=	nominalKp	Override value for Kp index if not using space weather file
type(sw_data_type),	private	::	sw_data			Space weather data type (from space_weather_module)

Type-Bound Procedures

procedure, public :: density => jacchia_roberts_density
procedure, public :: initialize => jr_init
procedure, public :: destroy => jr_cleanup
procedure, private :: exotherm
procedure, private :: rho_100
procedure, private :: rho_125
procedure, private :: rho_high
procedure, private :: parms_to_temp

Functions

private pure function parms_to_temp(me, f107_daily, f107a_avg, tkp) result(xtemp)

Compute exospheric temperature from F10.7 values using the Jacchia-Roberts formula

Arguments

Type	Intent		Name
class(jacchia_roberts_type),	intent(in)	::	me
real(kind=dp),	intent(in)	::	f107_daily	Daily F10.7 value
real(kind=dp),	intent(in)	::	f107a_avg	81-day average F10.7a
real(kind=dp),	intent(in)	::	tkp	Kp index

Return Value real(kind=dp)

Exospheric temperature (K)

private function jacchia_roberts_density(me, height, position, sun_vector, geo_lat, utc_mjd) result(density)

Compute atmospheric density using the Jacchia-Roberts model

Arguments

Type	Intent		Name
class(jacchia_roberts_type),	intent(inout)	::	me
real(kind=dp),	intent(in)	::	height	Spacecraft height above reference ellipsoid (km)
real(kind=dp),	intent(in)	::	position(3)	Spacecraft position vector (km, TOD/GCI)
real(kind=dp),	intent(in)	::	sun_vector(3)	Unit vector to Sun (TOD/GCI)
real(kind=dp),	intent(in)	::	geo_lat	Geodetic latitude (degrees)
real(kind=dp),	intent(in)	::	utc_mjd	UTC Modified Julian Date

Return Value real(kind=dp)

Atmospheric density (kg/m^3)

private function exotherm(me, position, sun_vector, geo, height, sun_dec, geo_lat) result(exotemp)

Compute the temperature of Earth's atmosphere and auxiliary temperature-related quantities at a given altitude

Arguments

Type	Intent		Name
class(jacchia_roberts_type),	intent(inout)	::	me
real(kind=dp),	intent(in)	::	position(3)	Spacecraft position (km)
real(kind=dp),	intent(in)	::	sun_vector(3)	Sun unit vector
type(geoparms_type),	intent(in)	::	geo	Geomagnetic parameters
real(kind=dp),	intent(in)	::	height	Spacecraft height (km)
real(kind=dp),	intent(in)	::	sun_dec	Sun declination (radians)
real(kind=dp),	intent(in)	::	geo_lat	Geodetic latitude (radians)

Return Value real(kind=dp)

Local exospheric temperature (K)

private function rho_100(me, height, temperature) result(density)

Compute density of the atmosphere between 90 and 100 km

Arguments

Type	Intent		Name
class(jacchia_roberts_type),	intent(inout)	::	me
real(kind=dp),	intent(in)	::	height	Spacecraft altitude (km)
real(kind=dp),	intent(in)	::	temperature	Exospheric temperature (K)

Return Value real(kind=dp)

Atmospheric density (g/cm^3)

private function rho_125(me, height, temperature) result(density)

Compute density of the atmosphere between 100 and 125 km

Arguments

Type	Intent		Name
class(jacchia_roberts_type),	intent(inout)	::	me
real(kind=dp),	intent(in)	::	height	Spacecraft altitude (km)
real(kind=dp),	intent(in)	::	temperature	Exospheric temperature (K)

Return Value real(kind=dp)

Atmospheric density (g/cm^3)

private function rho_high(me, height, temperature, t_500, sun_dec, geo_lat) result(density)

Compute density of the atmosphere between 125 and 2500 km

Arguments

Type	Intent		Name
class(jacchia_roberts_type),	intent(inout)	::	me
real(kind=dp),	intent(in)	::	height	Spacecraft altitude (km)
real(kind=dp),	intent(in)	::	temperature	Exospheric temperature at spacecraft altitude (K)
real(kind=dp),	intent(in)	::	t_500	Exospheric temperature at altitude of 500 km (K)
real(kind=dp),	intent(in)	::	sun_dec	Declination of the sun in TOD coordinates (radians)
real(kind=dp),	intent(in)	::	geo_lat	Geodetic latitude of spacecraft (radians)

Return Value real(kind=dp)

private pure function rho_cor(height, utc_mjd, geo_lat, geo) result(correction)

Calculates the density correction factor

Arguments

Type	Intent		Name
real(kind=dp),	intent(in)	::	height	Spacecraft altitude (km)
real(kind=dp),	intent(in)	::	utc_mjd	UTC Modified Julian Date
real(kind=dp),	intent(in)	::	geo_lat	Geodetic latitude of spacecraft (radians)
type(geoparms_type),	intent(in)	::	geo	Geomagnetic parameters

Return Value real(kind=dp)

Correction factor (multiplicative)

Subroutines

private subroutine jr_init(me, cb_polar_radius, status, filename, f107_override, f107a_override, kp_override)

Initialize the Jacchia-Roberts module with central body parameters.

Arguments

Type	Intent	Optional		Name
class(jacchia_roberts_type),	intent(inout)		::	me
real(kind=dp),	intent(in)		::	cb_polar_radius	Polar radius of central body (km)
integer(kind=ip),	intent(out)		::	status	Output status (0=success, non-zero=error)
character(len=*),	intent(in),	optional	::	filename	Path to CSSI space weather file
real(kind=dp),	intent(in),	optional	::	f107_override	Override value for F10.7 if not using space weather file
real(kind=dp),	intent(in),	optional	::	f107a_override	Override value for F10.7a if not using space weather file
real(kind=dp),	intent(in),	optional	::	kp_override	Override value for Kp index if not using space weather file

private subroutine jr_cleanup(me)

Clean up module resources

Arguments

Type	Intent	Optional	Attributes		Name
class(jacchia_roberts_type),	intent(inout)			::	me

public pure subroutine find_cstar_roots(c_star, tx, root1, root2, x_root, y_root)

Find the two real roots and the complex conjugate root pair of the degree-4 temperature-profile polynomial c_star(z).

Arguments

Type	Intent		Name
real(kind=dp),	intent(in)	::	c_star(5)	Quartic polynomial coefficients `c_star(1)..c_star(5)`
real(kind=dp),	intent(in)	::	tx	Temperature at 125 km (K), used for initial guesses
real(kind=dp),	intent(out)	::	root1	Larger real root (> 125 km)
real(kind=dp),	intent(out)	::	root2	Smaller real root (< 100 km)
real(kind=dp),	intent(out)	::	x_root	Real part of complex conjugate root pair
real(kind=dp),	intent(out)	::	y_root	Imaginary part of complex conjugate root pair

public pure subroutine find_cstar_roots_original(c_star, tx, root1, root2, x_root, y_root)

Original version of root-finding routine using the roots subroutine. This is based on the routine from GMAT.

Arguments

Type	Intent		Name
real(kind=dp),	intent(inout)	::	c_star(5)	Quartic polynomial coefficients `c_star(1)..c_star(5)` [note: modified in-place by deflation]
real(kind=dp),	intent(in)	::	tx	Temperature at 125 km (K) [not used here]
real(kind=dp),	intent(out)	::	root1	Larger real root (> 125 km)
real(kind=dp),	intent(out)	::	root2	Smaller real root (< 100 km)
real(kind=dp),	intent(out)	::	x_root	Real part of complex conjugate root pair
real(kind=dp),	intent(out)	::	y_root	Imaginary part of complex conjugate root pair

public subroutine prepare_flux_data(sw_data, flux_data, utc_mjd, kp_out, f107_out, f107a_out)

Prepare flux data for the current epoch, accounting for measurement timing. Kp selection matches the GMAT PrepareKpData historic-data branch:

Arguments

Type	Intent		Name
type(sw_data_type),	intent(inout)	::	sw_data	Space weather data manager
type(flux_data_type),	intent(in)	::	flux_data	Current day's flux data
real(kind=dp),	intent(in)	::	utc_mjd	Current epoch (MJD)
real(kind=dp),	intent(out)	::	kp_out	Selected Kp index for current epoch
real(kind=dp),	intent(out)	::	f107_out	Selected F10.7 daily value for current epoch
real(kind=dp),	intent(out)	::	f107a_out	Selected F10.7a average for current epoch

jacchia_roberts_module Module

Contents

Variables

Derived Types

Functions

Subroutines

Reference:

Uses

Variables

Derived Types

type, public :: geoparms_type

Components

type, public :: jacchia_roberts_type

Components

Type-Bound Procedures

Functions

private pure function parms_to_temp(me, f107_daily, f107a_avg, tkp) result(xtemp)

Arguments

Return Value real(kind=dp)

private function jacchia_roberts_density(me, height, position, sun_vector, geo_lat, utc_mjd) result(density)

Arguments

Return Value real(kind=dp)

private function exotherm(me, position, sun_vector, geo, height, sun_dec, geo_lat) result(exotemp)

Arguments

Return Value real(kind=dp)

private function rho_100(me, height, temperature) result(density)

Arguments

Return Value real(kind=dp)

private function rho_125(me, height, temperature) result(density)

Arguments

Return Value real(kind=dp)

private function rho_high(me, height, temperature, t_500, sun_dec, geo_lat) result(density)

Arguments

Return Value real(kind=dp)

private pure function rho_cor(height, utc_mjd, geo_lat, geo) result(correction)

Arguments

Return Value real(kind=dp)

Subroutines

private subroutine jr_init(me, cb_polar_radius, status, filename, f107_override, f107a_override, kp_override)

Arguments

private subroutine jr_cleanup(me)

Arguments

public pure subroutine find_cstar_roots(c_star, tx, root1, root2, x_root, y_root)

Arguments

public pure subroutine find_cstar_roots_original(c_star, tx, root1, root2, x_root, y_root)

Arguments

public subroutine prepare_flux_data(sw_data, flux_data, utc_mjd, kp_out, f107_out, f107a_out)

Arguments