!*==cm4field.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE CM4FIELD(Path,Unit,Load,Indx,Gmut,Cord,Pred,Curr,Coef, & Nhmf,Nlmf,Ut,Mut,Theta,Phi,Alt,Dst,F107,Bmdl, & Jmdl,Gmdl,Perr,Oerr,Cerr) * ******************************************************************************** * * 1. MODULE NAME - CM4FIELD * * 2. NAME DERIVATION - Comprehensive Model 4 magnetic FIELD evaluation * subroutine * * 3. FUNCTION - To evaluate the Comprehensive Model 4 (CM4) magnetic field * model at a given temporal and spatial location with the given * magnetic indices. It returns the local X, Y, Z (North, East, * Down) components of the B field vector from the main, primary * and induced magnetosphere, primary and induced ionosphere, and * toroidal field sources. Two evaluations of the main field are * accommodated per two given spherical harmonic degree ranges. * This is helpful when only sub-ranges of the full degree/order * range is of interest. In addition, the J field vector * associated with the primary and induced ionosphere and the * toroidal B field sources is also returned in local X, Y, Z * components and coefficients from the various sources are also * computed and returned in various forms. For details into the * nature of models of this type see Sabaka et al. (2002,2004). * * 4. LANGUAGE - Fortran 77 * * 5. ARGUMENTS - * * Name Type I/O Description * ---- ---- --- ----------- * * PATH(3) Chr I File absolute paths: * (1) Model coefficients. * (2) Dst index. * (3) F10.7 index. * * UNIT(3) Int I Logical file units: * (1) Model coefficients. * (2) Dst index. * (3) F10.7 index. * * LOAD(3) Log I/O File read flags: * = T, read file * = F, file previously read (return value). * (1) Model coefficients. * (2) Dst index. * (3) F10.7 index. * * INDX(2) Log I Index acquisition flags: * = T, from file. * = F, from argument. * (1) Dst index (DST). * (2) F10.7 index (F107). * * GMUT Log I Magnetic dipole universal time (MUT) acquisition * flag: * = T, compute from UT. * = F, from argument. * * CORD Log I Observation coordinate system flag: * = T, geodetic, i.e., local geographic on a * reference ellipsoid. * = F, geocentric, i.e., local geographic on a * reference sphere. * * PRED(6) Log I Model B magnetic field prediction source flags: * = T, compute. * = F, do not compute. * (1) Main field. * (2) Magnetospheric fields. * (3) Ionospheric fields. * (4) Toroidal field from: * (5) Magsat (= T) or Oersted (= F). * (6) Dawn (= T) or dusk (= F) for Magsat. * * CURR Log I Model J current field prediction flag: * = T, compute. * = F, do not compute. * * J currents are only computed for certain external * sources activated by PRED (see JMDL). * * COEF Log I Model coefficient generation flag: * = T, compute. * = F, do not compute. * * Coefficients are only computed for sources * activated by PRED (see GMDL). * * NHMF(2) Int I Maximum main field spherical harmonic * degree: * (1) for main field 1. * (2) for main field 2. * * These values must be < or = 65. * * NLMF(2) Int I Minimum main field spherical harmonic * degree: * (1) for main field 1. * (2) for main field 2. * * If NHMF(i) < NLMF(i), then this degree * range is not evaluated. * * UT Dble I Geographic universal time (years A.D.): * GMUT = T, used for secular, seasonal and diurnal * time scale variations. * GMUT = F, used only for secular and seasonal time * scale variations. * * MUT Dble I/O Magnetic dipole universal time (hours 0-24) used * for diurnal time scale variations: * GMUT = T, compute from UT and return. * GMUT = F, from argument. * * THETA Dble I Geodetic or geocentric colatitude (degrees). * * PHI Dble I Longitude (degrees). * * ALT Dble I Geodetic or geocentric altitude (km). * * DST Dble I/O Linearly interpolated hourly Dst magnetic index * value (nT): * INDX(1) = T, from file and return. * INDX(1) = F, from argument. * * F107 Dble I/O Linearly interpolated 3-monthly means of absolute * F10.7 solar radiation flux value * (10^(-22)W/m^2/Hz): * INDX(2) = T, from file and return. * INDX(2) = F, from argument. * * BMDL(3,7) Dble O Array storing computed B field vectors from various * sources (nT): * * Row label: * (1) X. * (2) Y. * (3) Z. * * Column label: * (1) Main field 1. * (2) Main field 2. * (3) Primary magnetospheric field. * (4) Induced magnetospheric field. * (5) Primary ionospheric field. * (6) Induced ionospheric field. * (7) Toroidal field. * * Note that toroidal B is produced by an insitu * poloidal J field through which a satellite has * sampled. Therefore, it is only considered accurate * near its associated RTAY shell radius. * * JMDL(3,4) Dble O Array storing computed J field vectors from certain * external sources: * * Row label: * (1) X. * (2) Y. * (3) Z or current function Psi. * * Column label: * (1) Induced magnetospheric field. * (2) Primary ionospheric field. * (3) Induced ionospheric field. * (4) Poloidal field. * * Note that the primary magnetospheric J field is a * culmination of many J field sources which extend * over a wide volume. Hence, an equivalent sheet * current representation is not considered useful * in this case. However, the induced magnetospheric * J are confined, more or less, to a thin outer * shell of the Earth. Therefore, the J computed * here is an equivalent current representation, which * is toroidal and considered to flow in a spherical * shell and RM. Hence, this J has no Z component. * Therefore, the equivalent current function Psi (a * stream function), along whose contours J flows, is * placed in JMDL(3,1) for the induced magnetospheric * J field. Both J and Psi are evaluated only on * their shells, i.e., ALT is ignored. For more * information on equivalent currents and current * functions see Chapman and Bartels (1940). * * Note that since primary and induced ionospheric J * are equivalent currents, they are toroidal and * considered to flow in spherical shells at RM+HION * and RM, respectively, and so have no Z component. * Therefore, the equivalent current function Psi (a * stream function), along whose contours J flows, is * placed in JMDL(3,2) and JMDL(3,3) for the primary * and induced ionospheric J fields, respectively. * Both J and Psi are evaluated only on their shells, * i.e., ALT is ignored. * * Note that poloidal J is an insitu field through * which a satellite has sampled and thus gives rise * to a toroidal B field. Therefore, they are only * considered accurate near their associated RTAY * shell radius. For more information on poloidal J * fields in thin sampling shells see Olsen (1997). * * Note all equivalent currents J are in A/m and all * current functions Psi are in kA. Poloidal J is in * nA/m^2. * * GMDL(L,*) Dble O Array storing coefficients from various sources, * where all coefficients are with respect to Schmidt- * normalized spherical harmonics (see Langel, 1987): * * Row length L is the maximum number of coefficients * N from the sources selected by PRED, NHMF and NLMF, * such that: * (1) Main field 1: N=NSHC(NHMF(1),NLMF(1))*(NXOR+1), * (2) Main field 2: N=NSHC(NHMF(2),NLMF(2))*(NXOR+1), * (3) Magnetospheric field: N=IXMG/2, * (4) Ionospheric field: N=IXSQ/2, * (5) Toroidal field: N=IXTO/2, * where NSHC(NMAX,NMIN) is the number of spherical * harmonic coefficients corresponding to NMAX and * NMIN, and NXOR is the maximum order of B-splines * describing secular variation. * * Column length is also a function of PRED, NHMF and * NLMF such that the order and number of columns are: * (1) Main field 1: 1 col if corresponding N>0, * (2) Main field 2: 1 col if corresponding N>0, * (3) Magnetospheric field: 2 cols (primary,induced), * (4) Ionospheric field: 2 cols (primary,induced), * (5) Toroidal field: 1 col. * * Note that columns are allocated only if the source * is selected, and in the case of the main field, * the number of coefficients are also greater than * zero. This allows for quasi-minimal memory * allocation for GMDL. If COEF = F, the calling * routine may use a dummy argument for GMDL so that * no additional memory is needed. * * The form and order of the coefficients is as * follows: * (1) Main fields 1 and 2: * * do j = 0, NXOR * do n = NLMF(), NHMF() * do m = 0, n * g(n,m,j) = jth time derivative of * cos(m*phi) term at UT in * geographic coordinates * if (m .gt. 0) then * h(n,m,j) = jth time derivative of * sin(m*phi) term at UT in * geographic coordinates * end if * end do * end do * end do * * (2) Primary and induced magnetospheric fields: * * do n = 1, NXMG * do m = 0, min(n, MXMG) * q(n,m) = cos(m*phi) term at UT, MUT * and DST in dipole coordinates * if (m .gt. 0) then * s(n,m) = sin(m*phi) term at UT, MUT * and DST in dipole coordinates * end if * end do * end do * * (3) Primary and induced ionospheric fields: * * do n = 1, NXSQ * do m = 0, min(n, MXSQ) * c(n,m) = cos(m*phi) term at UT, MUT, ALT * and F10.7 in dipole coordinates * if (m .gt. 0) then * d(n,m) = sin(m*phi) term at UT, MUT, * ALT and F10.7 in dipole * coordinates * end if * end do * end do * * Note that the ALT dependence here is simply * binary; when ALT is above (below) HION, the * primary coefficients are internal (external). * Note also that c(n,m) and d(n,m) correspond * to the usual spherical harmonic functions, NOT * quasi-dipole functions. * * (4) Toroidal field: * * do n = 1, NXTO * do m = 0, min(n, MXTO) * u(n,m) = cos(m*phi) term at UT, MUT and * ALT in dipole coordinates * if (m .gt. 0) then * v(n,m) = sin(m*phi) term at UT, MUT * and ALT in dipole coordinates * end if * end do * end do * * Note that u(n,m) and v(n,m) correspond to the * usual spherical harmonic functions, NOT quasi- * dipole functions. * * PERR Int I Error message print flag: * = 0, do not print. * = 1, print. * * OERR Int I Logical unit for error message print. * * CERR Int I/O Error return code: * = 0, normal. * > 0 and < 50, warning. * > 50, fatal. * * 6. FILE FORMATS - * * a) Magnetic field model coefficients: * * The magnetic field model file is an ASCII file composed of 25 * logical records read from logical unit UMDL=UNIT(1): * * READ(UMDL,1000) LSMF,LPOS,LCMF * READ(UMDL,1000) NMXI * READ(UMDL,1001) EPCH,RE,RP,RM * READ(UMDL,1002) (BORD(J),J=1,LSMF) * READ(UMDL,1002) (BKNO(J),J=1,LSMF) * READ(UMDL,1001) (BKPO(J),J=1,LPOS) * READ(UMDL,1001) (GAMF(J),J=1,LCMF) * READ(UMDL,1000) LCMG,LSMG * READ(UMDL,1000) PBMG,PEMG,PSMG,NXMG,MXMG,IXMG * READ(UMDL,1001) CNMP,ENMP,OMGS,OMGD,RE,RP,RM * READ(UMDL,1001) (((GPMG(I,J,K),I=1,LCMG),J=1,LSMG),K=1,2) * READ(UMDL,1001) (((GSMG(I,J,K),I=1,LCMG),J=1,LSMG),K=1,2) * READ(UMDL,1000) LCSQ,LSSQ * READ(UMDL,1000) PBSQ,PESQ,PSSQ,NXSQ,MXSQ,IXSQ * READ(UMDL,1001) CNMP,ENMP,OMGS,OMGD,RE,RP,RM,HION * READ(UMDL,1001) (((GPSQ(I,J,K),I=1,LCSQ),J=1,LSSQ),K=1,2) * READ(UMDL,1001) ((GSSQ(I,J),I=1,LCSQ),J=1,LSSQ) * READ(UMDL,1000) LCTO_MG,LSTO_MG,LRTO_MG * READ(UMDL,1000) PBTO_MG,PETO_MG,NTAY_MG,PSTO,NXTO,MXTO,IXTO * READ(UMDL,1001) CNMP,ENMP,OMGS,OMGD,RE,RP,RM,RTAY_DW,RTAY_DK * READ(UMDL,1001) ((((GCTO_MG(I,J,K,L),I=1,LCTO_MG),J=1,LSTO_MG), * K=1,LRTO_MG),L=1,2) * READ(UMDL,1000) LCTO_OR,LSTO_OR,LRTO_OR * READ(UMDL,1000) PBTO_OR,PETO_OR,NTAY_OR,PSTO,NXTO,MXTO,IXTO * READ(UMDL,1001) CNMP,ENMP,OMGS,OMGD,RE,RP,RM,RTAY_OR * READ(UMDL,1001) (((GCTO_OR(I,J,K),I=1,LCTO_OR),J=1,LSTO_OR), * K=1,LRTO_OR) * * with formats: * * 1000 FORMAT(10I8) * 1001 FORMAT(4E22.15) * 1002 FORMAT(22I4) * * where: * * Name Type Description * ---- ---- ----------- * * LSMF Int Total number of main field (MF) SHC's. * * LPOS Int Total number of B-spline knot positions * for all MF SHC's. * * LCMF Int Total number of MF model parameters. * * NMXI Int Maximum SH degree of MF model. * * EPCH Dble Epoch of MF model (years). * * RE Dble MF equatorial Earth radius (km). * * RP Dble MF polar Earth radius (km). * * RM Dble MF mean Earth radius (km). * * BORD(LSMF) Int Array storing the order of the B-spline * basis (4=cubic) for each MF SHC (note: * zero value indicates a static SHC). * * BKNO(LSMF) Int Array storing the interior knot number * of the B-spline basis for each MF SHC. * * BKPO(LPOS) Dble Array storing the knot positions of the * B-spline basis (including both exterior * knots) for each MF SHC. * * GAMF(LCMF) Dble Array storing the MF model parameters. * * LCMG Int Total number of spatial/diurnal modulated * magnetospheric field (MG) model parameters. * * LSMG Int Total number of MG fourier seasonal modes. * * PBMG Int Minimum MG diurnal wavenumber. * * PEMG Int Maximum MG diurnal wavenumber. * * PSMG Int Maximum MG seasonal wavenumber. * * NXMG Int Maximum SH degree of MG model. * * MXMG Int Maximum SH order of MG model. * * IXMG Int Total number of MG SHC's. * * CNMP Dble MG colatitude of North magnetic * dipole (degree). * * ENMP Dble MG east longitude of North magnetic * dipole (degree). * * OMGS Dble MG seasonal fundamental angular * frequency (rads/year). * * OMGD Dble MG diurnal fundamental angular * frequency (rads/hr). * * RE Dble MG equatorial Earth radius (km). * * RP Dble MG polar Earth radius (km). * * RM Dble MG mean Earth radius (km). * * GPMG(LCMG,LSMG,2) Dble Array storing the MG primary model * parameters, where GPMG(i,j,1) is the ith * spatial/diurnal modulated coefficient of * the jth seasonal mode having no Dst * dependence, and GPMG(i,j,2) is the same, * but with Dst dependence. * * GSMG(LCMG,LSMG,2) Dble Same as GPMG, except the MG induced model * parameters. * * LCSQ Int Total number of spatial/diurnal modulated * ionospheric field (SQ) model parameters. * * LSSQ Int Total number of SQ fourier seasonal modes. * * PBSQ Int Minimum SQ diurnal wavenumber. * * PESQ Int Maximum SQ diurnal wavenumber. * * PSSQ Int Maximum SQ seasonal wavenumber. * * NXSQ Int Maximum SH degree of SQ model. * * MXSQ Int Maximum SH order of SQ model. * * IXSQ Int Total number of SQ SHC's. * * CNMP Dble SQ colatitude of North magnetic * dipole (degree). * * ENMP Dble SQ east longitude of North magnetic * dipole (degree). * * OMGS Dble SQ seasonal fundamental angular * frequency (rads/year). * * OMGD Dble SQ diurnal fundamental angular * frequency (rads/hr). * * RE Dble SQ equatorial Earth radius (km). * * RP Dble SQ polar Earth radius (km). * * RM Dble SQ mean Earth radius (km). * * HION Dble Height of ionosphere where equivalent * currents conceptually flow (km). * * GPSQ(LCSQ,LSSQ,2) Dble Array storing the SQ primary model * parameters, where GPSQ(i,j,1) is the ith * spatial/diurnal modulated coefficient of * the jth seasonal mode below RM+HION, and * GPSQ(i,j,2) is the same, but for above. * * GSSQ(LCSQ,LSSQ) Dble Same as GPSQ, except the SQ induced model * parameters, all of which are for below * RM+HION. * * LCTO_MG Int Total number of spatial/diurnal modulated * toroidal field (TO) model parameters. * * LSTO_MG Int Total number of TO fourier seasonal modes. * * LRTO_MG Int Total number of TO radial expansion * parameters. * * PBTO_MG Int Minimum TO diurnal wavenumber. * * PETO_MG Int Maximum TO diurnal wavenumber. * * NTAY_MG Int Degree of radial Taylor series expansion. * * PSTO Int Maximum TO seasonal wavenumber. * * NXTO Int Maximum SH degree of TO model. * * MXTO Int Maximum SH order of TO model. * * IXTO Int Total number of TO SHC's. * * CNMP Dble TO colatitude of North magnetic * dipole (degree). * * ENMP Dble TO east longitude of North magnetic * dipole (degree). * * OMGS Dble TO seasonal fundamental angular * frequency (rads/year). * * OMGD Dble TO diurnal fundamental angular * frequency (rads/hr). * * RE Dble TO equatorial Earth radius (km). * * RP Dble TO polar Earth radius (km). * * RM Dble TO mean Earth radius (km). * * RTAY_DW Dble Radial Taylor series reference radius * for Magsat dawn TO. * * RTAY_DK Dble Radial Taylor series reference radius * for Magsat dusk TO. * * GCTO_MG(LCTO_MG,... Dble Array storing the Magsat TO model * LSTO_MG,... parameters,where GCTO(i,j,k,1) is the ith * LRTO_MG,2) spatial/diurnal modulated coefficient of the * jth seasonal mode of the kth radial * expansion coefficient for dawn local time, * and GCTO(i,j,k,2) is the same, but for * dusk local time (note: since this model * reduced only Magsat dawn and dusk data, * the diurnal dependence is not continuous * but is sampled at only two points (dawn * and dusk) which are broken out of the * first dimension and placed in the fourth * dimension of the array). * * (Note: The remaining 4 records hold TO information for Oersted which * is analogous to the 4 previous records for Magsat. However, * the Oersted TO is continuous in diurnal time.) * * b) Dst index: * * The Dst index file is a WDC format ASCII file composed of one logical * record per day read from logical unit UDST=UNIT(2): * * READ(UDST,1000) YEAR,MONTH,DAY,(JDST(I),I=24) * " * " * " * * with format: * * 1000 FORMAT(3X,2I2,1X,I2,10X,24I4) * * where: * * Name Type Description * ---- ---- ----------- * * YEAR Int Year since 1900. * * MONTH Int Month of year. * * DAY Int Day of month. * * JDST(24) Int Hourly means of Dst magnetic index values * (nT) for this day such that JDST(i) is mean * value over (i-1)th hour, thus values are * centered on the half-hour. * * These files types may be retrieved in one-year segments from: * * ftp://ftp.ngdc.noaa.gov/STP/GEOMAGNETIC_DATA/INDICES/DST/ * * as: * * dstXXXX * * where XXXX is the year A.D. Note that if Dst is to be calculated at * times spanning year boundaries, then these one-year segment files * should simply be appended. Of course, it is assumed that the * concatenated file is in increasing time order with no gaps. * * c) F10.7 index: * * The F10.7 index file is a MONTHPLT.ABS format ASCII file composed of * one logical record per month read from logical unit UFIN=UNIT(3): * * READ(UFIN,1001) YEAR,MONTH,JF107 * " * " * " * * with format: * * 1000 FORMAT(I4,1X,I2,1X,I4) * * where: * * Name Type Description * ---- ---- ----------- * * YEAR Int Year A.D. * * MONTH Int Month of year. * * JF107 Int (3-monthly means) x 10 of absolute F10.7 * solar radiation flux values * (10^(-22)W/m^2/Hz). * * This file type may be retrieved from: * * ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/ * * as: * * MONTHPLT.ABS * * 7. REFERENCES - * * Chapman, S. and J. Bartels, 1940. Geomagnetism, Clarendon Press, Oxford. * * Langel, R.A., 1987. Main Field (Chap. 4), in Geomagnetism, J. Jacobs * (ed.), pp. 249-512, Academic Press, New York. * * Olsen, N., 1997. Ionospheric F region currents at middle and low * latitudes estimated from Magsat data, J. Geophys. Res., 102, * 4563-4576. * * Sabaka, T.J., Olsen, N. and R.A. Langel, 2002. A comprehensive model of * the quiet-time, near-Earth magnetic field: phase 3, Geophys. J. Int., * 151, 32-68. * * Sabaka, T.J., Olsen, N. and M.E. Purucker, 2004. Extending comprehensive * models of the Earth's magnetic field with Oersted and CHAMP data, * Geophys. J. Int., doi: 10.1111/j.1365-246X.2004.02421.x. * * 8. HISTORY - * * Version Date Author Change Description * ------- ---- ------ ------------------ * * 1.0 7/26/02 T.J. Sabaka / RITSS Initial version. * 2.0 10/22/02 T.J. Sabaka / RITSS (1) Changed some variable names * for readability, * (2) Input DST and F107 variables * from call line or from files, * (3) Input MUT variable from call * line or compute from UT, * (4) Output model coefficients in * GMDL variable. * * 9. CONTACTS - * * Name: Terence J. Sabaka * * Address: Goddard Space Flight Center * Code 921 * Greenbelt, MD 20771 * * Phone: (301) 614-6493 * * Email: sabaka@geomag.gsfc.nasa.gov * ******************************************************************************** * IMPLICIT NONE !*--CM4FIELD686 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Alt , bc , bkpo , Bmdl , cdip , cego , cemo , & cemp , clat , cnmp , cosp , cpol , cpsi , ctgo , & cthe , ctmo , ctmp , doy , drad , Dst DOUBLE PRECISION dstt , dstx , ecto , edip , elon , enmp , epch , & epmg , epol , epsq , esmg , essq , F107 , f107x , & fdoy , frho , frto , fsrf , fyr , gamf DOUBLE PRECISION gcto_mg , gcto_or , Gmdl , gpmg , gpsq , gsmg , & gssq , hion , hq , ht , hymg , hysq , hyto , & omgd , omgs , Phi , pleg , psiz , rcur , re DOUBLE PRECISION rion , rlgm , rm , ro , rp , rrgt , rse , rt , & rtay , rtay_dk , rtay_dw , rtay_or , ru , sego , & semo , semp , sinp , spsi , stgo , sthe DOUBLE PRECISION stmo , stmp , taud , taus , tdmg , tdsq , tdto , & Theta , thetas , trig , tsmg , tssq , tsto , Ut , & wb , ws , xd , xg , yd , yg DOUBLE PRECISION zd , zg INTEGER i , I8SSUM , idom , idoy , igen , ILEG , imoh , imol , & imon , IXMF , IXMG , IXSQ , IXTO , iyr , iyrh , iyrl , j , & jaft , jdom , jdoy INTEGER jf107 , jmjd , jmon , js , jy , jyr , k , l , lcmf , & lcmg , lcsq , lcto , lpos , lrto , lsmf , lsmg , lssq , & lsto , lum1 , lum2 INTEGER lum3 , lum4 , lum5 , lum6 , lum7 , MAXYRDST , MAXYRF107 , & mjdh , mjdl , mjdy , mmdl , msec , mt , mu , MXMG , MXSQ , & MXTO , mz , nimf , nmax INTEGER nmin , nobo , noff , noga , nohq , nomn , nomx , nopo , & nout , NSHX , nsm1 , nsm2 , NSMG , NSSQ , nsto , NSTO_MG , & NSTO_OR , nt , ntay , NTAY_MG INTEGER NTAY_OR , nu , NXKN , NXMF , NXMG , NXOR , NXPO , NXSQ , & NXTO , nygo , NYMF , NYMG , NYSQ , nyto , NYTO_MG , & NYTO_OR , nz C*** End of declarations inserted by SPAG SAVE LOGICAL Curr , Coef , Cord , Gmut , omdl INTEGER*8 Perr , Oerr , Cerr INTEGER csys , PBMG , PEMG , PSMG , PBSQ , PESQ , PSSQ , pbto INTEGER PBTO_MG , PBTO_OR , peto , PETO_MG , PETO_OR , PSTO , p DOUBLE PRECISION Mut C======================================================================= C Main Field potential expansion parameters C======================================================================= PARAMETER (NXMF=65,NYMF=8840,NXOR=4,NXKN=19,NXPO=12415) C======================================================================= C Magnetospheric and coupled Induction potential expansion parameters C======================================================================= PARAMETER (PBMG=0,PEMG=5,PSMG=2,NXMG=11,MXMG=6,IXMG=226) C======================================================================= C Sq and coupled Induction potential expansion parameters C======================================================================= PARAMETER (PBSQ=0,PESQ=4,PSSQ=2,NXSQ=60,MXSQ=12,IXSQ=2736) C======================================================================= C Toroidal scalar or stream function expansion parameters C======================================================================= PARAMETER (PBTO_MG=0,PETO_MG=0,PBTO_OR=0,PETO_OR=4) PARAMETER (PSTO=2,NXTO=60,MXTO=12,IXTO=2736) PARAMETER (NTAY_MG=1,NTAY_OR=1) C======================================================================= C Set maximum number of years for index table limits C======================================================================= PARAMETER (MAXYRDST=100,MAXYRF107=100) C======================================================================= PARAMETER (IXMF=NXMF*(NXMF+2)) PARAMETER (ILEG=(NXMF+1)*(NXMF+2)) PARAMETER (NSMG=2*PSMG+1) PARAMETER (NSSQ=2*PSSQ+1) PARAMETER (NSTO_MG=PSTO+1) PARAMETER (NSTO_OR=2*PSTO+1) PARAMETER (NYMG=(PEMG-PBMG+1)*IXMG) PARAMETER (NYSQ=(PESQ-PBSQ+1)*IXSQ) PARAMETER (NYTO_MG=(PETO_MG-PBTO_MG+1)*IXTO) PARAMETER (NYTO_OR=(PETO_OR-PBTO_OR+1)*IXTO) C======================================================================= CHARACTER*(*) Path(3) CHARACTER*12 buff107 LOGICAL Pred(6) , Indx(2) , Load(3) INTEGER*8 Unit(3) INTEGER idim(12) INTEGER jdst(24) INTEGER bord(IXMF) , bkno(IXMF) INTEGER*8 Nhmf(2) , Nlmf(2) INTEGER us(IXMF) DOUBLE PRECISION Jmdl(3,4) DIMENSION dstx(MAXYRDST*366,24) DIMENSION f107x(MAXYRF107,12) DIMENSION bkpo(NXPO) DIMENSION Bmdl(3,7) DIMENSION bc(28) , rlgm(15) , rrgt(9) , rse(9) DIMENSION tsmg(2*(PSMG+1)) , tssq(2*(PSSQ+1)) , tsto(2*(PSTO+1)) DIMENSION tdmg(2*(PEMG+1)) , tdsq(2*(PESQ+1)) DIMENSION tdto(2*(PETO_MG+PETO_OR+1)) DIMENSION Gmdl(*) DIMENSION gamf(NYMF) DIMENSION gpmg(NYMG,NSMG,2) , gsmg(NYMG,NSMG,2) DIMENSION gpsq(NYSQ,NSSQ,2) , gssq(NYSQ,NSSQ) DIMENSION gcto_mg(NYTO_MG,NSTO_MG,NTAY_MG+1,2) DIMENSION gcto_or(NYTO_OR,NSTO_OR,NTAY_OR+1) DIMENSION epmg(NYMG) , esmg(NYMG) , epsq(NYSQ) , essq(NYSQ) DIMENSION ecto(NYTO_MG+NYTO_OR) DIMENSION hymg(NYMG,6) , hysq(NYSQ,6) DIMENSION hyto(3*(NYTO_MG+NYTO_OR)) DIMENSION pleg(ILEG) , rcur(2*ILEG+4*NXMF) , trig(2*(NXMF+1)) DIMENSION ht(2*NYMF) , hq(6*NYMF) , ws(IXMF) , & wb(2*(2*NXOR+NXKN+2)) C======================================================================= DATA drad/0.017453292519943D0/ , dstt/0.D0/ C======================================================================= IF ( Cerr.LT.50 ) THEN IF ( Load(1) ) THEN Load(1) = .FALSE. OPEN (Unit(1),FILE=Path(1)) REWIND Unit(1) READ (Unit(1),99004) lsmf , lpos , lcmf READ (Unit(1),99004) lum1 READ (Unit(1),99005) epch , re , rp , rm READ (Unit(1),99006) (bord(j),j=1,lsmf) READ (Unit(1),99006) (bkno(j),j=1,lsmf) READ (Unit(1),99005) (bkpo(j),j=1,lpos) READ (Unit(1),99005) (gamf(j),j=1,lcmf) READ (Unit(1),99004) lcmg , lsmg READ (Unit(1),99004) lum1 , lum2 , lum3 , lum4 , lum5 , lum6 READ (Unit(1),99005) cnmp , enmp , omgs , omgd , re , rp , & rm READ (Unit(1),99005) (((gpmg(i,j,k),i=1,lcmg),j=1,lsmg),k=1, & 2) READ (Unit(1),99005) (((gsmg(i,j,k),i=1,lcmg),j=1,lsmg),k=1, & 2) READ (Unit(1),99004) lcsq , lssq READ (Unit(1),99004) lum1 , lum2 , lum3 , lum4 , lum5 , lum6 READ (Unit(1),99005) cnmp , enmp , omgs , omgd , re , rp , & rm , hion READ (Unit(1),99005) (((gpsq(i,j,k),i=1,lcsq),j=1,lssq),k=1, & 2) READ (Unit(1),99005) ((gssq(i,j),i=1,lcsq),j=1,lssq) READ (Unit(1),99004) lcto , lsto , lrto READ (Unit(1),99004) lum1 , lum2 , lum3 , lum4 , lum5 , & lum6 , lum7 READ (Unit(1),99005) cnmp , enmp , omgs , omgd , re , rp , & rm , rtay_dw , rtay_dk READ (Unit(1),99005) ((((gcto_mg(i,j,k,l),i=1,lcto),j=1,lsto & ),k=1,lrto),l=1,2) READ (Unit(1),99004) lcto , lsto , lrto READ (Unit(1),99004) lum1 , lum2 , lum3 , lum4 , lum5 , & lum6 , lum7 READ (Unit(1),99005) cnmp , enmp , omgs , omgd , re , rp , & rm , rtay_or READ (Unit(1),99005) (((gcto_or(i,j,k),i=1,lcto),j=1,lsto), & k=1,lrto) cpol = cnmp*drad epol = enmp*drad ctmp = DCOS(cpol) stmp = DSIN(cpol) cemp = DCOS(epol) semp = DSIN(epol) rion = rm + hion ENDIF iyr = INT(Ut) fyr = Ut - DBLE(iyr) doy = fyr*DBLE(366-MIN(1,MOD(iyr,4))) idoy = INT(doy) fdoy = doy - DBLE(idoy) idoy = idoy + 1 msec = NINT(fdoy*86400000.D0) CALL YDTOMJDX(iyr,idoy,mjdy,imon,idom,idim) IF ( Gmut ) CALL GETMUT2(cnmp,enmp,iyr,idoy,msec,Mut) csys = 1 IF ( Cord ) csys = 0 IF ( Indx(1) ) THEN IF ( Load(2) ) THEN Load(2) = .FALSE. OPEN (Unit(2),FILE=Path(2)) REWIND Unit(2) jaft = 0 10 READ (Unit(2),99001,END=20) jyr , jmon , jdom , jdst 99001 FORMAT (3X,2I2,1X,I2,10X,24I4) jyr = jyr + 1900 IF ( jyr.LT.1957 ) jyr = jyr + 100 CALL YMDTOMJD(jyr,jmon,jdom,jmjd,jdoy) IF ( jaft.EQ.0 ) THEN jaft = 1 mjdl = jmjd ENDIF DO j = 1 , 24 dstx(jmjd-mjdl+1,j) = DBLE(jdst(j)) ENDDO GOTO 10 20 mjdh = jmjd ENDIF CALL INTDST(MAXYRDST,mjdl,mjdh,mjdy,msec,dstx,Dst,Perr,Oerr, & Cerr) IF ( Cerr.GT.49 ) GOTO 100 ENDIF IF ( Indx(2) ) THEN IF ( Load(3) ) THEN Load(3) = .FALSE. OPEN (Unit(3),FILE=Path(3)) REWIND Unit(3) jaft = 0 30 READ (Unit(3),99002,END=40) buff107 99002 FORMAT (A12) IF ( buff107(10:12).NE.'---' ) THEN READ (buff107,99003) jyr , jmon , jf107 99003 FORMAT (I4,1X,I2,1X,I4) IF ( jaft.EQ.0 ) THEN jaft = 1 iyrl = jyr imol = jmon ENDIF f107x(jyr-iyrl+1,jmon) = DBLE(jf107)/10.D0 ENDIF GOTO 30 40 iyrh = jyr imoh = jmon ENDIF CALL INTF107(MAXYRF107,iyrl,imol,iyrh,imoh,iyr,imon,idom, & idim,msec,f107x,F107,Perr,Oerr,Cerr) IF ( Cerr.GT.49 ) GOTO 100 ENDIF CALL R8VSET(1,21,0.D0,Bmdl) CALL R8VSET(1,12,0.D0,Jmdl) clat = Theta*drad elon = Phi*drad IF ( Coef ) THEN nout = 1 nygo = 0 IF ( Pred(2) ) nygo = MAX(nygo,IXMG/2) IF ( Pred(3) ) nygo = MAX(nygo,IXSQ/2) IF ( Pred(4) ) nygo = MAX(nygo,IXTO/2) ENDIF IF ( Pred(1) ) THEN igen = 1 nmax = MAX(Nhmf(1),Nhmf(2)) nmin = MIN(Nlmf(1),Nlmf(2)) nobo = NSHX(nmin-1,1,nmin-1,0) nopo = I8SSUM(1,nobo,bkno) + 2*nobo CALL BFIELD(igen,nmax,0,nmin,1,nmax,0,0,0,0,csys,3,2,0,epch, & re,rp,rm,Ut,clat,elon,Alt,Dst,dstt,rse,nz,mz,ro, & thetas,us,us,bord(nobo+1),bkno(nobo+1), & bkpo(nopo+1),us,us,us,us,ws,us,gamf,bc,gamf, & pleg,rcur,trig,us,ws,ht,hq,hq,Perr,Oerr,Cerr) IF ( Cerr.GT.49 ) GOTO 100 nomn = NSHX(int(Nlmf(1)-1),1,int(Nlmf(1)-1),0) nomx = NSHX(int(Nhmf(1)),1,int(Nhmf(1)),0) noff = nomn - nobo nsm1 = DIM(nomx,nomn) noga = I8SSUM(1,nomn,bord) + I8SSUM(1,nomn,bkno) + nomn nohq = I8SSUM(nobo+1,noff,bord) + I8SSUM(nobo+1,noff,bkno) & + noff nimf = I8SSUM(nomn+1,nsm1,bord) + I8SSUM(nomn+1,nsm1,bkno) & + nsm1 CALL BLSGEN(nimf,nz,3,Bmdl(1,1),gamf(noga+1),hq(nohq+1)) IF ( Coef ) THEN nopo = I8SSUM(1,nomn,bkno) + 2*nomn CALL GETGMF(NXOR,nsm1,epch,Ut,wb,gamf(noga+1),Gmdl(nout), & bkno(nomn+1),bord(nomn+1),bkpo(nopo+1),Perr, & Oerr,Cerr) IF ( Cerr.GT.49 ) GOTO 100 ENDIF nomn = NSHX(int(Nlmf(2)-1),1,int(Nlmf(2)-1),0) nomx = NSHX(int(Nhmf(2)),1,int(Nhmf(2)),0) noff = nomn - nobo nsm2 = DIM(nomx,nomn) noga = I8SSUM(1,nomn,bord) + I8SSUM(1,nomn,bkno) + nomn nohq = I8SSUM(nobo+1,noff,bord) + I8SSUM(nobo+1,noff,bkno) & + noff nimf = I8SSUM(nomn+1,nsm2,bord) + I8SSUM(nomn+1,nsm2,bkno) & + nsm2 CALL BLSGEN(nimf,nz,3,Bmdl(1,2),gamf(noga+1),hq(nohq+1)) IF ( Coef ) THEN nygo = MAX(nygo,nsm1*(NXOR+1)) nygo = MAX(nygo,nsm2*(NXOR+1)) nout = nout + nygo*MIN(1,nsm1) nopo = I8SSUM(1,nomn,bkno) + 2*nomn CALL GETGMF(NXOR,nsm2,epch,Ut,wb,gamf(noga+1),Gmdl(nout), & bkno(nomn+1),bord(nomn+1),bkpo(nopo+1),Perr, & Oerr,Cerr) IF ( Cerr.GT.49 ) GOTO 100 nout = nout + nygo*MIN(1,nsm2) ENDIF ENDIF IF ( Pred(2) .OR. Pred(3) .OR. Pred(4) ) THEN IF ( .NOT.Pred(1) ) CALL GEOCEN(csys,re,rp,rm,Alt,clat,ro, & thetas,sthe,cthe) psiz = thetas - clat cpsi = DCOS(psiz) spsi = DSIN(psiz) rlgm(1) = -cpsi rlgm(2) = 0.D0 rlgm(3) = -spsi rlgm(4) = 0.D0 rlgm(5) = 1.D0 rlgm(6) = 0.D0 rlgm(7) = spsi rlgm(8) = 0.D0 rlgm(9) = -cpsi ctgo = DCOS(thetas) stgo = DSIN(thetas) cego = DCOS(elon) sego = DSIN(elon) rrgt(1) = ctgo*cego rrgt(2) = ctgo*sego rrgt(3) = -stgo rrgt(4) = -sego rrgt(5) = cego rrgt(6) = 0.D0 rrgt(7) = stgo*cego rrgt(8) = stgo*sego rrgt(9) = ctgo CALL RMERGE(rlgm,rrgt) rrgt(1) = cemp*ctmp rrgt(2) = -semp rrgt(3) = cemp*stmp rrgt(4) = semp*ctmp rrgt(5) = cemp rrgt(6) = semp*stmp rrgt(7) = -stmp rrgt(8) = 0.D0 rrgt(9) = ctmp CALL RMERGE(rlgm,rrgt) xg = stgo*cego yg = stgo*sego zg = ctgo xd = rrgt(1)*xg + rrgt(4)*yg + rrgt(7)*zg yd = rrgt(2)*xg + rrgt(5)*yg + rrgt(8)*zg zd = rrgt(3)*xg + rrgt(6)*yg + rrgt(9)*zg cdip = DACOS(zd) edip = DATAN2(yd,xd) ctmo = zd stmo = DSIN(cdip) cemo = DCOS(edip) semo = DSIN(edip) rrgt(1) = -cemo*ctmo rrgt(2) = -semo rrgt(3) = -cemo*stmo rrgt(4) = -semo*ctmo rrgt(5) = cemo rrgt(6) = -semo*stmo rrgt(7) = stmo rrgt(8) = 0.D0 rrgt(9) = -ctmo CALL RMERGE(rlgm,rrgt) taus = omgs*Ut taud = omgd*Mut ENDIF IF ( Pred(2) ) THEN igen = 1 CALL BFIELD(igen,NXMG,NXMG,1,1,MXMG,MXMG,0,0,0,1,3,0,0,epch, & re,rp,rm,Ut,cdip,edip,Alt,Dst,dstt,rse,nu,mu,ru, & thetas,us,us,us,us,ws,us,us,us,us,ws,us,gsmg,bc, & gsmg,pleg,rcur,trig,us,ws,ht,hq,hq,Perr,Oerr, & Cerr) IF ( Cerr.GT.49 ) GOTO 100 js = nu/2 jy = 1 CALL TRIGMP(PSMG,taus,tsmg) CALL TRIGMP(PEMG,taud,tdmg) DO p = PBMG , PEMG cosp = tdmg(1+p) sinp = tdmg(2+p+PEMG) CALL MPOTENT(NXMG,MXMG,nu,NYMG,cosp,sinp,hq(1),hymg(jy,1) & ) CALL MPOTENT(NXMG,MXMG,nu,NYMG,cosp,sinp,hq(js+1), & hymg(jy,4)) jy = jy + IXMG ENDDO bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 CALL MSEASON(PSMG,NSMG,NYMG,Dst,tsmg,epmg,gpmg) CALL BLSGEN(NYMG,NYMG,3,bc,epmg,hymg(1,4)) CALL LTRANS(1,1,bc,rlgm,Bmdl(1,3)) bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 CALL MSEASON(PSMG,NSMG,NYMG,Dst,tsmg,esmg,gsmg) CALL BLSGEN(NYMG,NYMG,3,bc,esmg,hymg(1,1)) CALL LTRANS(1,1,bc,rlgm,Bmdl(1,4)) IF ( Curr ) THEN bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 CALL JTBELOW(PBMG,PEMG,NXMG,MXMG,ro,rm,NYMG,hymg(1,1)) CALL BLSGEN(NYMG,NYMG,3,bc,esmg,hymg(1,1)) CALL LTRANS(1,1,bc,rlgm,Jmdl(1,1)) ENDIF IF ( Coef ) THEN CALL GETGXF(PBMG,PEMG,NXMG,MXMG,js,epmg,Gmdl(nout),tdmg) nout = nout + nygo CALL GETGXF(PBMG,PEMG,NXMG,MXMG,js,esmg,Gmdl(nout),tdmg) nout = nout + nygo ENDIF ENDIF IF ( Pred(3) ) THEN igen = 1 fsrf = 1.D0 + 0.01485D0*F107 IF ( ro.LT.rion ) THEN CALL BFIELD(igen,NXSQ,NXSQ,1,1,MXSQ,MXSQ,0,0,0,1,3,0,0, & epch,re,rp,rm,Ut,cdip,edip,Alt,Dst,dstt,rse, & nu,mu,ru,thetas,us,us,us,us,ws,us,us,us,us, & ws,us,gssq,bc,gssq,pleg,rcur,trig,us,ws,ht, & hq,hq,Perr,Oerr,Cerr) IF ( Cerr.GT.49 ) GOTO 100 js = nu/2 jy = 1 CALL TRIGMP(PSSQ,taus,tssq) CALL TRIGMP(PESQ,taud,tdsq) DO p = PBSQ , PESQ cosp = tdsq(1+p) sinp = tdsq(2+p+PESQ) CALL MPOTENT(NXSQ,MXSQ,nu,NYSQ,cosp,sinp,hq(1), & hysq(jy,1)) CALL MPOTENT(NXSQ,MXSQ,nu,NYSQ,cosp,sinp,hq(js+1), & hysq(jy,4)) jy = jy + IXSQ ENDDO bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 CALL ISEASON(PSSQ,NSSQ,NYSQ,fsrf,tssq,epsq,gpsq(1,1,1)) CALL BLSGEN(NYSQ,NYSQ,3,bc,epsq,hysq(1,4)) CALL LTRANS(1,1,bc,rlgm,Bmdl(1,5)) bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 CALL ISEASON(PSSQ,NSSQ,NYSQ,fsrf,tssq,essq,gssq(1,1)) CALL BLSGEN(NYSQ,NYSQ,3,bc,essq,hysq(1,1)) CALL LTRANS(1,1,bc,rlgm,Bmdl(1,6)) IF ( Curr ) THEN bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 CALL JTABOVE(PBSQ,PESQ,NXSQ,MXSQ,ro,rion,NYSQ, & hysq(1,4)) CALL BLSGEN(NYSQ,NYSQ,3,bc,epsq,hysq(1,4)) CALL LTRANS(1,1,bc,rlgm,Jmdl(1,2)) bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 CALL JTBELOW(PBSQ,PESQ,NXSQ,MXSQ,ro,rm,NYSQ,hysq(1,1)) CALL BLSGEN(NYSQ,NYSQ,3,bc,essq,hysq(1,1)) CALL LTRANS(1,1,bc,rlgm,Jmdl(1,3)) ENDIF IF ( Coef ) THEN CALL GETGXF(PBSQ,PESQ,NXSQ,MXSQ,js,epsq,Gmdl(nout), & tdsq) nout = nout + nygo CALL GETGXF(PBSQ,PESQ,NXSQ,MXSQ,js,essq,Gmdl(nout), & tdsq) nout = nout + nygo ENDIF ELSE CALL BFIELD(igen,NXSQ,0,1,1,MXSQ,0,0,0,0,1,3,0,0,epch,re, & rp,rm,Ut,cdip,edip,Alt,Dst,dstt,rse,nu,mu,ru, & thetas,us,us,us,us,ws,us,us,us,us,ws,us,gssq, & bc,gssq,pleg,rcur,trig,us,ws,ht,hq,hq,Perr, & Oerr,Cerr) IF ( Cerr.GT.49 ) GOTO 100 jy = 1 CALL TRIGMP(PSSQ,taus,tssq) CALL TRIGMP(PESQ,taud,tdsq) DO p = PBSQ , PESQ cosp = tdsq(1+p) sinp = tdsq(2+p+PESQ) CALL MPOTENT(NXSQ,MXSQ,nu,NYSQ,cosp,sinp,hq(1), & hysq(jy,1)) jy = jy + IXSQ ENDDO bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 CALL ISEASON(PSSQ,NSSQ,NYSQ,fsrf,tssq,epsq,gpsq(1,1,2)) CALL BLSGEN(NYSQ,NYSQ,3,bc,epsq,hysq(1,1)) CALL LTRANS(1,1,bc,rlgm,Bmdl(1,5)) bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 CALL ISEASON(PSSQ,NSSQ,NYSQ,fsrf,tssq,essq,gssq(1,1)) CALL BLSGEN(NYSQ,NYSQ,3,bc,essq,hysq(1,1)) CALL LTRANS(1,1,bc,rlgm,Bmdl(1,6)) IF ( Curr ) THEN bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 CALL JTBELOW(PBSQ,PESQ,NXSQ,MXSQ,ro,rion,NYSQ, & hysq(1,1)) CALL BLSGEN(NYSQ,NYSQ,3,bc,epsq,hysq(1,1)) CALL LTRANS(1,1,bc,rlgm,Jmdl(1,2)) bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 CALL JTBCONT(PBSQ,PESQ,NXSQ,MXSQ,rion,rm,NYSQ, & hysq(1,1)) CALL BLSGEN(NYSQ,NYSQ,3,bc,essq,hysq(1,1)) CALL LTRANS(1,1,bc,rlgm,Jmdl(1,3)) ENDIF IF ( Coef ) THEN CALL GETGXF(PBSQ,PESQ,NXSQ,MXSQ,nu,epsq,Gmdl(nout), & tdsq) nout = nout + nygo CALL GETGXF(PBSQ,PESQ,NXSQ,MXSQ,nu,essq,Gmdl(nout), & tdsq) nout = nout + nygo ENDIF ENDIF ENDIF IF ( Pred(4) ) THEN IF ( .NOT.(Pred(5)) ) THEN pbto = PBTO_OR peto = PETO_OR nyto = NYTO_OR nsto = NSTO_OR ntay = NTAY_OR rtay = rtay_or omdl = .TRUE. mmdl = 1 ELSEIF ( Pred(6) ) THEN pbto = PBTO_MG peto = PETO_MG nyto = NYTO_MG nsto = NSTO_MG ntay = NTAY_MG rtay = rtay_dw omdl = .FALSE. mmdl = 1 ELSE pbto = PBTO_MG peto = PETO_MG nyto = NYTO_MG nsto = NSTO_MG ntay = NTAY_MG rtay = rtay_dk omdl = .FALSE. mmdl = 2 ENDIF igen = 1 CALL BFIELD(igen,NXTO,0,1,1,MXTO,0,0,0,0,1,3,0,0,epch,re,rp, & rm,Ut,cdip,edip,0.D0,Dst,dstt,rse,nt,mt,rt, & thetas,us,us,us,us,ws,us,us,us,us,ws,us,gcto_mg, & bc,gcto_mg,pleg,rcur,trig,us,ws,ht,hq,hq,Perr, & Oerr,Cerr) IF ( Cerr.LE.49 ) THEN frto = rm/ro frho = (ro-rtay)/rm jy = 1 CALL TRIGMP(PSTO,taus,tsto) CALL TRIGMP(peto,taud,tdto) DO p = pbto , peto cosp = tdto(1+p) sinp = tdto(2+p+peto) CALL MSTREAM(NXTO,MXTO,nt,nyto,cosp,sinp,frto,hq(1), & hyto(jy)) jy = jy + IXTO ENDDO bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 IF ( omdl ) THEN CALL TSEARAD(omdl,PSTO,ntay,nsto,nyto,frho,tsto,ecto, & gcto_or(1,1,1)) ELSE CALL TSEARAD(omdl,PSTO,ntay,nsto,nyto,frho,tsto,ecto, & gcto_mg(1,1,1,mmdl)) ENDIF CALL BLSGEN(nyto,nyto,2,bc,ecto,hyto(1)) CALL LTRANS(1,1,bc,rlgm,Bmdl(1,7)) IF ( Coef ) CALL GETGXF(pbto,peto,NXTO,MXTO,nt,ecto, & Gmdl(nout),tdto) IF ( Curr ) THEN bc(1) = 0.D0 bc(2) = 0.D0 bc(3) = 0.D0 CALL JPOLOID(pbto,peto,NXTO,MXTO,ro,rm,nt,nyto,tdto, & hq(1),hyto(1)) CALL BLSGEN(nyto,nyto,1,bc(3),ecto,hyto(2*nyto+1)) IF ( omdl ) THEN CALL TSEARDR(omdl,PSTO,ntay,nsto,nyto,frho,tsto, & ecto,gcto_or(1,1,1)) ELSE CALL TSEARDR(omdl,PSTO,ntay,nsto,nyto,frho,tsto, & ecto,gcto_mg(1,1,1,mmdl)) ENDIF CALL BLSGEN(nyto,nyto,2,bc(1),ecto,hyto(1)) CALL LTRANS(1,1,bc,rlgm,Jmdl(1,4)) ENDIF ENDIF ENDIF ENDIF 100 RETURN 99004 FORMAT (10I8) 99005 FORMAT (4E22.15) 99006 FORMAT (22I4) END !*==ymdtomjd.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE YMDTOMJD(Yearad,Month,Dayofmonth,Mjd,Dayofyear) IMPLICIT NONE !*--YMDTOMJD1289 C*** Start of declarations inserted by SPAG INTEGER Mjd , Month C*** End of declarations inserted by SPAG INTEGER Yearad , Dayofmonth , Dayofyear , remyear INTEGER daysuptomonth(12) DATA daysuptomonth/0 , 31 , 59 , 90 , 120 , 151 , 181 , 212 , & 243 , 273 , 304 , 334/ remyear = Yearad - 1900 IF ( remyear.GT.0 ) THEN remyear = remyear - 1 Mjd = 1461*(remyear/4) + 15384 remyear = MOD(remyear,4) Dayofyear = daysuptomonth(Month) + Dayofmonth IF ( Month.GT.2 ) Dayofyear = Dayofyear + DIM(remyear,2) Mjd = Mjd + 365*remyear + Dayofyear ELSE Dayofyear = daysuptomonth(Month) + Dayofmonth Mjd = 15019 + Dayofyear ENDIF END !*==ydtomjdx.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE YDTOMJDX(Yearad,Dayofyear,Mjd,Month,Dayofmonth, & Daysinmonth) IMPLICIT NONE !*--YDTOMJDX1314 C*** Start of declarations inserted by SPAG INTEGER j , leapcum , leapday , Mjd , Month C*** End of declarations inserted by SPAG INTEGER Yearad , Dayofmonth , Dayofyear , remyear INTEGER daysuptomonth(12) , Daysinmonth(12) DATA daysuptomonth/0 , 31 , 59 , 90 , 120 , 151 , 181 , 212 , & 243 , 273 , 304 , 334/ remyear = Yearad - 1900 IF ( remyear.GT.0 ) THEN remyear = remyear - 1 Mjd = 1461*(remyear/4) + 15384 remyear = MOD(remyear,4) Mjd = Mjd + 365*remyear + Dayofyear leapday = DIM(remyear,2) ELSE Mjd = 15019 + Dayofyear leapday = 0 ENDIF DO j = 12 , 1 , -1 leapcum = MIN(1,DIM(j,2))*leapday IF ( daysuptomonth(j)+leapcum.LT.Dayofyear ) THEN Month = j Dayofmonth = Dayofyear - daysuptomonth(j) - leapcum GOTO 100 ENDIF ENDDO 100 Daysinmonth(1) = 31 Daysinmonth(2) = 28 + leapday Daysinmonth(3) = 31 Daysinmonth(4) = 30 Daysinmonth(5) = 31 Daysinmonth(6) = 30 Daysinmonth(7) = 31 Daysinmonth(8) = 31 Daysinmonth(9) = 30 Daysinmonth(10) = 31 Daysinmonth(11) = 30 Daysinmonth(12) = 31 END !*==intdst.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE INTDST(Mxyr,Mjdl,Mjdh,Mjdy,Msec,Dstx,Dst,Perr,Oerr, & Cerr) IMPLICIT NONE !*--INTDST1358 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Dst , Dstx , ttop INTEGER jbot , jtop , Mjdh , Mjdl , mjdt , Mjdy , mobs , Msec , & Mxyr C*** End of declarations inserted by SPAG INTEGER*8 Perr , Oerr , Cerr INTEGER hour , htop , hbot DIMENSION Dstx(Mxyr*366,24) hour = Msec/3600000 mjdt = Mjdy + hour/24 hour = MOD(hour,24) + 1 mobs = MOD(Msec,3600000) IF ( mobs.LE.1800000 ) THEN ttop = DBLE(mobs+1800000)/3600000.D0 IF ( hour.GT.1 ) THEN jtop = mjdt jbot = mjdt htop = hour hbot = hour - 1 ELSE jtop = mjdt jbot = mjdt - 1 htop = 1 hbot = 24 ENDIF ELSE ttop = DBLE(mobs-1800000)/3600000.D0 IF ( hour.LT.24 ) THEN jtop = mjdt jbot = mjdt htop = hour + 1 hbot = hour ELSE jtop = mjdt + 1 jbot = mjdt htop = 1 hbot = 24 ENDIF ENDIF IF ( (jbot.LT.Mjdl) .OR. (jtop.GT.Mjdh) ) THEN Cerr = 50 IF ( Perr.NE.0 ) WRITE (Oerr,99001) Cerr 99001 FORMAT (' ',/,1X,'SUBROUTINE INTDST -- ERROR CODE ',I2, & ' -- T LIES OUTSIDE OF DST TABLE TIME SPAN -- ABORT',/, & ' ') Dst = -1.D12 ELSE Dst = ttop*Dstx(jtop-Mjdl+1,htop) + (1.D0-ttop) & *Dstx(jbot-Mjdl+1,hbot) ENDIF RETURN END !*==intf107.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE INTF107(Mxyr,Iyrl,Imol,Iyrh,Imoh,Iyr,Imon,Idom,Idim, & Msec,F107x,F107,Perr,Oerr,Cerr) IMPLICIT NONE !*--INTF1071415 C*** Start of declarations inserted by SPAG DOUBLE PRECISION F107 , F107x , tadd , thlf , tint , tobs , ttop INTEGER Idom , Imoh , Imol , Imon , Iyr , Iyrh , Iyrl , mbot , & Msec , mtop , Mxyr C*** End of declarations inserted by SPAG INTEGER*8 Perr , Oerr , Cerr INTEGER ytop , ybot INTEGER Idim(12) DIMENSION F107x(Mxyr,12) thlf = 0.5D0*DBLE(Idim(Imon)) tobs = DBLE(Idom-1) + DBLE(Msec)/86400000.D0 IF ( tobs.LE.thlf ) THEN IF ( Imon.GT.1 ) THEN tadd = 0.5D0*DBLE(Idim(Imon-1)) tobs = tobs + tadd tint = thlf + tadd ttop = tobs/tint ytop = Iyr ybot = Iyr mtop = Imon mbot = Imon - 1 ELSE tobs = tobs + 15.5D0 tint = thlf + 15.5D0 ttop = tobs/tint ytop = Iyr ybot = Iyr - 1 mtop = 1 mbot = 12 ENDIF ELSEIF ( Imon.LT.12 ) THEN tadd = 0.5D0*DBLE(Idim(Imon+1)) tobs = tobs - thlf tint = thlf + tadd ttop = tobs/tint ytop = Iyr ybot = Iyr mtop = Imon + 1 mbot = Imon ELSE tobs = tobs - 15.5D0 tint = thlf + 15.5D0 ttop = tobs/tint ytop = Iyr + 1 ybot = Iyr mtop = 1 mbot = 12 ENDIF IF ( (ybot.LT.Iyrl) .OR. (ytop.GT.Iyrh) .OR. & ((ybot.EQ.Iyrl) .AND. (mbot.LT.Imol)) .OR. & ((ytop.EQ.Iyrh) .AND. (mtop.GT.Imoh)) ) THEN Cerr = 50 IF ( Perr.NE.0 ) WRITE (Oerr,99001) Cerr 99001 FORMAT (' ',/,1X,'SUBROUTINE INTF107 -- ERROR CODE ',I2, & ' -- T LIES OUTSIDE OF F10.7 TABLE TIME SPAN -- ABORT', & /,' ') F107 = -1.D0 ELSE F107 = ttop*F107x(ytop-Iyrl+1,mtop) + (1.D0-ttop) & *F107x(ybot-Iyrl+1,mbot) ENDIF RETURN END !*==getmut2.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE GETMUT2(Thenmp,Phinmp,Iyear,Iday,Msec,Mut) IMPLICIT NONE !*--GETMUT21483 C*** Start of declarations inserted by SPAG DOUBLE PRECISION cphd , cth0 , cths , dr , eadj , eopp , gmts , & gst , phi0 , Phinmp , phis , sdec , slong , & sphd , srasn , sth0 , sths , the0 , Thenmp , thes INTEGER Iday , Iyear , Msec C*** End of declarations inserted by SPAG DOUBLE PRECISION Mut DATA dr/0.017453292519943D0/ the0 = Thenmp*dr phi0 = Phinmp*dr cth0 = DCOS(the0) sth0 = DSIN(the0) gmts = DBLE(Msec)/1000.D0 CALL SUN2(Iyear,Iday,gmts,gst,slong,srasn,sdec) thes = (90.D0-sdec)*dr cths = DCOS(thes) sths = DSIN(thes) phis = (srasn-gst)*dr cphd = DCOS(phis-phi0) sphd = DSIN(phis-phi0) eopp = sths*sphd eadj = cth0*sths*cphd - sth0*cths Mut = 12.D0 - DATAN2(eopp,eadj)/dr/15.D0 END !*==sun2.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE SUN2(Iyr,Iday,Secs,Gst,Slong,Srasn,Sdec) IMPLICIT NONE !*--SUN21511 C*** Start of declarations inserted by SPAG DOUBLE PRECISION cosined , dj , fday , g , Gst , obliq , rad , & Sdec , Secs , sined , Slong , slp , Srasn , t , & vl INTEGER Iday , Iyr C*** End of declarations inserted by SPAG DATA rad/0.572957795130833D+02/ IF ( (Iyr.LT.1901) .OR. (Iyr.GT.2099) ) THEN Gst = 0.D0 Slong = 0.D0 Srasn = 0.D0 Sdec = 0.D0 ELSE fday = Secs/86400.D0 dj = 365.D0*DBLE(Iyr-1900) + DBLE((Iyr-1901)/4) + DBLE(Iday) & + fday - 0.5D0 t = dj/36525.D0 vl = DMOD(279.696678D0+0.9856473354D0*dj,360.D0) Gst = DMOD(279.690983D0+0.9856473354D0*dj+360.D0*fday+180.D0, & 360.D0) g = DMOD(358.475845D0+0.985600267D0*dj,360.D0)/rad Slong = vl + (1.91946D0-0.004789D0*t)*DSIN(g) & + 0.020094D0*DSIN(2.D0*g) obliq = (23.45229D0-0.0130125D0*t)/rad slp = (Slong-0.005686D0)/rad sined = DSIN(obliq)*DSIN(slp) cosined = DSQRT(1.D0-sined**2) Sdec = rad*DATAN(sined/cosined) Srasn = 180.D0 - rad*DATAN2(1.D0/DTAN(obliq)*sined/cosined, & -COS(slp)/cosined) ENDIF END !*==rmerge.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE RMERGE(Rmrg,Rmlt) IMPLICIT NONE !*--RMERGE1547 C*** Start of declarations inserted by SPAG DOUBLE PRECISION r1 , r2 , r3 , Rmlt , Rmrg C*** End of declarations inserted by SPAG SAVE DIMENSION Rmrg(*) , Rmlt(*) r1 = Rmrg(1) r2 = Rmrg(2) r3 = Rmrg(3) Rmrg(1) = r1*Rmlt(1) + r2*Rmlt(4) + r3*Rmlt(7) Rmrg(2) = r1*Rmlt(2) + r2*Rmlt(5) + r3*Rmlt(8) Rmrg(3) = r1*Rmlt(3) + r2*Rmlt(6) + r3*Rmlt(9) r1 = Rmrg(4) r2 = Rmrg(5) r3 = Rmrg(6) Rmrg(4) = r1*Rmlt(1) + r2*Rmlt(4) + r3*Rmlt(7) Rmrg(5) = r1*Rmlt(2) + r2*Rmlt(5) + r3*Rmlt(8) Rmrg(6) = r1*Rmlt(3) + r2*Rmlt(6) + r3*Rmlt(9) r1 = Rmrg(7) r2 = Rmrg(8) r3 = Rmrg(9) Rmrg(7) = r1*Rmlt(1) + r2*Rmlt(4) + r3*Rmlt(7) Rmrg(8) = r1*Rmlt(2) + r2*Rmlt(5) + r3*Rmlt(8) Rmrg(9) = r1*Rmlt(3) + r2*Rmlt(6) + r3*Rmlt(9) END !*==tsearad.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE TSEARAD(Full,Ks,Kr,Ns,Ng,F,T,E,G) IMPLICIT NONE !*--TSEARAD1575 C*** Start of declarations inserted by SPAG DOUBLE PRECISION E , F , G , s , T , z INTEGER i , j , k , Kr , Ks , Ng , Ns C*** End of declarations inserted by SPAG SAVE LOGICAL Full DIMENSION E(Ng) , G(Ng,Ns,Kr+1) , T(*) CALL R8VSET(1,Ng,0.D0,E) j = 1 s = 1.D0 CALL R8VLINKT(1,1,Ng,s,G(1,j,1),E) DO i = 1 , Ks j = j + 1 s = T(i+1) CALL R8VLINKT(1,1,Ng,s,G(1,j,1),E) IF ( Full ) THEN j = j + 1 s = T(i+Ks+2) CALL R8VLINKT(1,1,Ng,s,G(1,j,1),E) ENDIF ENDDO z = 1.D0 DO k = 1 , Kr j = 1 z = z*F/DBLE(k) CALL R8VLINKT(1,1,Ng,z,G(1,j,k+1),E) DO i = 1 , Ks j = j + 1 s = T(i+1)*z CALL R8VLINKT(1,1,Ng,s,G(1,j,k+1),E) IF ( Full ) THEN j = j + 1 s = T(i+Ks+2)*z CALL R8VLINKT(1,1,Ng,s,G(1,j,k+1),E) ENDIF ENDDO ENDDO END !*==tseardr.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE TSEARDR(Full,Ks,Kr,Ns,Ng,F,T,E,G) IMPLICIT NONE !*--TSEARDR1617 C*** Start of declarations inserted by SPAG DOUBLE PRECISION E , F , G , s , T , z INTEGER i , j , k , Kr , Ks , Ng , Ns C*** End of declarations inserted by SPAG SAVE LOGICAL Full DIMENSION E(Ng) , G(Ng,Ns,Kr+1) , T(*) CALL R8VSET(1,Ng,0.D0,E) z = 1.D0 DO k = 1 , Kr j = 1 CALL R8VLINKT(1,1,Ng,z,G(1,j,k+1),E) DO i = 1 , Ks j = j + 1 s = T(i+1)*z CALL R8VLINKT(1,1,Ng,s,G(1,j,k+1),E) IF ( Full ) THEN j = j + 1 s = T(i+Ks+2)*z CALL R8VLINKT(1,1,Ng,s,G(1,j,k+1),E) ENDIF ENDDO z = z*F/DBLE(k) ENDDO END !*==mseason.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE MSEASON(Ks,Ns,Ng,D,T,E,G) IMPLICIT NONE !*--MSEASON1646 C*** Start of declarations inserted by SPAG DOUBLE PRECISION D , E , G , s , T INTEGER i , j , Ks , Ng , Ns C*** End of declarations inserted by SPAG SAVE DIMENSION E(Ng) , G(Ng,Ns,2) , T(*) CALL R8VSET(1,Ng,0.D0,E) j = 1 s = 1.D0 CALL R8VLINKT(1,1,Ng,s,G(1,j,1),E) CALL R8VLINKT(1,1,Ng,D,G(1,j,2),E) DO i = 1 , Ks j = j + 1 s = T(i+1) CALL R8VLINKT(1,1,Ng,s,G(1,j,1),E) s = s*D CALL R8VLINKT(1,1,Ng,s,G(1,j,2),E) j = j + 1 s = T(i+Ks+2) CALL R8VLINKT(1,1,Ng,s,G(1,j,1),E) s = s*D CALL R8VLINKT(1,1,Ng,s,G(1,j,2),E) ENDDO END !*==iseason.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE ISEASON(Ks,Ns,Ng,F,T,E,G) IMPLICIT NONE !*--ISEASON1674 C*** Start of declarations inserted by SPAG DOUBLE PRECISION E , F , G , s , T INTEGER i , j , Ks , Ng , Ns C*** End of declarations inserted by SPAG SAVE DIMENSION E(Ng) , G(Ng,Ns) , T(*) CALL R8VSET(1,Ng,0.D0,E) j = 1 CALL R8VLINKT(1,1,Ng,F,G(1,j),E) DO i = 1 , Ks j = j + 1 s = F*T(i+1) CALL R8VLINKT(1,1,Ng,s,G(1,j),E) j = j + 1 s = F*T(i+Ks+2) CALL R8VLINKT(1,1,Ng,s,G(1,j),E) ENDDO END !*==mpotent.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE MPOTENT(Nmax,Mmax,Nd,Nz,Cphi,Sphi,D,Z) IMPLICIT NONE !*--MPOTENT1696 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Cphi , D , Sphi , Z INTEGER id , iz , m , Mmax , n , Nd , nd2 , Nmax , Nz C*** End of declarations inserted by SPAG SAVE DIMENSION D(*) , Z(Nz,*) nd2 = Nd + Nd id = 0 iz = 0 DO n = 1 , Nmax id = id + 1 iz = iz + 1 Z(iz,1) = D(id)*Cphi Z(iz,2) = D(Nd+id)*Cphi Z(iz,3) = D(nd2+id)*Cphi iz = iz + 1 Z(iz,1) = D(id)*Sphi Z(iz,2) = D(Nd+id)*Sphi Z(iz,3) = D(nd2+id)*Sphi DO m = 1 , MIN(n,Mmax) id = id + 2 iz = iz + 1 Z(iz,1) = D(id-1)*Cphi + D(id)*Sphi Z(iz,2) = D(Nd+id-1)*Cphi + D(Nd+id)*Sphi Z(iz,3) = D(nd2+id-1)*Cphi + D(nd2+id)*Sphi iz = iz + 1 Z(iz,1) = D(id)*Cphi - D(id-1)*Sphi Z(iz,2) = D(Nd+id)*Cphi - D(Nd+id-1)*Sphi Z(iz,3) = D(nd2+id)*Cphi - D(nd2+id-1)*Sphi iz = iz + 1 Z(iz,1) = D(id-1)*Cphi - D(id)*Sphi Z(iz,2) = D(Nd+id-1)*Cphi - D(Nd+id)*Sphi Z(iz,3) = D(nd2+id-1)*Cphi - D(nd2+id)*Sphi iz = iz + 1 Z(iz,1) = D(id)*Cphi + D(id-1)*Sphi Z(iz,2) = D(Nd+id)*Cphi + D(Nd+id-1)*Sphi Z(iz,3) = D(nd2+id)*Cphi + D(nd2+id-1)*Sphi ENDDO ENDDO END !*==jtabove.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE JTABOVE(Pmin,Pmax,Nmax,Mmax,R,Rref,Nz,Z) IMPLICIT NONE !*--JTABOVE1740 C*** Start of declarations inserted by SPAG DOUBLE PRECISION fcur , ffac , fgeo , fpsi , frmu , frpw , R , & Rref , Z , ztemp INTEGER iz , m , Mmax , n , Nmax , Nz C*** End of declarations inserted by SPAG SAVE INTEGER Pmin , Pmax , p DIMENSION Z(Nz,*) DATA fgeo/0.795774715459478D-03/ frmu = Rref/R iz = 0 DO p = Pmin , Pmax frpw = fgeo DO n = 1 , Nmax ffac = frpw*DBLE(n+n+1) fcur = ffac/DBLE(n+1) fpsi = -Rref*ffac/DBLE(n*(n+1)) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = -fcur*Z(iz,2) Z(iz,2) = fcur*ztemp Z(iz,3) = fpsi*Z(iz,3) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = -fcur*Z(iz,2) Z(iz,2) = fcur*ztemp Z(iz,3) = fpsi*Z(iz,3) DO m = 1 , MIN(n,Mmax) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = -fcur*Z(iz,2) Z(iz,2) = fcur*ztemp Z(iz,3) = fpsi*Z(iz,3) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = -fcur*Z(iz,2) Z(iz,2) = fcur*ztemp Z(iz,3) = fpsi*Z(iz,3) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = -fcur*Z(iz,2) Z(iz,2) = fcur*ztemp Z(iz,3) = fpsi*Z(iz,3) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = -fcur*Z(iz,2) Z(iz,2) = fcur*ztemp Z(iz,3) = fpsi*Z(iz,3) ENDDO frpw = frpw*frmu ENDDO ENDDO END !*==jtbelow.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE JTBELOW(Pmin,Pmax,Nmax,Mmax,R,Rref,Nz,Z) IMPLICIT NONE !*--JTBELOW1797 C*** Start of declarations inserted by SPAG DOUBLE PRECISION fcur , ffac , fgeo , fpsi , frbg , frmu , frpw , & R , Rref , Z , ztemp INTEGER iz , m , Mmax , n , Nmax , Nz C*** End of declarations inserted by SPAG SAVE INTEGER Pmin , Pmax , p DIMENSION Z(Nz,*) DATA fgeo/0.795774715459478D-03/ frmu = R/Rref frbg = fgeo*frmu**3 iz = 0 DO p = Pmin , Pmax frpw = frbg DO n = 1 , Nmax ffac = frpw*DBLE(n+n+1) fcur = ffac/DBLE(n) fpsi = -Rref*ffac/DBLE(n*(n+1)) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = fcur*Z(iz,2) Z(iz,2) = -fcur*ztemp Z(iz,3) = fpsi*Z(iz,3) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = fcur*Z(iz,2) Z(iz,2) = -fcur*ztemp Z(iz,3) = fpsi*Z(iz,3) DO m = 1 , MIN(n,Mmax) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = fcur*Z(iz,2) Z(iz,2) = -fcur*ztemp Z(iz,3) = fpsi*Z(iz,3) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = fcur*Z(iz,2) Z(iz,2) = -fcur*ztemp Z(iz,3) = fpsi*Z(iz,3) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = fcur*Z(iz,2) Z(iz,2) = -fcur*ztemp Z(iz,3) = fpsi*Z(iz,3) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = fcur*Z(iz,2) Z(iz,2) = -fcur*ztemp Z(iz,3) = fpsi*Z(iz,3) ENDDO frpw = frpw*frmu ENDDO ENDDO END !*==jtbcont.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE JTBCONT(Pmin,Pmax,Nmax,Mmax,Rold,Rnew,Nz,Z) IMPLICIT NONE !*--JTBCONT1855 C*** Start of declarations inserted by SPAG DOUBLE PRECISION fcur , fpsi , frbg , frmu , Rnew , Rold , Z INTEGER iz , m , Mmax , n , Nmax , Nz C*** End of declarations inserted by SPAG SAVE INTEGER Pmin , Pmax , p DIMENSION Z(Nz,*) frmu = Rold/Rnew frbg = frmu**2 iz = 0 DO p = Pmin , Pmax fpsi = frbg DO n = 1 , Nmax fcur = fpsi*frmu iz = iz + 1 Z(iz,1) = fcur*Z(iz,1) Z(iz,2) = fcur*Z(iz,2) Z(iz,3) = fpsi*Z(iz,3) iz = iz + 1 Z(iz,1) = fcur*Z(iz,1) Z(iz,2) = fcur*Z(iz,2) Z(iz,3) = fpsi*Z(iz,3) DO m = 1 , MIN(n,Mmax) iz = iz + 1 Z(iz,1) = fcur*Z(iz,1) Z(iz,2) = fcur*Z(iz,2) Z(iz,3) = fpsi*Z(iz,3) iz = iz + 1 Z(iz,1) = fcur*Z(iz,1) Z(iz,2) = fcur*Z(iz,2) Z(iz,3) = fpsi*Z(iz,3) iz = iz + 1 Z(iz,1) = fcur*Z(iz,1) Z(iz,2) = fcur*Z(iz,2) Z(iz,3) = fpsi*Z(iz,3) iz = iz + 1 Z(iz,1) = fcur*Z(iz,1) Z(iz,2) = fcur*Z(iz,2) Z(iz,3) = fpsi*Z(iz,3) ENDDO fpsi = fpsi*frmu ENDDO ENDDO END !*==mstream.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE MSTREAM(Nmax,Mmax,Nd,Nz,Cphi,Sphi,Faor,D,Z) IMPLICIT NONE !*--MSTREAM1903 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Cphi , D , Faor , Sphi , Z INTEGER id , iz , m , Mmax , n , Nd , Nmax , Nz C*** End of declarations inserted by SPAG SAVE DIMENSION D(*) , Z(Nz,*) id = 0 iz = 0 DO n = 1 , Nmax id = id + 1 iz = iz + 1 Z(iz,1) = Faor*D(Nd+id)*Cphi Z(iz,2) = -Faor*D(id)*Cphi iz = iz + 1 Z(iz,1) = Faor*D(Nd+id)*Sphi Z(iz,2) = -Faor*D(id)*Sphi DO m = 1 , MIN(n,Mmax) id = id + 2 iz = iz + 1 Z(iz,1) = Faor*(D(Nd+id-1)*Cphi+D(Nd+id)*Sphi) Z(iz,2) = -Faor*(D(id-1)*Cphi+D(id)*Sphi) iz = iz + 1 Z(iz,1) = Faor*(D(Nd+id)*Cphi-D(Nd+id-1)*Sphi) Z(iz,2) = -Faor*(D(id)*Cphi-D(id-1)*Sphi) iz = iz + 1 Z(iz,1) = Faor*(D(Nd+id-1)*Cphi-D(Nd+id)*Sphi) Z(iz,2) = -Faor*(D(id-1)*Cphi-D(id)*Sphi) iz = iz + 1 Z(iz,1) = Faor*(D(Nd+id)*Cphi+D(Nd+id-1)*Sphi) Z(iz,2) = -Faor*(D(id)*Cphi+D(id-1)*Sphi) ENDDO ENDDO END !*==jpoloid.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE JPOLOID(Pmin,Pmax,Nmax,Mmax,R,Rm,Nd,Nz,T,D,Z) IMPLICIT NONE !*--JPOLOID1940 C*** Start of declarations inserted by SPAG DOUBLE PRECISION cosp , D , fdeg , fhtj , frtj , R , Rm , sinp , & T , u0 , Z , ztemp INTEGER id , iz , m , Mmax , n , Nd , nd2 , Nmax , Nz C*** End of declarations inserted by SPAG SAVE INTEGER Pmin , Pmax , p DIMENSION D(*) , T(*) , Z(Nz,*) DATA u0/0.12566370614359158D-02/ nd2 = Nd + Nd fhtj = 1.D0/Rm/u0 frtj = Rm/R**2/u0 iz = 0 DO p = Pmin , Pmax id = 0 cosp = T(1+p) sinp = T(2+p+Pmax) DO n = 1 , Nmax fdeg = frtj*DBLE(n) id = id + 1 iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = fhtj*Z(iz,2) Z(iz,2) = -fhtj*ztemp Z(iz,3) = fdeg*D(nd2+id)*cosp iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = fhtj*Z(iz,2) Z(iz,2) = -fhtj*ztemp Z(iz,3) = fdeg*D(nd2+id)*sinp DO m = 1 , MIN(n,Mmax) id = id + 2 iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = fhtj*Z(iz,2) Z(iz,2) = -fhtj*ztemp Z(iz,3) = fdeg*(D(nd2+id-1)*cosp+D(nd2+id)*sinp) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = fhtj*Z(iz,2) Z(iz,2) = -fhtj*ztemp Z(iz,3) = fdeg*(D(nd2+id)*cosp-D(nd2+id-1)*sinp) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = fhtj*Z(iz,2) Z(iz,2) = -fhtj*ztemp Z(iz,3) = fdeg*(D(nd2+id-1)*cosp-D(nd2+id)*sinp) iz = iz + 1 ztemp = Z(iz,1) Z(iz,1) = fhtj*Z(iz,2) Z(iz,2) = -fhtj*ztemp Z(iz,3) = fdeg*(D(nd2+id)*cosp+D(nd2+id-1)*sinp) ENDDO ENDDO ENDDO END !*==blsgen.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE BLSGEN(Nc,Nd,Ni,B,C,Dldc) IMPLICIT NONE !*--BLSGEN2000 C*** Start of declarations inserted by SPAG DOUBLE PRECISION B , C , Dldc , R8SDOT INTEGER j , Nc , Nd , Ni C*** End of declarations inserted by SPAG DIMENSION B(*) , C(*) , Dldc(Nd,*) DO j = 1 , Ni B(j) = B(j) + R8SDOT(1,1,Nc,Dldc(1,j),C) ENDDO END !*==getgmf.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE GETGMF(Nder,Ns,Ep,Tm,B,C,G,H,O,P,Perr,Oerr,Cerr) IMPLICIT NONE !*--GETGMF2013 C*** Start of declarations inserted by SPAG DOUBLE PRECISION B , bsum , C , Ep , G , gd , gm , P , Tm INTEGER i , ic , ig , ik , im , j , jc , jm , js , k , la , lb , & m , n , nb , Nder , ni , nk , no , Ns INTEGER null , nw C*** End of declarations inserted by SPAG INTEGER*8 Perr , Oerr , Cerr INTEGER d , y INTEGER H(*) , O(*) DIMENSION B(*) , C(*) , G(*) , P(*) DATA null/0/ ic = 1 ik = 1 DO i = 1 , Ns G(i) = C(ic) ig = i DO j = 1 , Nder ig = ig + Ns G(ig) = 0.D0 ENDDO n = O(i) IF ( n.GT.0 ) THEN k = H(i) IF ( (Tm.LT.P(ik)) .OR. (Tm.GT.P(ik+k+1)) ) THEN Cerr = 56 IF ( Perr.NE.0 ) WRITE (Oerr,99001) Cerr 99001 FORMAT (' ',/,1X,'SUBROUTINE GETGMF -- ERROR CODE ',I2, & ' -- T LIES OUTSIDE OF KNOT DOMAIN -- ABORT',/, & ' ') GOTO 100 ENDIF m = n + 1 nk = k + 2 nb = n + k ni = nb + 1 no = ni + 1 nw = no + no CALL R8SLT(1,nk,Ep,P(ik),y) la = MIN(nk,y+1) CALL R8SLT(1,nk,Tm,P(ik),y) lb = MIN(nk,y+1) CALL R8VSET(1,nw,0.D0,B) CALL DBSPLN(la,Ep,m,null,k,P(ik),B(la-1),B(nw+1)) CALL DBSPLN(lb,Tm,m,null,k,P(ik),B(no+lb-1),B(nw+1)) CALL R8VSUB(1,no+1,1,no,B,B,B) bsum = 0.D0 gm = 0.D0 jc = ic + nb DO j = ni , 2 , -1 bsum = bsum + B(j) im = MIN(j,nk) jm = MAX(j-n,1) gm = gm + (P(ik+im-1)-P(ik+jm-1))*bsum*C(jc) jc = jc - 1 ENDDO G(i) = G(i) + gm/DBLE(n) js = ic + lb - 1 ig = i DO d = 0 , Nder - 1 ig = ig + Ns CALL DBSPLN(lb,Tm,n,d,k,P(ik),B(1),B(nw+1)) gd = 0.D0 jc = js DO j = 1 , n gd = gd + B(j)*C(jc) jc = jc + 1 ENDDO G(ig) = gd ENDDO ic = ic + nb ik = ik + nk ENDIF ic = ic + 1 ENDDO 100 RETURN END !*==dbspln.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE DBSPLN(L,T,N,D,K,X,B,W) IMPLICIT NONE !*--DBSPLN2094 C*** Start of declarations inserted by SPAG DOUBLE PRECISION B , delx , difi , difj , fn , T , temp , W , X INTEGER i , ik , iw , j , jk , K , L , lenb , m , N C*** End of declarations inserted by SPAG INTEGER posb , posx , posw , D , q , y DIMENSION B(*) , X(*) , W(*) q = N - D IF ( q.EQ.1 ) THEN B(1) = 1.D0 ELSE m = L ik = MIN(m,K+2) jk = MAX(m-1,1) difi = X(ik) - T delx = X(ik) - X(jk) IF ( delx.EQ.0.D0 ) THEN B(q) = 0.D0 ELSE B(q) = 1.D0/delx ENDIF posb = q - 1 DO j = 2 , q jk = MAX(m-j,1) temp = difi*B(posb+1) delx = X(ik) - X(jk) IF ( delx.EQ.0.D0 ) THEN temp = 0.D0 ELSEIF ( j.LT.N ) THEN temp = temp/delx ENDIF B(posb) = temp posb = posb - 1 ENDDO B(q+1) = 0.D0 m = m + 1 DO i = 2 , q ik = MIN(m,K+2) difi = X(ik) - T posb = q DO j = i , q jk = MAX(m-j,1) difj = T - X(jk) temp = difj*B(posb) + difi*B(posb+1) delx = X(ik) - X(jk) IF ( delx.EQ.0.D0 ) THEN temp = 0.D0 ELSEIF ( j.LT.N ) THEN temp = temp/delx ENDIF B(posb) = temp posb = posb - 1 ENDDO m = m + 1 ENDDO ENDIF posx = L + N - 1 posw = D + N DO i = 1 , N lenb = MIN(q,posw-D) CALL R8VSET(1,posw,0.D0,W) CALL R8VGATHP(1,1,D+1,lenb,B,W) DO j = 1 , D ik = posx jk = posx - q - j iw = posw fn = DBLE(q+j-1) DO m = j , D IF ( j.LT.D ) THEN delx = X(MAX(1,MIN(ik,K+2))) - X(MAX(jk,1)) IF ( delx.EQ.0.D0 ) THEN W(iw) = 0.D0 ELSE W(iw) = fn*(W(iw-1)-W(iw))/delx ENDIF ELSE W(iw) = fn*(W(iw-1)-W(iw)) ENDIF ik = ik - 1 jk = jk - 1 iw = iw - 1 ENDDO ENDDO posx = posx - 1 posw = posw - 1 ENDDO CALL R8VGATHP(D+1,1,1,N,W,B) END !*==getgxf.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE GETGXF(Pmin,Pmax,Nmax,Mmax,Ng,E,G,T) IMPLICIT NONE !*--GETGXF2185 C*** Start of declarations inserted by SPAG DOUBLE PRECISION cosp , E , G , sinp , T INTEGER ie , ig , m , Mmax , n , Ng , Nmax C*** End of declarations inserted by SPAG SAVE INTEGER Pmin , Pmax , p DIMENSION E(*) , G(*) , T(*) CALL R8VSET(1,Ng,0.D0,G) ie = 0 DO p = Pmin , Pmax ig = 1 cosp = T(1+p) sinp = T(2+p+Pmax) DO n = 1 , Nmax G(ig) = G(ig) + E(ie+1)*cosp + E(ie+2)*sinp ig = ig + 1 ie = ie + 2 DO m = 1 , MIN(n,Mmax) G(ig) = G(ig) + (E(ie+1)+E(ie+3)) & *cosp + (E(ie+4)-E(ie+2))*sinp ig = ig + 1 G(ig) = G(ig) + (E(ie+4)+E(ie+2)) & *cosp + (E(ie+1)-E(ie+3))*sinp ig = ig + 1 ie = ie + 4 ENDDO ENDDO ENDDO END !*==bfield.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE BFIELD(Rgen,Nmxi,Nmxe,Nmni,Nmne,Mmxi,Mmxe,Mmni,Mmne, & Grad,Ctyp,Dtyp,Ityp,Etyp,Ep,Re,Rp,Rm,Tm,Clat, & Elon,H,Dst,Dstt,Rse,Nc,Na,Ro,Theta,Atyp,Dsti, & Bori,Bkni,Bkpi,Tdgi,Dste,Bore,Bkne,Bkpe,Tdge,A, & B,C,P,R,T,U,W,Dsdc,Dldc,Dlda,Perr,Oerr,Cerr) IMPLICIT NONE !*--BFIELD2222 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , B , Bkpe , Bkpi , C , Clat , Dlda , Dldc , & Dsdc , Dst , Dstt , Elon , Ep , H , P , phi , R , & Re , Rm , Ro DOUBLE PRECISION Rp , Rse , T , Theta , Tm , W INTEGER ia , ie , ii , Ityp , mmax , mmin , Mmne , Mmni , Mmxe , & Mmxi , Na , Nc , nce , nci , nmax , Nmne , Nmni , Nmxe , & Nmxi , np INTEGER ns , nse , nsi C*** End of declarations inserted by SPAG SAVE INTEGER*8 Perr , Oerr , Cerr Integer Rgen , Ctyp , Dtyp , Etyp , Grad INTEGER Atyp(*) , Bore(*) , Bori(*) , Bkne(*) , Bkni(*) , Dste(*) & , Dsti(*) INTEGER Tdge(*) , Tdgi(*) , U(*) DIMENSION Bkpe(*) , Bkpi(*) , Dlda(*) , Dldc(*) , Dsdc(*) , Rse(*) & , A(*) , B(*) DIMENSION C(*) , P(*) , R(*) , T(*) , W(*) IF ( Cerr.LE.49 ) THEN phi = Elon CALL PREBF(Rgen,Ityp,Etyp,Dtyp,Grad,Nmni,Nmxi,Nmne,Nmxe,Mmni, & Mmxi,Mmne,Mmxe,nmax,mmin,mmax,ns,nsi,nse,Nc,nci,nce, & Na,np,ii,ie,Atyp,Dsti,Bori,Bkni,Tdgi,Dste,Bore,Bkne, & Tdge,U,Perr,Oerr,Cerr) IF ( Cerr.LT.50 ) THEN IF ( Nc.GT.0 ) THEN CALL FDSDC(Rgen,Ityp,Etyp,nsi,nse,Nc,nci,Ep,Tm,Dst,Dstt, & Dsti,Bori,Bkni,Bkpi,Tdgi,Dste,Bore,Bkne,Bkpe, & Tdge,U,W,Dsdc,Perr,Oerr,Cerr) IF ( Cerr.GE.50 ) GOTO 99999 CALL FDLDS(Rgen,Grad,Ctyp,Clat,phi,H,Re,Rp,Rm,Ro,nsi,Nc, & nci,np,ii,ie,Nmni,Nmxi,Nmne,Nmxe,nmax,Mmni, & Mmxi,Mmne,Mmxe,mmin,mmax,Theta,P,R,T,U,W,Dldc, & Perr,Oerr,Cerr) IF ( Cerr.GE.50 ) GOTO 99999 IF ( Rgen.GT.0 ) THEN Rgen = 0 CALL R8VSET(1,28+36*Grad,0.D0,B) CALL FDLDC(Grad,Nc,Dsdc,Dldc) CALL BLGEN(Grad,Nc,B,C,Dldc) CALL BNGEN(B) ENDIF IF ( Dtyp.EQ.2 ) THEN CALL TEC(Grad,Atyp(1),Nc,Theta,phi,B,Dldc,W) CALL TSE(Grad,Atyp(2),Nc,Rse,B,Dldc,W) ENDIF ENDIF CALL R8VGATHP(1,1,15,14,B,B) IF ( Grad.EQ.1 ) CALL R8VGATHP(29,1,47,18,B,B) ia = 0 IF ( Na.GT.0 ) THEN CALL R8VSET(1,6*Na,0.D0,Dlda) IF ( Dtyp.EQ.1 ) THEN CALL TBI(Atyp(1),Na,ia,A,B,Dlda) ELSEIF ( Dtyp.EQ.2 ) THEN CALL TMS(Grad,Atyp(3),Nc,Na,ia,A,B,Dldc,Dlda,W) CALL TNM(Grad,Atyp(4),Nc,Na,ia,A,B,Dldc,Dlda,W) CALL TVN(Grad,Atyp(5),Nc,Na,ia,A,B,Dldc,Dlda,W) CALL TBI(Atyp(6),Na,ia,A,B,Dlda) ENDIF ENDIF ENDIF ENDIF 99999 END !*==prebf.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE PREBF(Rgen,Ityp,Etyp,Dtyp,Grad,Nmni,Nmxi,Nmne,Nmxe, & Mmni,Mmxi,Mmne,Mmxe,Nmax,Mmin,Mmax,Ns,Nsi,Nse,Nc, & Nci,Nce,Na,Np,Ii,Ie,Atyp,Dsti,Bori,Bkni,Tdgi, & Dste,Bore,Bkne,Tdge,U,Perr,Oerr,Cerr) IMPLICIT NONE !*--PREBF2293 C*** Start of declarations inserted by SPAG INTEGER I8SSUM , idst , Ie , Ii , isvr , Ityp , Mmax , Mmin , & Mmne , Mmni , Mmxe , Mmxi , Na , Nc , Nce , Nci , NLPX , & Nmax , Nmne , Nmni INTEGER Nmxe , Nmxi , Np , Ns , Nse , NSHX , Nsi , nx C*** End of declarations inserted by SPAG SAVE INTEGER*8 Perr , Oerr , Cerr Integer Rgen , Dtyp , Etyp , esvr , edst , & Grad INTEGER Atyp(*) , Bore(*) , Bori(*) , Bkne(*) , Bkni(*) , Dste(*) & , Dsti(*) INTEGER Tdge(*) , Tdgi(*) , U(*) DATA nx/0/ IF ( Rgen.EQ.1 ) THEN IF ( MIN(Nmni,Nmxi,Nmne,Nmxe).LT.0 ) THEN Cerr = 50 IF ( Perr.NE.0 ) WRITE (Oerr,99001) Cerr 99001 FORMAT (' ',/,1X,'SUBROUTINE BFIELD -- ERROR CODE ',I2, & ' -- NMNI, NMXI, NMNE, OR NMXE < 0 -- ABORT',/,' ') GOTO 100 ENDIF IF ( MIN(Mmni,Mmxi,Mmne,Mmxe).LT.0 ) THEN Cerr = 51 IF ( Perr.NE.0 ) WRITE (Oerr,99002) Cerr 99002 FORMAT (' ',/,1X,'SUBROUTINE BFIELD -- ERROR CODE ',I2, & ' -- MMNI, MMXI, MMNE, OR MMXE < 0 -- ABORT',/,' ') GOTO 100 ENDIF IF ( (Mmni.GT.Mmxi) .OR. (Mmne.GT.Mmxe) ) THEN Cerr = 52 IF ( Perr.NE.0 ) WRITE (Oerr,99003) Cerr 99003 FORMAT (' ',/,1X,'SUBROUTINE BFIELD -- ERROR CODE ',I2, & ' -- EITHER MMNI > MMXI OR MMNE > MMXE -- ABORT',/, & ' ') GOTO 100 ENDIF IF ( (Mmxi.GT.Nmxi) .OR. (Mmxe.GT.Nmxe) ) THEN Cerr = 53 IF ( Perr.NE.0 ) WRITE (Oerr,99004) Cerr 99004 FORMAT (' ',/,1X,'SUBROUTINE BFIELD -- ERROR CODE ',I2, & ' -- EITHER MMXI > NMXI OR MMXE > NMXE -- ABORT',/, & ' ') GOTO 100 ENDIF isvr = MOD(Ityp,3) idst = Ityp/3 esvr = MOD(Etyp,3) edst = Etyp/3 Nmax = MAX(Nmxi,Nmxe) Mmin = MIN(Mmni,Mmne) Mmax = MAX(Mmxi,Mmxe) Nsi = NSHX(Nmxi,Nmni,Mmxi,Mmni) Nse = NSHX(Nmxe,Nmne,Mmxe,Mmne) Ns = Nsi + Nse Np = NLPX(Nmax,Mmax,Mmin) Ii = NLPX(Nmni-1,Mmax,Mmin) Ie = NLPX(Nmne-1,Mmax,Mmin) Nci = 0 IF ( Nsi.GT.0 ) THEN CALL I8VSET(1,Nsi,1,U) IF ( isvr.EQ.1 ) THEN CALL I8VADD(1,1,1,Nsi,Tdgi,U,U) ELSEIF ( isvr.EQ.2 ) THEN CALL I8VADD(1,1,1,Nsi,Bori,U,U) CALL I8VADD(1,1,1,Nsi,Bkni,U,U) ENDIF IF ( idst.EQ.1 ) CALL I8VADD(1,1,1,Nsi,Dsti,U,U) Nci = I8SSUM(1,Nsi,U) ENDIF Nce = 0 IF ( Nse.GT.0 ) THEN CALL I8VSET(Nsi+1,Nse,1,U) IF ( esvr.EQ.1 ) THEN CALL I8VADD(1,Nsi+1,Nsi+1,Nse,Tdge,U,U) ELSEIF ( esvr.EQ.2 ) THEN CALL I8VADD(1,Nsi+1,Nsi+1,Nse,Bore,U,U) CALL I8VADD(1,Nsi+1,Nsi+1,Nse,Bkne,U,U) ENDIF IF ( edst.EQ.1 ) CALL I8VADD(1,Nsi+1,Nsi+1,Nse,Dste,U,U) Nce = I8SSUM(Nsi+1,Nse,U) ENDIF Nc = Nci + Nce Rgen = 7 ENDIF Rgen = Rgen + nx Na = 0 nx = 0 IF ( Dtyp.EQ.1 ) THEN Na = Na + MIN(1,Atyp(1))*3 ELSEIF ( Dtyp.EQ.2 ) THEN nx = nx + Atyp(1) nx = nx + Atyp(2) Na = Na + MIN(1,Atyp(3))*3 Na = Na + MIN(1,Atyp(4))*3 Na = Na + MIN(1,Atyp(5))*3 nx = nx + Na Na = Na + MIN(1,Atyp(6))*3 ENDIF nx = MIN(1,nx) 100 RETURN END !*==fdlds.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE FDLDS(Rgen,Grad,Ctyp,Clat,Phi,H,Re,Rp,Rm,Ro,Nsi,Nc,Nci, & Np,Ii,Ie,Nmni,Nmxi,Nmne,Nmxe,Nmax,Mmni,Mmxi,Mmne, & Mmxe,Mmin,Mmax,Theta,P,R,T,U,W,Dldc,Perr,Oerr, & Cerr) IMPLICIT NONE !*--FDLDS2401 C*** Start of declarations inserted by SPAG DOUBLE PRECISION ar , Clat , clato , costhe , cotthe , cscth2 , & cscthe , Dldc , fa , fc , fd , fdrr , fm , & fmfpst , fn , fnfp , fnm1 , fnp1 , fnp2 , fp DOUBLE PRECISION fprr , fr , fs , H , P , pbppp , pbppr , pbppt , & pbrpp , pbrpr , pbrpt , pbtpp , pbtpr , pbtpt , & Phi , phio , R , ra , Re , Rm DOUBLE PRECISION Ro , roo , Rp , sinthe , T , Theta , thetao , W INTEGER ic , id , Ie , Ii , ip , lend , m , Mmax , Mmin , Mmne , & Mmni , Mmxe , Mmxi , n , Nc , Nci , Nmax , Nmne , Nmni , & Nmxe INTEGER Nmxi , Np , Nsi C*** End of declarations inserted by SPAG SAVE INTEGER*8 Perr , Oerr , Cerr integer Ctyp , pgen , Rgen , tgen , Grad INTEGER U(*) DIMENSION Dldc(*) , P(*) , R(*) , T(*) , W(*) DATA roo/0.D0/ , phio/0.D0/ , thetao/0.D0/ , clato/0.D0/ CALL GEOCEN(Ctyp,Re,Rp,Rm,H,Clat,Ro,Theta,sinthe,costhe) pgen = Rgen - 6 tgen = pgen IF ( phio.NE.Phi ) THEN pgen = 1 Rgen = Rgen + 1 phio = Phi ENDIF IF ( thetao.NE.Theta ) THEN tgen = 1 Rgen = Rgen + 1 thetao = Theta ENDIF IF ( clato.NE.Clat ) THEN Rgen = Rgen + 1 clato = Clat ENDIF IF ( roo.NE.Ro ) THEN Rgen = Rgen + 1 roo = Ro ENDIF IF ( Rgen.GT.0 ) THEN IF ( sinthe.EQ.0.D0 ) THEN IF ( Grad.EQ.0 ) THEN Cerr = 1 IF ( Perr.NE.0 ) WRITE (Oerr,99001) Cerr 99001 FORMAT (' ',/,1X,'SUBROUTINE BFIELD -- ERROR CODE ',I2, &' -- GEOGRAPHIC POLAR POSITION DETECTED, B-PHI INDETERMINABLE -- W &ARNING',/,' ') ELSE Cerr = 2 IF ( Perr.NE.0 ) WRITE (Oerr,99002) Cerr 99002 FORMAT (' ',/,1X,'SUBROUTINE BFIELD -- ERROR CODE ',I2, & ' -- GEOGRAPHIC POLAR POSITION DETECTED, B-PHI AND' & ,/,36X, & '-- PHI-DERIVATIVE GRADIENT COMPONENTS INDETERMINABLE -- WARNING' & ,/,' ') ENDIF Phi = 0.D0 cscthe = 0.D0 cscth2 = 0.D0 cotthe = 0.D0 ELSE cscthe = 1.D0/sinthe cscth2 = cscthe*cscthe cotthe = costhe*cscthe ENDIF IF ( tgen.GT.0 ) CALL SCHMIT(Grad,Rgen,Nmax,Mmin,Mmax,sinthe, & costhe,P,R) IF ( pgen.GT.0 ) CALL TRIGMP(Mmax,Phi,T) CALL R8VSET(1,(6+18*Grad)*Nc,0.D0,Dldc) ic = 0 id = 1 ip = Ii ar = Rm/Ro fa = ar**(Nmni+1) DO n = Nmni , Nmxi fnp1 = DBLE(n+1) fnp2 = DBLE(n+2) fa = fa*ar DO m = Mmin , MIN(n,Mmax) ip = ip + 1 IF ( (m.GE.Mmni) .AND. (m.LE.Mmxi) ) THEN ic = ic + 1 fm = DBLE(m) fp = fa*P(ip) fd = -fa*P(Np+ip) fc = T(m+1) fs = T(Mmax+m+2) fnfp = fnp1*fp fmfpst = fm*fp*cscthe lend = U(ic) CALL R8VSET(id,lend,fd*fc,Dldc) CALL R8VSET(Nc+id,lend,fmfpst*fs,Dldc) CALL R8VSET(2*Nc+id,lend,fnfp*fc,Dldc) IF ( Grad.EQ.1 ) THEN fprr = fp/Ro fdrr = fd/Ro pbtpt = -fa*P(2*Np+ip)/Ro pbtpp = -fm*fdrr*cscthe pbtpr = -fnp2*fdrr pbppt = -fm*(fdrr+fprr*cotthe)*cscthe pbppp = fm*fm*fprr*cscth2 pbppr = -fnp2*fm*fprr*cscthe pbrpt = -fnp1*fdrr pbrpp = -fnp1*fm*fprr*cscthe pbrpr = -fnp1*fnp2*fprr CALL R8VSET(6*Nc+id,lend,pbtpt*fc,Dldc) CALL R8VSET(7*Nc+id,lend,pbtpp*fs,Dldc) CALL R8VSET(8*Nc+id,lend,pbtpr*fc,Dldc) CALL R8VSET(9*Nc+id,lend,pbppt*fs,Dldc) CALL R8VSET(10*Nc+id,lend,pbppp*fc,Dldc) CALL R8VSET(11*Nc+id,lend,pbppr*fs,Dldc) CALL R8VSET(12*Nc+id,lend,pbrpt*fc,Dldc) CALL R8VSET(13*Nc+id,lend,pbrpp*fs,Dldc) CALL R8VSET(14*Nc+id,lend,pbrpr*fc,Dldc) ENDIF id = id + lend IF ( m.GT.0 ) THEN ic = ic + 1 lend = U(ic) CALL R8VSET(id,lend,fd*fs,Dldc) CALL R8VSET(Nc+id,lend,-fmfpst*fc,Dldc) CALL R8VSET(2*Nc+id,lend,fnfp*fs,Dldc) IF ( Grad.EQ.1 ) THEN CALL R8VSET(6*Nc+id,lend,pbtpt*fs,Dldc) CALL R8VSET(7*Nc+id,lend,-pbtpp*fc,Dldc) CALL R8VSET(8*Nc+id,lend,pbtpr*fs,Dldc) CALL R8VSET(9*Nc+id,lend,-pbppt*fc,Dldc) CALL R8VSET(10*Nc+id,lend,pbppp*fs,Dldc) CALL R8VSET(11*Nc+id,lend,-pbppr*fc,Dldc) CALL R8VSET(12*Nc+id,lend,pbrpt*fs,Dldc) CALL R8VSET(13*Nc+id,lend,-pbrpp*fc,Dldc) CALL R8VSET(14*Nc+id,lend,pbrpr*fs,Dldc) ENDIF id = id + lend ENDIF ENDIF ENDDO ENDDO ip = Ie ra = Ro/Rm fr = ra**(Nmne-2) DO n = Nmne , Nmxe fnm1 = DBLE(n-1) fn = DBLE(n) fr = fr*ra DO m = Mmin , MIN(n,Mmax) ip = ip + 1 IF ( (m.GE.Mmne) .AND. (m.LE.Mmxe) ) THEN ic = ic + 1 fm = DBLE(m) fp = fr*P(ip) fd = -fr*P(Np+ip) fc = T(m+1) fs = T(Mmax+m+2) fnfp = -fn*fp fmfpst = fm*fp*cscthe lend = U(ic) CALL R8VSET(id,lend,fd*fc,Dldc) CALL R8VSET(Nc+id,lend,fmfpst*fs,Dldc) CALL R8VSET(2*Nc+id,lend,fnfp*fc,Dldc) IF ( Grad.EQ.1 ) THEN fprr = fp/Ro fdrr = fd/Ro pbtpt = -fr*P(2*Np+ip)/Ro pbtpp = -fm*fdrr*cscthe pbtpr = fnm1*fdrr pbppt = -fm*(fdrr+fprr*cotthe)*cscthe pbppp = fm*fm*fprr*cscth2 pbppr = fnm1*fm*fprr*cscthe pbrpt = fn*fdrr pbrpp = fn*fm*fprr*cscthe pbrpr = -fn*fnm1*fprr CALL R8VSET(6*Nc+id,lend,pbtpt*fc,Dldc) CALL R8VSET(7*Nc+id,lend,pbtpp*fs,Dldc) CALL R8VSET(8*Nc+id,lend,pbtpr*fc,Dldc) CALL R8VSET(9*Nc+id,lend,pbppt*fs,Dldc) CALL R8VSET(10*Nc+id,lend,pbppp*fc,Dldc) CALL R8VSET(11*Nc+id,lend,pbppr*fs,Dldc) CALL R8VSET(12*Nc+id,lend,pbrpt*fc,Dldc) CALL R8VSET(13*Nc+id,lend,pbrpp*fs,Dldc) CALL R8VSET(14*Nc+id,lend,pbrpr*fc,Dldc) ENDIF id = id + lend IF ( m.GT.0 ) THEN ic = ic + 1 lend = U(ic) CALL R8VSET(id,lend,fd*fs,Dldc) CALL R8VSET(Nc+id,lend,-fmfpst*fc,Dldc) CALL R8VSET(2*Nc+id,lend,fnfp*fs,Dldc) IF ( Grad.EQ.1 ) THEN CALL R8VSET(6*Nc+id,lend,pbtpt*fs,Dldc) CALL R8VSET(7*Nc+id,lend,-pbtpp*fc,Dldc) CALL R8VSET(8*Nc+id,lend,pbtpr*fs,Dldc) CALL R8VSET(9*Nc+id,lend,-pbppt*fc,Dldc) CALL R8VSET(10*Nc+id,lend,pbppp*fs,Dldc) CALL R8VSET(11*Nc+id,lend,-pbppr*fc,Dldc) CALL R8VSET(12*Nc+id,lend,pbrpt*fs,Dldc) CALL R8VSET(13*Nc+id,lend,-pbrpp*fc,Dldc) CALL R8VSET(14*Nc+id,lend,pbrpr*fs,Dldc) ENDIF id = id + lend ENDIF ENDIF ENDDO ENDDO CALL TDC(Grad,Nc,Clat,Theta,Dldc,W) ENDIF RETURN END !*==geocen.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE GEOCEN(Ctyp,Re,Rp,Rm,H,Clat,R,Theta,Sinthe,Costhe) IMPLICIT NONE !*--GEOCEN2614 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Clat , costh2 , Costhe , H , R , Re , re2 , & rhoew , rhons , Rm , Rp , rp2 , sinth2 , Sinthe , & Theta INTEGER ifirst C*** End of declarations inserted by SPAG SAVE INTEGER Ctyp DOUBLE PRECISION kappa DATA ifirst/1/ Theta = Clat R = Rm + H Sinthe = DSIN(Theta) Costhe = DCOS(Theta) IF ( Ctyp.EQ.0 ) THEN IF ( ifirst.EQ.1 ) THEN ifirst = 0 re2 = Re*Re rp2 = Rp*Rp ENDIF sinth2 = Sinthe*Sinthe costh2 = Costhe*Costhe kappa = DSQRT(re2*sinth2+rp2*costh2) rhoew = (re2/kappa) + H rhons = (rp2/kappa) + H Theta = DATAN2(rhoew*Sinthe,rhons*Costhe) R = DSQRT(rhoew**2*sinth2+rhons**2*costh2) Sinthe = DSIN(Theta) Costhe = DCOS(Theta) ENDIF END !*==schmit.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE SCHMIT(Grad,Rgen,Nmax,Mmin,Mmax,Sinthe,Costhe,P,R) IMPLICIT NONE !*--SCHMIT2649 C*** Start of declarations inserted by SPAG DOUBLE PRECISION costh2 , Costhe , cotthe , cscth2 , cscthe , & cthsth , P , R , root3 , sinth2 , Sinthe INTEGER kd0 , kd1 , kd2 , knnp , kp0 , kp1 , kp2 , kq0 , kq1 , & kq2 , ksec , ksm2 , ktes , ktm2 , l , Mmax , Mmin , n , & nd0sec , nd0tes INTEGER nd1tes , Nmax , np , np1sec , np1tes , np2tes , nq0msq , & nq0nnp C*** End of declarations inserted by SPAG SAVE INTEGER Rgen , Grad DIMENSION P(*) , R(*) DATA root3/1.732050807568877D0/ IF ( Rgen.GT.6 ) CALL SRECUR(Grad,Nmax,Mmin,Mmax,ksm2,ktm2,np, & np1sec,nd0sec,np1tes,np2tes,nd0tes, & nd1tes,nq0nnp,nq0msq,R) IF ( Nmax.GE.1 ) THEN kp0 = 1 kp2 = kp0 IF ( Mmin.LE.0 ) THEN P(kp0) = Costhe P(np+kp0) = -Sinthe IF ( Grad.EQ.1 ) P(2*np+kp0) = -Costhe kp0 = kp0 + 1 ENDIF IF ( Mmax.GE.1 ) THEN P(kp0) = Sinthe P(np+kp0) = Costhe IF ( Grad.EQ.1 ) P(2*np+kp0) = -Sinthe kp0 = kp0 + 1 ENDIF IF ( Nmax.GE.2 ) THEN kp1 = kp0 sinth2 = Sinthe*Sinthe costh2 = Costhe*Costhe cthsth = Costhe*Sinthe IF ( Mmin.LE.0 ) THEN P(kp0) = (3.D0*costh2-1.D0)/2.D0 P(np+kp0) = -3.D0*cthsth IF ( Grad.EQ.1 ) P(2*np+kp0) = -3.D0*(2.D0*costh2-1.D0) kp0 = kp0 + 1 ENDIF IF ( (Mmin.LE.1) .AND. (Mmax.GE.1) ) THEN P(kp0) = root3*cthsth P(np+kp0) = root3*(2.D0*costh2-1.D0) IF ( Grad.EQ.1 ) P(2*np+kp0) = -4.D0*root3*cthsth kp0 = kp0 + 1 ENDIF IF ( Mmax.GE.2 ) THEN P(kp0) = root3*sinth2/2.D0 P(np+kp0) = root3*cthsth IF ( Grad.EQ.1 ) P(2*np+kp0) = root3*(2.D0*costh2-1.D0) kp0 = kp0 + 1 ENDIF kd2 = np + kp2 kd1 = np + kp1 kd0 = np + kp0 kq2 = np + kd2 kq1 = np + kd1 kq0 = np + kd0 ksec = 1 ktes = 1 knnp = 1 IF ( Sinthe.EQ.0.D0 ) THEN DO n = 3 , Nmax l = MAX(0,MIN(n-1,Mmax)-Mmin+1) IF ( l.GT.0 ) THEN CALL R8VSET(kp0,l,0.D0,P) CALL R8VLINKQ(kp1,np1tes+ktes,kp0,l,Costhe,P,R,P) CALL R8VLINKQ(kp2,np2tes+ktes,kp0,l,-1.D0,P,R,P) CALL R8VSET(kd0,l,0.D0,P) CALL R8VLINKQ(kd1,np1tes+ktes,kd0,l,Costhe,P,R,P) CALL R8VLINKQ(kd2,np2tes+ktes,kd0,l,-1.D0,P,R,P) IF ( Grad.EQ.1 ) THEN CALL R8VSET(kq0,l,0.D0,P) CALL R8VLINKQ(kq1,np1tes+ktes,kq0,l,Costhe,P,R, & P) CALL R8VLINKQ(kp1,np1tes+ktes,kq0,l,-Costhe,P,R, & P) CALL R8VLINKQ(kq2,np2tes+ktes,kq0,l,-1.D0,P,R,P) ENDIF ktes = ktes + l ENDIF IF ( n.LE.Mmax ) THEN P(kp0+l) = 0.D0 P(kd0+l) = 0.D0 IF ( Grad.EQ.1 ) P(kq0+l) = 0.D0 l = l + 1 ENDIF kp2 = kp1 kp1 = kp0 kp0 = kp0 + l kd2 = kd1 kd1 = kd0 kd0 = kd0 + l kq2 = kq1 kq1 = kq0 kq0 = kq0 + l ENDDO ELSE cotthe = Costhe/Sinthe cscthe = 1.D0/Sinthe cscth2 = cscthe*cscthe DO n = 3 , Nmax l = MAX(0,MIN(n-1,Mmax)-Mmin+1) IF ( l.GT.0 ) THEN CALL R8VSET(kp0,l,0.D0,P) CALL R8VLINKQ(kp1,np1tes+ktes,kp0,l,Costhe,P,R,P) CALL R8VLINKQ(kp2,np2tes+ktes,kp0,l,-1.D0,P,R,P) CALL R8VSET(kd0,l,0.D0,P) CALL R8VLINKQ(kp0,nd0tes+ktes,kd0,l,cotthe,P,R,P) CALL R8VLINKQ(kp1,nd1tes+ktes,kd0,l,-cscthe,P,R,P) IF ( Grad.EQ.1 ) THEN CALL R8VSET(kq0,l,0.D0,P) CALL R8VLINKQ(kp0,nq0msq+ktm2,kq0,l,cscth2,P,R, & P) CALL R8VLINKT(kp0,kq0,l,-R(nq0nnp+knnp),P,P) CALL R8VLINKT(kd0,kq0,l,-cotthe,P,P) ENDIF ktes = ktes + l ENDIF IF ( n.LE.Mmax ) THEN P(kp0+l) = R(np1sec+ksec)*Sinthe*P(kp0-1) P(kd0+l) = R(nd0sec+ksec)*cotthe*P(kp0+l) IF ( Grad.EQ.1 ) P(kq0+l) & = (R(nq0msq-ksm2+n+1)*cscth2-R(nq0nnp+knnp)) & *P(kp0+l) - cotthe*P(kd0+l) ksec = ksec + 1 l = l + 1 ENDIF knnp = knnp + 1 kp2 = kp1 kp1 = kp0 kp0 = kp0 + l kd0 = kd0 + l kq0 = kq0 + l ENDDO ENDIF ENDIF ENDIF END !*==srecur.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE SRECUR(Grad,Nmax,Mmin,Mmax,Ksm2,Ktm2,Npall,Nad1,Nad2, & Nad3,Nad4,Nad5,Nad6,Nad7,Nad8,R) IMPLICIT NONE !*--SRECUR2795 C*** Start of declarations inserted by SPAG DOUBLE PRECISION f2n , f2nl1 , fdsqrt , fm , fmsq , fn , fnl1 , & fnl1sq , fnp1 , fnsq , R INTEGER kmax , kmin , kmsq , knnp , ksec , Ksm2 , ktes , Ktm2 , & lmax , m , Mmax , Mmin , n , Nad1 , Nad2 , Nad3 , Nad4 , & Nad5 , Nad6 , Nad7 INTEGER Nad8 , NLPX , Nmax , nnnp1 , Npall , nsect , ntess C*** End of declarations inserted by SPAG SAVE INTEGER Grad DIMENSION R(*) lmax = MIN(Nmax,2) kmax = MIN(Mmax,2) kmin = MIN(Mmin,2) Ksm2 = MIN(Mmin,3) Ktm2 = DIM(Mmin,3) + 1 Npall = NLPX(Nmax,Mmax,Mmin) ntess = Npall - NLPX(lmax,kmax,kmin) + kmax - Mmax nsect = MAX(0,Mmax-2) nnnp1 = MAX(0,Nmax-2) Nad1 = 0 Nad2 = nsect Nad3 = nsect + Nad2 Nad4 = ntess + Nad3 Nad5 = ntess + Nad4 Nad6 = ntess + Nad5 Nad7 = ntess + Nad6 Nad8 = nnnp1 + Nad7 ksec = 0 ktes = 0 knnp = 0 DO n = 3 , Nmax fnp1 = DBLE(n+1) fn = DBLE(n) fnl1 = DBLE(n-1) f2n = 2.D0*fn f2nl1 = f2n - 1.D0 fnsq = fn*fn fnl1sq = fnl1*fnl1 IF ( n.LE.Mmax ) THEN ksec = ksec + 1 R(Nad1+ksec) = DSQRT(f2nl1/f2n) R(Nad2+ksec) = fn ENDIF IF ( Grad.EQ.1 ) THEN knnp = knnp + 1 R(Nad7+knnp) = fn*fnp1 ENDIF DO m = Mmin , MIN(n-1,Mmax) ktes = ktes + 1 fm = DBLE(m) fmsq = fm*fm fdsqrt = DSQRT(fnsq-fmsq) R(Nad3+ktes) = f2nl1/fdsqrt R(Nad4+ktes) = DSQRT(fnl1sq-fmsq)/fdsqrt R(Nad5+ktes) = fn R(Nad6+ktes) = fdsqrt ENDDO ENDDO IF ( Grad.EQ.1 ) THEN kmsq = 0 DO m = Ksm2 , Mmax kmsq = kmsq + 1 fm = DBLE(m) R(Nad8+kmsq) = fm*fm ENDDO ENDIF END !*==trigmp.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE TRIGMP(Mmax,Phi,T) IMPLICIT NONE !*--TRIGMP2867 C*** Start of declarations inserted by SPAG INTEGER m , Mmax DOUBLE PRECISION Phi , T C*** End of declarations inserted by SPAG SAVE DIMENSION T(*) T(1) = 1.D0 T(Mmax+2) = 0.D0 IF ( Mmax.GT.0 ) THEN T(2) = DCOS(Phi) T(Mmax+3) = DSIN(Phi) DO m = 2 , Mmax T(m+1) = 2.D0*T(2)*T(m) - T(m-1) T(Mmax+m+2) = 2.D0*T(2)*T(Mmax+m+1) - T(Mmax+m) ENDDO ENDIF END !*==tdc.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE TDC(Grad,Nc,Clat,Theta,Dldc,R) IMPLICIT NONE !*--TDC2888 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Clat , cpsi , Dldc , psi , R , spsi , Theta INTEGER ione , Nc , null C*** End of declarations inserted by SPAG SAVE INTEGER Grad DIMENSION R(*) , Dldc(*) DATA null/0/ , ione/1/ psi = Theta - Clat spsi = DSIN(psi) cpsi = DCOS(psi) R(1) = -cpsi R(2) = 0.D0 R(3) = -spsi R(4) = 0.D0 R(5) = 1.D0 R(6) = 0.D0 R(7) = spsi R(8) = 0.D0 R(9) = -cpsi CALL LTRANV(ione,Nc,Nc,R,Dldc) IF ( Grad.EQ.1 ) THEN CALL LTRANV(null,3*Nc,3*Nc,R,Dldc(6*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(6*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(9*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(12*Nc+1)) ENDIF END !*==fdsdc.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE FDSDC(Rgen,Ityp,Etyp,Nsi,Nse,Nc,Nci,Ta,Tb,Dst,Dstt, & Dsti,Bori,Bkni,Bkpi,Tdgi,Dste,Bore,Bkne,Bkpe, & Tdge,U,W,Dsdc,Perr,Oerr,Cerr) IMPLICIT NONE !*--FDSDC2922 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Bkpe , Bkpi , Dsdc , Dst , Dstt , Ta , Tb , tbo , & W INTEGER idst , isvr , Ityp , Nc , Nci , Nse , Nsi C*** End of declarations inserted by SPAG SAVE INTEGER*8 Perr , Oerr , Cerr Integer Etyp , esvr , edst , Rgen , tgen INTEGER Bore(*) , Bori(*) , Bkne(*) , Bkni(*) , Dste(*) , Dsti(*) & , Tdge(*) INTEGER Tdgi(*) , U(*) DIMENSION Bkpi(*) , Bkpe(*) , Dsdc(*) , W(*) DATA tbo/0.D0/ tgen = MAX(0,Rgen-6) IF ( tbo.NE.Tb ) THEN tgen = MIN(1,tgen+Ityp+Etyp) Rgen = Rgen + tgen tbo = Tb ENDIF IF ( tgen.GT.0 ) THEN CALL R8VSET(1,2*Nc,0.D0,Dsdc) IF ( Nsi.GT.0 ) THEN isvr = MOD(Ityp,3) idst = Ityp/3 CALL I8VCUM(1,1,Nsi,U) CALL R8VSCATS(1,Nsi,1.D0,U,Dsdc(1)) CALL R8VSCATS(1,Nsi,0.D0,U,Dsdc(Nc+1)) CALL I8VADDS(1,1,Nsi,1,U,U) IF ( isvr.EQ.1 ) THEN CALL TAYLOR(Nc,Nsi,Ta,Tb,Tdgi,U(1),W,Dsdc(1)) ELSEIF ( isvr.EQ.2 ) THEN CALL BSPLYN(Nc,Nsi,Ta,Tb,Bori,Bkni,Bkpi,U(1),W,Dsdc(1), & Perr,Oerr,Cerr) IF ( Cerr.GE.50 ) GOTO 20 ENDIF IF ( idst.EQ.1 ) CALL DSTORM(Nc,Nsi,Dst,Dstt,Dsti,U(1), & Dsdc(1)) 20 CALL I8VDEL(1,1,Nsi,U) ENDIF IF ( Nse.GT.0 ) THEN esvr = MOD(Etyp,3) edst = Etyp/3 CALL I8VCUM(1,Nsi+1,Nse,U) CALL R8VSCATS(Nsi+1,Nse,1.D0,U,Dsdc(Nci+1)) CALL R8VSCATS(Nsi+1,Nse,0.D0,U,Dsdc(Nc+Nci+1)) CALL I8VADDS(Nsi+1,Nsi+1,Nse,1,U,U) IF ( esvr.EQ.1 ) THEN CALL TAYLOR(Nc,Nse,Ta,Tb,Tdge,U(Nsi+1),W,Dsdc(Nci+1)) ELSEIF ( esvr.EQ.2 ) THEN CALL BSPLYN(Nc,Nse,Ta,Tb,Bore,Bkne,Bkpe,U(Nsi+1),W, & Dsdc(Nci+1),Perr,Oerr,Cerr) IF ( Cerr.GE.50 ) GOTO 40 ENDIF IF ( edst.EQ.1 ) CALL DSTORM(Nc,Nse,Dst,Dstt,Dste,U(Nsi+1), & Dsdc(Nci+1)) 40 CALL I8VDEL(1,Nsi+1,Nse,U) ENDIF ENDIF END !*==taylor.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE TAYLOR(Nc,Ns,Ta,Tb,Tdeg,U,Dsdt,Dsdc) IMPLICIT NONE !*--TAYLOR2984 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Dsdc , Dsdt , dt , Ta , Tb INTEGER i , iu , j , n , Nc , Ns C*** End of declarations inserted by SPAG SAVE INTEGER Tdeg(*) , U(*) DIMENSION Dsdc(*) , Dsdt(*) dt = Tb - Ta DO i = 1 , Ns n = Tdeg(i) IF ( n.GT.0 ) THEN iu = U(i) Dsdt(1) = 1.D0 DO j = 1 , n Dsdt(j+1) = Dsdt(j)*dt/DBLE(j) ENDDO CALL R8VGATHP(2,1,iu,n,Dsdt,Dsdc) CALL R8VGATHP(1,1,Nc+iu,n,Dsdt,Dsdc) U(i) = U(i) + n ENDIF ENDDO END !*==bsplyn.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE BSPLYN(Nc,Ns,Ta,Tb,Bord,Bkno,Bkpo,U,Dtdb,Dsdc,Perr, & Oerr,Cerr) IMPLICIT NONE !*--BSPLYN3011 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Bkpo , Dsdc , Dtdb , Ta , Tb INTEGER i , ik , iu , k , l , m , n , Nc , Ns C*** End of declarations inserted by SPAG SAVE INTEGER*8 Perr , Oerr , Cerr Integer Bord(*) , Bkno(*) , U(*) DIMENSION Bkpo(*) , Dsdc(*) , Dtdb(*) ik = 1 DO i = 1 , Ns n = Bord(i) k = Bkno(i) l = n + k iu = U(i) IF ( n.GT.0 ) THEN m = n + 1 CALL SBSPLN(Ta,Tb,m,k,Bkpo(ik),Dtdb,Dsdc(iu),Perr,Oerr,Cerr) IF ( Cerr.GE.50 ) GOTO 99999 CALL R8VSET(1,l+1,0.D0,Dtdb) CALL TBSPLN(Tb,n,k,Bkpo(ik),Dtdb,Perr,Oerr,Cerr) IF ( Cerr.GE.50 ) GOTO 99999 CALL R8VGATHP(1,1,Nc+iu,l,Dtdb,Dsdc) ELSEIF ( k.GT.0 ) THEN CALL R8VSET(iu,l,0.D0,Dsdc) CALL R8VSET(Nc+iu,l,0.D0,Dsdc) ENDIF U(i) = U(i) + l ik = ik + k + 2 ENDDO 99999 END !*==sbspln.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE SBSPLN(Ta,Tb,N,K,Bkpo,Dtdb,Dsdc,Perr,Oerr,Cerr) IMPLICIT NONE !*--SBSPLN3044 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Bkpo , Dsdc , Dtdb , R8SSUM , rn , Ta , Tb INTEGER i , ik , ip , jk , K , N , nl , np , ns C*** End of declarations inserted by SPAG SAVE INTEGER*8 Perr , Oerr , Cerr DIMENSION Bkpo(*) , Dsdc(*) , Dtdb(*) IF ( N.GT.1 ) THEN nl = N + K + 1 CALL R8VSET(1,2*nl,0.D0,Dtdb) CALL TBSPLN(Tb,N,K,Bkpo,Dtdb(1),Perr,Oerr,Cerr) IF ( Cerr.LT.50 ) THEN CALL TBSPLN(Ta,N,K,Bkpo,Dtdb(nl+1),Perr,Oerr,Cerr) IF ( Cerr.LT.50 ) THEN CALL R8VSUB(nl+1,1,1,nl,Dtdb,Dtdb,Dtdb) rn = 1.D0/DBLE(N-1) np = N + K - 1 ns = np DO i = 1 , np ip = i + 1 ik = MIN(ip,K+2) jk = MAX(ip-N+1,1) Dsdc(i) = (Bkpo(ik)-Bkpo(jk))*R8SSUM(ip,ns,Dtdb) ns = ns - 1 ENDDO IF ( np.GT.0 ) CALL R8VSCALE(1,np,rn,Dsdc) ENDIF ENDIF ENDIF END !*==tbspln.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE TBSPLN(T,N,K,Bkpo,Dtdb,Perr,Oerr,Cerr) IMPLICIT NONE !*--TBSPLN3078 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Bkpo , delk , difi , difj , Dtdb , T , temp INTEGER i , ik , j , jk , K , l , m , N C*** End of declarations inserted by SPAG SAVE INTEGER*8 Perr , Oerr , Cerr Integer p , y , posb DIMENSION Bkpo(*) , Dtdb(*) IF ( (T.GE.Bkpo(1)) .AND. (T.LE.Bkpo(K+2)) ) THEN CALL R8SLT(1,K+2,T,Bkpo,y) l = MIN(K+2,y+1) p = l + N - 2 IF ( N.EQ.1 ) THEN Dtdb(p) = 1.D0 ELSEIF ( N.GT.1 ) THEN m = l ik = MIN(m,K+2) jk = MAX(m-1,1) difi = Bkpo(ik) - T delk = Bkpo(ik) - Bkpo(jk) IF ( delk.EQ.0.D0 ) THEN Dtdb(p) = 0.D0 ELSE Dtdb(p) = 1.D0/delk ENDIF posb = p - 1 DO j = 2 , N jk = MAX(m-j,1) temp = difi*Dtdb(posb+1) delk = Bkpo(ik) - Bkpo(jk) IF ( delk.EQ.0.D0 ) THEN temp = 0.D0 ELSEIF ( j.LT.N ) THEN temp = temp/delk ENDIF Dtdb(posb) = temp posb = posb - 1 ENDDO Dtdb(p+1) = 0.D0 m = m + 1 DO i = 2 , N ik = MIN(m,K+2) difi = Bkpo(ik) - T posb = p DO j = i , N jk = MAX(m-j,1) difj = T - Bkpo(jk) temp = difj*Dtdb(posb) + difi*Dtdb(posb+1) delk = Bkpo(ik) - Bkpo(jk) IF ( delk.EQ.0.D0 ) THEN temp = 0.D0 ELSEIF ( j.LT.N ) THEN temp = temp/delk ENDIF Dtdb(posb) = temp posb = posb - 1 ENDDO m = m + 1 ENDDO ENDIF ELSE Cerr = 54 IF ( Perr.NE.0 ) WRITE (Oerr,99001) Cerr 99001 FORMAT (' ',/,1X,'SUBROUTINE TBSPLN -- ERROR CODE ',I2, & ' -- T LIES OUTSIDE OF KNOT DOMAIN -- ABORT',/,' ') ENDIF RETURN END !*==dstorm.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE DSTORM(Nc,Ns,Dst,Dstt,Dstm,U,Dsdc) IMPLICIT NONE !*--DSTORM3149 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Dsdc , Dst , Dstt INTEGER i , iu , n , Nc , Ns C*** End of declarations inserted by SPAG SAVE INTEGER Dstm(*) , U(*) DIMENSION Dsdc(*) DO i = 1 , Ns n = Dstm(i) IF ( n.GT.0 ) THEN iu = U(i) CALL R8VSET(iu,n,Dst,Dsdc) CALL R8VSET(Nc+iu,n,Dstt,Dsdc) U(i) = U(i) + n ENDIF ENDDO END !*==fdldc.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE FDLDC(Grad,Nc,Dsdc,Dldc) IMPLICIT NONE !*--FDLDC3170 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Dldc , Dsdc INTEGER i , j , Nc C*** End of declarations inserted by SPAG SAVE INTEGER Grad DIMENSION Dsdc(*) , Dldc(*) i = 1 DO j = 1 , 3 CALL R8VMUL(Nc+1,i,3*Nc+i,Nc,Dsdc,Dldc,Dldc) i = i + Nc ENDDO i = 1 DO j = 1 , 3 CALL R8VMUL(1,i,i,Nc,Dsdc,Dldc,Dldc) i = i + Nc ENDDO IF ( Grad.EQ.1 ) THEN i = 1 DO j = 1 , 9 CALL R8VMUL(Nc+1,6*Nc+i,15*Nc+i,Nc,Dsdc,Dldc,Dldc) i = i + Nc ENDDO i = 1 DO j = 1 , 9 CALL R8VMUL(1,6*Nc+i,6*Nc+i,Nc,Dsdc,Dldc,Dldc) i = i + Nc ENDDO ENDIF END !*==blgen.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE BLGEN(Grad,Nc,B,C,Dldc) IMPLICIT NONE !*--BLGEN3204 C*** Start of declarations inserted by SPAG DOUBLE PRECISION B , C , Dldc , R8SDOT INTEGER i , j , Nc C*** End of declarations inserted by SPAG SAVE INTEGER Grad DIMENSION B(*) , C(*) , Dldc(*) i = 1 DO j = 1 , 6 B(j) = B(j) + R8SDOT(i,1,Nc,Dldc,C) i = i + Nc ENDDO IF ( Grad.EQ.1 ) THEN i = 1 DO j = 29 , 46 B(j) = B(j) + R8SDOT(6*Nc+i,1,Nc,Dldc,C) i = i + Nc ENDDO ENDIF END !*==bngen.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE BNGEN(B) IMPLICIT NONE !*--BNGEN3228 C*** Start of declarations inserted by SPAG DOUBLE PRECISION B , cd , cdt , cf , cft , ch , cht , ci , cit , & cx , cxt , cy , cyt , cz , czt C*** End of declarations inserted by SPAG SAVE DIMENSION B(*) cx = B(1) cy = B(2) cz = B(3) cxt = B(4) cyt = B(5) czt = B(6) ch = DSQRT(cx*cx+cy*cy) cf = DSQRT(ch*ch+cz*cz) IF ( ch.NE.0.D0 ) THEN cd = 2.D0*DATAN(cy/(ch+cx)) cdt = (cx*cyt-cy*cxt)/ch/ch cht = (cx*cxt+cy*cyt)/ch ELSE cd = 0.D0 cdt = 0.D0 cht = 0.D0 ENDIF IF ( cf.NE.0.D0 ) THEN ci = 2.D0*DATAN(cz/(cf+ch)) cit = (ch*czt-cz*cht)/cf/cf cft = (ch*cht+cz*czt)/cf ELSE ci = 0.D0 cit = 0.D0 cft = 0.D0 ENDIF B(7) = cd B(8) = ci B(9) = ch B(10) = cf B(11) = cdt B(12) = cit B(13) = cht B(14) = cft END !*==tec.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE TEC(Grad,K,Nc,Theta,Phi,B,Dldc,R) IMPLICIT NONE !*--TEC3273 C*** Start of declarations inserted by SPAG DOUBLE PRECISION B , costhe , cphi , Dldc , Phi , R , sinthe , & sphi , Theta INTEGER ione , K , Nc , null C*** End of declarations inserted by SPAG SAVE INTEGER Grad DIMENSION Dldc(*) , B(*) , R(*) DATA null/0/ , ione/1/ IF ( K.GE.1 ) THEN sinthe = DSIN(Theta) costhe = DCOS(Theta) sphi = DSIN(Phi) cphi = DCOS(Phi) R(1) = -costhe*cphi R(2) = -sphi R(3) = -sinthe*cphi R(4) = -costhe*sphi R(5) = cphi R(6) = -sinthe*sphi R(7) = sinthe R(8) = 0.D0 R(9) = -costhe CALL LTRANS(ione,ione,B(1),R,B(1)) CALL LTRANS(ione,ione,B(4),R,B(4)) CALL LTRANV(ione,Nc,Nc,R,Dldc(1)) CALL LTRANV(null,Nc,Nc,R,Dldc(3*Nc+1)) CALL BNGEN(B) IF ( Grad.EQ.1 ) THEN CALL LTRANV(null,3,3,R,B(29)) CALL LTRANV(null,3,3,R,B(38)) CALL LTRANS(ione,ione,B(29),R,B(29)) CALL LTRANS(ione,ione,B(32),R,B(32)) CALL LTRANS(ione,ione,B(35),R,B(35)) CALL LTRANS(ione,ione,B(38),R,B(38)) CALL LTRANS(ione,ione,B(41),R,B(41)) CALL LTRANS(ione,ione,B(44),R,B(44)) CALL LTRANV(null,3*Nc,3*Nc,R,Dldc(6*Nc+1)) CALL LTRANV(null,3*Nc,3*Nc,R,Dldc(15*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(6*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(9*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(12*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(15*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(18*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(21*Nc+1)) ENDIF ENDIF END !*==tse.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE TSE(Grad,K,Nc,Rse,B,Dldc,R) IMPLICIT NONE !*--TSE3325 C*** Start of declarations inserted by SPAG DOUBLE PRECISION B , Dldc , R , Rse INTEGER ione , K , Nc , null C*** End of declarations inserted by SPAG SAVE INTEGER Grad DIMENSION Dldc(*) , Rse(*) , B(*) , R(*) DATA null/0/ , ione/1/ IF ( K.GE.1 ) THEN CALL R8VGATHP(1,1,1,9,Rse,R) CALL LTRANS(ione,ione,B(1),R,B(1)) CALL LTRANS(ione,ione,B(4),R,B(4)) CALL LTRANV(ione,Nc,Nc,R,Dldc(1)) CALL LTRANV(null,Nc,Nc,R,Dldc(3*Nc+1)) CALL BNGEN(B) IF ( Grad.EQ.1 ) THEN CALL LTRANV(null,3,3,R,B(29)) CALL LTRANV(null,3,3,R,B(38)) CALL LTRANS(ione,ione,B(29),R,B(29)) CALL LTRANS(ione,ione,B(32),R,B(32)) CALL LTRANS(ione,ione,B(35),R,B(35)) CALL LTRANS(ione,ione,B(38),R,B(38)) CALL LTRANS(ione,ione,B(41),R,B(41)) CALL LTRANS(ione,ione,B(44),R,B(44)) CALL LTRANV(null,3*Nc,3*Nc,R,Dldc(6*Nc+1)) CALL LTRANV(null,3*Nc,3*Nc,R,Dldc(15*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(6*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(9*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(12*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(15*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(18*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(21*Nc+1)) ENDIF ENDIF END !*==tms.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE TMS(Grad,K,Nc,Na,Ia,A,B,Dldc,Dlda,R) IMPLICIT NONE !*--TMS3364 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , B , ceulx , ceuly , ceulz , Dlda , Dldc , & eulx , euly , eulz , R , seulx , seuly , seulz INTEGER Ia , ione , K , Na , Nc , null C*** End of declarations inserted by SPAG SAVE INTEGER Grad DIMENSION Dlda(*) , Dldc(*) , A(*) , B(*) , R(*) DATA null/0/ , ione/1/ IF ( K.GE.1 ) THEN eulx = A(Ia+1) euly = A(Ia+2) eulz = A(Ia+3) seulx = DSIN(eulx) ceulx = DCOS(eulx) seuly = DSIN(euly) ceuly = DCOS(euly) seulz = DSIN(eulz) ceulz = DCOS(eulz) CALL FDLDEU(K,Na,Ia,seulx,ceulx,seuly,ceuly,seulz,ceulz,R,B, & Dlda) R(1) = ceuly*ceulz R(2) = -seulz R(3) = -seuly*ceulz R(4) = ceuly*ceulx*seulz + seuly*seulx R(5) = ceulx*ceulz R(6) = ceuly*seulx - seuly*seulz*ceulx R(7) = -ceuly*seulx*seulz + seuly*ceulx R(8) = -seulx*ceulz R(9) = ceuly*ceulx + seuly*seulx*seulz CALL LTRANS(ione,ione,B(1),R,B(1)) CALL LTRANS(ione,ione,B(4),R,B(4)) CALL LTRANV(ione,Nc,Nc,R,Dldc(1)) CALL LTRANV(null,Nc,Nc,R,Dldc(3*Nc+1)) CALL LTRANV(null,Na,Ia,R,Dlda(1)) CALL LTRANV(null,Na,Ia,R,Dlda(3*Na+1)) CALL BNGEN(B) IF ( Grad.EQ.1 ) THEN CALL LTRANV(null,3,3,R,B(29)) CALL LTRANV(null,3,3,R,B(38)) CALL LTRANS(ione,ione,B(29),R,B(29)) CALL LTRANS(ione,ione,B(32),R,B(32)) CALL LTRANS(ione,ione,B(35),R,B(35)) CALL LTRANS(ione,ione,B(38),R,B(38)) CALL LTRANS(ione,ione,B(41),R,B(41)) CALL LTRANS(ione,ione,B(44),R,B(44)) CALL LTRANV(null,3*Nc,3*Nc,R,Dldc(6*Nc+1)) CALL LTRANV(null,3*Nc,3*Nc,R,Dldc(15*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(6*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(9*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(12*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(15*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(18*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(21*Nc+1)) ENDIF Ia = Ia + 3 ENDIF END !*==fdldeu.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE FDLDEU(K,Na,Ia,Seulx,Ceulx,Seuly,Ceuly,Seulz,Ceulz,R,B, & Dlda) IMPLICIT NONE !*--FDLDEU3427 C*** Start of declarations inserted by SPAG DOUBLE PRECISION B , Ceulx , Ceuly , Ceulz , Dlda , R , Seulx , & Seuly , Seulz INTEGER i , Ia , ione , j , K , Na C*** End of declarations inserted by SPAG SAVE DIMENSION Dlda(*) , B(*) , R(*) DATA ione/1/ IF ( K.EQ.1 ) THEN i = Ia DO j = 1 , 6 Dlda(i+1) = 0.D0 Dlda(i+2) = 0.D0 Dlda(i+3) = 0.D0 i = i + Na ENDDO ELSE R(1) = 0.D0 R(2) = 0.D0 R(3) = 0.D0 R(4) = -Ceuly*Seulx*Seulz + Seuly*Ceulx R(5) = -Seulx*Ceulz R(6) = Ceuly*Ceulx + Seuly*Seulz*Seulx R(7) = -Ceuly*Ceulx*Seulz - Seuly*Seulx R(8) = -Ceulx*Ceulz R(9) = -Ceuly*Seulx + Seuly*Ceulx*Seulz CALL LTRANS(Na,ione,B(1),R,Dlda(Ia+1)) CALL LTRANS(Na,ione,B(4),R,Dlda(3*Na+Ia+1)) R(1) = -Seuly*Ceulz R(2) = 0.D0 R(3) = -Ceuly*Ceulz R(4) = -Seuly*Ceulx*Seulz + Ceuly*Seulx R(5) = 0.D0 R(6) = -Seuly*Seulx - Ceuly*Seulz*Ceulx R(7) = Seuly*Seulx*Seulz + Ceuly*Ceulx R(8) = 0.D0 R(9) = -Seuly*Ceulx + Ceuly*Seulx*Seulz CALL LTRANS(Na,ione,B(1),R,Dlda(Ia+2)) CALL LTRANS(Na,ione,B(4),R,Dlda(3*Na+Ia+2)) R(1) = -Ceuly*Seulz R(2) = -Ceulz R(3) = Seuly*Seulz R(4) = Ceuly*Ceulx*Ceulz R(5) = -Ceulx*Seulz R(6) = -Seuly*Ceulz*Ceulx R(7) = -Ceuly*Seulx*Ceulz R(8) = Seulx*Seulz R(9) = Seuly*Seulx*Ceulz CALL LTRANS(Na,ione,B(1),R,Dlda(Ia+3)) CALL LTRANS(Na,ione,B(4),R,Dlda(3*Na+Ia+3)) ENDIF END !*==tnm.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE TNM(Grad,K,Nc,Na,Ia,A,B,Dldc,Dlda,R) IMPLICIT NONE !*--TNM3483 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , B , cchix , cchiy , cchiz , chix , chiy , & chiz , Dlda , Dldc , R , schix , schiy , schiz INTEGER Ia , ione , K , Na , Nc , null C*** End of declarations inserted by SPAG SAVE INTEGER Grad DIMENSION Dlda(*) , Dldc(*) , A(*) , B(*) , R(*) DATA null/0/ , ione/1/ IF ( K.GE.1 ) THEN chix = A(Ia+1) chiy = A(Ia+2) chiz = A(Ia+3) schix = DSIN(chix) cchix = DCOS(chix) schiy = DSIN(chiy) cchiy = DCOS(chiy) schiz = DSIN(chiz) cchiz = DCOS(chiz) CALL FDLDNO(K,Na,Ia,schix,cchix,schiy,cchiy,schiz,cchiz,R,B, & Dlda) R(1) = 1.D0 R(2) = 0.D0 R(3) = 0.D0 R(4) = schix R(5) = cchix R(6) = 0.D0 R(7) = schiy*cchiz R(8) = schiy*schiz R(9) = cchiy CALL LTRANS(ione,ione,B(1),R,B(1)) CALL LTRANS(ione,ione,B(4),R,B(4)) CALL LTRANV(ione,Nc,Nc,R,Dldc(1)) CALL LTRANV(null,Nc,Nc,R,Dldc(3*Nc+1)) CALL LTRANV(null,Na,Ia,R,Dlda(1)) CALL LTRANV(null,Na,Ia,R,Dlda(3*Na+1)) CALL BNGEN(B) IF ( Grad.EQ.1 ) THEN CALL LTRANV(null,3,3,R,B(29)) CALL LTRANV(null,3,3,R,B(38)) CALL LTRANS(ione,ione,B(29),R,B(29)) CALL LTRANS(ione,ione,B(32),R,B(32)) CALL LTRANS(ione,ione,B(35),R,B(35)) CALL LTRANS(ione,ione,B(38),R,B(38)) CALL LTRANS(ione,ione,B(41),R,B(41)) CALL LTRANS(ione,ione,B(44),R,B(44)) CALL LTRANV(null,3*Nc,3*Nc,R,Dldc(6*Nc+1)) CALL LTRANV(null,3*Nc,3*Nc,R,Dldc(15*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(6*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(9*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(12*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(15*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(18*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(21*Nc+1)) ENDIF Ia = Ia + 3 ENDIF END !*==fdldno.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE FDLDNO(K,Na,Ia,Schix,Cchix,Schiy,Cchiy,Schiz,Cchiz,R,B, & Dlda) IMPLICIT NONE !*--FDLDNO3546 C*** Start of declarations inserted by SPAG DOUBLE PRECISION B , Cchix , Cchiy , Cchiz , Dlda , R , Schix , & Schiy , Schiz INTEGER i , Ia , ione , j , K , Na C*** End of declarations inserted by SPAG SAVE DIMENSION Dlda(*) , B(*) , R(*) DATA ione/1/ IF ( K.EQ.1 ) THEN i = Ia DO j = 1 , 6 Dlda(i+1) = 0.D0 Dlda(i+2) = 0.D0 Dlda(i+3) = 0.D0 i = i + Na ENDDO ELSE R(1) = 0.D0 R(2) = 0.D0 R(3) = 0.D0 R(4) = Cchix R(5) = -Schix R(6) = 0.D0 R(7) = 0.D0 R(8) = 0.D0 R(9) = 0.D0 CALL LTRANS(Na,ione,B(1),R,Dlda(Ia+1)) CALL LTRANS(Na,ione,B(4),R,Dlda(3*Na+Ia+1)) R(1) = 0.D0 R(2) = 0.D0 R(3) = 0.D0 R(4) = 0.D0 R(5) = 0.D0 R(6) = 0.D0 R(7) = Cchiy*Cchiz R(8) = Cchiy*Schiz R(9) = -Schiy CALL LTRANS(Na,ione,B(1),R,Dlda(Ia+2)) CALL LTRANS(Na,ione,B(4),R,Dlda(3*Na+Ia+2)) R(1) = 0.D0 R(2) = 0.D0 R(3) = 0.D0 R(4) = 0.D0 R(5) = 0.D0 R(6) = 0.D0 R(7) = -Schiy*Schiz R(8) = Schiy*Cchiz R(9) = 0.D0 CALL LTRANS(Na,ione,B(1),R,Dlda(Ia+3)) CALL LTRANS(Na,ione,B(4),R,Dlda(3*Na+Ia+3)) ENDIF END !*==tvn.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE TVN(Grad,K,Nc,Na,Ia,A,B,Dldc,Dlda,R) IMPLICIT NONE !*--TVN3602 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , B , Dlda , Dldc , R , slpx , slpy , slpz INTEGER Ia , ione , K , Na , Nc , null C*** End of declarations inserted by SPAG SAVE INTEGER Grad DIMENSION Dlda(*) , Dldc(*) , A(*) , B(*) , R(*) DATA null/0/ , ione/1/ IF ( K.GE.1 ) THEN slpx = A(Ia+1) slpy = A(Ia+2) slpz = A(Ia+3) CALL FDLDSL(K,Na,Ia,B,Dlda) R(1) = slpx R(2) = 0.D0 R(3) = 0.D0 R(4) = 0.D0 R(5) = slpy R(6) = 0.D0 R(7) = 0.D0 R(8) = 0.D0 R(9) = slpz CALL LTRANS(ione,ione,B(1),R,B(1)) CALL LTRANS(ione,ione,B(4),R,B(4)) CALL LTRANV(ione,Nc,Nc,R,Dldc(1)) CALL LTRANV(null,Nc,Nc,R,Dldc(3*Nc+1)) CALL LTRANV(null,Na,Ia,R,Dlda(1)) CALL LTRANV(null,Na,Ia,R,Dlda(3*Na+1)) CALL BNGEN(B) IF ( Grad.EQ.1 ) THEN CALL LTRANV(null,3,3,R,B(29)) CALL LTRANV(null,3,3,R,B(38)) CALL LTRANS(ione,ione,B(29),R,B(29)) CALL LTRANS(ione,ione,B(32),R,B(32)) CALL LTRANS(ione,ione,B(35),R,B(35)) CALL LTRANS(ione,ione,B(38),R,B(38)) CALL LTRANS(ione,ione,B(41),R,B(41)) CALL LTRANS(ione,ione,B(44),R,B(44)) CALL LTRANV(null,3*Nc,3*Nc,R,Dldc(6*Nc+1)) CALL LTRANV(null,3*Nc,3*Nc,R,Dldc(15*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(6*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(9*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(12*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(15*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(18*Nc+1)) CALL LTRANV(null,Nc,Nc,R,Dldc(21*Nc+1)) ENDIF Ia = Ia + 3 ENDIF END !*==fdldsl.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE FDLDSL(K,Na,Ia,B,Dlda) IMPLICIT NONE !*--FDLDSL3656 C*** Start of declarations inserted by SPAG DOUBLE PRECISION B , Dlda INTEGER i , Ia , j , K , Na C*** End of declarations inserted by SPAG SAVE DIMENSION Dlda(*) , B(*) i = Ia DO j = 1 , 6 Dlda(i+1) = 0.D0 Dlda(i+2) = 0.D0 Dlda(i+3) = 0.D0 i = i + Na ENDDO IF ( K.GT.1 ) THEN i = Ia Dlda(i+1) = B(1) i = i + Na Dlda(i+2) = B(2) i = i + Na Dlda(i+3) = B(3) i = i + Na Dlda(i+1) = B(4) i = i + Na Dlda(i+2) = B(5) i = i + Na Dlda(i+3) = B(6) ENDIF END !*==tbi.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE TBI(K,Na,Ia,A,B,Dlda) IMPLICIT NONE !*--TBI3688 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , B , biax , biay , biaz , Dlda INTEGER Ia , K , Na C*** End of declarations inserted by SPAG SAVE DIMENSION Dlda(*) , A(*) , B(*) IF ( K.GE.1 ) THEN biax = A(Ia+1) biay = A(Ia+2) biaz = A(Ia+3) CALL FDLDBI(K,Na,Ia,Dlda) B(1) = B(1) + biax B(2) = B(2) + biay B(3) = B(3) + biaz CALL BNGEN(B) Ia = Ia + 3 ENDIF END !*==fdldbi.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE FDLDBI(K,Na,Ia,Dlda) IMPLICIT NONE !*--FDLDBI3710 C*** Start of declarations inserted by SPAG DOUBLE PRECISION Dlda INTEGER i , Ia , j , K , Na C*** End of declarations inserted by SPAG SAVE DIMENSION Dlda(*) i = Ia DO j = 1 , 6 Dlda(i+1) = 0.D0 Dlda(i+2) = 0.D0 Dlda(i+3) = 0.D0 i = i + Na ENDDO IF ( K.GT.1 ) THEN i = Ia Dlda(i+1) = 1.D0 i = i + Na Dlda(i+2) = 1.D0 i = i + Na Dlda(i+3) = 1.D0 ENDIF END !*==ltrans.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE LTRANS(N,M,Q,R,S) IMPLICIT NONE !*--LTRANS3736 C*** Start of declarations inserted by SPAG INTEGER M , N DOUBLE PRECISION Q , q1 , q2 , q3 , R , S C*** End of declarations inserted by SPAG SAVE DIMENSION Q(*) , R(*) , S(*) q1 = Q(1) q2 = Q(M+1) q3 = Q(M+M+1) S(1) = q1*R(1) + q2*R(2) + q3*R(3) S(N+1) = q1*R(4) + q2*R(5) + q3*R(6) S(N+N+1) = q1*R(7) + q2*R(8) + q3*R(9) END !*==ltranv.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE LTRANV(Rfac,N,M,R,V) IMPLICIT NONE !*--LTRANV3753 C*** Start of declarations inserted by SPAG INTEGER M , N DOUBLE PRECISION R , V C*** End of declarations inserted by SPAG SAVE INTEGER Rfac DIMENSION R(*) , V(*) IF ( M.GT.0 ) THEN IF ( Rfac.EQ.1 ) THEN R(10) = R(4)/R(1) R(11) = R(5) - R(2)*R(10) R(12) = R(6) - R(3)*R(10) R(13) = R(7)/R(1) R(14) = (R(8)-R(2)*R(13))/R(11) R(15) = (R(9)-R(3)*R(13)) - R(12)*R(14) R(13) = R(13) - R(10)*R(14) ENDIF CALL R8VSCALE(1,M,R(1),V) CALL R8VLINKT(N+1,1,M,R(2),V,V) CALL R8VLINKT(N+N+1,1,M,R(3),V,V) CALL R8VSCALE(N+1,M,R(11),V) CALL R8VLINKT(1,N+1,M,R(10),V,V) CALL R8VLINKT(N+N+1,N+1,M,R(12),V,V) CALL R8VSCALE(N+N+1,M,R(15),V) CALL R8VLINKT(1,N+N+1,M,R(13),V,V) CALL R8VLINKT(N+1,N+N+1,M,R(14),V,V) ENDIF END !*==nshx.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 FUNCTION NSHX(Nmax,Nmin,Mmax,Mmin) IMPLICIT NONE !*--NSHX3785 C*** Start of declarations inserted by SPAG INTEGER kmax , Mmax , Mmin , Nmax , Nmin , NSHX C*** End of declarations inserted by SPAG kmax = Mmax + 1 NSHX = MAX(0,kmax**2-Mmin**2+MIN(Nmin,Mmin)**2-MIN(Nmin,kmax) & **2+(Nmax-Mmax-DIM(Nmin,kmax))*(2*Mmax+1) & +(DIM(Nmin,Mmin)-Nmax+Mmin-1)*MAX(0,2*Mmin-1)) END !*==nlpx.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 FUNCTION NLPX(Nmax,Mmax,Mmin) IMPLICIT NONE !*--NLPX3797 C*** Start of declarations inserted by SPAG INTEGER mdif , Mmax , Mmin , NLPX , Nmax C*** End of declarations inserted by SPAG mdif = MAX(0,MIN(Nmax,Mmax)-Mmin+1) NLPX = mdif*(mdif+1)/2 + mdif*DIM(Nmax,Mmax) + Mmin - 1 END !*==i8ssum.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 FUNCTION I8SSUM(Abeg,Alen,A) IMPLICIT NONE !*--I8SSUM3807 C*** Start of declarations inserted by SPAG INTEGER i , I8SSUM C*** End of declarations inserted by SPAG INTEGER aadr , Abeg , Alen INTEGER A(*) I8SSUM = 0 aadr = Abeg DO i = 1 , Alen I8SSUM = I8SSUM + A(aadr) aadr = aadr + 1 ENDDO END !*==i8vset.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE I8VSET(Abeg,Alen,S,A) IMPLICIT NONE !*--I8VSET3823 C*** Start of declarations inserted by SPAG INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , Alen , S INTEGER A(*) aadr = Abeg DO i = 1 , Alen A(aadr) = S aadr = aadr + 1 ENDDO END !*==i8vadd.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE I8VADD(Abeg,Bbeg,Cbeg,Vlen,A,B,C) IMPLICIT NONE !*--I8VADD3839 C*** Start of declarations inserted by SPAG INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , badr , Bbeg , cadr , Cbeg , Vlen INTEGER A(*) , B(*) , C(*) aadr = Abeg badr = Bbeg cadr = Cbeg DO i = 1 , Vlen C(cadr) = B(badr) + A(aadr) aadr = aadr + 1 badr = badr + 1 cadr = cadr + 1 ENDDO END !*==i8vadds.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE I8VADDS(Abeg,Bbeg,Vlen,S,A,B) IMPLICIT NONE !*--I8VADDS3859 C*** Start of declarations inserted by SPAG INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , badr , Bbeg , Vlen , S INTEGER A(*) , B(*) aadr = Abeg badr = Bbeg DO i = 1 , Vlen B(badr) = A(aadr) + S aadr = aadr + 1 badr = badr + 1 ENDDO END !*==i8vcum.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE I8VCUM(Abas,Abeg,Alen,A) IMPLICIT NONE !*--I8VCUM3877 C*** Start of declarations inserted by SPAG INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , Alen , Abas , aprv , acur INTEGER A(*) aprv = A(Abeg) A(Abeg) = Abas aadr = Abeg + 1 DO i = 1 , Alen - 1 acur = A(aadr) A(aadr) = A(aadr-1) + aprv aprv = acur aadr = aadr + 1 ENDDO END !*==i8vdel.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE I8VDEL(Abas,Abeg,Alen,A) IMPLICIT NONE !*--I8VDEL3897 C*** Start of declarations inserted by SPAG INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , Alen , Abas , aprv , acur INTEGER A(*) aprv = Abas aadr = Abeg DO i = 1 , Alen acur = A(aadr) A(aadr) = acur - aprv aprv = acur aadr = aadr + 1 ENDDO END !*==r8vset.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE R8VSET(Abeg,Alen,S,A) IMPLICIT NONE !*--R8VSET3916 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , S INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , Alen DIMENSION A(*) aadr = Abeg DO i = 1 , Alen A(aadr) = S aadr = aadr + 1 ENDDO END !*==r8sdot.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 FUNCTION R8SDOT(Abeg,Bbeg,Vlen,A,B) IMPLICIT NONE !*--R8SDOT3933 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , B , R8SDOT INTEGER i C*** End of declarations inserted by SPAG INTEGER aadr , Abeg , badr , Bbeg , Vlen DIMENSION A(*) , B(*) R8SDOT = 0.D0 aadr = Abeg badr = Bbeg DO i = 1 , Vlen R8SDOT = R8SDOT + A(aadr)*B(badr) aadr = aadr + 1 badr = badr + 1 ENDDO END !*==r8ssum.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 FUNCTION R8SSUM(Abeg,Alen,A) IMPLICIT NONE !*--R8SSUM3952 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , R8SSUM INTEGER i C*** End of declarations inserted by SPAG INTEGER aadr , Abeg , Alen DIMENSION A(*) R8SSUM = 0.D0 aadr = Abeg DO i = 1 , Alen R8SSUM = R8SSUM + A(aadr) aadr = aadr + 1 ENDDO END !*==r8slt.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE R8SLT(Abeg,Alen,S,A,J) IMPLICIT NONE !*--R8SLT3969 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , S INTEGER i , J C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , Alen DIMENSION A(*) aadr = Abeg DO i = 1 , Alen IF ( S.LT.A(aadr) ) THEN J = i - 1 GOTO 99999 ENDIF aadr = aadr + 1 ENDDO J = Alen 99999 END !*==r8vsub.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE R8VSUB(Abeg,Bbeg,Cbeg,Vlen,A,B,C) IMPLICIT NONE !*--R8VSUB3990 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , B , C INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , badr , Bbeg , cadr , Cbeg , Vlen DIMENSION A(*) , B(*) , C(*) aadr = Abeg badr = Bbeg cadr = Cbeg DO i = 1 , Vlen C(cadr) = B(badr) - A(aadr) aadr = aadr + 1 badr = badr + 1 cadr = cadr + 1 ENDDO END !*==r8vmul.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE R8VMUL(Abeg,Bbeg,Cbeg,Vlen,A,B,C) IMPLICIT NONE !*--R8VMUL4011 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , B , C INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , badr , Bbeg , cadr , Cbeg , Vlen DIMENSION A(*) , B(*) , C(*) aadr = Abeg badr = Bbeg cadr = Cbeg DO i = 1 , Vlen C(cadr) = B(badr)*A(aadr) aadr = aadr + 1 badr = badr + 1 cadr = cadr + 1 ENDDO END !*==r8vscale.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE R8VSCALE(Abeg,Alen,S,A) IMPLICIT NONE !*--R8VSCALE4032 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , S INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , Alen DIMENSION A(*) aadr = Abeg DO i = 1 , Alen A(aadr) = S*A(aadr) aadr = aadr + 1 ENDDO END !*==r8vscats.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE R8VSCATS(Qbeg,Qlen,S,Q,A) IMPLICIT NONE !*--R8VSCATS4049 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , S INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER qadr , Qbeg , Qlen INTEGER Q(*) DIMENSION A(*) qadr = Qbeg DO i = 1 , Qlen A(Q(qadr)) = S qadr = qadr + 1 ENDDO END !*==r8vlinkt.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE R8VLINKT(Abeg,Bbeg,Vlen,S,A,B) IMPLICIT NONE !*--R8VLINKT4067 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , B , S INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , badr , Bbeg , Vlen DIMENSION A(*) , B(*) aadr = Abeg badr = Bbeg DO i = 1 , Vlen B(badr) = B(badr) + S*A(aadr) aadr = aadr + 1 badr = badr + 1 ENDDO END !*==r8vlinkq.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE R8VLINKQ(Abeg,Bbeg,Cbeg,Vlen,S,A,B,C) IMPLICIT NONE !*--R8VLINKQ4086 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , B , C , S INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , badr , Bbeg , cadr , Cbeg , Vlen DIMENSION A(*) , B(*) , C(*) aadr = Abeg badr = Bbeg cadr = Cbeg DO i = 1 , Vlen C(cadr) = C(cadr) + S*A(aadr)*B(badr) aadr = aadr + 1 badr = badr + 1 cadr = cadr + 1 ENDDO END !*==r8vgathp.spg processed by SPAG 6.50Rc at 13:25 on 30 Jul 2010 SUBROUTINE R8VGATHP(Abeg,Ainc,Bbeg,Blen,A,B) IMPLICIT NONE !*--R8VGATHP4107 C*** Start of declarations inserted by SPAG DOUBLE PRECISION A , B INTEGER i C*** End of declarations inserted by SPAG SAVE INTEGER aadr , Abeg , Ainc , badr , Bbeg , Blen DIMENSION A(*) , B(*) aadr = Abeg badr = Bbeg DO i = 1 , Blen B(badr) = A(aadr) aadr = aadr + Ainc badr = badr + 1 ENDDO END