PROGRAM GETNDCS C This program will retrieve Dave Evan's Hemispheric Power Input C values from I.S. database file EHP781102A, the 3 hour values C of Kp and ap and daily values of AP, F10.7 flux and sunspot C number from I.S. database file GPI600101A, and the hourly C Interplanetary Magnetic Field Solar Wind data from NCAR I.S. C database IMF631127A. In addition, it will write out those data for C those times that can be used (i.e. when data were recorded) for C Hemispheric Power, Crosstail Potential, and the IMF. The format C of the output will allow lexical reads to be made--the first C number is the UT date, the second number is the UT hour, the C third the UT minute, and the last either POWER, CTPOTEN, or IMF. C It will also plot these data versus time. C C This version has been installed on a Cray, and uses the version C of rdblk6 which calls rdtape; it does not require the longer C work array; see MXLWRK below. C C File linkages: C Done as part of the jcl of the batch job: C stdin = fort.5 - job control = JCTLUN C stdout = fort.6 - print output and diagnostics = MSGUN C fort.61 = fort.61 - lexical write output file (not C explicitly opened). C gmeta = fort.2 - graphics metacode plot output C Input data files (NCAR I.S. binary format with Cray block envelope) C NOTE: The fortran unit number file name linkage must be done prior C to execution because the OPNBLK entry in the Cray version of C RDBLK6 does not attempt to open the named files. C EHP781102A ( = fort.11 ) = IUNEHP C GPI600101A ( = fort.12 ) = IUNGPI C IMF631127A ( = fort.13 ) = IUNIMF C AEI780101A ( = fort.14 ) = IUNAE C DST570101A ( = fort.15 ) = IUNDST C Additional source files referenced but not included here (all C files should be in directory ~/pgms/is/src for this version): C rdblk6.f cbfrd.f timdif.f hours.f newtim.f cvti2r.f getehp.f C gethpd.f getgpd.f getimf.f getae.f getdst.f gtmtic.f gaxtic.f C gnumb.f pltyvt.f parsb.f C Define unit number to read job control parameters: PARAMETER (JCTLUN=5) C Declarations for RDBLK6 which are the same for all get__ subs PARAMETER (NEOFS=9999 , MSGUN=6 , IOCKSM=0 , + LMWD=64 , MXVALS=10000 , + MXWDS=MXVALS*16/LMWD , MXLWRK= 17 ) CUNIX+ MXWDS=MXVALS*16/LMWD , MXLWRK= 17+20 + 2*512*64/LMWD) CUNIX version w/o CBFRD does not need extra work array space C Declarations for getehp: PARAMETER (MXNPWR = 2000) C MXNPWR = maximum number of values to be returned by GETEHP DIMENSION ITMEHP(4,MXNPWR) , HPEHP(MXNPWR) , NDXEHP(MXNPWR) , + IQLEHP(MXNPWR) C Declarations for rdblk6 specifically for getehp CHARACTER*80 NAMEHP PARAMETER (IUNEHP=11) DIMENSION LREHP(MXVALS) , IUPEHP(MXVALS) , IBKEHP(MXWDS) , + IWKEHP(MXLWRK) C Declarations for getgpi: C PARAMETER (MXNGPI = 40 , MXNG3H = MXNGPI*8) PARAMETER (MXNGPI = 2192 , MXNG3H = MXNGPI*8) C MXNGPI = maximum number of values to be retrieved by GETGPI DIMENSION ITM24G(4,MXNGPI) , APDG(MXNGPI) , F107G(MXNGPI) , + F107AG(MXNGPI) , SSNG(MXNGPI) , ISTGPI(MXNGPI) , + ITM03G(4,MXNG3H) , XKPG(MXNG3H) , APG(MXNG3H) , + IKP(8) , IAP(8) , I3HT(8) , RF107(MXNGPI+80) DATA I3HT/0130 , 0430 , 0730 , 1030 , 1330 , 1630 , 1930 , 2230/ C Declarations for rdblk6 specifically for getgpi CHARACTER*80 NAMGPI,toplab,xlab dimension xdays(366),xdays8(366*8) PARAMETER (IUNGPI=12) DIMENSION LRGPI(MXVALS) , IUPGPI(MXVALS) , IBKGPI(MXWDS) , + IWKGPI(MXLWRK) C Declarations for getimf: PARAMETER (MXNBY = 750) C MXNBY = maximum number of hourly IMF's to be retreived C by GETIMF DIMENSION ITMIMF(4,MXNBY) , + BX(MXNBY) , BY(MXNBY) , BZ(MXNBY) , + BYGSE(MXNBY) , BZGSE(MXNBY) , IQLIMF(MXNBY) , + SWD(MXNBY) , SWV(MXNBY) , ISTIMF(MXNBY) C Declarations for rdblk6 specifically for getimf CHARACTER*80 NAMIMF PARAMETER (IUNIMF=13) DIMENSION LRIMF(MXVALS) , IUPIMF(MXVALS) , IBKIMF(MXWDS) , + IWKIMF(MXLWRK) C Declarations for getae: PARAMETER (MXNAE = 750) DIMENSION ITMAE(4,MXNAE) , AE(MXNAE), ISTAE(MXNAE) C Declarations for rdblk6 specifically for getae CHARACTER*80 NAMAE PARAMETER (IUNAE=14) DIMENSION LRAE(MXVALS) , IUPAE(MXVALS) , IBKAE(MXWDS) , + IWKAE(MXLWRK) C Declarations for getdst: PARAMETER (MXNDST = 750) DIMENSION ITMDST(4,MXNDST) , DST(MXNDST) , ISTDST(MXNDST) C Declarations for rdblk6 specifically for getdst CHARACTER*80 NAMDST PARAMETER (IUNDST=15) DIMENSION LRDST(MXVALS) , IUPDST(MXVALS) , IBKDST(MXWDS) , + IWKDST(MXLWRK) C Miscellaneous declarations DIMENSION IBEGT(4) , IENDT(4) , ITMP1(4) , ITMP2(4) , + IBTHR(4) , IETHR(4) , IBTG(4) , IETG(4) DOUBLE PRECISION ICSDIF , IDNHCS , I1HR , I1DA , M40DA PARAMETER (MISS=-32767 , RMISS=-32767. , I1HR=360000D0 , + TMISS=-32766.9 , NOPRT=99999 , I1DA=24D0*I1HR ) C Declarations for parameters found from AE DIMENSION CPFAE(MXNAE,2) , NOACP(2) , TIMAE(MXNAE), AENM(MXNAE), | HPFAE(MXNAE) C Declarations for parameters found from Hp (P) DIMENSION TIMHP(MXNPWR) , CPFHP(MXNPWR,2), NOHCP(2) C Declarations for parameters found from Kp DIMENSION HPFKP(MXNG3H) , CPFKP(MXNG3H,3), NOKCP(3), | TIMKP(MXNG3H) C Declarations for parameters found from IMF (and some others) DIMENSION TMNOG(MXNBY) , SWNOG(MXNBY) , SDNOG(MXNBY) , | POT(MXNBY,3) , POTF(MXNBY,3) , NOICP(3), NOICPF(3) do i=1,366 xdays(i) = float(i) do ii=1,8 xdays8((i-1)*8+ii) = float(i)+float(ii-1)/8. enddo enddo C Get and print desired options and local input file names READ (JCTLUN,*) IGETHP,IGETGPI,IGETIMF,IGETAE,IGETDST,ICALCT C IGETHP = 1,0 if do,not get HP data C IGETGPI = 1,0 if do,not get GPI data (SS#,F107,Kp,Ap,..) C IGETIMF = 1,0 if do,not get IMF data (IMF,VSW,DSW,CPOTEN) C IGETAE = 1,0 if do,not get AE hourly data C IGETDST = 1,0 if do,not get Dst C ICALCT = 1,0 if do,not calculate crostail potential and lexical write WRITE (MSGUN,'(''1 PROGRAM GETNDCS - OPTIONS SELECTED:'',/, + '' IGETHP IGETGPI IGETIMF IGETAE IGETDST ICALCT'',/,6I8)') + IGETHP,IGETGPI,IGETIMF, IGETAE,IGETDST, ICALCT READ (JCTLUN,*) NAMEHP,NAMGPI,NAMIMF,NAMAE,NAMDST IF (IGETHP .EQ. 1) THEN CALL PARSB(NAMEHP,IBH,IEH) WRITE(MSGUN,'('' Hemispheric power data are read from '',A) + ') NAMEHP(IBH:IEH) ENDIF IF (IGETGPI .EQ. 1) THEN CALL PARSB(NAMGPI,IBG,IEG) WRITE(MSGUN,'('' GPI (SS#,F107,Kp,Ap,..) are read from '',A) + ') NAMGPI(IBG:IEG) ENDIF IF (IGETIMF .EQ. 1) THEN CALL PARSB(NAMIMF,IBI,IEI) WRITE(MSGUN,'('' Hourly IMF data are read from '',A) + ') NAMIMF(IBI:IEI) ENDIF IF (IGETAE .EQ. 1) THEN CALL PARSB(NAMAE,IBA,IEA) WRITE(MSGUN,'('' Hourly ae indices are read from '',A) + ') NAMAE(IBA:IEA) ENDIF IF (IGETDST .EQ. 1) THEN CALL PARSB(NAMDST,IBD,IED) WRITE(MSGUN,'('' Hourly DST data are read from '',A) + ') NAMDST(IBD:IED) ENDIF C Open GKS for plotting CALL OPNGKS C Loop on each time interval requested 2 CONTINUE C Get and print desired time interval or terminate IF (ICALCT .EQ. 1) THEN READ (JCTLUN,*,END=1000) (IBEGT(I),I=1,4),(IENDT(J),J=1,4) + , LEXBDAY WRITE (MSGUN,'( + '' Begin time (year moday hrmin centisec) ='',4I5,/, + '' End time ='',4I5,/, + '' and lexical day (LEXBDAY) ='',I3, + '' matches begin day'',2I5)') + (IBEGT(I),I=1,4),(IENDT(J),J=1,4),LEXBDAY,(IBEGT(K),K=1,2) ELSE READ (JCTLUN,*,END=1000) (IBEGT(I),I=1,4),(IENDT(J),J=1,4) WRITE (MSGUN,'( + '' Begin time (year moday hrmin centisec) ='',4I5,/, + '' End time ='',4I5)') + (IBEGT(I),I=1,4),(IENDT(J),J=1,4) ENDIF C Check input time NOK = 0 IF(IBEGT(1) .LT. 1940 .OR. IBEGT(1) .GT. 2100)NOK = 1 IF(IBEGT(2) .LT. 0101 .OR. IBEGT(2) .GT. 1231)NOK = 1 IF(IBEGT(3) .LT. 0 .OR. IBEGT(3) .GT. 2400)NOK = 1 IF(IBEGT(4) .LT. 0 .OR. IBEGT(4) .GT. 5999)NOK = 1 IF(IENDT(1) .LT. 1940 .OR. IENDT(1) .GT. 2100)NOK = 1 IF(IENDT(2) .LT. 0101 .OR. IENDT(2) .GT. 1231)NOK = 1 IF(IENDT(3) .LT. 0 .OR. IENDT(3) .GT. 2400)NOK = 1 IF(IENDT(4) .LT. 0 .OR. IENDT(4) .GT. 5999)NOK = 1 IF(NOK .GT. 0)GO TO 9000 C Determine NDAS the number of days indicies requested C and NHRTOT the number of hours difference CALL TIMDIF(IBEGT,IENDT,ICSDIF) IF (ICSDIF .LT. 0D0) GO TO 9050 NHRTOT = IDINT( (ICSDIF-1D0) / I1HR ) + 1 DO 3 I=1,2 IBTG(I) = IBEGT(I) 3 IETG(I) = IENDT(I) IBTG(3) = 0 IBTG(4) = 0 IETG(3) = 2400 IETG(4) = 0 CALL TIMDIF(IBTG,IETG,ICSDIF) NDAS = IDINT( (ICSDIF-1D0) / I1DA ) + 1 C Get Ae indicies IF(IGETAE .EQ. 0) GO TO 100 C Open NCAR I.S. format (binary) for read CALL OPNBLK(IUNAE,NAMAE(IBA:IEA),NEOFS,IOCKSM,MSGUN,LMWD,MXVALS, + MXWDS,IWKAE,MXLWRK) WRITE(6,'('' --------- AE INDEX ----------'')') IF(NHRTOT .GT. MXNAE)GO TO 9250 NOKAE = 0 ITMP1(1) = IBEGT(1) ITMP1(2) = IBEGT(2) ITMP1(3) = 100*(IBEGT(3)/100) + 30 ITMP1(4) = 0 ITMP2(4) = 0 DO 10 I=1,NHRTOT ITMAE(1,I) = ITMP1(1) ITMAE(2,I) = ITMP1(2) ITMAE(3,I) = ITMP1(3) ITMAE(4,I) = 0 ISTAT = NOPRT CALL GETAE(MSGUN,ITMP1,324,AE(I) , ISTAT , + IWKAE,IBKAE,IUPAE,LRAE) ISTAE(I) = ISTAT IF (AE(I) .GT. TMISS) NOKAE = NOKAE + 1 C Increment date ITMP2(1) = ITMP1(1) ITMP2(2) = ITMP1(2) ITMP2(3) = ITMP1(3) 10 CALL NEWTIM(ITMP2,ITMP1,I1HR,MSGUN) C Print Ae values except when entire date is missing: WRITE(6,"(' Time (mid-point) Ae(hourly)' + ,/, ' year mmdd hrmn nT')") LY = 0 LMD = 0 DO 20 I=1,NHRTOT IF(ISTAE(I) .EQ. 6)THEN WRITE(6,9600)(ITMAE(J,I),J=1,3) GO TO 20 ENDIF IF(ISTAE(I) .NE. 0 .AND. + ITMAE(1,I) .EQ. LY .AND. ITMAE(2,I) .EQ. LMD)GO TO 20 LY = ITMAE(1,I) LMD = ITMAE(2,I) IF(ISTAE(I) .EQ. 2)WRITE(6,9300)(ITMAE(J,I),J=1,2) IF(ISTAE(I) .EQ. 3)WRITE(6,9400)(ITMAE(J,I),J=1,2) IF(ISTAE(I) .EQ. 0)WRITE(6,'(3I6,F9.0)')(ITMAE(J,I),J=1,3),AE(I) 20 CONTINUE WRITE(6,'('' '',/,'' '')') 100 CONTINUE C Get Hemispheric Power Data NVEHP = 0 IF(IGETHP .EQ. 0) GO TO 200 C Open NCAR I.S. format (binary) for read CALL OPNBLK(IUNEHP,NAMEHP(IBH:IEH),NEOFS,IOCKSM,MSGUN,LMWD,MXVALS, + MXWDS,IWKEHP,MXLWRK) WRITE(6,'('' ------ ESTIMATED HEMISPHERIC POWER INPUT -------'')') CALL GETEHP(MSGUN,IBEGT,IENDT,MXNPWR,NVEHP, + ITMEHP,HPEHP,NDXEHP,IQLEHP,ISTEHP, + IWKEHP,IBKEHP,IUPEHP,LREHP) WRITE(6,'('' Requested Time Interval # Values''/, + '' -------Begin------- --------End-------- Found '',/, + 8I5,I8)') (IBEGT(I),I=1,4),(IENDT(J),J=1,4),NVEHP IF(NVEHP .GT. 0) + WRITE (6,'('' UT Date/Time Power(GW) PwrIndex PwrQualifier'', + /,(3I5,F8.1,I7,I13))') ((ITMEHP(I,J),I=1,3), + HPEHP(J),NDXEHP(J),IQLEHP(J),J=1,NVEHP) WRITE(6,'('' '',/,'' '')') 200 CONTINUE C Get GPI (Kp, ap, AP, F10.7, Sunspot Number) IF (IGETGPI .EQ. 0) GO TO 300 C Open NCAR I.S. format (binary) for read CALL OPNBLK(IUNGPI,NAMGPI(IBG:IEG),NEOFS,IOCKSM,MSGUN,LMWD,MXVALS, + MXWDS,IWKGPI,MXLWRK) WRITE(6,'('' ------------- GEOPHYSICAL INDICIES ------------'')') IF(NDAS .GT. MXNGPI)GO TO 9100 C Since 81 day average F10.7 is desired, compute it here using C GETGPD to fetch a days worth of values for each call rather C than using GETGPI which would force a rewind for each new C 81 day average. NTOTD = NDAS + 80 IBSAV = 41 IESAV = NDAS + 40 M40DA = -I1DA* 40D0 CALL NEWTIM(IBEGT,ITMP1,M40DA,MSGUN) ITMP1(3) = 0 ITMP1(4) = 0 ITMP2(3) = 0 ITMP2(4) = 0 ISAVD = 0 NOKAPD = 0 NOKF10 = 0 NOKF1A = 0 NOKSSN = 0 NOKKP = 0 NOKAP3 = 0 DO 220 I=1,NTOTD ISTAT = NOPRT CALL GETGPD(MSGUN,ITMP1,IKP,IAP,IAPD,IF107,IFQFR,ISSN, + ISTAT,IWKGPI,IBKGPI,IUPGPI,LRGPI) C Always save F10.7 values for 81 day averages RF107(I) = RMISS IF(IF107 .NE. MISS)RF107(I) = FLOAT(IF107) / 10. IF(I .GE. IBSAV .AND. I .LE. IESAV)THEN C Save daily average values; hr-min is fixed at 1200, C representing the mean time of the daily average values. ISAVD = ISAVD + 1 ISTGPI(ISAVD) = ISTAT ITM24G(1,ISAVD) = ITMP1(1) ITM24G(2,ISAVD) = ITMP1(2) ITM24G(3,ISAVD) = 1200 ITM24G(4,ISAVD) = 0000 APDG (ISAVD) = FLOAT(IAPD) F107G (ISAVD) = RF107(I) SSNG (ISAVD) = FLOAT(ISSN) IF(IAPD .NE. MISS) NOKAPD = NOKAPD + 1 IF(IF107 .NE. MISS) NOKF10 = NOKF10 + 1 IF(ISSN .NE. MISS) NOKSSN = NOKSSN + 1 C Get 3hr average indicies; hr-min is fixed to 0130,0430...2230 C representing the mean time of day for each index value. IB3 = (ISAVD-1)*8 + 1 IE3 = IB3 + 7 J = 0 DO 210 I3H=IB3,IE3 J = J + 1 ITM03G(1,I3H) = ITMP1(1) ITM03G(2,I3H) = ITMP1(2) ITM03G(3,I3H) = I3HT(J) ITM03G(4,I3H) = 0000 APG (I3H) = FLOAT(IAP(J)) XKPG(I3H) = RMISS IF(IKP(J) .NE. MISS)THEN XKPG(I3H) = FLOAT(IKP(J)) / 10. NOKKP = NOKKP + 1 ENDIF IF(IAP(J) .NE. MISS) NOKAP3 = NOKAP3 + 1 210 CONTINUE ENDIF C Increment date ITMP2(1) = ITMP1(1) ITMP2(2) = ITMP1(2) 220 CALL NEWTIM(ITMP2,ITMP1,I1DA,MSGUN) C Compute F10.7 solar flux 81 day average from saved values DO 240 I=1,NDAS SUM = 0. CNT = 0. DO 230 J=I,I+80 IF(RF107(J) .LT. TMISS)GO TO 230 SUM = SUM + RF107(J) CNT = CNT + 1. 230 CONTINUE F107AG(I) = RMISS IF(CNT .GT. 40.)THEN F107AG(I) = SUM / CNT NOKF1A = NOKF1A + 1 ENDIF 240 CONTINUE C Establish the number of 3hr GPI assigned N3HRG = NDAS*8 C Print GPI values except when entire date is missing: c WRITE(6,'('' Time (mid-point) Daily 81day Daily Daily'',/, c + '' year mmdd hrmn SS# F10.7 F10.7 AP'')') c DO 250 I=1,NDAS c IF(ISTGPI(I) .EQ. 2)WRITE(6,9300)(ITM24G(J,I),J=1,2) c IF(ISTGPI(I) .EQ. 3)WRITE(6,9400)(ITM24G(J,I),J=1,2) c IF(ISTGPI(I) .EQ. 6)WRITE(6,9600)(ITM24G(J,I),J=1,2) c IF(ISTGPI(I) .EQ. 0)WRITE(6,'(3I6,F7.0,2F8.1,F7.0)') c + (ITM24G(J,I),J=1,3),SSNG(I),F107AG(I),F107G(I),APDG(I) c 250 CONTINUE c WRITE(6,'('' '',/, c + '' Time (mid-point) 3hr 3hr'',/, c + '' year mmdd hrmn Kp ap'')') c DO 260 I=1,NDAS c I1 = (I-1)*8 + 1 c IF(ISTGPI(I) .EQ. 2)WRITE(6,9300)(ITM03G(J,I1),J=1,2) c IF(ISTGPI(I) .EQ. 3)WRITE(6,9400)(ITM03G(J,I1),J=1,2) c IF(ISTGPI(I) .EQ. 6)WRITE(6,9600)(ITM03G(J,I1),J=1,2) c IF(ISTGPI(I) .NE. 0)GO TO 260 c WRITE(6,'(3I6,F9.1,F7.0)')(ITM03G(J,I1),J=1,3),XKPG(I1),APG(I1) c I2 = I1+1 c I8 = I1+7 c WRITE(6,'(12X,I6,F9.1,F7.0)')(ITM03G(3,J),XKPG(J),APG(J),J=I2,I8) c 260 CONTINUE c WRITE(6,'('' '',/,'' '')') write(6,"('IYD F107D F107A KP')") ludat = 90 do i=1,ndas imo = itm24g(2,i)/100 ida = itm24g(2,i)-imo*100 iyr = itm24g(1,i)-1900 iyd = idate2iyd(imo,ida,iyr) I1 = (I-1)*8 + 1 I2 = I1+1 I8 = I1+7 write(ludat,"(i5,2f8.1,8f5.1)") iyd,f107g(i),f107ag(i), + (xkpg(j),j=(i-1)*8+1,(i-1)*8+8) enddo 300 CONTINUE C Get IMF and solar wind data IF (IGETIMF .EQ. 0) GO TO 400 C Open NCAR I.S. format (binary) for read CALL OPNBLK(IUNIMF,NAMIMF(IBI:IEI),NEOFS,IOCKSM,MSGUN,LMWD,MXVALS, + MXWDS,IWKIMF,MXLWRK) WRITE(6,'('' ----------------------------------- IMF AND SOLAR'', + '' WIND ------------------------------------'')') IF(NHRTOT .GT. MXNBY)GO TO 9200 NOKSD = 0 NOKSV = 0 NOKBX = 0 NOKBY = 0 NOKBZ = 0 SDMIN = 0. SDMAX = 0. SVMIN = 1.E6 SVMAX = 0. BMIN = 1.E6 BMAX =-1.E6 ITMP1(1) = IBEGT(1) ITMP1(2) = IBEGT(2) ITMP1(3) = 100*(IBEGT(3)/100) + 30 ITMP1(4) = 0 ITMP2(4) = 0 DO 310 I=1,NHRTOT ITMIMF(1,I) = ITMP1(1) ITMIMF(2,I) = ITMP1(2) ITMIMF(3,I) = ITMP1(3) ITMIMF(4,I) = 0 ISTAT = NOPRT CALL GETIMF(MSGUN,ITMP1, + BX(I),BY(I),BZ(I),BYGSE(I),BZGSE(I),SWDEN,SWSPD, + IQLIMF(I),ISTAT,IWKIMF,IBKIMF,IUPIMF,LRIMF) ISTIMF(I)= ISTAT SWD(I) = SWDEN/1000000. SWV(I) = SWSPD/1000. IF(SWDEN .LT. TMISS)SWD(I) = RMISS IF(SWSPD .LT. TMISS)SWV(I) = RMISS C Determine plot stuff: # ok values and domain IF(SWD(I) .GT. TMISS)THEN NOKSD = NOKSD + 1 SDMAX = AMAX1(SWD(I),SDMAX) ENDIF IF(SWV(I) .GT. TMISS)THEN NOKSV = NOKSV + 1 SVMIN = AMIN1(SWV(I),SVMIN) SVMAX = AMAX1(SWV(I),SVMAX) ENDIF IF(BX(I) .GT. TMISS)THEN NOKBX = NOKBX + 1 IF(BX(I) .LT. BMIN) BMIN = BX(I) IF(BX(I) .GT. BMAX) BMAX = BX(I) ENDIF IF(BY(I) .GT. TMISS)THEN NOKBY = NOKBY + 1 IF(BY(I) .LT. BMIN) BMIN = BY(I) IF(BY(I) .GT. BMAX) BMAX = BY(I) ENDIF IF(BZ(I) .GT. TMISS)THEN NOKBZ = NOKBZ + 1 IF(BZ(I) .LT. BMIN) BMIN = BZ(I) IF(BZ(I) .GT. BMAX) BMAX = BZ(I) ENDIF C Increment date ITMP2(1) = ITMP1(1) ITMP2(2) = ITMP1(2) ITMP2(3) = ITMP1(3) 310 CALL NEWTIM(ITMP2,ITMP1,I1HR,MSGUN) C Print IMF values except when entire date is missing: WRITE(6,'( +'' Mean ------- GSM [nT] ------- --- GSE [nT] --- '', + '' Solar Wind'',/, +'' ---- Time ---- Bx By Bz By Bz '', + '' [#/cm3] [km/s] Qualifier'')') LY = 0 LMD = 0 DO 320 I=1,NHRTOT IF(ISTIMF(I) .NE. 0 .AND. + ITMIMF(1,I) .EQ. LY .AND. ITMIMF(2,I) .EQ. LMD)GO TO 320 LY = ITMIMF(1,I) LMD = ITMIMF(2,I) IF(ISTIMF(I) .EQ. 2)WRITE(6,9300)(ITMIMF(J,I),J=1,2) IF(ISTIMF(I) .EQ. 3)WRITE(6,9400)(ITMIMF(J,I),J=1,2) IF(ISTIMF(I) .EQ. 4)WRITE(6,9500)(ITMIMF(J,I),J=1,2) IF(ISTIMF(I) .EQ. 0)WRITE(6,'(3I5,5F9.2,F12.1,F9.1,I7)') + (ITMIMF(J,I),J=1,3) , BX(I),BY(I),BZ(I), + BYGSE(I),BZGSE(I),SWD(I),SWV(I),IQLIMF(I) 320 CONTINUE WRITE(6,'('' '',/,'' '')') 400 CONTINUE C Get DST indicies IF(IGETDST .EQ. 0)GO TO 500 C Open NCAR I.S. format (binary) for read CALL OPNBLK(IUNDST,NAMDST(IBD:IED),NEOFS,IOCKSM,MSGUN,LMWD,MXVALS, + MXWDS,IWKDST,MXLWRK) WRITE(6,'('' --------- DST INDEX ---------'')') IF(NHRTOT .GT. MXNDST)GO TO 9260 NOKDST = 0 ITMP1(1) = IBEGT(1) ITMP1(2) = IBEGT(2) ITMP1(3) = 100*(IBEGT(3)/100) + 30 ITMP1(4) = 0 ITMP2(4) = 0 DO 410 I=1,NHRTOT ITMDST(1,I) = ITMP1(1) ITMDST(2,I) = ITMP1(2) ITMDST(3,I) = ITMP1(3) ITMDST(4,I) = 0 ISTAT = NOPRT CALL GETDST(MSGUN,ITMP1,DST(I), ISTAT , + IWKDST,IBKDST,IUPDST,LRDST) ISTDST(I) = ISTAT IF(DST(I) .GT. TMISS) NOKDST = NOKDST + 1 C Increment date ITMP2(1) = ITMP1(1) ITMP2(2) = ITMP1(2) ITMP2(3) = ITMP1(3) 410 CALL NEWTIM(ITMP2,ITMP1,I1HR,MSGUN) C Print DST values except when entire date is missing: WRITE(6,'('' Time (mid-point) DST(hourly)'',/, + '' year mmdd hrmn nT '')') LY = 0 LMD = 0 DO 420 I=1,NHRTOT IF(ISTDST(I) .EQ. 6)THEN WRITE(6,9600)(ITMDST(J,I),J=1,3) GO TO 420 ENDIF IF(ISTDST(I) .NE. 0 .AND. + ITMDST(1,I) .EQ. LY .AND. ITMDST(2,I) .EQ. LMD)GO TO 420 LY = ITMDST(1,I) LMD = ITMDST(2,I) IF(ISTDST(I) .EQ. 2)WRITE(6,9300)(ITMDST(J,I),J=1,2) IF(ISTDST(I) .EQ. 3)WRITE(6,9400)(ITMDST(J,I),J=1,2) IF(ISTDST(I) .EQ. 0)WRITE(6,'(3I6,F9.0)') + (ITMDST(J,I),J=1,3),DST(I) 420 CONTINUE WRITE(6,'('' '',/,'' '')') 500 CONTINUE C Done reading indicies C If requested calculate other parameters from the indicies IF (ICALCT .EQ. 0) GO TO 599 IF (IGETAE .EQ. 1) THEN NOKAE = 0 DO 510 I=1,NHRTOT HPFAE(I) = RMISS CPFAE(I,1) = RMISS CPFAE(I,2) = RMISS IF(AE(I) .GT. TMISS) THEN C Calculate hemispheric power from AE using Maeda's formula for 1 hr AE HPFAE(I) = 4.28 + 0.0624 * AE(I) C Calculate CP from P using CP=max(25,-48.05+35.82*ln(P)) CPFAE(I,1) = AMAX1( -48.05 + 35.82*ALOG(HPFAE(I)), 25. ) C Calculate CP from AE using Reiff, 1981, JGR, p 7639 (corr coeff 55%) CPFAE(I,2) = 41. + 0.11 * AE(I) C Refind NOKAE so as to define AENM and TIMAE NOKAE = NOKAE + 1 AENM(NOKAE) = AE(I) TIMAE(NOKAE) = MOD(ITMAE(2,I),100)*24 + ITMAE(3,I)/100 + | FLOAT( MOD(ITMAE(3,I),100)) / 60. ENDIF 510 CONTINUE ENDIF IF (IGETHP .EQ. 1 .AND. NVEHP .GT. 0) THEN DO 520 I=1,NVEHP TIMHP(I) = MOD(ITMEHP(2,I),100)*24 + ITMEHP(3,I)/100 + | FLOAT( MOD(ITMEHP(3,I),100)) / 60. CPFHP(I,1) = 16. + 1.3*HPEHP(I) CPFHP(I,2) = AMAX1( 25., -48.05+35.82*ALOG(HPEHP(I)) ) 520 CONTINUE ENDIF IF (IGETGPI .EQ. 1 .AND. NOKKP .GT. 0) THEN DO 530 I=1,N3HRG TIMKP(I) = MOD(ITM03G(2,I),100)*24 + ITM03G(3,I)/100 + + FLOAT( MOD(ITM03G(3,I),100)) / 60. C Median value of P for Kp=0 from Maeda plots is about 5 GW. HPFKP(I) = AMAX1(3.0, -2.78+9.33*XKPG(I) ) CPFKP(I,1) = AMAX1(-4.21+22.47*ALOG(HPFKP(I)), 10.) CPFKP(I,2) = AMAX1(-48.05+35.82*ALOG(HPFKP(I)), 25.) C From Reiff (71% correlation coeff), Reiff, JGR, 1981, p7639 (refers C to correlation coeff but not to formula CPFKP(I,3) = 29. + 11. * XKPG(I) 530 CONTINUE ENDIF C Must do this last, because will fill in IMF gaps with other C parameters found previously, i.e., CP(AE,P,Kp) NUMSW = 0 DO 535 I=1,NHRTOT DO 535 K=1,3 POT (I,K) = RMISS 535 POTF(I,K) = RMISS IF (IGETIMF .EQ. 1 .AND. NOKBZ .GT. 0) THEN NUMSW = 0 DO 540 I=1,NHRTOT DO 538 K=1,3 POT (I,K) = RMISS 538 POTF(I,K) = RMISS C Save all non-missing solar wind speed and density separately IF(SWV(I) .GT. TMISS)THEN NUMSW = NUMSW + 1 TMNOG(NUMSW) = MOD(ITMIMF(2,I),100)*24 + ITMIMF(3,I)/100 + + FLOAT( MOD(ITMIMF(3,I),100)) / 60. SWNOG(NUMSW) = SWV(I) SDNOG(NUMSW) = SWD(I) ENDIF 540 CONTINUE DO 550 I=1,NHRTOT TM = MOD(ITMIMF(2,I),100)*24 + ITMIMF(3,I)/100 + + FLOAT( MOD(ITMIMF(3,I),100)) / 60. SWIND = SWV(I) SWDEN = SWD(I) IF (SWIND.LT.TMISS .AND. BZ(I).GT.TMISS .AND. NUMSW.GT.0) THEN C Interpolate Vsw and Nsw values if within 24 hours SWIND = TERP(TMNOG,SWNOG,TM,NUMSW,24.,IFLGWSP) SWDEN = TERP(TMNOG,SDNOG,TM,NUMSW,24.,IFLGDSP) IF (IFLGWSP .EQ. 1) SWIND = SWV(I) IF (IFLGDSP .EQ. 1) SWDEN = SWD(I) ENDIF FHP = RMISS FAE = RMISS FKP = RMISS C Find interpolated P value if within 1 hour IF (NVEHP .GT. 0) FHP = TERP(TIMHP,HPEHP,TM,NVEHP,1.0,IFLGSP) C Find interpolated AE value if within 2 hours IF (IGETAE .EQ. 1 .AND. NOKAE .GT. 0) + FAE = TERP(TIMAE,AENM ,TM,NOKAE,2. ,IFLGAE) C Find interpolated Kp value if within 3 hours IF (IGETGPI .EQ. 1 .AND. NOKKP .GT. 0) + FKP = TERP(TIMKP,XKPG ,TM,N3HRG,3. ,IFLGSP) C FINDPOT calculates the crosstail potential CALL FINDPOT(BX(I),BY(I),BZ(I),SWDEN,SWIND,FHP,FAE,FKP, + POT(I,1),POT(I,2),POT(I,3)) 550 CONTINUE C Subroutine FILPOT filters the data for the crosstail potential CALL FILPOT(POT,POTF,NHRTOT,0,POT) CALL FILPOT(POT(1,2),POTF(1,2),NHRTOT,0,POT) CALL FILPOT(POT(1,3),POTF(1,3),NHRTOT,1,POT) ENDIF C Print out indicies and calculated values on 6 and lexical reads on 61 C WRITE (6,'(''1'')') WRITE (61,'(''C lexical day '',I3,'' = year'',I5,'' month'',I3, + '' day'',I3)') LEXBDAY,IBEGT(1),IBEGT(2)/100,MOD(IBEGT(2),100) IMODA1 = MOD(IBEGT(2),100) LEXADD = LEXBDAY - IMODA1 IF (IGETDST .EQ. 1 .AND. NOKDST .GT. 0) THEN CALL LEXWRT (' ','DST =', 0, DST, ITMDST, NHRTOT, 4, LEXADD) ENDIF IF (IGETAE .EQ. 1 .AND. NOKAE .GT. 0) THEN WRITE(6,"(' year mmdd hrmn Ae(hourly) calc P CP (P)',3X, + 'CP (AE)')") DO 560 I=1,NHRTOT IF(AE(I) .GT. TMISS)WRITE(6,'(3I6,F9.0,3F9.1)')(ITMAE(J,I),J=1,3) | ,AE(I),HPFAE(I),CPFAE(I,1),CPFAE(I,2) 560 CONTINUE CALL LEXWRT ('C Hourly AE (nT)', 'AE =', 0, AE, ITMAE, NHRTOT, | 4, LEXADD) CALL LEXWRT ('C P from AE (nT), P=4.28+0.0624*AE', 'CALC P =', | 1, HPFAE, ITMAE, NHRTOT, 4, LEXADD) CALL LEXWRT ('C CP from P(AE), CP=max(25,-48.05+35.82*ln(P)', | 'CALC CP =', 1, CPFAE(1,1), ITMAE, NHRTOT, 4, LEXADD) CALL LEXWRT ('C CP from AE, CP=41+0.11*AE', | 'CALC CP =', 1, CPFAE(1,2), ITMAE, NHRTOT, 4, LEXADD) ENDIF IF (IGETHP .EQ. 1 .AND. NVEHP .GT. 0) THEN WRITE (6,"(1X/1X,' year mmdd hrmn P(GW) CP(P) CP(lnP)')") DO 570 I=1,NVEHP WRITE (6,"(1X,3I5,3F8.1)") (ITMEHP(J,I),J=1,3),HPEHP(I), | CPFHP(I,1),CPFHP(I,2) 570 CONTINUE CALL LEXWRT ('C Estimated hemispheric power (GW)', 'POWER =', 1, | HPEHP, ITMEHP, NVEHP, 4, LEXADD) CALL LEXWRT ('C CP from P, CP = 16 + 1.3*P', 'CALC CP =', 1, | CPFHP(1,1), ITMEHP, NVEHP, 4, LEXADD) CALL LEXWRT ('C CP from P, CP = MAX(25.,-48.05+35.82*LN(P))', | 'CALC CP =', 1, CPFHP(1,2), ITMEHP, NVEHP, 4, LEXADD) ENDIF IF (IGETGPI .EQ. 1 .AND. NOKKP .GT. 0) THEN WRITE (6,"(1X/1X,' year mmdd hrmn Kp P(Kp) CP lnPs', | ' CP lnP CP (Kp)')") DO 580 I=1,N3HRG WRITE (6,"(1X,3I5,5F8.1)") (ITM03G(J,I),J=1,3),XKPG(I), | HPFKP(I),(CPFKP(I,K),K=1,3) 580 CONTINUE CALL LEXWRT (' ','SUNSPOTS =', 0, SSNG, ITM24G, NDAS, 4,LEXADD) CALL LEXWRT (' ','F107A =', 1, F107AG, ITM24G, NDAS, 4, LEXADD) CALL LEXWRT (' ','F107 =', 1, F107G, ITM24G, NDAS, 4, LEXADD) CALL LEXWRT (' ','APD =', 0, APDG, ITM24G, NDAS, 4, LEXADD) CALL LEXWRT (' ','KP =', 1, XKPG, ITM03G, N3HRG, 4, LEXADD) CALL LEXWRT ('C P from Kp, P = max(0.1,-2.78+9.33*Kp)', | 'CALC P =', 1, HPFKP, ITM03G, N3HRG, 4, LEXADD) CALL LEXWRT ('C CP from P(Kp), CP = max(10.,-4.21+22.47*LN(P))', | 'CALC CP =', 1, CPFKP(1,1), ITM03G, N3HRG, 4, LEXADD) CALL LEXWRT ('C CP from P(Kp), CP = max(25.,-48.05+35.82*LN(P))' | ,'CALC CP =', 1, CPFKP(1,2), ITM03G, N3HRG, 4, LEXADD) CALL LEXWRT ('C CP from Kp, CP = 29 + 11 * Kp', | 'CALC CP =', 1, CPFKP(1,3), ITM03G, N3HRG, 4, LEXADD) ENDIF IF (IGETIMF .EQ. 1 .AND. NOKBZ .GT. 0) THEN WRITE (6,"(1X/1X,' year mmdd hrmn Bx By Bz',5X, | 'Vsw Vden CPuVBQ1 CPus3ng CP vBs3 CP VBQ1 CP s3ng')") DO 590 I=1,NHRTOT IF (BZ(I) .LT. TMISS) GO TO 590 WRITE (6,"(1X,3I5,10F8.1)") (ITMIMF(J,I),J=1,3),BX(I),BY(I), | BZ(I),SWV(I),SWD(I),(POT(I,J),J=1,2),(POTF(I,K),K=1,3) 590 CONTINUE CALL LEXWRT (' ','BXIMF =', 1, BX, ITMIMF, NHRTOT, 4, LEXADD) CALL LEXWRT (' ','BYIMF =', 1, BY, ITMIMF, NHRTOT, 4, LEXADD) CALL LEXWRT (' ','BZIMF =', 1, BZ, ITMIMF, NHRTOT, 4, LEXADD) CALL LEXWRT (' ','VSW =', 0, SWV, ITMIMF, NHRTOT, 4, LEXADD) CALL LEXWRT (' ','DENSW =', 1, SWD, ITMIMF, NHRTOT, 4, LEXADD) CALL LEXWRT ('C Filtered CP from vBsin3 IMF with gaps', | 'CALC CP =', 1, POTF(1,1), ITMIMF, NHRTOT, 4, LEXADD) CALL LEXWRT ('C Filtered CP from vBQ1 IMF with gaps', | 'CALC CP =', 1, POTF(1,2), ITMIMF, NHRTOT, 4, LEXADD) CALL LEXWRT ('C Filtered CP from vBsin3 IMF with gaps filled wit |h CP(AE,P,Kp)', 'CALC CP =', 1, POTF(1,3), ITMIMF,NHRTOT,4,LEXADD) ENDIF 599 CONTINUE C Build plot time axis begin (and end) time which is rounded C down (and up) to next even hour. IBTHR(1) = IBEGT(1) IBTHR(2) = IBEGT(2) IBTHR(3) = 100*(IBEGT(3)/100) IBTHR(4) = 0 IDNHCS = DBLE(NHRTOT) * I1HR c CALL NEWTIM(IBTHR,IETHR,IDNHCS,MSGUN) C Plot Ae IF(IGETAE .EQ. 0)GO TO 600 IF(NOKAE .EQ. 0)GO TO 600 AEMIN = 0. AEMAX = RMISS CALL PLTYVT('Hourly AE Index','[nT]',AE,ITMAE,MXNAE,NHRTOT, + AEMIN,AEMAX,IBTHR,IETHR) 600 CONTINUE C Plot DST IF(IGETDST .EQ. 0)GO TO 700 DSMIN = RMISS DSMAX = RMISS IF(NOKDST .GT. 0)CALL PLTYVT('Hourly Dst Index','[nT]', + DST,ITMDST,MXNDST,NHRTOT,DSMIN,DSMAX,IBTHR,IETHR) 700 CONTINUE C Plot IMF Bx By and Bz IF(IGETIMF .EQ. 0)GO TO 800 IF(NOKBX .GT. 0)CALL PLTYVT('Hourly IMF - Bx','[nT]', + BX,ITMIMF,MXNIMF,NHRTOT,BMIN,BMAX,IBTHR,IETHR) IF(NOKBY .GT. 0)CALL PLTYVT('Hourly IMF - By','[nT]', + BY,ITMIMF,MXNIMF,NHRTOT,BMIN,BMAX,IBTHR,IETHR) IF(NOKBZ .GT. 0)CALL PLTYVT('Hourly IMF - Bz','[nT]', + BZ,ITMIMF,MXNIMF,NHRTOT,BMIN,BMAX,IBTHR,IETHR) IF(NOKSV .GT. 0)CALL PLTYVT('Hourly solar wind speed','[km/s]', + SWV,ITMIMF,MXNIMF,NHRTOT,SVMIN,SVMAX,IBTHR,IETHR) IF(NOKSD .GT. 0)CALL PLTYVT('Hourly solar wind density','[#/cm3]', + SWD,ITMIMF,MXNIMF,NHRTOT,SDMIN,SDMAX,IBTHR,IETHR) 800 CONTINUE C Lenhart tape indicies Kp, ap, Ap, F10.7, Sunspot Number c IF(IGETGPI .EQ. 0)GO TO 900 goto 900 C Build plot time axis begin and end times for GPI parameters C which extend backwards to 0 UT and forwards to 24 UT IBTG(1) = ITM24G(1,1) IBTG(2) = ITM24G(2,1) IBTG(3) = 0 IBTG(4) = 0 IETG(1) = ITM24G(1,NDAS) IETG(2) = ITM24G(2,NDAS) IETG(3) = 2400 IETG(4) = 0 ADMIN = 0. ADMAX = RMISS IF(NOKAPD .GT. 0)CALL PLTYVT('Daily Ap Index',' ', + APDG,ITM24G,MXNGPI,NDAS,ADMIN,ADMAX,IBTG,IETG) F1MIN = RMISS F1MAX = RMISS IF(NOKF10 .GT. 0)CALL PLTYVT('Daily 10.7 cm Solar Flux', + '[W/m2/Hz*10-23]', + F107G,ITM24G,MXNGPI,NDAS,F1MIN,F1MAX,IBTG,IETG) FAMIN = RMISS FAMAX = RMISS IF(NOKF1A .GT. 0)CALL PLTYVT('81 Day 10.7 cm Solar Flux', + '[W/m2/Hz*10-23]', + F107AG,ITM24G,MXNGPI,NDAS,FAMIN,FAMAX,IBTG,IETG) SSMIN = RMISS SSMAX = RMISS IF(NOKSSN .GT. 0)CALL PLTYVT('Daily Sunspot Number',' ', + SSNG,ITM24G,MXNGPI,NDAS,SSMIN,SSMAX,IBTG,IETG) XKMIN = 0. XKMAX = 9. IF(NOKKP .GT. 0)CALL PLTYVT('3-hour Kp',' ', + XKPG,ITM03G,MXNG3H,N3HRG,XKMIN,XKMAX,IBTG,IETG) A3MIN = 0. A3MAX = RMISS IF(NOKAP3 .GT. 0)CALL PLTYVT('3-hour ap',' ', + APG,ITM03G,MXNG3H,N3HRG,A3MIN,A3MAX,IBTG,IETG) 900 CONTINUE c c 5/94: Skip above plots and just make simple line plots of fluxes c and kp (this code being run 1 year at a time): c subroutine pltline(toplab,xlab,x,mxx,nx,y,mxy,ny) c write(xlab,"('DAYS of ',i4)") iyd/1000+1900 call pltline('Daily 10.7 cm Solar Flux',xlab,xdays,366,ndas, + f107g,mxngpi,ndas) call pltline('81 Day 10.7 cm Solar Flux',xlab,xdays,366,ndas, + f107ag,mxngpi,ndas) call pltline('3-hour Kp',xlab,xdays8,366*8,ndas*8, + xkpg,mxng3h,n3hrg) C Determine maximum power (GW) for consistent plot axis labeling C for requested power parameters: Evans HP, HP from Kp GWMIN = 0. GWMAX = RMISS NOKEHP = NVEHP IF (ICALCT .EQ. 1) THEN IF (IGETHP .EQ. 1) THEN NOKEHP = 0 DO 910 I=1,NVEHP IF(HPEHP(I) .LT. TMISS)GO TO 910 NOKEHP = NOKEHP + 1 GWMAX = AMAX1(HPEHP(I),GWMAX) 910 CONTINUE ENDIF IF (IGETGPI .EQ. 1) THEN NOKCHP = 0 IF (NOKKP .GT. 0) THEN DO 920 I=1,N3HRG IF (HPFKP(I) .GT. TMISS) THEN NOKCHP = NOKCHP + 1 GWMAX = AMAX1(HPFKP(I),GWMAX) ENDIF 920 CONTINUE ENDIF ENDIF IF (IGETAE .EQ. 1) THEN NOKAHP = 0 DO 930 I=1,NHRTOT IF (HPFAE(I) .GT. TMISS) THEN NOKAHP = NOACHP + 1 GWMAX = AMAX1(HPFAE(I),GWMAX) ENDIF 930 CONTINUE ENDIF ENDIF C Plot the "GW" power units parameters IF(IGETHP .EQ. 1 .AND. NOKEHP .GT. 0) + CALL PLTYVT('Estimated Hemispheric Power','[GW]', + HPEHP,ITMEHP,MXNPWR,NVEHP,GWMIN,GWMAX,IBTHR,IETHR) C Plot the calculated parameters if requested IF (ICALCT .EQ. 0) GO TO 980 IF(IGETGPI .EQ. 1 .AND. NOKCHP .GT. 0) + CALL PLTYVT('H. Power = max(0.1,-2.78+9.33Kp)','[GW]', + HPFKP,ITM03G,MXNG3H,N3HRG,GWMIN,GWMAX,IBTHR,IETHR) IF(IGETAE .EQ. 1 .AND. NOKAHP .GT. 0) + CALL PLTYVT('H. Power = 4.28+0.0624AE','[GW]', + HPFAE,ITMAE,MXNAE,NHRTOT,GWMIN,GWMAX,IBTHR,IETHR) C Determine maximum potential (KV) for consistent plot axis C labeling for requested parameters: 3 versions of CALPOT VKMIN = 0. VKMAX = RMISS NOKCP(1) = 0 NOKCP(2) = 0 NOKCP(3) = 0 IF (IGETGPI .EQ. 1 .AND. NOKKP .GT. 0) THEN DO 940 K=1,3 DO 940 I=1,N3HRG IF (CPFKP(I,K) .GT. TMISS) THEN NOKCP(K) = NOKCP(K) + 1 VKMAX = AMAX1(CPFKP(I,K),VKMAX) ENDIF 940 CONTINUE ENDIF NOACP(1) = 0 NOACP(2) = 0 IF (IGETAE .EQ. 1) THEN DO 950 K=1,2 DO 950 I=1,NHRTOT IF(CPFAE(I,K) .GT. TMISS)THEN NOACP(K) = NOACP(K) + 1 VKMAX = AMAX1(CPFAE(I,K),VKMAX) ENDIF 950 CONTINUE ENDIF NOHCP(1) = 0 NOHCP(2) = 0 IF (IGETHP .EQ. 1) THEN DO 960 K=1,2 DO 960 I=1,NVEHP IF(CPFHP(I,K) .GT. TMISS)THEN NOHCP(K) = NOHCP(K) + 1 VKMAX = AMAX1(CPFHP(I,K),VKMAX) ENDIF 960 CONTINUE ENDIF NOICP(1) = 0 NOICP(2) = 0 NOICP(3) = 0 NOICPF(1) = 0 NOICPF(2) = 0 NOICPF(3) = 0 IF (IGETIMF .EQ. 1) THEN DO 970 K=1,3 DO 970 I=1,NHRTOT IF(POT(I,K) .GT. TMISS)THEN NOICP(K) = NOICP(K) + 1 VKMAX = AMAX1(POT(I,K),VKMAX) ENDIF IF(POTF(I,K) .GT. TMISS)THEN NOICPF(K) = NOICPF(K) + 1 VKMAX = AMAX1(POTF(I,K),VKMAX) ENDIF 970 CONTINUE ENDIF C Plot the "KV" potential units parameters IF(IGETGPI .EQ. 1)THEN IF(NOKCP(1).GT. 0)CALL PLTYVT('P(Kp) CP = max{10,-4+22Ln(P)}', + '[KV]',CPFKP,ITM03G,MXNGPI,N3HRG,VKMIN,VKMAX, + IBTHR,IETHR) IF(NOKCP(2).GT.0)CALL PLTYVT('P(Kp) CP = max{25,-48+36Ln(P)}', + '[KV]',CPFKP(1,2),ITM03G,MXNGPI,N3HRG,VKMIN,VKMAX, + IBTHR,IETHR) IF(NOKCP(3) .GT. 0)CALL PLTYVT('Crosstail = 29+11Kp', + '[KV]',CPFKP(1,3),ITM03G,MXNGPI,N3HRG,VKMIN,VKMAX, + IBTHR,IETHR) ENDIF IF(IGETHP .EQ. 1)THEN IF(NOHCP(1) .GT. 0)CALL PLTYVT('Crosstail = 16+1.3P', + '[KV]',CPFHP(1,1),ITMEHP,MXNPWR,NVEHP,VKMIN,VKMAX,IBTHR,IETHR) IF(NOHCP(2).GT.0)CALL PLTYVT('Crosstail = max{25,-48+36Ln(P)}', + '[KV]',CPFHP(1,2),ITMEHP,MXNPWR,NVEHP,VKMIN,VKMAX,IBTHR,IETHR) ENDIF IF(IGETAE .EQ. 1)THEN IF(NOACP(2) .GT. 0)CALL PLTYVT('Crosstail = 41+.11AE', + '[KV]',CPFAE(1,2),ITMAE,MXNAE,NHRTOT,VKMIN,VKMAX,IBTHR,IETHR) IF(NOACP(1).GT.0)CALL PLTYVT('P(AE) CP = max{25,-48+36Ln(P)}', + '[KV]',CPFAE(1,1),ITMAE,MXNAE,NHRTOT,VKMIN,VKMAX,IBTHR,IETHR) ENDIF IF(IGETIMF .EQ. 1)THEN IF(NOICPF(1) .GT. 0)CALL PLTYVT('Crosstail = vBsin3 filtered', + '[KV]',POTF(1,1),ITMIMF,MXNIMF,NHRTOT,VKMIN,VKMAX,IBTHR,IETHR) IF(NOICPF(2) .GT. 0)CALL PLTYVT('Crosstail = vBQ1 filtered', + '[KV]',POTF(1,2),ITMIMF,MXNIMF,NHRTOT,VKMIN,VKMAX,IBTHR,IETHR) IF(NOICPF(3) .GT. 0)CALL PLTYVT('Crosstail = vBsin3 filtered gap +s filled with CP(AE,P,Kp)', + '[KV]',POTF(1,3),ITMIMF,MXNIMF,NHRTOT,VKMIN,VKMAX,IBTHR,IETHR) ENDIF 980 CONTINUE GO TO 2 1000 CONTINUE C Close GKS plotting unit CALL CLSGKS GO TO 9999 C Error trap diagnostics 9000 WRITE(MSGUN,'('' Bad time input. Begin date: '',4I5,/, + '' End date: '',4I5)')IBEGT,IENDT STOP 9050 WRITE(MSGUN,'('' Beginning date later than ending date input'',/, + '' Begin date: '',4I5,/, + '' End date: '',4I5)')IBEGT,IENDT STOP 9100 WRITE(MSGUN,'('' Dimension exceeded. Increase MXNGPI to'',I8)') + NDAS STOP 9200 WRITE(MSGUN,'('' Dimension exceeded. Increase MXNBY to'',I8)') + NHRTOT STOP 9250 WRITE(MSGUN,'('' Dimension exceeded. Increase MXNAE to'',I8)') + NHRTOT 9260 WRITE(MSGUN,'('' Dimension exceeded. Increase MXNDST to'',I8)') + NHRTOT STOP 9300 FORMAT(' Requested date',2I5,' is before first available') 9400 FORMAT(' Requested date',2I5,' is after last available') 9500 FORMAT(' No data available for date',2I5) 9600 FORMAT(' Indicies dataset corrupted? Data should be', + ' continuous yet',/, + ' somehow read past requested time',3I5) C Normal termination point 9999 CONTINUE END