Program Code for Ion Production Rates in the Upper Atmosphere
/******************************************************************************
* This PL/1 subroutine calculates the ion production rates due to
* solar EUV radiation for the constituents H, He, N, O, N2, O2, NO.
* The data for the solar flux and ionization cross sections are largely
* based on publications by H.E. Hinteregger.
*
* Input parameters are:
* NC: number of constituents (=7)
* NSP: number of spectral regions (=38)
* SPEC: =1: solar minimum; =2: solar maximum
* QSPEC: multiplicator for solar spectrum
* OQST: ionisation level for oxygen (='OALL': all levels;
* or 'O4S', 'O2DO', 'O2PO, 'O4PE', 'O2PE')
* TAUAB: optical depth for absorption of solar radiation (vector of length NSP)
* NEUDEN: neutral densities (in cm-3)
* (vector of length NC, in the order N2,O,He,N,O2,H,NO)
*
* Output parameter is:
* QP: ion production rates (in cm-3sec-1)
* (vector of length NC, in the order N2,O,He,N,O2,H,NO)
*
*
* This program has worked well and correct for me in the past (using IBM PL/1 compilers),
* but I can't guarantee this for all compilers, especially if the code is translated
* into a different language.
* If there are any problems, please let me know and I mention the issue here.
* (for contact details see http://www.plasmaphysics.org.uk/feedback.htm ).
*
* Thomas Smid, April 2007
*********************************************************************************/
IONPROD: PROCEDURE(NC,NSP,SPEC,QSPEC,OQST,TAUAB,NEUDEN,QP) ;
DEFAULT RANGE(A:Z) FLOAT DECIMAL VALUE (DEC FLOAT(16)) ;
DCL
SOLSPEC1(38) INIT(0.44,0.17,1.87,1.40,0.51,0.08,1.36,0.60,7.76,
0.87,0.74,0.21,0.41,0.33,0.31,0.51,0.80,1.59,0.48,0.63,1.84,0.40,0.26,
0.39,0.17,0.20,0.24,0.78,0.87,1.93,4.43,3.70,4.84,1.71,4.37,1.94,2.48,
297.0) DEC FLOAT(16) ;
DCL
SOLSPEC2(38) INIT(1.37,0.47,5.70,7.14,1.08,5.73,12.16,4.69,14.40,6.83,
1.53,2.54,1.53,0.74,1.82,1.65,1.52,4.30,1.05,2.48,3.87,1.37,0.54,0.75,
0.43,0.44,1.19,1.51,2.45,4.85,12.22,9.85,10.22,4.08,11.85,6.10,6.09,
806.0) DEC FLOAT(16) ;
DCL
IQH(38) INIT(0.0065,0.036,0.10,0.20,0.33,0.46,0.40,0.55,0.55,0.90,
1.20,1.30,1.70,2.00,2.20,2.80,3.20,3.50,3.40,3.80,4.00,4.00,4.80,5.10,
5.20,5.60,5.60,5.60,5.60,5.50,4.00,0,0,0,0,0,0,0),
IQHE(38) INIT(0.21,0.44,1.00,1.67,2.16,2.67,2.53,3.05,3.05,3.59,4.35,
4.20,5.65,6.53,6.76,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
DEC FLOAT(16) ;
DCL
IQN(38) INIT(0.064,0.27,0.69,1.40,2.00,2.68,2.45,3.21,3.21,3.91,2.60,
2.70,3.30,4.00,4.50,5.20,5.50,5.90,5.75,2.75,2.85,2.80,3.00,3.10,3.20,
3.25,3.30,3.30,3.30,3.30,0,0,0,0,0,0,0,0) DEC FLOAT(16) ;
DCL
IQNO(38) INIT(0.38,2.22,4.05,5.88,7.02,8.04,7.71,8.74,8.74,9.55,11.12,
11.38,13.21,14.68,15.04,16.87,18.75,19.50,19.50,18.94,17.81,18.19,
15.38,12.75,10.69,7.69,7.50,8.06,7.50,9.38,14.63,18.75,12.94,12.94,
8.44,7.88,8.44,2.10) DEC FLOAT(16) ;
DCL
IQN2(38) INIT(0.60,2.32,5.40,8.15,9.65,10.60,10.08,11.58,11.60,14.60,
18.00,17.51,21.07,21.80,21.85,24.53,24.69,23.20,22.38,23.10,23.20,
23.22,25.06,23.00,23.20,23.77,18.39,10.18,16.75,0,0,0,0,0,0,0,0,0)
DEC FLOAT(16) ;
DCL
AQN2(38) INIT(0.60,2.32,5.40,8.15,9.65,10.60,10.08,11.58,11.60,
14.60,18.00,17.51,21.07,21.80,21.85,24.53,24.69,23.20,22.38,23.10,
23.20,23.22,29.75,26.30,30.94,35.46,26.88,19.26,30.71,15.05,46.63,
16.99,0.70,36.16,0,0,0,0) DEC FLOAT(16) ;
DCL
IQO4S(38) INIT(0.32,1.03,1.62,1.95,2.15,2.33,2.23,2.45,2.45,2.61,2.81,
2.77,2.99,3.15,3.28,3.39,3.50,3.58,3.46,3.67,3.74,3.73,4.04,4.91,4.20,
4.18,4.18,4.28,4.23,4.38,4.18,2.12,0,0,0,0,0,0) DEC FLOAT(16) ;
DCL
IQO2DO(38) INIT(0.34,1.14,2.00,2.62,3.02,3.39,3.18,3.62,3.63,3.98,
4.37,4.31,4.75,5.04,5.23,5.36,5.47,5.49,5.30,5.51,5.50,5.50,5.52,6.44,
3.80,0,0,0,0,0,0,0,0,0,0,0,0,0) DEC FLOAT(16) ;
DCL
IQO2PO(38) INIT(0.22,0.75,1.30,1.70,1.95,2.17,2.04,2.32,2.32,2.52,
2.74,2.70,2.93,3.06,3.13,3.15,3.16,3.10,3.02,3.05,2.98,2.97,0.47,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) DEC FLOAT(16) ;
DCL
IQO4PE(38) INIT(0.10,0.34,0.58,0.73,0.82,0.89,0.85,0.91,0.91,0.93,
0.92,0.92,0.55,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
DEC FLOAT(16) ;
DCL
IQO2PE(38) INIT(0.03,0.27,0.46,0.54,0.56,0.49,0.52,0.41,0.41,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
DEC FLOAT(16) ;
DCL
IQO2(38) INIT(1.19,3.00,7.09,10.3,12.8,14.8,14.2,16.0,16.0,17.2,
18.4,18.1,19.5,20.4,20.9,23.5,25.6,22.0,25.1,26.1,25.8,26.0,23.3,25.0,
24.3,8.5,11.0,9.3,9.2,7.4,4.4,12.1,5.0,14.6,1.0,0,0.7,0)
DEC FLOAT(16) ;
DCL
AQO2(38) INIT(1.19,3.00,7.09,10.3,12.8,14.8,14.2,16.0,16.0,17.2,18.4,
18.1,19.5,20.4,20.9,23.5,25.6,22.0,25.1,26.1,25.8,26.0,25.6,27.0,30.9,
22.5,18.0,30.7,27.5,28.6,10.6,17.8,7.5,19.8,1.6,1.0,1.0,0.01)
REAL DEC FLOAT(16) ;
DCL NEUDEN(*),TAUAB(*),QP(*) CTL DEC FLOAT(16) ;
DCL(I,NC,NSP,SPEC) DEC FIXED(7,0) ;
DCL OQST CHARACTER(4) ;
ON UNDERFLOW ;
/* if arrays not correctly dimensioned return 1E50 as production rate */
IF DIM(IQH,1)=NSP&NC=7 THEN GO TO CONT ;
ELSE DO ;
QP=1E50 ;
GO TO ENDE ;
END ;
CONT:
IF SPEC=2 THEN SOLSPEC1=SOLSPEC2 ;
SOLSPEC1=SOLSPEC1*QSPEC ;
DO I=1 TO NSP ;
SOLSPEC1(I)=SOLSPEC1(I)*EXP(-TAUAB(I)) ;
END ;
IQN2=IQN2*SOLSPEC1 ;
IF OQST='OALL' THEN
IQO4S=(IQO4S+IQO2DO+IQO2PO+IQO4PE+IQO2PE)*SOLSPEC1 ;
IF OQST='O4S ' THEN IQO4S=IQO4S*SOLSPEC1 ;
IF OQST='O2DO' THEN IQO4S=IQO2DO*SOLSPEC1 ;
IF OQST='O2PO' THEN IQO4S=IQO2PO*SOLSPEC1 ;
IF OQST='O4PE' THEN IQO4S=IQO4PE*SOLSPEC1 ;
IF OQST='O2PE' THEN IQO4S=IQO2PE*SOLSPEC1 ;
IQHE=IQHE*SOLSPEC1 ;
IQN=IQN*SOLSPEC1 ;
IQO2=IQO2*SOLSPEC1 ;
IQH=IQH*SOLSPEC1 ;
IQNO=IQNO*SOLSPEC1 ;
QP(1)=SUM(IQN2)*NEUDEN(1) ;
QP(2)=SUM(IQO4S)*NEUDEN(2) ;
QP(3)=SUM(IQHE)*NEUDEN(3) ;
QP(4)=SUM(IQN)*NEUDEN(4) ;
QP(5)=SUM(IQO2)*NEUDEN(5) ;
QP(6)=SUM(IQH)*NEUDEN(6) ;
QP(7)=SUM(IQNO)*NEUDEN(7) ;
QP=QP*1E-09 ;
ENDE: END IONPROD ;
|