!====================================================================== !====================================================================== ! main lem_v5.f !====================================================================== ! lem_v5: ! Galileo EPD Web Tool (LEMMS channels only) ! ! Features: ! (1). Assumes that a user-generated, formatted, parameter ! input file 'lem_v3' exists ! (2). Assumes a Galileo NAV file with pre-determined contents ! See comments below. ! (3). Assumes a EPD data file with pre-determined contents ! (4). Restored (from lem_v1.f) useful diagnostic variables and ! pitch cosine in comoving frame ! See comments below. ! (5). Code added to perform user-defined plasma convection ! velocity options jtvr=2 ! ! Fortran code written for Linux systems. Compile with: ! g77 -o lem_v5 -O -W lem_v5.f ! or ! gfortran -o lem_v5 -O -W lem_v5.f ! ! History: ! Created: 23-Feb-07 rbd ! Changed: 06-Mar-07 rbd ! v3 has TP, TO, and TS channels, 6/23/08, cp ! Changed: 11-Jul-08 rbd (see point (4) above) ! Changed: 24-Jul-08 rbd (flux -> rate_bgc, v4 ->v5) ! Changed: 28-Jul-08 rbd (see (5) above) !====================================================================== !====================================================================== program lem_v5 character*80 dfnav, dfdat, dfout integer*4 mdate(3), intim_nav(5), intim_dat(5), ms_rec(16) real*8 dsh_n, dsh_d, dsh_s, dsh_e real*4 tld(7), ctld(7), stld(7) real*4 pld(16), cpld(16), spld(16) real*4 xld_mn(7,16), yld_mn(7,16), zld_mn(7,16) real*4 gal_jse(3), corot_jse(3), galv_sc(3) real*4 corot_sc(3), jup_sc(3), jup_sc_unit(3), sun_sc(3) real*4 r_hat(3), c_hat(3), t_hat(3) real*4 pamag(16), bxd(16), byd(16), bzd(16) ! Auxilliary 16-element vector to hold input 16-sector ... ! ... background-corrected rate data (cnts/sec) and ... real*4 A0_rate(16),A1_rate(16),A2_rate(16) real*4 A3_rate(16),A4_rate(16),A5_rate(16) real*4 TP1_rate(16),TP2_rate(16),TP3_rate(16) real*4 TO2_rate(16),TO3_rate(16),TO4_rate(16) real*4 TS1_rate(16),TS2_rate(16),TS3_rate(16) ! ... raw rate data (cnts/sec): real*4 A0_rrate(16),A1_rrate(16),A2_rrate(16) real*4 A3_rrate(16),A4_rrate(16),A5_rrate(16) real*4 TP1_rrate(16),TP2_rrate(16),TP3_rrate(16) real*4 TO2_rrate(16),TO3_rrate(16),TO4_rrate(16) real*4 TS1_rrate(16),TS2_rrate(16),TS3_rrate(16) !====================================================================== ! LEMMS: !----------------------- ! r_lem(i:15, k:16): ! EPD LEMMS background-corrected rate channel i, sector k ! rr_lem(i:15, j: 4, k:16): ! EPD LEMMS raw rate channel i, species j, sector k !---------------------------------------------------------------------- ! eg A0 is i=1, A1 is 2, A2 is 3, etc. !---------------------------------------------------------------------- real*4 r_lem(15,16), rr_lem(15,16) ! 16-sect/15-ch (input) real*4 rav_lem(15) ! sect avg., 15 ch. (derived) real*4 rtf(15), egm(15), flx(15) ! re-derived flux integer*4 nr_lem(15), nf_lem(15) ! number good pnts. ch. real*4 esp_ind(15) ! energy spectral indices !---------------------------------------------------------------------- ! Geometric factors gf_lem(i:15) of LEMMS channel i: ! Follow Mauk et al., JGR, 109, A09S12, 2004 and base TOF channels ! on his formula in the appendix divided by 2. If you divide formula ! by 2.0, then formula and plot match. !---------------------------------------------------------------------- real*4 gf_lem(15) data gf_lem/0.0052,0.0056,0.006, ! A0, A1, A2 * 0.006,0.006,0.006, ! A3, A4, A5 * 0.0016,0.0014,0.0007, ! TP1-TP3 * 0.0028,0.0034,0.0037, ! TO2-TO4 * 0.0022,0.003,0.0035/ ! TS1-TS3 !---------------------------------------------------------------------- ! Energy passbands elo_lem(i:15, j:4) and ehi_lem(i:15, j:4) (keV) ! and efficiencies of LEMMS channel i, species j (2=H, 2=He, 3=O, 4=S) !---------------------------------------------------------------------- real*4 elo_lem(15,4), ehi_lem(15,4), eff_lem(15,4) real*4 atn(4), atm(4) data atn/ 1.0, 2.0, 8.0, 16.0/ ! atomic number data atm/ 1.0, 4.0, 16.0, 32.0/ ! mass number !----------------------- ! Species H (j=1): !----------------------- data elo_lem(1,1), ehi_lem(1,1), eff_lem(1,1)/ ! A0/H * 22.0, 42.0, 1.00/, * elo_lem(2,1), ehi_lem(2,1), eff_lem(2,1)/ ! A1/H * 42.0, 65.0, 1.00/, * elo_lem(3,1), ehi_lem(3,1), eff_lem(3,1)/ ! A2/H * 65.0, 120.0, 1.00/, * elo_lem(4,1), ehi_lem(4,1), eff_lem(4,1)/ ! A3/H * 120.0, 280.0, 1.00/, * elo_lem(5,1), ehi_lem(5,1), eff_lem(5,1)/ ! A4/H * 280.0, 515.0, 1.00/, * elo_lem(6,1), ehi_lem(6,1), eff_lem(6,1)/ ! A5/H * 515.0, 825.0, 1.00/, * elo_lem(7,1), ehi_lem(7,1), eff_lem(7,1)/ ! TP1/H * 130.0, 220.0, 1.0/ * elo_lem(8,1), ehi_lem(8,1), eff_lem(8,1)/ ! TP2/H * 220.0, 540.0, 1.0/ * elo_lem(9,1), ehi_lem(9,1), eff_lem(9,1)/ ! TP3/H * 540.0, 1140.0, 1.0/ * elo_lem(10,1), ehi_lem(10,1), eff_lem(10,1)/ ! TO2/H * 0.0, 0.0, 0.0/ * elo_lem(11,1), ehi_lem(11,1), eff_lem(11,1)/ ! TO3/H * 0.0, 0.0, 0.0/ * elo_lem(12,1), ehi_lem(12,1), eff_lem(12,1)/ ! TO4/H * 0.0, 0.0, 0.0/ * elo_lem(13,1), ehi_lem(13,1), eff_lem(13,1)/ ! TS1/H * 0.0, 0.0, 0.0/ * elo_lem(14,1), ehi_lem(14,1), eff_lem(14,1)/ ! TS2/H * 0.0, 0.0, 0.0/ * elo_lem(15,1), ehi_lem(15,1), eff_lem(15,1)/ ! TS3/H * 0.0, 0.0, 0.0/ !----------------------- !----------------------- ! Species He (j=2): !----------------------- data elo_lem(1,2), ehi_lem(1,2), eff_lem(1,2)/ ! A0/He * 32.0, 50.0, 1.00/, * elo_lem(2,2), ehi_lem(2,2), eff_lem(2,2)/ ! A1/He * 50.0, 71.0, 1.00/, * elo_lem(3,2), ehi_lem(3,2), eff_lem(3,2)/ ! A2/He * 71.0, 129.0, 1.00/, * elo_lem(4,2), ehi_lem(4,2), eff_lem(4,2)/ ! A3/He * 129.0, 299.0, 1.00/, * elo_lem(5,2), ehi_lem(5,2), eff_lem(5,2)/ ! A4/He * 299.0, 560.0, 1.00/, * elo_lem(6,2), ehi_lem(6,2), eff_lem(6,2)/ ! A5/He * 560.0, 845.0, 1.00/, * elo_lem(7,2), ehi_lem(7,2), eff_lem(7,2)/ ! TP1/He * 0.0, 0.0, 0.0/ * elo_lem(8,2), ehi_lem(8,2), eff_lem(8,2)/ ! TP2/He * 0.0, 0.0, 0.0/ * elo_lem(9,2), ehi_lem(9,2), eff_lem(9,2)/ ! TP3/He * 0.0, 0.0, 0.0/ * elo_lem(10,2), ehi_lem(10,2), eff_lem(10,2)/ ! TO2/He * 0.0, 0.0, 0.0/ * elo_lem(11,2), ehi_lem(11,2), eff_lem(11,2)/ ! TO3/He * 0.0, 0.0, 0.0/ * elo_lem(12,2), ehi_lem(12,2), eff_lem(12,2)/ ! TO4/He * 0.0, 0.0, 0.0/ * elo_lem(13,2), ehi_lem(13,2), eff_lem(13,2)/ ! TS1/He * 0.0, 0.0, 0.0/ * elo_lem(14,2), ehi_lem(14,2), eff_lem(14,2)/ ! TS2/He * 0.0, 0.0, 0.0/ * elo_lem(15,2), ehi_lem(15,2), eff_lem(15,2)/ ! TS3/He * 0.0, 0.0, 0.0/ !----------------------- !----------------------- ! Species O (j=3): !----------------------- data elo_lem(1,3), ehi_lem(1,3), eff_lem(1,3)/ ! A0/O * 45.0, 107.0, 0.70/, * elo_lem(2,3), ehi_lem(2,3), eff_lem(2,3)/ ! A1/O * 93.0, 133.0, 0.80/, * elo_lem(3,3), ehi_lem(3,3), eff_lem(3,3)/ ! A2/O * 127.0, 210.0, 0.90/, * elo_lem(4,3), ehi_lem(4,3), eff_lem(4,3)/ ! A3/O * 200.0, 410.0, 1.00/, * elo_lem(5,3), ehi_lem(5,3), eff_lem(5,3)/ ! A4/O * 415.0, 650.0, 1.00/, * elo_lem(6,3), ehi_lem(6,3), eff_lem(6,3)/ ! A5/O * 650.0, 1000.0, 1.00/, * elo_lem(7,3), ehi_lem(7,3), eff_lem(7,3)/ ! TP1/O * 0.0, 0.0, 0.0/ * elo_lem(8,3), ehi_lem(8,3), eff_lem(8,3)/ ! TP2/O * 0.0, 0.0, 0.0/ * elo_lem(9,3), ehi_lem(9,3), eff_lem(9,3)/ ! TP3/O * 0.0, 0.0, 0.0/ * elo_lem(10,3), ehi_lem(10,3), eff_lem(10,3)/ ! TO2/O * 416.0, 816.0, 1.0/ * elo_lem(11,3), ehi_lem(11,3), eff_lem(11,3)/ ! TO3/O * 816.0, 1792.0, 1.0/ * elo_lem(12,3), ehi_lem(12,3), eff_lem(12,3)/ ! TO4/O * 1792.0, 8992.0, 1.0/ * elo_lem(13,3), ehi_lem(13,3), eff_lem(13,3)/ ! TS1/O * 0.0, 0.0, 0.0/ * elo_lem(14,3), ehi_lem(14,3), eff_lem(14,3)/ ! TS2/O * 0.0, 0.0, 0.0/ * elo_lem(15,3), ehi_lem(15,3), eff_lem(15,3)/ ! TS3/O * 0.0, 0.0, 0.0/ !----------------------- !----------------------- ! Species S (j=4): !----------------------- data elo_lem(1,4), ehi_lem(1,4), eff_lem(1,4)/ ! A0/S * 82.0, 157.0, 0.55/, * elo_lem(2,4), ehi_lem(2,4), eff_lem(2,4)/ ! A1/S * 137.0, 180.0, 0.72/, * elo_lem(3,4), ehi_lem(3,4), eff_lem(3,4)/ ! A2/S * 167.0, 270.0, 0.85/, * elo_lem(4,4), ehi_lem(4,4), eff_lem(4,4)/ ! A3/S * 250.0, 480.0, 1.00/, * elo_lem(5,4), ehi_lem(5,4), eff_lem(5,4)/ ! A4/S * 480.0, 750.0, 1.00/, * elo_lem(6,4), ehi_lem(6,4), eff_lem(6,4)/ ! A5/S * 750.0, 1150.0, 1.00/, * elo_lem(7,4), ehi_lem(7,4), eff_lem(7,4)/ ! TP1/S * 0.0, 0.0, 0.0/ * elo_lem(8,4), ehi_lem(8,4), eff_lem(8,4)/ ! TP2/S * 0.0, 0.0, 0.0/ * elo_lem(9,4), ehi_lem(9,4), eff_lem(9,4)/ ! TP3/S * 0.0, 0.0, 0.0/ * elo_lem(10,4), ehi_lem(10,4), eff_lem(10,4)/ ! TO2/S * 0.0, 0.0, 0.0/ * elo_lem(11,4), ehi_lem(11,4), eff_lem(11,4)/ ! TO3/S * 0.0, 0.0, 0.0/ * elo_lem(12,4), ehi_lem(12,4), eff_lem(12,4)/ ! TO4/S * 0.0, 0.0, 0.0/ * elo_lem(13,4), ehi_lem(13,4), eff_lem(13,4)/ ! TS1/S * 512.0, 960.0, 1.0/ * elo_lem(14,4), ehi_lem(14,4), eff_lem(14,4)/ ! TS2/S * 960.0, 1984.0, 1.0/ * elo_lem(15,4), ehi_lem(15,4), eff_lem(15,4)/ ! TS3/S * 1984.0, 9920.0, 1.0/ !----------------------- !====================================================================== data dtr, rtd/ 1.745329e-02, 5.729578e+01/ data pm, qe, sol/ 1.67e-24, 4.8e-10, 3.0e+10/ data c_ergtkev, c_kevterg/ 6.2422e+08, 1.602e-09/ data omega_J, rad_J/ 1.7585e-04, 7.14e+04/ ! rad/sec, km ! npts1=0 avg1=0.0 sumsq1=0.0 npts2=0. avg2=0. sumsq2=0. ! ! Current date and time call idate (mdate) !---------------------------------------------------------------------- ! Open parameter input file (user-generated) !---------------------------------------------------------------------- open (unit=1, file='lem_v5.inp', form='formatted', status='old') !---------------------------------------------------------------------- ! Open Galileo navigation file: ! nhnav = number of header info. records prior to first nav record ! dfnav = name of formatted nav file !---------------------------------------------------------------------- do i=1,4 read (1,*) enddo read (1,*) nhnav read (1,10) dfnav 10 format (a80) open (unit=2, file=dfnav, form='formatted', status='old') ! Read and discard nav file header records do i=1,nhnav read (2,*) enddo !---------------------------------------------------------------------- ! Open EPD data file: ! nhdat = number of header info. records prior to first data record ! dfdat = name of formatted data file !---------------------------------------------------------------------- read (1,*) read (1,*) read (1,*) nhdat read (1,10) dfdat open (unit=4, file=dfdat, form='formatted', status='old') ! Read and discard EPD data file header records do i=1,nhdat read (4,*) enddo !---------------------------------------------------------------------- ! Enter name of output disk file for listing results: ! ndo = 0: do not write results to output file ! 1: write results to output file ! dfout = name of formatted output data file !---------------------------------------------------------------------- read (1,*) read (1,*) read (1,*) ndo read (1,10) dfout if (ndo .eq. 1) then open (unit=8, file=dfout, form='formatted', status='new') endif !---------------------------------------------------------------------- ! Enter start and end times to process (all integers): ! start time: ksy, ksd, ksh, ksm, kss ! end time: key, ked, keh, kem, kes !---------------------------------------------------------------------- do i=1,3 read (1,*) enddo read (1,*) ksy, ksd, ksh, ksm, kss read (1,*) read (1,*) key, ked, keh, kem, kes ! Get decimal sequential hour form integer time values call c_seqt (1, dsh_s, idsh_s, ksy, ksd, ksh, ksm, kss) call c_seqt (1, dsh_e, idsh_e, key, ked, keh, kem, kes) !---------------------------------------------------------------------- ! Enter EPD channel (jcht) to tranform and ions species (jion): ! jcht: e.g., 1=A0, 2=A1, see input file for choices lem_v3.inp ! jion: 1=H, 2=He, 3=O, 4=S !---------------------------------------------------------------------- do i=1,7 read (1,*) enddo read (1,*) jcht, jion !---------------------------------------------------------------------- ! Open file to write output from this run !---------------------------------------------------------------------- ! open (unit=27, file='ratout.txt', form='formatted', status='new') ! open (unit=28, file='gamout.txt', form='formatted', status='new') ! open (unit=29, file='traout.txt', form='formatted', status='new') kk=0 ! initialize for write if ((jcht .lt. 1) .or. (jcht .gt. 15)) then print *, 'Invalid channel definition - stop.' stop endif if ((jion .lt. 1) .or. (jion .gt. 4)) then print *, 'Invalid species definition - stop.' stop endif ! Check if channel 'jcht' has nonzero efficiency for species 'jion' if (eff_lem(jcht,jion) .lt. 0.10) then print *, ' Species not measured by this channel - stop' stop endif !---------------------------------------------------------------------- ! Get mean ion speed for this channel/species combination !---------------------------------------------------------------------- Tch = Sqrt(elo_lem(jcht,jion)*ehi_lem(jcht,jion)) ! mean KE (keV) vch = Sqrt(2.0*Tch*c_kevterg/(atm(jion)*pm))/1.0e+05 ! mean v, km/s ! write(6,*) 'Channel number and kinetic energy:', jion, tch !---------------------------------------------------------------------- ! Energy Spectum Index Options: ! d_sp_ind = default index to use when not calculable from EPD data !---------------------------------------------------------------------- do i=1,4 read (1,*) enddo read (1,*) d_esp_ind !---------------------------------------------------------------------- ! Plasma Convection Velocity Options (all speeds in km/s): ! jvtr : 0 = Rigid corotation direction and magnitude[*cfact=0.0-1.0] ! : 1 = Rigid corotation direction, user-defined magnitude vc_cus ! : 2 = User-defined convection velocity (vt_mag, vt_clt, vt_azm) ! ! For option 2: ! A local cartesian coordinate system centered on and moving with ! Galileo is defined with: ! c_hat[1:3] = unit vector in local corotation direction ! r_hat[1:3] = unit radial vector from Jupiter to Galileo ! t_hat[1:3] = unit vector c_hat x r_hat ! c_hat axis is polar axis from which polar angle vt_clt is defined ! (vt_clt = [0,180] deg) ! r_hat-t_hat plane is plane in which azimuth is measured, with ! vt_azm increasing from r_hat toward t_hat in right-hand ! sense along polar axis (vt_azm = [0,360)) !---------------------------------------------------------------------- do i=1,6 read (1,*) enddo read (1,*) jvtr, cfact, vc_cus, vt_mag, vt_clt, vt_azm close (unit=1) ! ! Now we have finished reading the input file !---------------------------------------------------------------------- ! Calculate look angles of EPD motor positions 1-7 ! relative to spacecraft-frame cartesian system K[x,y,z]. !---------------------------------------------------------------------- do m=1,7 tld(m) = Float(m-1)*30.6 ctld(m) = Cos(tld(m)*dtr) stld(m) = Sin(tld(m)*dtr) enddo !---------------------------------------------------------------------- ! Calculate look angles of EPD spin sectors 0-15 ! relative to spacecraft-frame cartesian system K[x,y,z] !---------------------------------------------------------------------- do n=1,16 pld(n) = Float(2*n-1)*11.25 ! 0-360 deg cpld(n) = Cos(pld(n)*dtr) spld(n) = Sin(pld(n)*dtr) enddo !---------------------------------------------------------------------- ! Calculate K[x,y,z] comps. of look-direction unit vector ! ld_mn = [xld_mn, yld_nm, zld_mn]. {Note: Unit vector of ! incident particle is then -ld_mn.} !---------------------------------------------------------------------- do m=1,7 do n=1,16 xld_mn(m,n) = stld(m)*cpld(n) yld_mn(m,n) = stld(m)*spld(n) zld_mn(m,n) = ctld(m) enddo enddo close (unit=1) ! close parameter input file n_nav = 0 ! initialize nav record counter n_dat = 0 ! initialize data record counter !====================================================================== ! Nav Info. and EPD Data Input Loop !====================================================================== 1000 continue !---------------------------------------------------------------------- ! Read a NAV file record (one record per time tag) !---------------------------------------------------------------------- ! intim_nav(1) = year (4-digit) [i] ! intim_nav(2) = doy (1-365 or 366) [i] ! intim_nav(3) = hour (0-23) [i] ! intim_nav(4) = min (0-59) [i] ! intim_nav(5) = sec (0-59) [i] ! fmpcn = current motor position (0-7) [r] ! fmpln = motor position of previous spin (0-7) [r] ! fnsbn = begin sector (0-63) [r] ! fnsbn = end sector (0-63) [r] ! gal_jse[1:3] = position vector of Galileo wrt Jupiter in ! JSE system (km) ! corot_jse[1:3] = corotation velocity at Galileo in JSE ! system (km/s) ! galv_sc[1:3] = Galileo velocity in S/C system (km/s) ! corot_sc[1:3] = corotation unit vector direction at Galileo ! in S/C system ! jup_sc[1:3] = Position vector of Jupiter wrt Galileo in ! S/C system (km) ! sun_sc[1:3] = Position vector of Sun wrt Galileo in ! S/C system (km) !---------------------------------------------------------------------- 1100 read (2,*, end=2000) intim_nav, fmpcn, fmpln, * fnsbn, fnsen, gal_jse, corot_jse, galv_sc, * corot_sc, jup_sc, sun_sc mpcn = Ifix(fmpcn) ! current motor position [i] mpln = Ifix(fmpln) ! previous motor position [i] nsbn = Ifix(fnsbn) ! start sector (0-63) [i] nsen = Ifix(fnsen) ! end sector (0-63) [i] n_nav = n_nav + 1 ! NAV record counter if (nsbn .ne. nsen) then write (6,1110) n_nav 1110 format (' Start & end sectors unequal: NAV record ', i6/) goto 1100 endif kny = intim_nav(1) knd = intim_nav(2) knh = intim_nav(3) knm = intim_nav(4) kns = intim_nav(5) call c_seqt (1, dsh_n, idsh_n, kny, knd, knh, knm, kns) !---------------------------------------------------------------------- ! Read EPD data file records (multiple records per time tag) !---------------------------------------------------------------------- ! intim_dat(1) = year [i] ! intim_dat(2) = doy [i] ! intim_dat(3) = hour [i] ! intim_dat(4) = min [i] ! intim_dat(5) = sec [i] ! mpcd = current motor positon (0-7) [i] ! nsbd = begin sector of spin (0-63) [i] ! nebd = end sector of spin (0-63) [i] ! pamag(16) = pitch angles of the 16 macro-sectors (0-180 deg) [r] ! *NOTE: The 16 pang elememts are always listed in order ! increasing macro-sector number 0-15 ! bxd(16) = x-component of the B (nT) wrt S/C system K[x,y,z] ! byd(16) = y-component of the B (nT) wrt S/C system K[x,y,z] ! bzd(16) = z-component of the B (nT) wrt S/C system K[x,y,z] !---------------------------------------------------------------------- 1200 read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, pamag read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, bxd read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, byd read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, bzd read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A0_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A1_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A2_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A3_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A4_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A5_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TP1_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TP2_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TP3_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TO2_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TO3_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TO4_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TS1_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TS2_rate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TS3_rate ! read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A0_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A1_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A2_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A3_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A4_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, A5_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TP1_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TP2_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TP3_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TO2_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TO3_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TO4_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TS1_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TS2_rrate read (4,*, end=3000) intim_dat, mpcd, nsbd, nsed, TS3_rrate kdy = intim_dat(1) kdd = intim_dat(2) kdh = intim_dat(3) kdm = intim_dat(4) kds = intim_dat(5) call c_seqt (1, dsh_d, idsh_d, kdy, kdd, kdh, kdm, kds) n_dat = n_dat + 1 ! EPD data record counter !---------------------------------------------------------------------- ! Test for nav file and data file time synchonicity. ! Check for equivalence of integer year, doy, horr, min, and sec ! between the current NAV an EPD data record. !---------------------------------------------------------------------- if (dsh_n .lt. dsh_d) goto 1100 if (dsh_n .gt. dsh_d) goto 1200 !---------------------------------------------------------------------- ! Check nav file times against chosen start and end times. !---------------------------------------------------------------------- if (dsh_n .lt. dsh_s) goto 1000 if (dsh_n .gt. dsh_e) goto 4000 !////////////////////////// SAMPLE OUTPUT ///////////////////////////// write (6,1400) intim_nav, n_nav, * intim_dat, n_dat 1400 format (/' nav: yrn,doyn,hhn,mmn,ssn,n_nav: ', 2i5,3i3,i5/ * ' dat: yrd,doyd,hhd,mmd,ssd,n_dat: ', 2i5,3i3,i5) gjr = Sqrt(gal_jse(1)**2 + gal_jse(2)**2 + gal_jse(3)**2) crj = Sqrt(corot_jse(1)**2 + corot_jse(2)**2 + corot_jse(3)**2) crs = Sqrt(corot_sc(1)**2 + corot_sc(2)**2 + corot_sc(3)**2) gav = Sqrt(galv_sc(1)**2 + galv_sc(2)**2 + galv_sc(3)**2) jsc = Sqrt(jup_sc(1)**2 + jup_sc(2)**2 + jup_sc(3)**2) ! unit vector to Jupiter from Galileo = -unit radial from J to G jup_sc_unit(1) = jup_sc(1)/jsc jup_sc_unit(2) = jup_sc(2)/jsc jup_sc_unit(3) = jup_sc(3)/jsc ! test that unit vector from Galileo to Jupiter dotted with unit ! corotation vector at Glileo vanishes (within single-precion) csc_d_jscu = corot_sc(1)*jup_sc_unit(1) * + corot_sc(2)*jup_sc_unit(2) * + corot_sc(3)*jup_sc_unit(3) !---------------------------------------------------------------------- ! rigid corotation speed in spin equatorial plane at radius gjr crv = omega_J*gjr*rad_J !---------------------------------------------------------------------- write (6,1410) mpcn, mpln, nsbn, nsen, * gal_jse(1)*rad_J, * gal_jse(2)*rad_J, * gal_jse(3)*rad_J, * gjr*rad_J, galv_sc, gav, * corot_jse, crj, crv, * corot_sc, crs, jup_sc_unit, * csc_d_jscu, sun_sc, * bxd, byd, bzd, A1_rate 1410 format (' ', 4i5/ * ' gal_jse: ', 4(1pe11.3)/ * ' galv_sc: ', 4(1pe11.3)/ * ' corot_jse: ', 5(1pe11.3)/ * ' corot_sc: ', 4(1pe11.3)/ * 'jup_sc_unit: ', 3(1pe11.3)/ * ' csc_d_jscu: ', 1(1pe13.5)/ * ' sun_sc: ', 3(1pe11.3)/ * ' Bx: ', 4(1pe11.3)/ * ' ', 4(1pe11.3)/ * ' ', 4(1pe11.3)/ * ' ', 4(1pe11.3)/ * ' By: ', 4(1pe11.3)/ * ' ', 4(1pe11.3)/ * ' ', 4(1pe11.3)/ * ' ', 4(1pe11.3)/ * ' Bz: ', 4(1pe11.3)/ * ' ', 4(1pe11.3)/ * ' ', 4(1pe11.3)/ * ' ', 4(1pe11.3)/ * ' A1_rate: ', 4(1pe11.3)/ * ' ', 4(1pe11.3)/ * ' ', 4(1pe11.3)/ * ' ', 4(1pe11.3)) !////////////////////////////////////////////////////////////////////// !---------------------------------------------------------------------- ! Process current record !---------------------------------------------------------------------- ! Ignore data in motor position 0. Read in a new record. if (mpcn .eq. 0) goto 1000 !---------------------------------------------------------------------- ! Fill in 2D arrays r_lem & rr_lem with vector (sectored) channel data !---------------------------------------------------------------------- do j=1,16 ! Sectored rates r_lem(1,j) = A0_rate(j) r_lem(2,j) = A1_rate(j) r_lem(3,j) = A2_rate(j) r_lem(4,j) = A3_rate(j) r_lem(5,j) = A4_rate(j) r_lem(6,j) = A5_rate(j) r_lem(7,j) = TP1_rate(j) r_lem(8,j) = TP2_rate(j) r_lem(9,j) = TP3_rate(j) r_lem(10,j) = TO2_rate(j) r_lem(11,j) = TO3_rate(j) r_lem(12,j) = TO4_rate(j) r_lem(13,j) = TS1_rate(j) r_lem(14,j) = TS2_rate(j) r_lem(15,j) = TS3_rate(j) ! Sectored fluxes rr_lem(1,j) = A0_rrate(j) rr_lem(2,j) = A1_rrate(j) rr_lem(3,j) = A2_rrate(j) rr_lem(4,j) = A3_rrate(j) rr_lem(5,j) = A4_rrate(j) rr_lem(6,j) = A5_rrate(j) rr_lem(7,j) = TP1_rrate(j) rr_lem(8,j) = TP2_rrate(j) rr_lem(9,j) = TP3_rrate(j) rr_lem(10,j) = TO2_rrate(j) rr_lem(11,j) = TO3_rrate(j) rr_lem(12,j) = TO4_rrate(j) rr_lem(13,j) = TS1_rrate(j) rr_lem(14,j) = TS2_rrate(j) rr_lem(15,j) = TS3_rrate(j) enddo ! ! Write the unprocessed intensity of the channel in question ! to the output file number 27, named above ! We do our own conversion from rate to intensity here fact1=ehi_lem(jcht,jion) - elo_lem(jcht,jion) fact2=gf_lem(jcht)*eff_lem(jcht,jion)*fact1 do jj=1,16 if (r_lem(jcht,jj).gt.0) then aj=r_lem(jcht,jj)/fact2 ! This is the unprocessed intensity, aj write(27,*) n_dat,kdy,kdd,kdh,kdm,kds,aj npts1=npts1+1 ! we will average the logs for our metric so that ! no point gets too much weight aja=alog10(aj) avg1=avg1+aja sumsq1=sumsq1+(aja*aja) endif enddo ! !---------------------------------------------------------------------- ! Calculate sector-averaged rates and fluxes !---------------------------------------------------------------------- i1 = 1 ! Have the capacity to do all channels i2 = 15 ! if wrong species will get NAN for tof channels ! Initialize accumulation buffers do i=i1,i2 rav_lem(i) = 0.0 nr_lem(i) = 0 nr_lem(i) = 0 enddo do i=i1,i2 do j=1,16 ! Accumulate rates rr = r_lem(i,j) if ((rr .gt. 1.e-10) .and. (rr .lt. 1.e+10)) then rav_lem(i) = rav_lem(i) + rr nr_lem(i) = nr_lem(i) + 1 endif enddo ! Sector-averaged rate if (nr_lem(i) .gt. 0) then rav_lem(i) = rav_lem(i)/Float(nr_lem(i)) else rav_lem(i) = -99.00 endif ! Derive separate sector-averaged fluxes for chosen ion species (jion) rtf(i) = gf_lem(i)*eff_lem(i,jion)* ! rate->flux * (ehi_lem(i,jion) - elo_lem(i,jion)) egm(i) = Sqrt(ehi_lem(i,jion)*elo_lem(i,jion)) ! log-mean if (rav_lem(i) .gt. 0.0) then flx(i) = rav_lem(i)/rtf(i) else flx(i) = -99.0 endif enddo !====================================================================== ! Get estimates for energy spectral indices !---------------------------------------------------------------------- ! Chans. A0-A5 (6 chans, 5 indices) ! esp_ind(1) = Index using A0 & A1 [i=1 & i=2] ! esp_ind(2) = Index using A1 & A2 [i=2 & i=3] ! esp_ind(3) = Index using A2 & A3 [i=3 & i=4] ! esp_ind(4) = Index using A3 & A4 [i=4 & i=5] ! esp_ind(5) = Index using A4 & A5 [i=5 & i=6] ! temp esp_ind(6) = defn.: esp_ind(5) !---------------------------------------------------------------------- do i=1,5 if ((flx(i) .gt. 0.0) .and. (flx(i+1) .gt. 0.0)) then esp_ind(i) = -Log(flx(i)/flx(i+1))/Log(egm(i)/egm(i+1)) endif enddo esp_ind(6) = esp_ind(5) ! temp.: assign A5 the A4 spec. ind. !---------------------------------------------------------------------- ! Chans. TP1-TP3 (3 chans., 2 indices) ! esp_ind(7) = Index using TP1 & TP2 [i=7 & i=8] ! esp_ind(8) = Index using TP2 & TP3 [i=8 & i=9] ! temp esp_ind(9) = defn.: esp_ind(8) !---------------------------------------------------------------------- do i=7,8 if ((flx(i) .gt. 0.0) .and. (flx(i+1) .gt. 0.0)) then esp_ind(i) = -Log(flx(i)/flx(i+1))/Log(egm(i)/egm(i+1)) endif enddo esp_ind(9) = esp_ind(8) ! temp.: assign TP3 the TP2 spec. ind. !---------------------------------------------------------------------- ! Chans. TO2-TO4 (3 chans., 2 indices) ! esp_ind(10) = Index using TO2 & TO3 [i=10 & i=11] ! esp_ind(11) = Index using TO3 & TO4 [i=11 & i=12] ! temp esp_ind(12) = defn.: esp_ind(11) !---------------------------------------------------------------------- do i=10,11 if ((flx(i) .gt. 0.0) .and. (flx(i+1) .gt. 0.0)) then esp_ind(i) = -Log(flx(i)/flx(i+1))/Log(egm(i)/egm(i+1)) endif enddo esp_ind(12) = esp_ind(11) ! temp.: assign TO4 the TO3 spec. ind. !---------------------------------------------------------------------- ! Chans. TS1-TS3 (3 chans., 2 indices) ! esp_ind(13) = Index using TS1 & TS2 [i=13 & i=14] ! esp_ind(14) = Index using TS2 & TS3 [i=14 & i=15] ! temp esp_ind(15) = defn.: esp_ind(14) !---------------------------------------------------------------------- do i=13,14 if ((flx(i) .gt. 0.0) .and. (flx(i+1) .gt. 0.0)) then esp_ind(i) = -Log(flx(i)/flx(i+1))/Log(egm(i)/egm(i+1)) endif enddo esp_ind(15) = esp_ind(14) ! temp.: assign TS3 the TS2 spec. ind. ! write gamma values used to file 28, named above ! if not species in question will get NAN for egms ! ! print *, ' esp_ind: ', (esp_ind(i), i=i1,i2) !====================================================================== !====================================================================== ! Frame Transformation Procedure !====================================================================== !---------------------------------------------------------------------- ! Determine Plasma Convection Velocity !---------------------------------------------------------------------- Vcon_x = 0.0 Vcon_y = 0.0 Vcon_z = 0.0 !---------------------------------------------------------------------- ! jvtr: 0 = Rigid corotation direction and magnitude[*cfact=0.0-1.0] !---------------------------------------------------------------------- if (jvtr .eq. 0) then Vcon_x = corot_sc(1)*crv*cfact - galv_sc(1) Vcon_y = corot_sc(2)*crv*cfact - galv_sc(2) Vcon_z = corot_sc(3)*crv*cfact - galv_sc(3) !---------------------------------------------------------------------- ! jvtr: 1 = Rigid corotation direction, user-defined magnitude vc_cus !---------------------------------------------------------------------- else if (jvtr .eq. 1) then Vcon_x = corot_sc(1)*vc_cus - galv_sc(1) Vcon_y = corot_sc(2)*vc_cus - galv_sc(2) Vcon_z = corot_sc(3)*vc_cus - galv_sc(3) !---------------------------------------------------------------------- ! jvtr: 2 = User-defined convection velocity (vt_mag, vt_clt, vt_azm) !---------------------------------------------------------------------- ! Construct local cartensian system centered on and moving with S/C: ! r_hat[1:3] = unit radial vector from Jupiter to Galileo ! c_hat[1:3] = unit vector in direction of corotation ! t_hat[1:3] = unit vector c_hat x r_hat ! Let c_hat axis be polar axis from which polar angle (colatitude) ! angle vt_clt [0,180] is defined. ! Let R_hat-t_hat plane be equatorial plane in which azimuth angle ! vt_azm [0,360) is defined, with vt_azm measured from r_hat ! CCW toward t_hat [note: wrt usual spherical polars, identify ! z -> c-hat, x -> r_hat, y -> t_hat] !---------------------------------------------------------------------- else if (jvtr .eq. 2) then do i=1,3 c_hat(i) = corot_sc(i) ! unit corotation vector r_hat(i) = -jup_sc_unit(i) ! unit radial vector J->G enddo ! construct unit vector t_hat (=c_hat x r_hat) normal to c_hat-r_hat plane t_hat(1) = c_hat(2)*r_hat(3) - c_hat(3)*r_hat(2) t_hat(2) = c_hat(3)*r_hat(1) - c_hat(1)*r_hat(3) t_hat(3) = c_hat(1)*r_hat(2) - c_hat(2)*r_hat(1) ! get components of convection velocity along local axes [c_hat,r_hat,t_hat] V_c_hat = vt_mag*Cos(vt_clt*dtr) V_r_hat = vt_mag*Sin(vt_clt*dtr)*Cos(vt_azm*dtr) V_t_hat = vt_mag*Sin(vt_clt*dtr)*Sin(vt_azm*dtr) ! get components of convection velocity along S/C axes [xsc,ysc,zsc] Vcx = V_c_hat*c_hat(1) + V_r_hat*r_hat(1) + V_t_hat*t_hat(1) Vcy = V_c_hat*c_hat(2) + V_r_hat*r_hat(2) + V_t_hat*t_hat(2) Vcz = V_c_hat*c_hat(3) + V_r_hat*r_hat(3) + V_t_hat*t_hat(3) ! subtract off S/C velocity Vcon_x = Vcx - galv_sc(1) Vcon_y = Vcy - galv_sc(2) Vcon_z = Vcz - galv_sc(3) endif Vcon = Sqrt(Vcon_x**2 + Vcon_y**2 + Vcon_z**2) Vcon_ux = Vcon_x/Vcon ! x-comp. of unit vector Vcon_u Vcon_uy = Vcon_y/Vcon ! y-comp. of unit vector Vcon_u Vcon_uz = Vcon_z/Vcon ! z-comp. of unit vector Vcon_u print *, 'Vcon_x,Vcon_y,Vcon_z,Vcon: ', * Vcon_x,Vcon_y,Vcon_z,Vcon ! ! Write a header for 28, gamma file if (kk.eq.0) write(28,*) jcht, jion, crv, vcon kk=1 write(28,*) n_dat,kdy,kdd,kdh,kdm,kds,esp_ind(jcht) !---------------------------------------------------------------------- ! Get begin sector ns_beg (0-15) from micro-sector nsbd (0-63) by ! doing integer divide by 4 (i.e., 4 micro-sectors per macro-sector) !---------------------------------------------------------------------- ns_beg = nsbd/4 ! start sector (0-15) print *, 'nsbd, ns_beg:', nsbd, ns_beg !---------------------------------------------------------------------- ! Now calculate the 16 consecutive sectors in current record !---------------------------------------------------------------------- do n=1,16 mtmp = ns_beg + (n-1) ms_rec(n) = Mod(mtmp,16) ! must be [1,16] print *, 'n, ms_rec(n):', n, ms_rec(n) enddo !---------------------------------------------------------------------- ! For current record, determine look directions of the 16 sector ! positions at current motor position. Look direction for current ! record, sector n [1-16], is unit vector ld = (xld,xld,xld), and ! vector velocity of entering ion is v = (vxld,vyld,vzld). !---------------------------------------------------------------------- ! Current motor positon [1-7] mp = mpcn print *, 'mp,tld(mp):', mp, tld(mp) ! !---------------------------------------------------------------------- ! Look directions 16 sectors of current record ! ms_rec(n:16) contains numbers of the 16 macro-sectors [0,15] !---------------------------------------------------------------------- do n=1,16 ns = ms_rec(n) + 1 ! index [1-16] of sector [0-15] !---------------------------------------------------------------------- ! EPD look direction unit vector: ! LEMMS in 180 deg end: multiply components by -1 !---------------------------------------------------------------------- xld = -xld_mn(mp,ns) ! x-comp. wrt SC sys. K[x,y,z] yld = -yld_mn(mp,ns) ! y-comp. wrt SC sys. K[x,y,z] zld = -zld_mn(mp,ns) ! z-comp. wrt SC sys. K[x,y,z] !---------------------------------------------------------------------- ! Unit vector of particle entry velocity in SC frame K[x,y,z]: ! t_v_mn = [t_vx_mn, t_vy_mn, t_vz_mn] !---------------------------------------------------------------------- t_vx_mn = -xld ! vx-comp. wrt SC sys. K[x,y,z] t_vy_mn = -yld ! vy-comp. wrt SC sys. K[x,y,z] t_vz_mn = -zld ! vz-comp. wrt SC sys. K[x,y,z] !---------------------------------------------------------------------- ! Cos(t_v_mn,Vcon_u): Cos[angle between (unit vector of particle ! entry velocity in SC frame K) and (unit vector of convection ! velocity Vcon_u)] !---------------------------------------------------------------------- Cos_tv_Vcon = t_vx_mn*Vcon_ux + * t_vy_mn*Vcon_uy + * t_vz_mn*Vcon_uz Sin2_tv_Vcon = 1.0 - Cos_tv_Vcon**2 !---------------------------------------------------------------------- ! Derived particle speed in SC frame K[x,y,z] at this direction (m,n) ! given the look-direction-independent constant speed vch defined ! in comoving frame K'[x',y',z']: !---------------------------------------------------------------------- v_mn = vch*((Vcon/vch)*Cos_tv_Vcon + * Sqrt(1.0 - ((Vcon/vch)**2)*Sin2_tv_Vcon)) c v_mn = vch + Vcon*Cos_tv_Vcon ! approx. to O(Vcon/vch) !---------------------------------------------------------------------- ! Derived components of unit vector to look-direction-independent ! constant speed vch defined in comoving frame K'[x',y',z']: ! t_vp_mn = [t_vxp_mn, t_vyp_mn, t_vzp_mn] ! NOTE: Needed to calc. pitch cosine in comoving frame K'[x',y',z'] !---------------------------------------------------------------------- t_vxp_mn = (t_vx_mn*v_mn - Vcon_x)/vch t_vyp_mn = (t_vy_mn*v_mn - Vcon_y)/vch t_vzp_mn = (t_vz_mn*v_mn - Vcon_z)/vch !---------------------------------------------------------------------- ! Dot product of ld=(xld,yld,zld) with B=(bxd,byd,bzd) gives the ! pitch cosine of the sector's look direction (sld). To get the pitch ! cosine of a particle that enters that sector, change the sign ! of the sector's look direction pitch cosine. !---------------------------------------------------------------------- Bx = bxd(n) By = byd(n) Bz = bzd(n) Bm = Sqrt(Bx**2 + By**2 + Bz**2) b_ux = Bx/Bm b_uy = By/Bm b_uz = Bz/Bm !---------------------------------------------------------------------- ! Pitch cosine and pitch angle of particle velocity ! in SC frame K[x,y,z]: !---------------------------------------------------------------------- pcos_v = t_vx_mn*b_ux + t_vy_mn*b_uy + t_vz_mn*b_uz pang_v = ACos(pcos_v)*rtd !---------------------------------------------------------------------- ! Pitch cosine and pitch angle of particle velocity ! in comoving frame K'[x',y',z']: !---------------------------------------------------------------------- pcos_vp = t_vxp_mn*b_ux + t_vyp_mn*b_uy + t_vzp_mn*b_uz pang_vp = ACos(pcos_vp)*rtd !................................................................... ! Note: Because the NAV file always lists the magnetic pitch-angles ! of the 16 macro-sectors in order, from sector 0 through 15, ! to compare my results pang to those on the NAV records use ! the index ns to increment pamag, i.e., as pamag(ns). !................................................................... c print *, ' n,ms_rec(n),ns,pld(ns),pang_v,pamag(ns):', c * n,ms_rec(n),ns,pld(ns),pang_v,pamag(ns) rate = r_lem(jcht,n) ! sector rate from file flxd = rr_lem(jcht,n) ! sector flux from file flux = rate/rtf(jcht) ! rederived flux !---------------------------------------------------------------------- ! Velocity components & K.E. of incident particle in SC frame K !---------------------------------------------------------------------- vx_mn = -t_vx_mn*vch ! x-comp. vy_mn = -t_vy_mn*vch ! y-comp. vz_mn = -t_vz_mn*vch ! z-comp. T_mn = vx_mn**2 + vy_mn**2 + vz_mn**2 v_mn = Sqrt(T_mn) !---------------------------------------------------------------------- ! Velocity components & K.E. of incident particle in comoving frame K' !---------------------------------------------------------------------- vxp_mn = vx_mn - Vcon_x ! x'-comp vyp_mn = vy_mn - Vcon_y ! y'-comp vzp_mn = vz_mn - Vcon_z ! z'-comp Tp_mn = vxp_mn**2 + vyp_mn**2 + vzp_mn**2 vp_mn = Sqrt(Tp_mn) !---------------------------------------------------------------------- ! Get energy spetral index for this channel. ! Try indices eps_ind(jcht) and esp_ind(jcht+1) first; ! use default d_esp_ind if these are unsatisfactory. ! Note: ==> For now, channels A1-A5 only, i1=2, i2=6 !---------------------------------------------------------------------- sp_ind = 0.0 ! for the compiler if ((jcht .ge. i1) .and. (jcht .lt. i2)) then if (Abs(esp_ind(jcht)) .lt. 6.0) then sp_ind = esp_ind(jcht) c elseif (Abs(esp_ind(jcht+1)) .lt. 6.0) then c sp_ind = esp_ind(jcht+1) else sp_ind = d_esp_ind endif endif cgf = (Tp_mn/T_mn)**(sp_ind + 1) ! C-G factor rate_p = rate*cgf flux_p = flux*cgf ! print *, ! * 'n,ms_rec(n),pamag(ns),pang_v,v_mn,vp_mn,', ! * 'sp_ind,rate,rate_p:', ! * n,ms_rec(n),pamag(ns),pang_v,v_mn,vp_mn, ! * sp_ind,rate,rate_p c This is what Rob had: c if (ndo .gt. 0) then c write (8,*) intim_dat, n, ms_rec(n), pamag(ns), pang_v, c * v_mn, vp_mn, rate, rate_p c endif if (ndo .gt. 0) then write (8,*) intim_dat(2),intim_dat(3),intim_dat(4), * intim_dat(5),pang_v, pang_vp, rate, rate_p endif ! ! write the transformed intensity to file 29, named above ! If the gamma is 1-3, there is little reason ! to preserve pitch angle and do the spectral slopes separately ! if (flux_p.gt.0.) then write(29,*) n_dat,kdy,kdd,kdh,kdm,kds,flux_p npts2=npts2+1 ! Average the log of the new transformed intensity fpa=alog10(flux_p) avg2=avg2+fpa sumsq2=sumsq2+(fpa*fpa) endif ! enddo !====================================================================== goto 1000 2000 write (6, 2010) 2010 format (' EOF on Nav data file.') goto 5000 3000 write (6, 3010) 3010 format (' EOF on EPD data file.') goto 5000 4000 write (6, 4010) 4010 format (' User-defined end time encountered.') 5000 close (unit = 2) close (unit = 4) if (npts1.gt.0) then avg1=avg1/float(npts1) fact3=sumsq1-(npts1*avg1*avg1) sdev1=sqrt(fact3/float(npts1-1)) write(6,*) 0.,npts1,avg1,sdev1,sdev1/avg1 endif if (npts2.gt.0) then avg2=avg2/float(npts2) fact4=sumsq2-(npts2*avg2*avg2) sdev2=sqrt(fact4/float(npts2-1)) write(6,*) 1.,npts2,avg2,sdev2,sdev2/avg2 endif write(6,*) 'jcht, jion, crv, vcon' write(6,*) jcht, jion, crv, vcon stop end !==================================================================== ! subroutine c_seqt: !-------------------- ! Sequential time in hours or days is elapsed time from ! 1960/001/00:00:00 ! ! k=-1: Convert decimal sequential HOUR (seqt) to integer: ! sequential HOUR (iseqt), year (iyr), doy (idoy), ! hour (ihr), minute (imn), and second (isc). ! ! k=-2: Convert decimal sequential DAY (seqt) to integer: ! sequential DAY (iseqt), year (iyr), doy (idoy), ! hour (ihr), minute (imn), and second (isc). ! ! k= 1: Convert integer year (iyr), doy (idoy), hour (ihr), ! minute (imn), and second (isc) to integer sequential ! HOUR (iseqt) and decimal sequential HOUR (seqt). ! ! k= 2: Convert integer year (iyr), doy (idoy), hour (ihr), ! minute (imn), and second (isc) to integer sequential ! DAY (iseqt) and decimal sequential DAY (seqt). ! ! notes: ! (1) uses 4-digit year, iyr (1960 .ge. iyr .le. 2022) ! (2) idoy = 1 to 365 or 366 ! (3) ihr = 0 to 23 ! (4) imn = 0 to 59 ! (5) isc = 0 to 59 ! (6) seqt must be real*8 in calling program ! ! created: 04-jun-97 RBDecker ! changed: 28-jun-00 RBDecker !==================================================================== subroutine c_seqt (k,seqt,iseqt,iyr,idoy,ihr,imn,isc) integer*4 jsd(64) real*8 seqt, sqd, f, fhr, fmn, fsc, eps ! doy 0, 1960 through 2022 (add 1 to get doy 1) data jsd/ 0, 366, 731, 1096, ! 000/1960 to 000/1963 * 1461, 1827, 2192, 2557, ! 000/1964 to 000/1967 * 2922, 3288, 3653, 4018, ! 000/1967 to 000/1971 * 4383, 4749, 5114, 5479, ! 000/1972 to 000/1975 * 5844, 6210, 6575, 6940, ! 000/1976 to 000/1979 * 7305, 7671, 8036, 8401, ! 000/1980 to 000/1983 * 8766, 9132, 9497, 9862, ! 000/1984 to 000/1987 * 10227, 10593, 10958, 11323, ! 000/1988 to 000/1991 * 11688, 12054, 12419, 12784, ! 000/1992 to 000/1995 * 13149, 13515, 13880, 14245, ! 000/1996 to 000/1999 * 14610, 14976, 15341, 15706, ! 000/2000 to 000/2003 * 16071, 16437, 16802, 17167, ! 000/2004 to 000/2007 * 17532, 17898, 18263, 18628, ! 000/2008 to 000/2011 * 18993, 19359, 19724, 20089, ! 000/2012 to 000/2015 * 20454, 20820, 21185, 21550, ! 000/2016 to 000/2019 * 21915, 22281, 22646, 23011/ ! 000/2020 to 000/2023 data eps/ 1.0d-04/ ! Compiler wants next two statements sqd = 0.0d0 isqd = 0 !-------------------------------------------------------------------- ! k=-1 or -2: seqt -> [iseqt=Int(seqt), iyr, idoy, ihr, imn, isc] !-------------------------------------------------------------------- if (k .lt. 0) then iseqt = seqt ! integer sequential time if (k .eq. -1) then sqd = seqt/24.d0 ! decimal sequential day, k=-1 isqd = sqd ! integer sequential day, k=-1 else if (k .eq. -2) then sqd = seqt ! decimal sequential day, k=-2 isqd = sqd ! integer sequential day, k=-2 end if nn = 1 do n=1,63 if ((isqd .gt. jsd(n)) .and. (isqd .le. jsd(n+1))) then nn = n go to 10 end if end do 10 iyr = nn + 1959 idoy = isqd - jsd(nn) fhr = (sqd - Float(isqd))*24.d0 ihr = fhr + eps ! integer hour of doy fmn = (fhr - Float(ihr))*60.d0 imn = fmn + eps ! integer minute of hour fsc = (fmn - Float(imn))*60.d0 isc = fsc + eps ! integer second of minute ! test for isc=60, imn=60, ihr=24 if (isc .lt. 60) go to 20 isc = 0 imn = imn + 1 20 if (imn .lt. 60) go to 25 imn = 0 ihr = ihr + 1 25 if (ihr .lt. 24) return ihr = 0 isqd = isqd + 1 end if !-------------------------------------------------------------------- ! k=1 or 2: [iyr, idoy, ihr, imn, isc] -> [seqt, iseqt=Int(seqt)] !-------------------------------------------------------------------- if (k .gt. 0) then jyr = iyr if (jyr .lt. 100) jyr = jyr + 1900 n = jyr - 1959 ! no. years since 1960 isqd = jsd(n) + idoy ! integer sequential doy f = Float(ihr) * + Float(imn)/60.d0 * + Float(isc)/3600.d0 ! decimal hour of day if (k .eq. 1) then khr = isqd*24 iseqt = khr + ihr ! integer sequential hour seqt = Float(khr) + f ! decimal sequential hour else if (k .eq. 2) then iseqt = isqd ! integer sequential day seqt = Float(iseqt) + f/24.d0 ! decimal sequential day end if end if return end