      module define_atcfunix_rec
        type atcfunix_card    ! Define a new type for an atcfunix card
          sequence
          character*2   atx_basin       ! Hurricane Center Acronym
          character*2   atx_storm_num   ! Storm number (03, etc)
          character*10  atx_ymdh        ! Fcst start date (yyyymmddhh)
          character*2   atx_technum     ! Always equal to 03 for models
          character*4   atx_model       ! Name of forecast model
          integer       atx_fcsthr      ! Fcst hour (i3.3)
          integer       atx_lat         ! Storm Lat (*10), always >0
          character*1   atx_latns       ! 'N' or 'S'
          integer       atx_lon         ! Storm Lon (*10), always >0
          character*1   atx_lonew       ! 'E' or 'W'
          integer       atx_vmax        ! max sfc wind speed (knots)
          integer       atx_pcen        ! Min central pressure (mb)
c          character*2   atx_stype       ! Storm type (= 'XX' for now)
c          integer       atx_wthresh     ! Wind threshold (34, 50, 64 kt)
c          character*3   atx_quad_id     ! NEQ = start with NE quadrant
c          integer       atx_ne_quad     ! Wind radius in NE quad (nm)
c          integer       atx_se_quad     ! Wind radius in SE quad (nm)
c          integer       atx_sw_quad     ! Wind radius in SW quad (nm)
c          integer       atx_nw_quad     ! Wind radius in NW quad (nm)
        end type atcfunix_card
      end module define_atcfunix_rec

      module maxparms
        integer, parameter :: ncmaxmem = 10 ! max # ncep ensemble perts
        integer, parameter :: nrmaxmem = 10 ! max # Tom ensemble perts
        integer, parameter :: n0maxmem = 20 ! max # ncep_bc ensemble perts
        integer, parameter :: h0maxmem = 20 ! max # hwrf ensemble
        integer, parameter :: y0maxmem = 10 ! max # COAMPS ensemble
        integer, parameter :: d0maxmem = 10 ! max # GFDL ensemble
        integer, parameter :: e0maxmem = 40 ! max # HWRF+GFDL+COAMPS  ensemble
        integer, parameter :: e1maxmem = 30 ! max # HWRF+GFDL
        integer, parameter :: e2maxmem = 30 ! max # HWRF+NAVY
        integer, parameter :: e3maxmem = 20 ! max # GFDL+NAVY

        integer, parameter :: ukmaxmem = 23 ! max # UKMET ensemble perts
        integer, parameter :: frmaxmem = 34 ! max # FRANCE ensemble perts

cJ.Peng-------2010-09-29--------------------------------------
        integer, parameter :: gcmaxmem = 20 ! max # gv ensemble perts
        integer, parameter :: g0maxmem = 20 ! max # FNMOC_bc ensemble perts

        integer, parameter :: ccmaxmem = 20 ! max # cmc ensemble perts
        integer, parameter :: c0maxmem = 20 ! max # cmc_bc ensemble perts

        integer, parameter :: ecmaxmem = 50 ! max # ecmwf ensemble perts
        integer, parameter :: ezmaxmem = 20 ! max # ecmwf ensemble perts

        integer, parameter :: srmaxmem = 21 ! max # sref ensemble perts
cJ.Peng-------2010-10-29------------2010-11-02-----------------
        integer, parameter :: c2maxmem = 40 ! max # NCEP+CMC ens pert
        integer, parameter :: f2maxmem = 40 ! max # NCEP+FNMOC ens pert

        integer, parameter :: c3maxmem = 90 ! max # NCEP+CMC+EC ens pert
        integer, parameter :: n2maxmem = 70 ! max # NCEP+EC ens pert
        integer, parameter :: m2maxmem = 70 ! max # NCEP+EC ens pert

        integer, parameter :: c4maxmem = 110 ! max # NCEP+CMC+EC+FNMOC
        integer, parameter :: c5maxmem = 133 ! max # NCEP+CMC+EC+FNMOC+UK

        integer, parameter :: c6maxmem = 153 ! max # NCEP+CMC+EC+FNMOC+UK+

        integer, parameter :: n3maxmem = 93 ! max # NCEP+EC+UK

        integer, parameter :: f3maxmem = 60 ! max # NCEP+CMC+FNMOC

        integer, parameter :: s3maxmem =  2 ! max # eemn+2emn ens pert
        integer, parameter :: p3maxmem =  3 ! max # ee/ae/ce+mn ens pert

        integer, parameter :: d3maxmem =  3 ! max # AVNO+CMC+EMX  pert
        integer, parameter :: d6maxmem =  6 ! max # AVNO+CMC+EMX+AEMN+..

        integer, parameter :: x6maxmem =  6 ! max # EEMN+AEMN+CEMN+FEMN+
        integer, parameter :: d4maxmem =  4 ! max # AVNO+CMC+EMX+3EMN
        integer, parameter :: d2maxmem =  2 ! max # 3EMN+3DET..

        integer, parameter :: ncminmem =  4 ! min # of ncep ensemble 
                         ! perts needed at a given fcst hr to get a mean
        integer, parameter :: nrminmem =  4 ! min # of Tom ensemble
        integer, parameter :: n0minmem =  8 ! min # of ncep_bc ensemble
        integer, parameter :: h0minmem =  8 ! min # of hwrf ensemble
        integer, parameter :: y0minmem =  4 ! min # of coamps ensemble
        integer, parameter :: d0minmem =  4 ! min # of GFDL ensemble
        integer, parameter :: e0minmem = 16 ! min # of GFDL+HWRF+COAMPS  ensemble
        integer, parameter :: e1minmem = 12 ! min # of GFDL+HWRF
        integer, parameter :: e2minmem = 12 ! min # of HWRF+NAVY
        integer, parameter :: e3minmem = 8 ! min # of GFDL+NAVY

        integer, parameter :: ukminmem =  9 ! min # of UKMET ensemble
        integer, parameter :: frminmem = 14 ! min # of FRANCE ensemble

cJ.Peng----2010-09-29---------------------------------------------
        integer, parameter :: gcminmem =  8 ! min # of gv ensemble
        integer, parameter :: g0minmem =  8 ! min # of gv ensemble

        integer, parameter :: ccminmem =  8 ! min # of cmc ensemble
        integer, parameter :: c0minmem =  8 ! min # of cmc ensemble

        integer, parameter :: ecminmem = 20 ! min # of ecmwf ensemble
                         ! perts needed at a given fcst hr to get a mean

        integer, parameter :: ezminmem =  8 ! min # of ecmwf ensemble

        integer, parameter :: srminmem =  8 ! min # of sref ensemble 
                         ! perts needed at a given fcst hr to get a mean
cJ.Peng-------2010-10-29-------------2010-11-02--------------------------
c        integer, parameter :: c2minmem = 20 ! min # of NCEP+CMC ensemble
c                         ! perts needed at a given fcst hr to get a mean
cJ.Peng------------2010-11-19--------------------------
        integer, parameter :: c2minmem = 22 ! min # of NCEP+CMC ensemble
                         ! perts needed at a given fcst hr to get a mean
        integer, parameter :: f2minmem = 22 ! min # of NCEP+FNMOC ensemble

        integer, parameter :: c3minmem = 50 ! min # of NCEP+CMC+EC ensemble
        integer, parameter :: n2minmem = 39 ! min # of NCEP+EC ensemble
        integer, parameter :: m2minmem = 39 ! min # of NCEP+EC ensemble

        integer, parameter :: c4minmem = 61 ! min # of NCEP+CMC+EC+FNMOC
        integer, parameter :: c5minmem = 74 ! min # of NCEP+CMC+EC+FNMOC+UK

        integer, parameter :: c6minmem = 85 ! min # of NCEP+CMC+EC+FNMOC+UK
        integer, parameter :: n3minmem = 52 ! min # of NCEP+EC+UK

        integer, parameter :: f3minmem = 33 ! min # of NCEP+CMC+FNMOC

        integer, parameter :: s3minmem =  2 ! min # of NCEP+CMC+EC ensemble
        integer, parameter :: p3minmem =  3 ! min # of NCEP+CMC+EC ensemble

        integer, parameter :: d3minmem =  3 ! min # of AVNO+CMC+EMX ensemble
        integer, parameter :: d6minmem =  6 ! min # of AVNO+CMC+EMX+AEMN+.. 

        integer, parameter :: x6minmem =  6 ! min # of
        integer, parameter :: d4minmem =  4 ! min # of AVNO+CMC+EMX+3EMN
        integer, parameter :: d2minmem =  2 ! min # of 3DET+3EMN..

        integer, parameter :: maxstorms  = 15 ! max # storms at a time
        integer, parameter :: maxtimes   = 65 ! max # fcst times (0-384
                                              ! every 6 hrs = 65).
c        integer, parameter :: max_accum_prob_hour = 384 ! Max hour to 
c        integer, parameter :: max_accum_prob_hour = 120 ! Max hour to 
c                              ! which we'll carry out the accum probs
        real, parameter :: grdspc = 1.00      ! grid spacing of our grid
      end module maxparms

      module ensinfo
        USE maxparms
c        character*4 :: ncperts(ncmaxmem) = (/'AP01','AP02','AP03'
c     &          ,'AP04','AP05','AP06','AP07','AP08','AP09','AP10'
c     &          ,'AP11','AP12','AP13','AP14','AP15','AP16','AP17'
c     &          ,'AP18','AP19','AP20'/)

        character*4 :: ncperts(ncmaxmem) = (/'AP01','AP02','AP03'
     &       ,'AP04','AP05','AP06','AP07','AP08','AP09','AP10'/)

        character*4 :: nrperts(nrmaxmem) = (/'GR01','GR02','GR03'
     &          ,'GR04','GR05','GR06','GR07','GR08','GR09','GR10'
     &          /)

c        character*4 :: ncperts(ncmaxmem) = (/'GE01','GE02','GE03'
c     &          ,'GE04','GE05','GE06','GE07','GE08','GE09','GE10'
c     &          ,'GE11','GE12','GE13','GE14','GE15','GE16','GE17'
c     &          ,'GE18','GE19','GE20'/)

        character*4 :: n0perts(n0maxmem) = (/'AP01','AP02','AP03'
     &          ,'AP04','AP05','AP06','AP07','AP08','AP09','AP10'
     &          ,'AP11','AP12','AP13','AP14','AP15','AP16','AP17'
     &          ,'AP18','AP19','AP20'/)

        character*4 :: h0perts(h0maxmem) = (/'HW01','HW02','HW03'
     &          ,'HW04','HW05','HW06','HW07','HW08','HW09','HW10'
     &          ,'HW11','HW12','HW13','HW14','HW15','HW16','HW17'
     &          ,'HW18','HW19','HW20'/)

        character*4 :: y0perts(y0maxmem) = (/'C01C','C02C','C03C'
     &          ,'C04C','C05C','C06C','C07C','C08C','C09C','C10C'
     &          /)

        character*4 :: d0perts(d0maxmem) = (/'GP00','GP01','GP02'
     &          ,'GP03','GP04','GP05','GP06','GP07','GP08','GP09'
     &          /)

        character*4 :: e0perts(e0maxmem) = (/'GP00','GP01','GP02'
     &          ,'GP03','GP04','GP05','GP06','GP07','GP08','GP09'
     &          ,'C01C','C02C','C03C'
     &          ,'C04C','C05C','C06C','C07C','C08C','C09C','C10C'
     &          ,'HW01','HW02','HW03'
     &          ,'HW04','HW05','HW06','HW07','HW08','HW09','HW10'
     &          ,'HW11','HW12','HW13','HW14','HW15','HW16','HW17'
     &          ,'HW18','HW19','HW20'/)

        character*4 :: e1perts(e1maxmem) = (/'GP00','GP01','GP02'
     &          ,'GP03','GP04','GP05','GP06','GP07','GP08','GP09'
     &          ,'HW01','HW02','HW03'
     &          ,'HW04','HW05','HW06','HW07','HW08','HW09','HW10'
     &          ,'HW11','HW12','HW13','HW14','HW15','HW16','HW17'
     &          ,'HW18','HW19','HW20'/)

        character*4 :: e2perts(e2maxmem) = (/'C01C','C02C','C03C'
     &          ,'C04C','C05C','C06C','C07C','C08C','C09C','C10C'
     &          ,'HW01','HW02','HW03'
     &          ,'HW04','HW05','HW06','HW07','HW08','HW09','HW10'
     &          ,'HW11','HW12','HW13','HW14','HW15','HW16','HW17'
     &          ,'HW18','HW19','HW20'/)

        character*4 :: e3perts(e3maxmem) = (/'GP00','GP01','GP02'
     &          ,'GP03','GP04','GP05','GP06','GP07','GP08','GP09'
     &          ,'C01C','C02C','C03C'
     &          ,'C04C','C05C','C06C','C07C','C08C','C09C','C10C'
     &          /)

        character*4 :: ukperts(ukmaxmem) = (/'UK01','UK02'
     &   ,'UK03','UK04','UK05','UK06','UK07','UK08','UK09','UK10'
     &          ,'UK11','UK12','UK13','UK14','UK15','UK16','UK17'
     &          ,'UK18','UK19','UK20','UK21','UK22','UK23'/)

        character*4 :: frperts(frmaxmem) = (/'FR01','FR02','FR03'
     &          ,'FR04','FR05','FR06','FR07','FR08','FR09','FR10'
     &          ,'FR11','FR12','FR13','FR14','FR15','FR16','FR17'
     &          ,'FR18','FR19','FR20','FR21','FR22','FR23','FR24'
     &          ,'FR25','FR26','FR27','FR28','FR29','FR30','FR31'
     &          ,'FR32','FR33','FR34'/)

c#-------J.Peng----2010-09-29---------------------------------------------
        character*4 :: gcperts(gcmaxmem) = (/'FP01','FP02','FP03'
     &          ,'FP04','FP05','FP06','FP07','FP08','FP09','FP10'
     &          ,'FP11','FP12','FP13','FP14','FP15','FP16','FP17'
     &          ,'FP18','FP19','FP20'/)
        character*4 :: g0perts(g0maxmem) = (/'FP01','FP02','FP03'
     &          ,'FP04','FP05','FP06','FP07','FP08','FP09','FP10'
     &          ,'FP11','FP12','FP13','FP14','FP15','FP16','FP17'
     &          ,'FP18','FP19','FP20'/)


        character*4 :: ccperts(ccmaxmem) = (/'CP01','CP02','CP03'
     &          ,'CP04','CP05','CP06','CP07','CP08','CP09','CP10'
     &          ,'CP11','CP12','CP13','CP14','CP15','CP16','CP17'
     &          ,'CP18','CP19','CP20'/)
        character*4 :: c0perts(c0maxmem) = (/'CP01','CP02','CP03'
     &          ,'CP04','CP05','CP06','CP07','CP08','CP09','CP10'
     &          ,'CP11','CP12','CP13','CP14','CP15','CP16','CP17'
     &          ,'CP18','CP19','CP20'/)

cJ.Peng-------2010-10-29-------------------------------------
        character*4 :: c2perts(c2maxmem) = (/'AP01','AP02','AP03'
     &          ,'AP04','AP05','AP06','AP07','AP08','AP09','AP10'
     &          ,'AP11','AP12','AP13','AP14','AP15','AP16','AP17'
     &          ,'AP18','AP19','AP20'
     &          ,'CP01','CP02','CP03'
     &          ,'CP04','CP05','CP06','CP07','CP08','CP09','CP10'
     &          ,'CP11','CP12','CP13','CP14','CP15','CP16','CP17'
     &          ,'CP18','CP19','CP20'/)

cJ.Peng-------2012-04-27-------------------------------------
        character*4 :: f2perts(f2maxmem) = (/'AP01','AP02','AP03'
     &          ,'AP04','AP05','AP06','AP07','AP08','AP09','AP10'
     &          ,'AP11','AP12','AP13','AP14','AP15','AP16','AP17'
     &          ,'AP18','AP19','AP20'
     &          ,'FP01','FP02','FP03'
     &          ,'FP04','FP05','FP06','FP07','FP08','FP09','FP10'
     &          ,'FP11','FP12','FP13','FP14','FP15','FP16','FP17'
     &          ,'FP18','FP19','FP20'/)

cJ.Peng-------2010-11-02-------------------------------------------
        character*4 :: c3perts(c3maxmem) = (/'AP01','AP02','AP03'
     &          ,'AP04','AP05','AP06','AP07','AP08','AP09','AP10'
     &          ,'AP11','AP12','AP13','AP14','AP15','AP16','AP17'
     &          ,'AP18','AP19','AP20'
     &          ,'CP01','CP02','CP03'
     &          ,'CP04','CP05','CP06','CP07','CP08','CP09','CP10'
     &          ,'CP11','CP12','CP13','CP14','CP15','CP16','CP17'
     &          ,'CP18','CP19','CP20'
     &          ,'EN01','EN02','EN03'
     &          ,'EN04','EN05','EN06','EN07','EN08','EN09','EN10'
     &          ,'EN11','EN12','EN13','EN14','EN15','EN16','EN17'
     &          ,'EN18','EN19','EN20','EN21','EN22','EN23','EN24'
     &          ,'EN25','EP01','EP02','EP03','EP04','EP05','EP06'
     &          ,'EP07','EP08','EP09','EP10','EP11','EP12','EP13'
     &          ,'EP14','EP15','EP16','EP17','EP18','EP19','EP20'
     &          ,'EP21','EP22','EP23','EP24','EP25'/)

        character*4 :: n2perts(n2maxmem) = (/'AP01','AP02','AP03'
     &          ,'AP04','AP05','AP06','AP07','AP08','AP09','AP10'
     &          ,'AP11','AP12','AP13','AP14','AP15','AP16','AP17'
     &          ,'AP18','AP19','AP20'
     &          ,'EN01','EN02','EN03'
     &          ,'EN04','EN05','EN06','EN07','EN08','EN09','EN10'
     &          ,'EN11','EN12','EN13','EN14','EN15','EN16','EN17'
     &          ,'EN18','EN19','EN20','EN21','EN22','EN23','EN24'
     &          ,'EN25','EP01','EP02','EP03','EP04','EP05','EP06'
     &          ,'EP07','EP08','EP09','EP10','EP11','EP12','EP13'
     &          ,'EP14','EP15','EP16','EP17','EP18','EP19','EP20'
     &          ,'EP21','EP22','EP23','EP24','EP25'/)

        character*4 :: m2perts(m2maxmem) = (/'01AI','02AI','03AI'
     &          ,'04AI','05AI','06AI','07AI','08AI','09AI','10AI'
     &          ,'11AI','12AI','13AI','14AI','15AI','16AI','17AI'
     &          ,'18AI','19AI','20AI'
     &          ,'01N2','02N2','03N2'
     &          ,'04N2','05N2','06N2','07N2','08N2','09N2','10N2'
     &          ,'11N2','12N2','13N2','14N2','15N2','16N2','17N2'
     &          ,'18N2','19N2','20N2','21N2','22N2','23N2','24N2'
     &          ,'25N2'
     &          ,'01P2','02P2','03P2','04P2','05P2','06P2'
     &          ,'07P2','08P2','09P2','10P2','11P2','12P2','13P2'
     &          ,'14P2','15P2','16P2','17P2','18P2','19P2','20P2'
     &          ,'21P2','22P2','23P2','24P2','25P2'/)

        character*4 :: f3perts(f3maxmem) = (/'AP01','AP02','AP03'
     &          ,'AP04','AP05','AP06','AP07','AP08','AP09','AP10'
     &          ,'AP11','AP12','AP13','AP14','AP15','AP16','AP17'
     &          ,'AP18','AP19','AP20'
     &          ,'CP01','CP02','CP03'
     &          ,'CP04','CP05','CP06','CP07','CP08','CP09','CP10'
     &          ,'CP11','CP12','CP13','CP14','CP15','CP16','CP17'
     &          ,'CP18','CP19','CP20'
     &          ,'FP01','FP02','FP03'
     &          ,'FP04','FP05','FP06','FP07','FP08','FP09','FP10'
     &          ,'FP11','FP12','FP13','FP14','FP15','FP16','FP17'
     &          ,'FP18','FP19','FP20'/)


        character*4 :: c4perts(c4maxmem) = (/'AP01','AP02','AP03'
     &          ,'AP04','AP05','AP06','AP07','AP08','AP09','AP10'
     &          ,'AP11','AP12','AP13','AP14','AP15','AP16','AP17'
     &          ,'AP18','AP19','AP20'
     &          ,'CP01','CP02','CP03'
     &          ,'CP04','CP05','CP06','CP07','CP08','CP09','CP10'
     &          ,'CP11','CP12','CP13','CP14','CP15','CP16','CP17'
     &          ,'CP18','CP19','CP20'
     &          ,'EN01','EN02','EN03'
     &          ,'EN04','EN05','EN06','EN07','EN08','EN09','EN10'
     &          ,'EN11','EN12','EN13','EN14','EN15','EN16','EN17'
     &          ,'EN18','EN19','EN20','EN21','EN22','EN23','EN24'
     &          ,'EN25','EP01','EP02','EP03','EP04','EP05','EP06'
     &          ,'EP07','EP08','EP09','EP10','EP11','EP12','EP13'
     &          ,'EP14','EP15','EP16','EP17','EP18','EP19','EP20'
     &          ,'EP21','EP22','EP23','EP24','EP25'
     &          ,'FP01','FP02','FP03'
     &          ,'FP04','FP05','FP06','FP07','FP08','FP09','FP10'
     &          ,'FP11','FP12','FP13','FP14','FP15','FP16','FP17'
     &          ,'FP18','FP19','FP20'/)

        character*4 :: c5perts(c5maxmem) = (/'AP01','AP02','AP03'
     &          ,'AP04','AP05','AP06','AP07','AP08','AP09','AP10'
     &          ,'AP11','AP12','AP13','AP14','AP15','AP16','AP17'
     &          ,'AP18','AP19','AP20'
     &          ,'CP01','CP02','CP03'
     &          ,'CP04','CP05','CP06','CP07','CP08','CP09','CP10'
     &          ,'CP11','CP12','CP13','CP14','CP15','CP16','CP17'
     &          ,'CP18','CP19','CP20'
     &          ,'EN01','EN02','EN03'
     &          ,'EN04','EN05','EN06','EN07','EN08','EN09','EN10'
     &          ,'EN11','EN12','EN13','EN14','EN15','EN16','EN17'
     &          ,'EN18','EN19','EN20','EN21','EN22','EN23','EN24'
     &          ,'EN25','EP01','EP02','EP03','EP04','EP05','EP06'
     &          ,'EP07','EP08','EP09','EP10','EP11','EP12','EP13'
     &          ,'EP14','EP15','EP16','EP17','EP18','EP19','EP20'
     &          ,'EP21','EP22','EP23','EP24','EP25'
     &          ,'FP01','FP02','FP03'
     &          ,'FP04','FP05','FP06','FP07','FP08','FP09','FP10'
     &          ,'FP11','FP12','FP13','FP14','FP15','FP16','FP17'
     &          ,'FP18','FP19','FP20'
     &          ,'UK01','UK02','UK03'
     &          ,'UK04','UK05','UK06','UK07','UK08','UK09','UK10'
     &          ,'UK11','UK12','UK13','UK14','UK15','UK16','UK17'
     &          ,'UK18','UK19','UK20','UK21','UK22','UK23'/)

        character*4 :: c6perts(c6maxmem) = (/'AP01','AP02','AP03'
     &          ,'AP04','AP05','AP06','AP07','AP08','AP09','AP10'
     &          ,'AP11','AP12','AP13','AP14','AP15','AP16','AP17'
     &          ,'AP18','AP19','AP20'
     &          ,'CP01','CP02','CP03'
     &          ,'CP04','CP05','CP06','CP07','CP08','CP09','CP10'
     &          ,'CP11','CP12','CP13','CP14','CP15','CP16','CP17'
     &          ,'CP18','CP19','CP20'
     &          ,'EN01','EN02','EN03'
     &          ,'EN04','EN05','EN06','EN07','EN08','EN09','EN10'
     &          ,'EN11','EN12','EN13','EN14','EN15','EN16','EN17'
     &          ,'EN18','EN19','EN20','EN21','EN22','EN23','EN24'
     &          ,'EN25','EP01','EP02','EP03','EP04','EP05','EP06'
     &          ,'EP07','EP08','EP09','EP10','EP11','EP12','EP13'
     &          ,'EP14','EP15','EP16','EP17','EP18','EP19','EP20'
     &          ,'EP21','EP22','EP23','EP24','EP25'
     &          ,'FP01','FP02','FP03'
     &          ,'FP04','FP05','FP06','FP07','FP08','FP09','FP10'
     &          ,'FP11','FP12','FP13','FP14','FP15','FP16','FP17'
     &          ,'FP18','FP19','FP20'
     &          ,'UK01','UK02','UK03'
     &          ,'UK04','UK05','UK06','UK07','UK08','UK09','UK10'
     &          ,'UK11','UK12','UK13','UK14','UK15','UK16','UK17'
     &          ,'UK18','UK19','UK20','UK21','UK22','UK23'
     &          ,'GE01','GE02','GE03'
     &          ,'GE04','GE05','GE06','GE07','GE08','GE09','GE10'
     &          ,'GE11','GE12','GE13','GE14','GE15','GE16','GE17'
     &          ,'GE18','GE19','GE20'/)


        character*4 :: n3perts(n3maxmem) = (/'AP01','AP02','AP03'
     &          ,'AP04','AP05','AP06','AP07','AP08','AP09','AP10'
     &          ,'AP11','AP12','AP13','AP14','AP15','AP16','AP17'
     &          ,'AP18','AP19','AP20'
     &          ,'EN01','EN02','EN03'
     &          ,'EN04','EN05','EN06','EN07','EN08','EN09','EN10'
     &          ,'EN11','EN12','EN13','EN14','EN15','EN16','EN17'
     &          ,'EN18','EN19','EN20','EN21','EN22','EN23','EN24'
     &          ,'EN25','EP01','EP02','EP03','EP04','EP05','EP06'
     &          ,'EP07','EP08','EP09','EP10','EP11','EP12','EP13'
     &          ,'EP14','EP15','EP16','EP17','EP18','EP19','EP20'
     &          ,'EP21','EP22','EP23','EP24','EP25'
     &          ,'UK01','UK02','UK03'
     &          ,'UK04','UK05','UK06','UK07','UK08','UK09','UK10'
     &          ,'UK11','UK12','UK13','UK14','UK15','UK16','UK17'
     &          ,'UK18','UK19','UK20','UK21','UK22','UK23'/)

        character*4 :: s3perts(s3maxmem) = (/'EEMN','2EMN'/)
        character*4 :: p3perts(p3maxmem) = (/'EEMN','AEMN','CEMN'/)

        character*4 :: d3perts(d3maxmem) = (/'AVNO',' EMX',' CMC'/)

        character*4 :: d6perts(d6maxmem) = (/'AEMN','CEMN','EEMN'
     &                                      ,'AVNO',' CMC',' EMX'/)

        character*4 :: x6perts(x6maxmem) = (/'AEMN','CEMN','EEMN'
     &                                      ,'FEMN','UKMN','HPMN'/)

        character*4 :: d4perts(d4maxmem) = (/'AVNO',' EMX',' CMC'
     &                                      ,'3EMN'/)

        character*4 :: d2perts(d2maxmem) = (/'3EMN','3DET'/)


        character*4 :: ecperts(ecmaxmem) = (/'EN01','EN02','EN03'
     &          ,'EN04','EN05','EN06','EN07','EN08','EN09','EN10'
     &          ,'EN11','EN12','EN13','EN14','EN15','EN16','EN17'
     &          ,'EN18','EN19','EN20','EN21','EN22','EN23','EN24'
     &          ,'EN25','EP01','EP02','EP03','EP04','EP05','EP06'
     &          ,'EP07','EP08','EP09','EP10','EP11','EP12','EP13'
     &          ,'EP14','EP15','EP16','EP17','EP18','EP19','EP20'
     &          ,'EP21','EP22','EP23','EP24','EP25'/)

        character*4 :: ezperts(ezmaxmem) = (/'EN01','EN02','EN03'
     &          ,'EN04','EN05','EN06','EN07','EN08','EN09','EN10' 
     &          ,'EP01','EP02','EP03','EP04','EP05','EP06'
     &          ,'EP07','EP08','EP09','EP10'/)

        character*4 :: srperts(srmaxmem) = (/'EMC1','EMN1','EMN2'
     &          ,'EMN3','EMP1','EMP2','EMP3','NBC1','NBN1','NBN2'
     &          ,'NBN3','NBP1','NBP2','NBP3','NMC1','NMN1','NMN2'
     &          ,'NMN3','NMP1','NMP2','NMP3'/)

        integer  :: ncfcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: nrfcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: n0fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: h0fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: y0fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: d0fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: e0fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: e1fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: e2fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: e3fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: ukfcsthrs(maxtimes) = (/0,12,24,36,48,60,72,84,96
     &               ,108,120,132
     &               ,144,156,168,180,192,204,216,228,240,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: frfcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

c#-------J.Peng----2010-09-29---------------------------------------------
        integer  :: gcfcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)
        integer  :: g0fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: ccfcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)
        integer  :: c0fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)


cJ.Peng-------2010-10-29-------------------------------------
        integer  :: c2fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: f2fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

cJ.Peng-------2010-11-02------------------------------------------
c        integer  :: c3fcsthrs(maxtimes) = (/0,12,24,36,48,60,72,84,96
c     &               ,108,120,132,144,156,168,180,192,204,216,228,240
c     &               ,999,999,999,999,999,999,999,999,999,999,999,999
c     &               ,999,999,999,999,999,999,999,999,999,999,999,999
c     &               ,999,999,999,999,999,999,999,999,999,999,999,999
c     &               ,999,999,999,999,999,999,999,999/)

        integer  :: c3fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: n2fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: m2fcsthrs(maxtimes) = (/0,12,24,36,48,60,72,84,96
     &               ,108,120,132
     &               ,144,156,168,180,192,204,216,228,240,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: f3fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: c4fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: c5fcsthrs(maxtimes) = (/0,12,24,36,48,60,72,84,96
     &               ,108,120,132
     &               ,144,156,168,180,192,204,216,228,240,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: c6fcsthrs(maxtimes) = (/0,12,24,36,48,60,72,84,96
     &               ,108,120,132
     &               ,144,156,168,180,192,204,216,228,240,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: n3fcsthrs(maxtimes) = (/0,12,24,36,48,60,72,84,96
     &               ,108,120,132
     &               ,144,156,168,180,192,204,216,228,240,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: s3fcsthrs(maxtimes) = (/0,12,24,36,48,60,72,84,96
     &               ,108,120,132,144,156,168,180,192,204,216,228,240
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999/)
        integer  :: p3fcsthrs(maxtimes) = (/0,12,24,36,48,60,72,84,96
     &               ,108,120,132,144,156,168,180,192,204,216,228,240
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999/)

        integer  :: d3fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: d6fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: x6fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: d4fcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/) 

        integer  :: d2fcsthrs(maxtimes) = (/0,12,24,36,48,60,72,84,96
     &               ,108,120,132,144,156,168,180,192,204,216,228,240
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999/)

c------J.Peng------2011-01-03----------------------------------------
c        integer  :: ecfcsthrs(maxtimes) = (/0,12,24,36,48,60,72,84,96
c     &               ,108,120,132,144,156,168,180,192,204,216,228,240
c     &               ,999,999,999,999,999,999,999,999,999,999,999,999
c     &               ,999,999,999,999,999,999,999,999,999,999,999,999
c     &               ,999,999,999,999,999,999,999,999,999,999,999,999
c     &               ,999,999,999,999,999,999,999,999/)

        integer  :: ecfcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)

        integer  :: ezfcsthrs(maxtimes) = (/0,6,12,18,24,30,36,42,48,54
     &               ,60,66,72,78,84,90,96,102,108,114,120,126,132,138
     &               ,144,150,156,162,168,174,180,186,192,198,204,210
     &               ,216,222,228,234,240,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)


        integer  :: srfcsthrs(maxtimes) = (/0,3,6,9,12,15,18,21,24,27
     &               ,30,33,36,39,42,45,48,51,54,57,60,63,66,69
     &               ,72,75,78,81,84,87,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999,999,999,999,999,999,999,999
     &               ,999,999,999,999,999/)
        character*2, save :: run_basin
      end module ensinfo

      module grid_bounds
        integer, parameter :: maxbasins = 11
c---J.Peng-----2015-01-09---------------------------------------------
        character*2 :: basin(maxbasins) = (/'AL','EP','WP','CP'
     &                 ,'AA','BB','SI','SP','SC','EC','AU'/)
        real :: xmaxlat(maxbasins) = (/ 60.0, 60.0, 60.0, 60.0, 60.0
     &                        , 60.0, -1.0, -1.0,999.0,999.0,999.0/)

        real :: xminlat(maxbasins) = (/  1.0,  1.0,  1.0,  1.0,  1.0
     &                        ,  1.0,-60.0,-60.0,999.0,999.0,999.0/)

        real :: xmaxlon(maxbasins) = (/345.0,285.0,185.0,255.0,125.0
     &                        ,145.0,125.0,185.0,999.0,999.0,999.0/)

        real :: xminlon(maxbasins) = (/250.0,190.0, 90.0,160.0, 30.0
     &                        , 50.0, 30.0, 90.0,999.0,999.0,999.0/)
      end module grid_bounds

      module trig_vals
        real, save :: pi, dtr
        real, save :: dtk = 111.1949     ! Dist (km) over 1 deg lat
                                         ! using erad=6371.0e+3
        real, save :: erad = 6371.0e+3   ! Earth's radius (m)
        real, save :: ecircum = 40030.2  ! Earth's circumference
                                         ! (km) using erad=6371.e3
        real, save :: omega = 7.292e-5
      end module trig_vals

c
c--------------------------------------------------------------------
c
      program trakave
c
c$$$  MAIN PROGRAM DOCUMENTATION BLOCK
c
c Main Program: TRAKAVE      Get average of ensemble tracks
C   PRGMMR: MARCHOK          ORG: NP22        DATE: 2005-08-01
c
c ABSTRACT: This program reads in the tropical cyclone track files 
c   (in "atcfunix" format) from the various ensemble perturbation
c   members and produces an ensemble mean track.  Currently, the
c   program can handle members from the NCEP, ECMWF and CANADIAN
c   ensembles.
c
c Program history log:
c   2005-08-01  Marchok - Original operational version.
c   2006-08-02  Marchok - Changed pertubation IDs for Canadian ens
c                         (they no longer have "negative" perts).
c   2007-08-08  Marchok - Added code for averaging SREF ensemble tracks.
c
c Input files:
c   unit   11    Text file containing the concatenated atcfunix tracks
c                with all storms for all ensemble members at all 
c                forecast hours.
c
c Output files:
c   unit   52    Opened and closed internally, 1 file for each storm,
c                this series of output files contains the gridded 
c                strike probabilities, both for the individual forecast
c                hours as well as the accumulated probabilities.
c                pt (1,1) is the far NW point on this output grid, 
c                pt (imax,jmax) is the far SE point on this output grid.
c   unit   53    Output atcfunix file with mean ensemble track forecast
c                positions for each storm at each forecast hour.
c   unit   54    Output atcf file with mean ensemble track forecast
c                positions for each storm forecast hours 12-72.
c   unit   55    Output track file with mean ensemble track forecast
c                positions for each storm forecast hours 0-72.
c   unit   56    Output text file with ensemble track forecast spread
c                values for each storm at each forecast hour.
c   unit   57    Output text file with ensemble track forecast mode
c                values for each storm at each forecast hour.
c
c Subprograms called:
c   read_nlists  Read input namelists for input date & storm number info
c
c Attributes:
c   Language: Fortran90
c
c Notes:   
c   maxmem     This is the max number of members expected for a 
c              given Center's ensemble forecast
c   minmem     This is the minimum number of member tracks that need to 
c              exist at a given forecast hour in order to produce a mean
c              
c
c$$$
c
c-------
c
      USE define_atcfunix_rec; USE maxparms; USE trig_vals
      USE ensinfo

      type (atcfunix_card) saverec

      real,allocatable :: traklat(:,:,:)
      real,allocatable :: traklon(:,:,:)
cJ.Peng-04-15-2013
      real,allocatable :: tcvmx(:,:,:)
      real,allocatable :: tcpce(:,:,:)

      real    xmeanlat(maxstorms,maxtimes)
      real    xmeanlon(maxstorms,maxtimes)
cJ.Peng-04-15-2013
      real    xmeanvmx(maxstorms,maxtimes)
      real    xmeanpce(maxstorms,maxtimes)

      real    xspread(maxstorms,maxtimes)
      real    dthresh
      integer  :: fcsthrs(maxtimes)
      integer minmem,maxmem,fcst_interval
      character*4,allocatable :: perts(:)
      character*1 :: do_probs = 'y'
      character*1 :: modeflag = 'y'
      character*4 :: stormarr(maxstorms)
      character*4 :: catcf,catcf_lc
      character*5 :: cmodel
c
      pi = 4. * atan(1.)   ! Both pi and dtr were declared in module
      dtr = pi/180.0       ! trig_vals, but were not yet defined.

      xmeanlat = -9999.0; xmeanlon = -9999.0
cJ.Peng-04-15-2013
      xmeanvmx = -9999.0; xmeanpce = -9999.0

      xspread = -999.0

      call read_nlists (dthresh,cmodel)

      select case (cmodel)
        case ('ens');  maxmem = ncmaxmem; fcsthrs(:) = ncfcsthrs(:)
        case ('ref');  maxmem = nrmaxmem; fcsthrs(:) = nrfcsthrs(:)

        case ('ensb');  maxmem = n0maxmem; fcsthrs(:) = n0fcsthrs(:)
        case ('hf21');  maxmem = h0maxmem; fcsthrs(:) = h0fcsthrs(:)
        case ('ny10');  maxmem = y0maxmem; fcsthrs(:) = y0fcsthrs(:)
        case ('gd10');  maxmem = d0maxmem; fcsthrs(:) = d0fcsthrs(:)
        case ('egn3');  maxmem = e0maxmem; fcsthrs(:) = e0fcsthrs(:)
        case ('egt2');  maxmem = e1maxmem; fcsthrs(:) = e1fcsthrs(:)
        case ('ent2');  maxmem = e2maxmem; fcsthrs(:) = e2fcsthrs(:)
        case ('gnt2');  maxmem = e3maxmem; fcsthrs(:) = e3fcsthrs(:)

        case ('ukes');  maxmem = ukmaxmem; fcsthrs(:) = ukfcsthrs(:)

        case ('fres');  maxmem = frmaxmem; fcsthrs(:) = frfcsthrs(:)

c----------J.Peng----------2010-09-29-------------------------------
c        case ('gve');  maxmem = gcmaxmem; fcsthrs(:) = gcfcsthrs(:)
        case ('fens');  maxmem = gcmaxmem; fcsthrs(:) = gcfcsthrs(:)
        case ('fenb');  maxmem = g0maxmem; fcsthrs(:) = g0fcsthrs(:)

        case ('cens'); maxmem = ccmaxmem; fcsthrs(:) = ccfcsthrs(:)
        case ('cenb'); maxmem = c0maxmem; fcsthrs(:) = c0fcsthrs(:)


cJ.Peng-------2010-10-29--------------2010-11-02---------------
        case ('g2c'); maxmem = c2maxmem; fcsthrs(:) = c2fcsthrs(:)

        case ('naf2'); maxmem = f2maxmem; fcsthrs(:) = f2fcsthrs(:)

        case ('gce3'); maxmem = c3maxmem; fcsthrs(:) = c3fcsthrs(:)
        case ('nae2'); maxmem = n2maxmem; fcsthrs(:) = n2fcsthrs(:)
        case ('inae'); maxmem = m2maxmem; fcsthrs(:) = m2fcsthrs(:)

        case ('ncf3'); maxmem = f3maxmem; fcsthrs(:) = f3fcsthrs(:)

        case ('gce4'); maxmem = c4maxmem; fcsthrs(:) = c4fcsthrs(:)
        case ('gce5'); maxmem = c5maxmem; fcsthrs(:) = c5fcsthrs(:)

        case ('gce6'); maxmem = c6maxmem; fcsthrs(:) = c6fcsthrs(:)
        case ('neu3'); maxmem = n3maxmem; fcsthrs(:) = n3fcsthrs(:)


        case ('gce3a'); maxmem = s3maxmem; fcsthrs(:) = s3fcsthrs(:)
        case ('gce3b'); maxmem = p3maxmem; fcsthrs(:) = p3fcsthrs(:)

        case ('sing'); maxmem = d3maxmem; fcsthrs(:) = d3fcsthrs(:)
        case ('tot6'); maxmem = d6maxmem; fcsthrs(:) = d6fcsthrs(:)

        case ('six6'); maxmem = x6maxmem; fcsthrs(:) = x6fcsthrs(:)
        case ('tot4'); maxmem = d4maxmem; fcsthrs(:) = d4fcsthrs(:)

        case ('tot2'); maxmem = d2maxmem; fcsthrs(:) = d2fcsthrs(:)

        case ('eens');  maxmem = ecmaxmem; fcsthrs(:) = ecfcsthrs(:)
        case ('eefs');  maxmem = ezmaxmem; fcsthrs(:) = ezfcsthrs(:)

        case ('sref'); maxmem = srmaxmem; fcsthrs(:) = srfcsthrs(:)
      end select
         
      allocate (traklat(maxstorms,maxmem,maxtimes),stat=itla)
      allocate (traklon(maxstorms,maxmem,maxtimes),stat=itlo)
cJ.Peng-04-15-2013
      allocate (tcvmx(maxstorms,maxmem,maxtimes),stat=itlp)
      allocate (tcpce(maxstorms,maxmem,maxtimes),stat=itlq)

      allocate (perts(maxmem),stat=ipa)
      if (itla/=0.or.itlo/=0.or.ipa/=0.or.itlp/=0.or.itlq/=0) then
        print *,' '
        print *,'!!! ERROR in trakave allocating arrays.'
        print *,'!!! itla = ',itla,' itlo= ',itlo,' ipa= ',ipa
        stop 94
      endif

      traklat = -9999.0;  traklon = -9999.0
cJ.Peng-04-15-2013
      tcvmx = -9999.0
      tcpce = -9999.0

      select case (cmodel)
        case ('ens');  perts(:) = ncperts(:)
                       minmem   = ncminmem  
                       catcf    = 'AEMN'  
                       catcf_lc = 'aemn'  

        case ('ref');  perts(:) = nrperts(:)
                       minmem   = nrminmem
                       catcf    = 'GRMN'
                       catcf_lc = 'grmn'

        case ('ensb');  perts(:) = n0perts(:)
                       minmem   = n0minmem
                       catcf    = 'AEMB'
                       catcf_lc = 'aemb'

        case ('hf21');  perts(:) = h0perts(:)
                       minmem   = h0minmem
                       catcf    = 'HEMN'
                       catcf_lc = 'hemn'

        case ('ny10');  perts(:) = y0perts(:)
                       minmem   = y0minmem
                       catcf    = 'YEMN'
                       catcf_lc = 'yemn'

        case ('gd10');  perts(:) = d0perts(:)
                       minmem   = d0minmem
                       catcf    = 'GEMN'
                       catcf_lc = 'gemn'

        case ('egn3');  perts(:) = e0perts(:)
                       minmem   = e0minmem
                       catcf    = 'MEMN'
                       catcf_lc = 'memn'

        case ('egt2');  perts(:) = e1perts(:)
                       minmem   = e1minmem
                       catcf    = 'GHMN'
                       catcf_lc = 'ghmn'

        case ('ent2');  perts(:) = e2perts(:)
                       minmem   = e2minmem
                       catcf    = 'HYMN'
                       catcf_lc = 'hymn'

        case ('gnt2');  perts(:) = e3perts(:)
                       minmem   = e3minmem
                       catcf    = 'GYMN'
                       catcf_lc = 'gymn'

        case ('ukes');  perts(:) = ukperts(:)
                       minmem   = ukminmem
                       catcf    = 'UKMN'
                       catcf_lc = 'ukmn'

        case ('fres');  perts(:) = frperts(:)
                       minmem   = frminmem
                       catcf    = 'FRMN'
                       catcf_lc = 'frmn'

c        case ('gve');  perts(:) = gcperts(:)
c                       minmem   = gcminmem  
c                       catcf    = 'GVMN'  
c                       catcf_lc = 'gvmn'  
c-------J.Peng---------2010-09-29--------------------------
        case ('fens');  perts(:) = gcperts(:)
                       minmem   = gcminmem
                       catcf    = 'FEMN'
                       catcf_lc = 'femn'
        case ('fenb');  perts(:) = g0perts(:)
                       minmem   = g0minmem
                       catcf    = 'FEMB'
                       catcf_lc = 'femb'


        case ('cens'); perts(:) = ccperts(:)
                       minmem   = ccminmem  
                       catcf    = 'CEMN'  
                       catcf_lc = 'cemn'  

        case ('cenb'); perts(:) = c0perts(:)
                       minmem   = c0minmem
                       catcf    = 'CEMB'
                       catcf_lc = 'cemb'

cJ.Peng-------2010-10-29-------------------------------------
        case ('g2c'); perts(:) = c2perts(:)
                       minmem   = c2minmem 
                       catcf    = '2EMN' 
                       catcf_lc = '2emn'

        case ('naf2'); perts(:) = f2perts(:)
                       minmem   = f2minmem
                       catcf    = '2NAF'
                       catcf_lc = '2naf'

cJ.Peng-------2010-11-02-----------------------------
        case ('gce3'); perts(:) = c3perts(:)
                       minmem   = c3minmem
                       catcf    = '3EMN'
                       catcf_lc = '3emn'

        case ('nae2'); perts(:) = n2perts(:)
                       minmem   = n2minmem
                       catcf    = '2NAE'
                       catcf_lc = '2nae'

        case ('inae'); perts(:) = m2perts(:)
                       minmem   = m2minmem
                       catcf    = 'NAEI'
                       catcf_lc = 'naei'

        case ('ncf3'); perts(:) = f3perts(:)
                       minmem   = f3minmem
                       catcf    = '3NCF'
                       catcf_lc = '3ncf'

        case ('gce4'); perts(:) = c4perts(:)
                       minmem   = c4minmem
                       catcf    = '4EMN'
                       catcf_lc = '4emn'

        case ('gce5'); perts(:) = c5perts(:)
                       minmem   = c5minmem
                       catcf    = '5EMN'
                       catcf_lc = '5emn'

        case ('gce6'); perts(:) = c6perts(:)
	               minmem   = c6minmem
		       catcf    = 'MIX6'
		       catcf_lc = 'mix6'

        case ('neu3'); perts(:) = n3perts(:)
                       minmem   = n3minmem
                       catcf    = '3NEU'
                       catcf_lc = '3neu'

        case ('gce3a'); perts(:) = s3perts(:)
                       minmem   = s3minmem
                       catcf    = 'SEMN'
                       catcf_lc = 'semn'
        case ('gce3b'); perts(:) = p3perts(:)
                       minmem   = p3minmem
                       catcf    = 'PEMN'
                       catcf_lc = 'pemn'

        case ('sing'); perts(:) = d3perts(:)
                       minmem   = d3minmem
                       catcf    = '3DET'
                       catcf_lc = '3det'
        case ('tot6'); perts(:) = d6perts(:)
                       minmem   = d6minmem
                       catcf    = '6EMN'
                       catcf_lc = '6emn'

        case ('six6'); perts(:) = x6perts(:)
                       minmem   = x6minmem
                       catcf    = 'SIXM'
                       catcf_lc = 'sixm'

        case ('tot4'); perts(:) = d4perts(:)
                       minmem   = d4minmem
                       catcf    = 'T4MN'
                       catcf_lc = 't4mn'

        case ('tot2'); perts(:) = d2perts(:)
                       minmem   = d2minmem
                       catcf    = 'S6MN'
                       catcf_lc = 's6mn'


        case ('eens');  perts(:) = ecperts(:)
                       minmem   = ecminmem  
                       catcf    = 'EEMN'  
                       catcf_lc = 'eemn'  

        case ('eefs');  perts(:) = ezperts(:)
                       minmem   = ezminmem
                       catcf    = 'TEMN'
                       catcf_lc = 'temn'

        case ('sref'); perts(:) = srperts(:)
                       minmem   = srminmem  
                       catcf    = 'SFMN'  
                       catcf_lc = 'sfmn'  
      end select
 
      fcst_interval = fcsthrs(2) - fcsthrs(1)

cJ.Peng-04-15-2013
      call read_atcfunix (traklat,traklon,tcvmx,tcpce,
     & saverec,stormarr,maxmem,fcst_interval,perts)

      do i = 1,maxstorms
        print *,'i= ',i,'  stormarr(i)= ',stormarr(i)
      enddo

cJ.Peng-04-15-2013
      call get_mean_track (traklat,traklon,xmeanlat,xmeanlon
     &            ,tcvmx,tcpce,xmeanvmx,xmeanpce
     &            ,stormarr,saverec,maxmem,minmem,catcf,fcsthrs)

      call get_spread (traklat,traklon,xmeanlat,xmeanlon,xspread
     &                ,stormarr,saverec,maxmem,minmem,fcsthrs)

      if (do_probs == 'y') then
        call get_probs (traklat,traklon,xmeanlat,xmeanlon
     &             ,stormarr,saverec%atx_ymdh,dthresh,'indiv'
     &             ,modeflag,catcf,maxmem,fcsthrs,catcf_lc)
        call get_probs (traklat,traklon,xmeanlat,xmeanlon
     &             ,stormarr,saverec%atx_ymdh,dthresh,'accum'
     &             ,modeflag,catcf,maxmem,fcsthrs,catcf_lc)
      endif

      print *,' '
      print *,'NORMAL END TO TRAKAVE....'
      print *,' '

      deallocate (traklat,stat=ictla)
      deallocate (traklon,stat=ictlo)

cJ.Peng-04-15-2013
      deallocate (tcvmx,stat=ictlp)
      deallocate (tcpce,stat=ictlq)

      deallocate (perts,stat=icpa)

      if (ictla/=0.or.ictlo/=0.or.icpa/=0.or.ictlp/=0.or.ictlq/=0) then
        print *,' '
        print *,'!!! ERROR in trakave deallocating arrays.'
        print *,'!!! ictla = ',ictla,' ictlo= ',ictlo,' icpa= ',icpa
        stop 94
      endif

      stop 
      end

c
c----------------------------------------------------------------------
c
c----------------------------------------------------------------------
      subroutine get_probs (traklat,traklon,xmeanlat,xmeanlon
     &                     ,stormarr,cymdh,dthresh,prob_type
     &                     ,modeflag,catcf,maxmem,fcsthrs,catcf_lc)
c
c     ABSTRACT: This subroutine calculates the probability of a
c     given location (gridpoint) being within a specified distance
c     of an ensemble member trackpoint.  For the domain, we'll sort
c     through all the track points and get the most extreme
c     N, E, S and W points, then add a certain amount onto that on
c     each side.  Then, for all the gridpoints that are included,
c     we'll go one by one to each gridpoint and calculate the 
c     distances from each forecast track point at that hour to 
c     each gridpoint.  If that distance is less than the specified
c     distance threshold, we add it to the total, and then divide
c     that total count by the total number of members to get the 
c     probability at that gridopoint for that forecast hour.
c     
c     NOTE: This particular subroutine can calculate the 
c     probabilities individually for each forecast hour or for
c     the total accumulated probability over the forecast, 
c     depending on the value of the input variable prob_type.
c     If we are doing this for individual cases, then we will
c     also call  get_prob_mode in order to get the ensemble 
c     probability mode and print it out in a separate atcfunix
c     record.


      USE maxparms; USE grid_bounds; USE trig_vals; USE ensinfo

      real  traklat(maxstorms,maxmem,maxtimes)
      real  traklon(maxstorms,maxmem,maxtimes)
      real  xmeanlat(maxstorms,maxtimes)
      real  xmeanlon(maxstorms,maxtimes)
      real, allocatable :: xprob(:,:)
      real  :: xnm_per_degree = 60.0    
      real  xdist,dthresh,xlondiff,xlatdiff,xinterplon,xinterplat
      character*1 :: modeflag,hold_input_modeflag
      character*2 :: cbasin
      character*4 :: stormarr(maxstorms),cstorm,catcf,catcf_lc
      character*5 :: prob_type
      character*10   cymdh
      character*12   cacc
      character*46   outgrfile
      character*33   outtxtfile
      character      cthresh1*1,cthresh2*2,cthresh3*3,cthresh4*4
      integer :: fcsthrs(maxtimes)
      integer  basinix,ifh_interval,ibct
      integer, allocatable :: accum_ct(:,:)
      integer  max_accum_prob_hour
      logical(1),allocatable :: gridpoint_already_hit(:,:,:)
c
c      outtxtfile = 'enstrkprob.'//cymdh//'.ctlinfo.txt'
c      open (unit=51,file=outtxtfile,status='unknown')

      if (catcf_lc == 'sfmn') then
        max_accum_prob_hour = 87   ! SREF max fcst hour = 87 in Aug 2007
      else
c--------J.Peng----2011-01-11---------------------------------------
        max_accum_prob_hour = 168  ! For all others, go out to 120h
      endif

      hold_input_modeflag = modeflag

      stormloop: do ist = 1,maxstorms

        modeflag = hold_input_modeflag

        ! As of June, 2002, we will only compute these probabilities
        ! for storms in the AL, EP and WP basins.  This next code
        ! checks for those basin ID's and will cycle the stormloop
        ! if it's not one of those basins....

        cbasin = stormarr(ist)(1:2)
c--------J.Peng----2015-01-09---------------------------------
        if (cbasin == 'AL' .or. cbasin == 'EP' .or. cbasin == 'WP' 
     & .or. cbasin == 'CP' .or. cbasin == 'AA' .or. cbasin == 'BB'
     & .or. cbasin == 'SI' .or. cbasin == 'SP' )then

          continue
        else
          print *,' '
          print *,'+++ NOTE: In get_probs, we are NOT calculating'
          print *,'+++       the probabilities for storm '
     &                      ,stormarr(ist)
          print *,'+++       It is not in a basin we do probs for.'
          print *,' '
          cycle stormloop
        endif

        ! Get the index for the basin we're in.

        basinix = -99
        ibct = 1
        basinloop: do while (ibct <= maxbasins)
          if (cbasin == basin(ibct)) then
            basinix = ibct
            exit basinloop
          endif
          ibct = ibct + 1
        enddo basinloop
        if (basinix == -99) then
          print *,' '
          print *,'!!! ERROR IN GET_PROBS.  BASIN IS NOT '
          print *,'!!! RECOGNIZED.  BASIN = ',cbasin
          print *,'!!! EXITING....'
          stop 99
        endif

        ! In an older version of this program, we narrowed the area 
        ! over which we computed the probabilities to keep the file 
        ! size smaller, but that meant that each output file had 
        ! different starting and ending points.  In this new version
        ! we just calculate it over the whole basin.

        maxlat = int(xmaxlat(basinix))
        maxlon = int(xmaxlon(basinix))
        minlat = int(xminlat(basinix))
        minlon = int(xminlon(basinix))

        numlat = abs(maxlat-minlat) + 1
        numlon = abs(maxlon-minlon) + 1

        itemp = int(dthresh)
        if (itemp < 10) then
          write (cthresh1,'(i1)') itemp
          outgrfile = catcf_lc//'.trkprob.'//stormarr(ist)//'.'//
     &              cthresh1//'nm.'//cymdh//'.'//prob_type//'.ieee'
        else if (itemp >= 10 .and. itemp < 100) then
          write (cthresh2,'(i2)') itemp
          outgrfile = catcf_lc//'.trkprob.'//stormarr(ist)//'.'//
     &              cthresh2//'nm.'//cymdh//'.'//prob_type//'.ieee'
        else if (itemp >= 100 .and. itemp < 1000) then
          write (cthresh3,'(i3)') itemp
          outgrfile = catcf_lc//'.trkprob.'//stormarr(ist)//'.'//
     &              cthresh3//'nm.'//cymdh//'.'//prob_type//'.ieee'
        else if (itemp >= 1000) then
          write (cthresh4,'(i4)') itemp
          outgrfile = catcf_lc//'.trkprob.'//stormarr(ist)//'.'//
     &              cthresh4//'nm.'//cymdh//'.'//prob_type//'.ieee'
        endif

        print *,' '
        print *,'Before open 52, stormarr(ist) = ',stormarr(ist)
        print *,' '

        if (xmeanlon(ist,1) > 0.0) then
          print *,' '
          print *,'!!! OPEN:  '
          print *,'   open for prob_type = ',prob_type
          print *,'   ist = ',ist,'  xmeanlon = ',xmeanlon(ist,1)
          print *,'   file name = ',outgrfile
          print *,'   dthresh = ',dthresh
          inquire (52,access=cacc) 
          print *,'   cacc 1 = ',cacc
          open (unit=52,file=outgrfile,access='direct'
     &         ,form='unformatted',iostat=ios
     &         ,status='unknown',recl=numlat*numlon*4)
          print *,'   ios = ',ios
          inquire (52,access=cacc) 
          print *,'   cacc 2 = ',cacc
        else
          cycle stormloop
        endif

        cstorm = stormarr(ist)

        if (prob_type == 'indiv') then
          write (51,82) stormarr(ist),cymdh,numlon,minlon,numlat,minlat
        endif

   82   format (a4,1x,a10,' nx= ',i4,' beglon= ',i3
     &                   ,' ny= ',i4,' beglat= ',i3)

        ! Allocate the xprob array.  Pt (1,1) is in the far NW part 
        ! of the subgrid.

        if (.not.allocated(xprob)) then
          allocate (xprob(numlon,numlat),stat=ipa)
          if (ipa /= 0) then
            print *,' '
            print *,'!!! ERROR in sub get_probs allocating xprob array.'
            print *,'!!! ipa = ',ipa
            igpret = 94
            return
          endif
        endif

        if (.not.allocated(gridpoint_already_hit)) then
          allocate (gridpoint_already_hit(numlon,numlat,maxmem)
     &             ,stat=iga)
          if (iga /= 0) then
            print *,' '
            print *,'!!! ERROR in sub get_probs allocating '
            print *,'!!! gridpoint_already_hit array.'
            print *,'!!! iga = ',iga
            igpret = 94
            return
          endif
        endif

        if (.not.allocated(accum_ct)) then
          allocate (accum_ct(numlon,numlat),stat=iaa)
          if (iaa /= 0) then
            print *,' '
            print *,'!!! ERROR in sub get_probs allocating '
            print *,'!!! accum_ct array.'
            print *,'!!! iaa = ',iaa
            igpret = 94
            return
          endif
        endif

        ! Now go through each forecast time and get the
        ! probabilities at each grid point.  We check the xmeanlon
        ! at each forecast time to make sure that we still have a
        ! mean track (i.e., that the forecast isn't over or that
        ! the storm didn't dissipate).
        !
        ! Note below that the indexing of the xprob array will go
        ! so that pt (1,1) is in the upper left (most NW) corner
        ! of the subgrid, and point (maxlon,maxlat) is in the 
        ! lower right (most SE) corner of the subgrid.

        print *,' '
        print *,' '
        print *,'++++++++++++++++++++++++ '
        print *,' +++ NEW STORM +++'
        print *,'++++++++++++++++++++++++ '
        print *,'ist= ',ist,' maxlat= ',maxlat,' minlat= ',minlat
        print *,' '

c--------J.Peng-----2010-10-26---------------------
        gridpoint_already_hit = .false.
        accum_ct = 0

        ifhloop: do ifh = 1,maxtimes

          ifhour = fcsthrs(ifh) 
          print *,' '
          print *,'++++++++++++++++++++++++ '
          print *,' +++ NEW TIME +++'
          print *,'++++++++++++++++++++++++ '
          print *,'ifh= ',ifh,' ist= ',ist,'  forecast hour = ',ifhour

          if (ifhour > max_accum_prob_hour .and. prob_type == 'accum')
     &    then
            print *,' '
            print *,'!!! NOT CALCULATING ACCUM PROBS PAST '
     &             ,max_accum_prob_hour,'hours'
            exit ifhloop
          endif

          if (xmeanlon(ist,ifh) > 0.0) then    ! valid mean fcst at ifh

            xprob = 0.0
            latloop: do jlat = maxlat,minlat,-1

c              print *,'jlat= ',jlat

              jix = maxlat - jlat + 1
              ylat = float(jlat)

              lonloop: do ilon = minlon,maxlon

c                print *,'ilon= ',ilon

                iix = ilon - minlon + 1
                xlon = float(ilon)

                ict = 0
                ipt = 0
                memberloop: do im = 1,maxmem

                  ict = ict + 1
                  if (traklon(ist,im,ifh) > 0.0) then
                    call calcdist (xlon,ylat,traklon(ist,im,ifh)
     &                            ,traklat(ist,im,ifh),dtr,xdist)
                    if (xdist <= dthresh) then
                      if (prob_type == 'accum') then
c--------J.Peng-----2010-10-26---------------------
                        if (gridpoint_already_hit(iix,jix,im)) then
                          continue
                        else
                          accum_ct(iix,jix) = accum_ct(iix,jix) + 1
                          gridpoint_already_hit(iix,jix,im) = .true.
                        endif
c                          accum_ct(iix,jix) = accum_ct(iix,jix) + 1
                      else
                        ipt = ipt + 1
c                        print *,'++YES, dist= ',xdist,' ipt= ',ipt
                      endif
                    endif 
                  endif

                  ! * * * * * * * * * * * * * * * * * * * * * * * * * *
                  ! INTERPOLATE TRACKS FOR "ACCUMULATED" PROBABILITIES:
                  ! For accum probs, interpolate the tracks in between
                  ! time periods to get a smooth track and not miss 
                  ! any "in-between" points.  Don't need to do this for
                  ! the individual, instantaneous track probs at each
                  ! forecast hour.  Currently, since TPC only does their
                  ! accumulated probabilities out to 120h, we will only
                  ! do ours to 120h as well, so don't interpolate if the
                  ! the next forecast hour would be > 120.
                  ! * * * * * * * * * * * * * * * * * * * * * * * * * *

                  if (prob_type == 'accum' .and. ifh < maxtimes .and.
     &                fcsthrs(ifh+1) <= max_accum_prob_hour) then
                    if (traklon(ist,im,ifh)   > 0.0 .and. 
     &                  traklon(ist,im,ifh+1) > 0.0) then
                      ifh_interval = fcsthrs(ifh+1) - fcsthrs(ifh)
                      xlondiff = traklon(ist,im,ifh+1) - 
     &                           traklon(ist,im,ifh)
                      xlatdiff = traklat(ist,im,ifh+1) - 
     &                           traklat(ist,im,ifh)
                      interploop: do iint = 1,ifh_interval-1
                        xinterplon = traklon(ist,im,ifh) + 
     &                               (float(iint)/float(ifh_interval) *
     &                                xlondiff)
                        xinterplat = traklat(ist,im,ifh) + 
     &                               (float(iint)/float(ifh_interval) *
     &                                xlatdiff)
                        call calcdist (xlon,ylat,xinterplon
     &                                ,xinterplat,dtr,xdist)
                        if (xdist <= dthresh) then
c--------J.Peng-----2010-10-26---------------------
                          if (gridpoint_already_hit(iix,jix,im)) then
                            continue
                          else
                            accum_ct(iix,jix) = accum_ct(iix,jix) + 1
                            gridpoint_already_hit(iix,jix,im) = .true.
                          endif
c                            accum_ct(iix,jix) = accum_ct(iix,jix) + 1
                        endif
                      enddo interploop
                    endif
                  endif

                enddo memberloop

                if (maxmem > 0) then
c                  print *,'xprob: ipt= ',ipt,' ict= ',ict,' prob= '
c     &                     ,float(ipt)/float(maxmem)
                  xprob(iix,jix) = float(ipt)/float(maxmem)
                endif

              enddo lonloop

            enddo latloop

c           If the probability type is for an individual forecast hour,
c           then write out the data right now for this hour, and then
c           if the modeflag is set to find the probability mode, call
c           the subroutine  get_prob_mode.

            if (prob_type == 'indiv') then
c              do iii = 1,numlon
c                do jjj = 1,numlat
c                  if (xprob(iii,jjj) > 0.0) then
c                    print *,'HEY, iii= ',iii,' jjj= ',jjj,' xprob= '
c     &                     ,xprob(iii,jjj)
c                  endif
c                enddo
c              enddo 
cccccc-------J.Peng.------2010-09-30--------------------------------
              print *,'JPENG_ifh=',ifh 
              write (52,rec=ifh) ((xprob(i,j),i=1,numlon),j=1,numlat)
              if (modeflag == 'y') then
                call get_prob_mode (numlon,numlat,minlon,maxlon,minlat
     &               ,maxlat,xprob,xmodelon,xmodelat,xmodeval,ifhour
     &               ,cstorm,cymdh,ifh,dthresh,xmeanlon,xmeanlat,ist
     &               ,basinix,catcf,fcsthrs,modeflag)
              endif
            endif

          else

            cycle ifhloop

          endif
              
        enddo ifhloop

        if (prob_type == 'accum') then
          jloop: do jlat = maxlat,minlat,-1
            jix = maxlat - jlat + 1
            ylat = float(jlat)
            iloop: do ilon = minlon,maxlon
              iix = ilon - minlon + 1
              xlon = float(ilon)
              xprob(iix,jix) = float(accum_ct(iix,jix)) /
     &                         float(maxmem)
c              print *,'jix= ',jix,' iix= ',iix,' xprob= ',xprob(iix,jix)
c              write (61,120) jix,iix,float(jlat),float(ilon)
c     &                       ,xprob(iix,jix)
            enddo iloop
          enddo jloop
          irec=1
          write (52,rec=irec) ((xprob(i,j),i=1,numlon),j=1,numlat)
        endif

        deallocate (xprob,stat=icxa)
        if (icxa /= 0) then
          print *,' '
          print *,'!!! ERROR in trakave deallocating xprob array.'
          print *,'!!! icxa = ',icxa
          stop 94
        endif

        close (52)

      enddo stormloop

 120  format (1x,'jix= ',i3,' iix= ',i3,' jlat= ',f6.1,' ilon= '
     &          ,f6.1,' xprob= ',f6.2)

      icgaha = 0
      icaca  = 0

      if (allocated(gridpoint_already_hit)) then
        deallocate (gridpoint_already_hit,stat=icgaha)
      endif

      if (allocated(accum_ct)) then
        deallocate (accum_ct,stat=icaca)
      endif

      if (icgaha /= 0 .or. icaca /= 0) then
        print *,' '
        print *,'!!! ERROR in trakave deallocating gridpoint'
        print *,'!!! or accum_ct arrays.'
        print *,'!!! icgaha = ',icgaha,' icaca= ',icaca
        stop 94
      endif
c
      return
      end
       
c
c----------------------------------------------------------------
c
c----------------------------------------------------------------
      subroutine get_mean_track (traklat,traklon,xmeanlat,xmeanlon
     &            ,tcvmx,tcpce,xmeanvmx,xmeanpce
     &              ,stormarr,saverec,maxmem,minmem,catcf,fcsthrs)

c     ABSTRACT: This subroutine calculates the mean lat and lon 
c     positions for each storm at each forecast hour.  maxmem is
c     the maximum number of members for a given Center's forecast,
c     while minmem is the minimum number of member tracks that must
c     exist at a given forecast hour in order to produce a forecast
c     track.  Currently, I have this set at ~40% of maxmem, so for
c     NCEP (maxmem=10,minmem=4), CMC (maxmem=16,minmem=6),
c     ECMWF (maxmem=50,minmem=20)

      USE define_atcfunix_rec; USE maxparms
      USE trig_vals

      type (atcfunix_card) saverec

      real  traklat(maxstorms,maxmem,maxtimes)
      real  traklon(maxstorms,maxmem,maxtimes)
      real  xmeanlat(maxstorms,maxtimes)
      real  xmeanlon(maxstorms,maxtimes) 
cJ.Peng-04-15-2013----
      real  tcvmx(maxstorms,maxmem,maxtimes)
      real  tcpce(maxstorms,maxmem,maxtimes)
      real  xmeanvmx(maxstorms,maxtimes)
      real  xmeanpce(maxstorms,maxtimes)
      real  stdvmx(maxstorms,maxtimes)
      real  stdpce(maxstorms,maxtimes)
      real  yspread(maxstorms,maxtimes)

      integer :: fcsthrs(maxtimes)
      integer pairct
      character*4 :: stormarr(maxstorms)
      character*4 :: catcf
      yspread =-999.0 
c
      print *,' '
      print *,'************************** '
      print *,'*** In get_mean_track.... '
      print *,'************************** '
      print *,' '

      stormloop: do ist = 1,maxstorms

        print *,' '
        print *,'-------------------'
        print *,'storm number = ',ist
        print *,'storm name   = ',stormarr(ist)
        print *,' '

        timeloop: do ifh = 1,maxtimes

          ict = 0
          xlatsum = 0.0
          xlonsum = 0.0

cJ.Peng-04-15-2013-----
	  xvmxsum = 0.0
	  xpcesum = 0.0 

          print *,' '
          print *,' time index (ifh) = ',ifh
          print *,' '

          do imem = 1,maxmem
            if (traklon(ist,imem,ifh) > 0.0) then
              ict = ict + 1
              xlatsum = xlatsum + traklat(ist,imem,ifh)
cJ.Peng------03/27/2015-----------------------------------
c              xlonsum = xlonsum + traklon(ist,imem,ifh)
              if(traklon(ist,imem,ifh) < 30.0) then
                xlonsum = xlonsum + traklon(ist,imem,ifh)+360.0
              else
                xlonsum = xlonsum + traklon(ist,imem,ifh)
              endif
cJ.Peng------03/27/2015-----------------------------------

cJ.Peng-04-15-2013-----
              xvmxsum = xvmxsum + tcvmx(ist,imem,ifh)
   	      xpcesum = xpcesum + tcpce(ist,imem,ifh)

              print *,'mem= ',imem,' lat= ',traklat(ist,imem,ifh)
     &                            ,' lon= ',traklon(ist,imem,ifh)
            else
              print *,'mem= ',imem,' lat= ',traklat(ist,imem,ifh)
     &                            ,' lon= ',traklon(ist,imem,ifh)
            endif
          enddo


c          pairct = 0
c          if (traklon(ist,1,ifh) > 0.0 .and. traklon(ist,6,ifh) > 0.0)
c     &    then
c            pairct = pairct + 1
c          endif
c          if (traklon(ist,2,ifh) > 0.0 .and. traklon(ist,7,ifh) > 0.0)
c     &    then
c            pairct = pairct + 1
c          endif
c          if (traklon(ist,3,ifh) > 0.0 .and. traklon(ist,8,ifh) > 0.0)
c     &    then
c            pairct = pairct + 1
c          endif
c          if (traklon(ist,4,ifh) > 0.0 .and. traklon(ist,9,ifh) > 0.0)
c     &    then
c            pairct = pairct + 1
c          endif
c          if (traklon(ist,5,ifh) > 0.0 .and. traklon(ist,10,ifh) > 0.0)
c     &    then
c            pairct = pairct + 1
c          endif
c          if (ict > 2 .and. pairct >= 2) then


          if (ict >= minmem) then
            xmeanlat(ist,ifh) = xlatsum / ict
            xmeanlon(ist,ifh) = xlonsum / ict
            print *,'MEAN:  LAT= ',xmeanlat(ist,ifh)
     &                   ,' LON= ',xmeanlon(ist,ifh)

cJ.Peng-04-15-2013-----------
	    xmeanvmx(ist,ifh) = xvmxsum / ict
	    xmeanpce(ist,ifh) = xpcesum / ict
cJ.Peng-04-15-2013---spread for intensity----wind and SLP
            jct = 0
            svmxsum = 0.0
            spcesum = 0.0
            do imem = 1,maxmem
              if (traklon(ist,imem,ifh) > 0.0) then
		jct = jct + 1
		svmxsum = svmxsum +
     &               (tcvmx(ist,imem,ifh)-xmeanvmx(ist,ifh))**2
                spcesum = spcesum + 
     &               (tcpce(ist,imem,ifh)-xmeanpce(ist,ifh))**2
	      endif 
            enddo 
            if (jct >= minmem) then
	      stdvmx(ist,ifh) = sqrt(svmxsum/jct)
	      stdpce(ist,ifh) = sqrt(spcesum/jct)
            endif
cJ.Peng-04-15-2013---spread for track------------------
            xdistsum = 0.0
	    imemct   =   0
	    do imem = 1,maxmem
              if (traklon(ist,imem,ifh) > 0.0) then
		call calcdist (xmeanlon(ist,ifh),xmeanlat(ist,ifh)
     &              ,traklon(ist,imem,ifh),traklat(ist,imem,ifh)
     &              ,dtr,xdist)
		xdistsum = xdistsum + xdist
		imemct   = imemct + 1
              endif
            enddo
	    if (imemct >= minmem) then
	      yspread(ist,ifh) = xdistsum / float(imemct)
            endif
c---------------------------------------------------------------
c            call output_atcfunix (stormarr(ist),saverec%atx_ymdh
c     &           ,xmeanlat(ist,ifh),xmeanlon(ist,ifh),ifh,'mean'
c     &           ,-99.0,-99.0,catcf,fcsthrs)
cJ.Peng-04-15-2013---------
            call output_atcfunix1 (stormarr(ist),saverec%atx_ymdh
     &           ,xmeanlat(ist,ifh),xmeanlon(ist,ifh)
     &           ,stdvmx(ist,ifh),stdpce(ist,ifh),yspread(ist,ifh)
     &           ,xmeanvmx(ist,ifh),xmeanpce(ist,ifh),ifh,'mean'
     &           ,-99.0,-99.0,catcf,fcsthrs)

c          else if (ict > 0 .and. ict <= 2) then
          else if (ict > 0 .and. ict < minmem) then
            xmlat = xlatsum / ict
            xmlon = xlonsum / ict
            print *,' '
            print *,'!!! NOT ENOUGH MEMBERS FOR MEAN FOR STORM #',ist
            print *,'!!! STORM NAME = ',stormarr(ist)
            print *,'!!! MODEL NAME = ',catcf
            print *,'!!! NUMBER OF MEMBERS = ',ict,'  TIME INDEX = ',ifh
            print *,'!!! MINIMUM NUMBER OF MEMBERS NEEDED= ',minmem
            print *,'!!! MEAN POSITION WOULD HAVE BEEN: '
            print *,'!!! LAT: ',xmlat,'  LON: ',xmlon
          else
            print *,'!!! UNDEFINED MEAN POSITION !!!'
          endif

        enddo timeloop
c------J.Peng---------2010-09-29------------------------------------
c        if (xmeanlon(ist,1) > 0.0) then
c          call output_atcf (stormarr(ist),xmeanlon,xmeanlat
c     &                      ,saverec%atx_ymdh,ist,catcf)
c          call output_all (stormarr(ist),xmeanlon,xmeanlat
c     &                      ,saverec%atx_ymdh,ist,catcf)
c     &                    
c        endif

      enddo stormloop

      return
      end

c
c------------------------------------------------------------------
c
c------------------------------------------------------------------
      subroutine get_spread (traklat,traklon,xmeanlat,xmeanlon
     &          ,xspread,stormarr,saverec,maxmem,minmem,fcsthrs)
c
c     ABSTRACT: This subroutine calculates the spread of the 
c     ensemble track forecasts, defined as the average distance
c     of the ensemble members to the mean track.

      USE trig_vals; USE maxparms; USE define_atcfunix_rec
      USE ensinfo

      type (atcfunix_card) saverec

      real      traklat(maxstorms,maxmem,maxtimes)
      real      traklon(maxstorms,maxmem,maxtimes)
      real      xmeanlat(maxstorms,maxtimes)
      real      xmeanlon(maxstorms,maxtimes)
      real      xspread(maxstorms,maxtimes)
      real      xdistsum,xdist
      integer :: fcsthrs(maxtimes)
      integer   imemct,ict,ifh_interval,ifhour
      character*4 :: stormarr(maxstorms)
c
      print *,' '
      print *,'************************** '
      print *,'*** In get_spread.... '
      print *,'************************** '
      print *,' '

      ifh_interval = fcsthrs(2) - fcsthrs(1)

      stormloop: do ist = 1,maxstorms

        print *,' '
        print *,'-------------------'
        print *,'Spread calculation '
        print *,'storm number = ',ist
        print *,'storm name   = ',stormarr(ist)
        print *,' '

        timeloop: do ifh = 1,maxtimes

          ifhour = (ifh - 1) * ifh_interval

          print *,' '
          print *,' time index (ifh) = ',ifh,' fhour= '
     &           ,ifhour
          print *,' '

          if (xmeanlat(ist,ifh) > -999.0) then

            ! We have a valid mean for this storm and tau

            xdistsum = 0.0
            imemct   =   0

            write (*,81) 'MEAN:  LAT= ',xmeanlat(ist,ifh)
     &                        ,' LON= ',xmeanlon(ist,ifh)

            do imem = 1,maxmem
              if (traklon(ist,imem,ifh) > 0.0) then
                call calcdist (xmeanlon(ist,ifh),xmeanlat(ist,ifh)
     &              ,traklon(ist,imem,ifh),traklat(ist,imem,ifh)
     &              ,dtr,xdist)
                write (*,83) 'pert:  LAT= ',traklat(ist,imem,ifh)
     &              ,' LON= ',traklon(ist,imem,ifh),'  pert= ',imem
     &              ,'  xdist= ',xdist,' nm'
                xdistsum = xdistsum + xdist
                imemct   = imemct + 1
              endif
            enddo

            if (imemct >= minmem) then
              xspread(ist,ifh) = xdistsum / float(imemct)
              write (*,85) 'SPREAD: ',xspread(ist,ifh),'  nmem= '
     &                    ,imemct,'  ifh= ',ifh
              write (56,95) stormarr(ist),saverec%atx_ymdh,ifhour
     &             ,'SPREAD: ',xspread(ist,ifh),'  nmem= ',imemct
     &             ,'  ifh= ',ifh
            else
              xspread(ist,ifh) = -999.0
              write (*,87) 'SPREAD: NOT ENOUGH MEMBERS, nmem= ',imemct
              write (56,95) stormarr(ist),saverec%atx_ymdh,ifhour
     &                     ,'SPREAD: ',-999.0,'  nmem= ',imemct
     &                     ,'  ifh= ',ifh
            endif
       
          else
 
            write (*,89) 'SPREAD: No mean track for this time'
            write (56,95) stormarr(ist),saverec%atx_ymdh,ifhour
     &                   ,'SPREAD: ',-999.0,'  nmem= ',0
     &                   ,'  ifh= ',ifh
 
          endif

        enddo timeloop

        if (xmeanlat(ist,1) > -999.0) then
          ! We have at least a valid initial time for this storm, so we
          ! should print out the spread array....
          write (65,91) saverec%atx_ymdh,stormarr(ist)
     &                 ,(xspread(ist,ifh),ifh=1,8)
        endif

      enddo stormloop

  81  format (1x,a12,f7.2,2x,a6,f7.2)
  83  format (1x,a12,f7.2,2x,a6,f7.2,2x,a8,i2,2x,a9,f7.2,a3)
  85  format (1x,a8,f7.2,2x,a8,i2,a7,i2)
  87  format (1x,a34,i2)
  89  format (1x,a35)
  91  format (1x,a10,3x,a4,8(3x,f7.2))
  95  format (1x,a4,2x,a10,2x,'fhour= ',i3,3x,a8,f8.2,2x,a8,i2,a7,i2)
      
      return
      end
c
c------------------------------------------------------------------
c
c------------------------------------------------------------------
      subroutine get_prob_mode (numlon,numlat,minlon,maxlon,minlat
     &        ,maxlat,xprob,xmodelon,xmodelat,xmodeval,ifhour
     &        ,cstorm,cymdh,ifh,dthresh,xmeanlon,xmeanlat,ist
     &        ,basinix,catcf,fcsthrs,modeflag)
c
c     ABSTRACT: This subroutine will scan the array of probabilities
c     in order to find the highest probability, or mode.  Since there
c     is a very good chance that you may have bimodal or multimodal
c     cases in which there is more than one probability mode, we need
c     a way to select only one mode.  A barnes analysis, used for 
c     conducting the search, effectively applies a gaussian filter
c     to provide a sort of area-averaged value.  This will help to
c     highlight one area among any other similar areas.
c
c     re      input e-folding radius for barnes analysis
c     ri      input influence radius for searching for min/max

      USE maxparms; USE trig_vals; USE grid_bounds

      real        xprob(numlon,numlat)
      real        rlon(numlon),rlat(numlat)
      real        xmeanlon(maxstorms,maxtimes)
      real        xmeanlat(maxstorms,maxtimes)
      real        re,ri,blat,blon,bval,ctlon,ctlat,dell
      real        tlat,tlon,fixlon,fixlat,dthresh
      real        xmodelon,xmodelat,xmodeval
      integer     fcsthrs(maxtimes)
      integer     numlon,numlat,minlon,minlat,maxlon,maxlat,basinix
     
      character*1 modeflag
      character*4 cstorm,catcf
      character(*) cymdh

c      re   = 150.0
c      ri   = 300.0
      re   = 120.0
      ri   = 240.0
      fmax = -1.0e+10
      dell = grdspc
      nhalf = 3

c     First and foremost we need to check to see if our storm is out
c     of our grid bounds that we have prescribed for this probability
c     mode analysis.  If in fact it is outside those bounds, then
c     quit this analysis and return (this should only happen on very
c     rare occasions, such as if TPC keeps a storm under their control
c     that passes into the central Pacific instead of handing it off
c     to CPHC).

      if (xmeanlon(ist,ifh) > 0.0                 .and.
     &    xmeanlon(ist,ifh) >= xminlon(basinix)   .and.
     &    xmeanlon(ist,ifh) <= xmaxlon(basinix)   .and.
     &    xmeanlat(ist,ifh) >= xminlat(basinix)   .and.
     &    xmeanlat(ist,ifh) <= xmaxlat(basinix)) then
        continue
      else
        print *,' '
        print *,'!!! WARNING: in get_prob_mode, the current storm'
        print *,'!!! mean position is outside the bounds for this'
        print *,'!!! specified ocean domain.  Therefore, we cannot'
        print *,'!!! calculate the mode at this hour.'
        print *,'!!! RETURNING.....'
        print *,' '
        print *,'    Storm info:'
        print *,'       Storm index=         ',ist
        print *,'       Forecast hour index= ',ifh
        print *,'       Mean longitude=      ',xmeanlon(ist,ifh)
        print *,'       Mean latitude=       ',xmeanlat(ist,ifh)
        print *,' '
        print *,'    Grid info:'
        print *,'       basinix= ',basinix
        print *,'       min longitude= ',xminlon(basinix)
        print *,'       max longitude= ',xmaxlon(basinix)
        print *,'       min latitude= ',xminlat(basinix)
        print *,'       max latitude= ',xmaxlat(basinix)
        print *,' '
        return
      endif


c     First we need to be able to narrow our scan.  We will do this
c     by first scanning the array to find the northernmost, southern-
c     most, easternmost and westernmost extent of the probabilities
c     which are > 0.

      mini =  9999
      minj =  9999
      maxi = -9999
      maxj = -9999

      dmax = -1.0e+10
      do i = 1,numlon
        do j = 1,numlat
          if (xprob(i,j) > dmax) then
            dmax = xprob(i,j)
            id   = i
            jd   = j
          endif
          if (xprob(i,j) > 0.0) then
            if (i < mini) mini = i
            if (i > maxi) maxi = i
            if (j < minj) minj = j
            if (j > maxj) maxj = j
          endif
        enddo
      enddo

      print *,' '
      print *,'dmax:  at beg of get_prob_mode, dmax= ',dmax 
      print *,'i= ',id,' j= ',jd
      print *,'mini= ',mini
      print *,'maxi= ',maxi
      print *,'minj= ',minj
      print *,'maxj= ',maxj

      if (mini ==  9999 .or. minj ==  9999 .or. 
     &    maxi == -9999 .or. maxj == -9999) then
        print *,' '
        print *,'!!! ERROR in get_prob_mode, one or more of the max or'
        print *,'!!! min indeces for the xprob array could not '
        print *,'!!! be found.  The xprob array must = 0.'
        print *,'!!! EXITING.... '
        print *,' '
        STOP 99
      endif

c     Load all the lat and lon values into arrays.  Remember that
c     the xprob array goes from upper left (NW) to lower right (SE).

      do i = 1,numlon
        rlon(i) = float(minlon + i - 1)
      enddo

      do j = numlat,1,-1
        jix = numlat - j + 1
        rlat(jix) = float(maxlat - jix + 1)
      enddo

c     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c              First pass through barnes analysis....
c                   
c     For the first round through the  barnes analysis, we use the 
c     full grid resolution, and at each grid point, we look at an
c     area roughly 2x(ri).

      npts = ceiling(ri/(dtk*grdspc)) * 2

      ibct=0
      ibarnes_loopct = 0

c      print *,' '
c      print *,'before first barnes, mini= ',mini,' maxi= ',maxi
c      print *,'before first barnes, minj= ',minj,' maxj= ',maxj
c      print *,'npts= ',npts

      do i = mini,maxi
        do j = minj,maxj

          blon = rlon(i)        
          blat = rlat(j)        

c          write (*,*) ' '
c          write (*,*) ' - - - - - - - - - - - -'
c          write (*,*) ' '
c          write (*,59) cstorm,ifhour
c          write (*,61) 'before 1st barnes, i= ',i,j,blon,360.-blon
c     &                 ,blat,-99.0

          jjbeg = j - npts
          if (jjbeg < 1) jjbeg = 1
          jjend = j + npts
          if (jjend > numlat) jjend = numlat
          iibeg = i - npts
          if (iibeg < 1) iibeg = 1
          iiend = i + npts
          if (iiend > numlon) iiend = numlon

          ibct = ibct + 1
          call barnes (blon,blat,rlon,rlat,numlon,numlat,iibeg,jjbeg
     &                ,iiend,jjend,xprob,re,ri,bval,icount)

c          write (*,61) 'after  1st barnes, i= ',i,j,blon,360.-blon
c     &                 ,blat,bval

          ibarnes_loopct = ibarnes_loopct + icount

          if (bval > fmax) then
            fmax  = bval
            ctlon = blon
            ctlat = blat
          endif

        enddo
      enddo

  59  format (1x,'Storm= ',a4,'  forecast hour= ',i4)
  61  format (1x,a22,i3,' j= ',i3,' blon= ',f7.2,'E (',f7.2,'W) blat= '
     &       ,f7.2,' bval= ',f10.6)

      print *,' '
      print *,'++++ FIRST BARNES LOOP COMPLETED ++++'
      print *,' '
      print *,'After first loop, # calls to barnes = ',ibct
      print *,'Total # of internal barnes loop iterations = '
     &       ,ibarnes_loopct

      write (*,59) cstorm,ifhour
      write (6,63) ctlon,360.-ctlon,ctlat,fmax
  63  format (' After first run, ctlon= ',f8.3,'E (',f8.3,'W) ctlat= '
     &       ,f8.3,' barnval = ',f10.6)
      print *,' '

      if (ctlon <= xminlon(basinix) .or. ctlon >= xmaxlon(basinix) .or.
     &    ctlat <= xminlat(basinix) .or. ctlat >= xmaxlat(basinix)) then
        print *,' '
        print *,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        print *,'!!! EXIT: In get_prob_mode, the mode position that '
        print *,'!!! was found after 1 iteration of the barnes loop is'
        print *,'!!! equal to or greater than 1 of the boundaries, '
        print *,'!!! meaning that we are pushing up against the'
        print *,'!!! boundary with our search and it is possible that'
        print *,'!!! we are not finding the true mode.'
        print *,'!!! '
        print *,'!!! ctlon= ',ctlon,'  ctlat= ',ctlat
        print *,'!!! grid minlon = ',xminlon(basinix)
        print *,'!!! grid maxlon = ',xmaxlon(basinix)
        print *,'!!! grid minlat = ',xminlat(basinix)
        print *,'!!! grid maxlat = ',xmaxlat(basinix)
        print *,'!!! '
        print *,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        modeflag = 'n'
        return
      endif

c     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c              Second pass through barnes analysis....
c                   
c     Now that we've gotten a first estimate of the position, we 
c     go through this 3 more times, cutting the grid resolution
c     in half each time through in order to sharpen the estimate.
c     Cut the number of points (npts) in half from what we looked
c     at in the first barnes iteration so that our first position 
c     estimate is restricted from moving too much.

      npts = ceiling(ri/(dtk*grdspc))

      ibct = 0
      ibarnes_loopct = 0
      kloop: do k = 1,nhalf

        dell = 0.5*dell
        tlon = ctlon
        tlat = ctlat
        fmax = -1.0e+15

        print *,' '
        write (*,67) k,ctlon,360.-ctlon,ctlat
 67     format (1x,'In second barnes loop, k= ',i3,' ctlon= ',f8.3
     &         ,'E (',f8.3,'W) ctlat= ',f8.3)
        print *,' '

        jloop: do j = -npts,npts

          blat = tlat + dell*float(j)
         
          if (blat > float(maxlat) .or. blat < float(minlat)) then
            print *,' '
            print *,'!!! WARNING: In get_prob_mode 2nd barnes loop for'
            print *,'!!! nhalf = ',nhalf,', blat is outside of the '
            print *,'!!! grid lat bounds.  blat= ',blat
            print *,'!!! maxlat= ',maxlat,' minlat= ',minlat
            cycle jloop
          endif

          iloop: do i = -npts,npts

            blon = tlon + dell*float(i)

            if (blon > float(maxlon) .or. blon < float(minlon)) then
              print *,' '
              print *,'!!! WARNING: In get_prob_mode 2nd barnes loop' 
              print *,'!!! for nhalf = ',nhalf,', blon is outside of'
              print *,'!!! the grid lon bounds.  blon= ',blon
              print *,'!!! maxlon= ',maxlon,' minlon= ',minlon
              cycle iloop
            endif
      
            ifix = int(blon - float(minlon)) + 1
            jfix = int(float(maxlat) - blat) + 1

c            write (*,*) ' '
c            write (*,*) ' - - - - - - - - - - - -'
c            write (*,*) ' '
c            write (*,59) cstorm,ifhour
c            write (*,71) 'before 2nd barnes, ifix= ',ifix,jfix,k,blon
c     &                   ,360.-blon,blat,-99.0

            ! We cut npts in half so that the center position is
            ! restricted from moving too much in this refinement,
            ! but we need to double it here to make sure we 
            ! include enough points in the  barnes analysis.

            jjbeg = jfix - (npts*2)
            if (jjbeg < 1) jjbeg = 1
            jjend = jfix + (npts*2)
            if (jjend > numlat) jjend = numlat
            iibeg = ifix - (npts*2)
            if (iibeg < 1) iibeg = 1
            iiend = ifix + (npts*2)
            if (iiend > numlon) iiend = numlon

            ibct = ibct + 1
            call barnes (blon,blat,rlon,rlat,numlon,numlat,iibeg,jjbeg
     &                  ,iiend,jjend,xprob,re,ri,bval,icount)

c            write (*,71) 'after  2nd barnes, ifix= ',ifix,jfix,k,blon
c     &                   ,360.-blon,blat,bval

            ibarnes_loopct = ibarnes_loopct + icount

            if (bval > fmax) then
              fmax  = bval
              ctlon = blon
              ctlat = blat
            endif

          enddo iloop

        enddo jloop

        write (6,73) k,ctlon,360.-ctlon,ctlat,fmax

        if (ctlon <= xminlon(basinix) .or. ctlon >= xmaxlon(basinix) 
     &      .or. ctlat <= xminlat(basinix) 
     &      .or. ctlat >= xmaxlat(basinix)) then
          print *,' '
          print *,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          print *,'!!! EXIT: In get_prob_mode, the mode position that '
          print *,'!!! was found after iteration nhalf= ',k,' of the '
          print *,'!!! barnes loop is equal to or greater than 1 of the'
          print *,'!!! boundaries, meaning that we are pushing up '
          print *,'!!! against the boundary with our search and it is'
          print *,'!!! possible that we are not finding the true mode.'
          print *,'!!! '
          print *,'!!! ctlon= ',ctlon,'  ctlat= ',ctlat
          print *,'!!! grid minlon = ',xminlon(basinix)
          print *,'!!! grid maxlon = ',xmaxlon(basinix)
          print *,'!!! grid minlat = ',xminlat(basinix)
          print *,'!!! grid maxlat = ',xmaxlat(basinix)
          print *,'!!! '
          print *,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          modeflag = 'n'
          return
        endif

      enddo kloop

 71   format (1x,a25,i3,' jfix= ',i3,' k= ',i2,' blon= ',f7.2,'E ('
     &       ,f7.2,'W) blat= ',f7.2,' bval= ',f10.6)

 73   format (' nhalf barnes, k= ',i2,' ctlon= ',f8.3,'E (',f8.3
     &       ,'W)  ctlat= ',f8.3,' barnval = ',f10.6)

      print *,' '
      print *,'++++ SECOND BARNES LOOP COMPLETED ++++'
      print *,' '
      print *,'After second loop, # calls to barnes = ',ibct
      print *,'Total # of internal barnes loop iterations = '
     &       ,ibarnes_loopct

      write (*,59) cstorm,ifhour
      write (6,77) ctlon,360.-ctlon,ctlat,fmax
 77   format (' After second run, ctlon= ',f8.3,'E (',f8.3,'W) ctlat= '
     &       ,f8.3,' barnval = ',f10.6)
      print *,' '
      print *,'+ + + + + + + + + + + + + + + + + + + '

      xmodelon = ctlon
      xmodelat = ctlat
      bmaxval  = fmax


c     We have now gotten the position of the max probability, or
c     the mode, using the  barnes analysis.  This barnes analysis,
c     however, gives us essentially an area-averaged value.  But
c     what we really need is the actual raw probability value 
c     calculated previously in this program.  So we will look at
c     the eight gridpoints surrounding the lat/lon value obtained
c     by the  barnes analysis and pick out the max raw value.
c     For the  output atcfunix mode record, however, we will use
c     the lat/lon position obtained through the barnes analysis.

      print *,' '
      print *,' Getting (i,j) location of max prob mode....'
      print *,' numlon= ',numlon,'  numlat= ',numlat
      print *,' ifix= ',ifix,' jfix= ',jfix
      print *,' '

      ifix = int(ctlon - float(minlon)) + 1
      jfix = int(float(maxlat) - ctlat) + 1

      fmax = -1.0e+15

      if (xprob(ifix,jfix) > fmax) then
        ihold = ifix
        jhold = jfix
        fmax  = xprob(ifix,jfix)
      endif

      if (ifix < numlon) then
        ! Check the point just *east* of our 
        ! found mode position....
        if (xprob(ifix+1,jfix) > fmax) then
          ihold = ifix+1
          jhold = jfix
          fmax  = xprob(ifix+1,jfix)
        endif
      endif

      if (ifix > 1) then
        ! Check the point just *west* of our
        ! found mode position....
        if (xprob(ifix-1,jfix) > fmax) then
          ihold = ifix-1        
          jhold = jfix 
          fmax  = xprob(ifix-1,jfix)
        endif
      endif

      if (jfix < numlat) then
        ! Check the point just *south* of our 
        ! found mode position (remember: j=1 is
        ! the northernmost point....)
        if (xprob(ifix,jfix+1) > fmax) then
          ihold = ifix
          jhold = jfix+1
          fmax  = xprob(ifix,jfix+1)
        endif
      endif

      if (jfix > 1) then
        ! Check the point just *north* of our
        ! found mode position (remember: j=1 is
        ! the northernmost point....)
        if (xprob(ifix,jfix-1) > fmax) then
          ihold = ifix
          jhold = jfix-1
          fmax  = xprob(ifix,jfix-1)
        endif
      endif

      if (ifix < numlon .and. jfix > 1) then
        ! Check the point just *northeast* of our
        ! found mode position (remember: j=1 is
        ! the northernmost point....)
        if (xprob(ifix+1,jfix-1) > fmax) then
          ihold = ifix+1
          jhold = jfix-1
          fmax  = xprob(ifix+1,jfix-1)
        endif
      endif

      if (ifix < numlon .and. jfix < numlat) then
        ! Check the point just *southeast* of our
        ! found mode position (remember: j=1 is
        ! the northernmost point....)   
        if (xprob(ifix+1,jfix+1) > fmax) then
          ihold = ifix+1   
          jhold = jfix+1
          fmax  = xprob(ifix+1,jfix+1)
        endif
      endif

      if (ifix > 1 .and. jfix < numlat) then
        ! Check the point just *southwest* of our
        ! found mode position (remember: j=1 is
        ! the northernmost point....)   
        if (xprob(ifix-1,jfix+1) > fmax) then
          ihold = ifix-1   
          jhold = jfix+1
          fmax  = xprob(ifix-1,jfix+1)
        endif
      endif
        
      if (ifix > 1 .and. jfix > 1) then
        ! Check the point just *northwest* of our
        ! found mode position (remember: j=1 is
        ! the northernmost point....)
        if (xprob(ifix-1,jfix-1) > fmax) then
          ihold = ifix-1
          jhold = jfix-1
          fmax  = xprob(ifix-1,jfix-1)
        endif
      endif

      fixlon   = float(minlon + ihold - 1)
      fixlat   = float(maxlat - jhold + 1)
      xmodeval = fmax

      print *,' '
      print *,' Location (i,j) of max mode: '
      print *,'   i-location: ',ihold,'  fixlon: ',fixlon
      print *,'   j-location: ',jhold,'  fixlat: ',fixlat
      print *,'   fmax = xmodeval = ',xmodeval
      print *,' '

      print *,' '
      print *,'At end of get_prob_mode, '
      write (*,59) cstorm,ifhour
      write (*,81) ctlat,xmodelat
      write (*,83) ctlon,360.-ctlon,xmodelon,360.-xmodelon
      write (*,85) bmaxval,xmodeval
 81   format (1x,'Barnes lat: ',f7.2,'              fixlat: ',f7.2)
 83   format (1x,'Barnes lon: ',f7.2,'E (',f7.2,'W)  fixlon: '
     &       ,f7.2,'E (',f7.2,'W)')
 85   format (1x,'Barnes val: ',f10.6,'  fixval: ',f10.6)
      print *,'ifix= ',ihold,' jfix= ',jhold


c     Now make a call to output_atcfunix to create an output
c     atcfunix record for this mode value....

      call output_atcfunix (cstorm,cymdh,xmodelat,xmodelon
     &               ,ifh,'mode',xmodeval,dthresh,catcf,fcsthrs)

c
      return
      end

c
c----------------------------------------------------------------
c
c----------------------------------------------------------------
      subroutine barnes(flon,flat,rlon,rlat,imax,jmax,iibeg,jjbeg
     &                 ,iiend,jjend,fxy,re,ri,favg,icount)
c
c     ABSTRACT: This routine performs a single-pass barnes anaylsis
c     of fxy at the point (flon,flat). The e-folding radius (km)
c     and influence radius (km) are re and ri, respectively.
c
c     NOTE:  The input grid that is searched in this subroutine is most
c     likely NOT the model's full, original grid.  Instead, a smaller
c     subgrid of the original grid is searched, which is determined by
c     the npts value passed into this subroutine. 
c
c     INPUT:
c     flon    Lon value for center point about which barnes anl is done
c     flat    Lat value for center point about which barnes anl is done
c     rlon    Array of lon values for each grid point
c     rlat    Array of lat values for each grid point
c     imax    Max number of pts in x-direction on input grid
c     jmax    Max number of pts in y-direction on input grid
c     iibeg   i index for grid point to start barnes anlysis (upp left)
c     jjbeg   j index for grid point to start barnes anlysis (upp left)
c     iiend   i index for last grid point in barnes anlysis (low right)
c     jjend   j index for last grid point in barnes anlysis (low right)
c     fxy     Real array of data on which to perform barnes analysis
c     re      input e-folding radius for barnes analysis
c     ri      input influence radius for searching for min/max
c
c     OUTPUT:
c     favg    Average value about the point (flon,flat)
c
      real      fxy(imax,jmax), rlon(imax), rlat(jmax)
      real      flon,flat,favg,res,re,ri,wts,dist,wt

c     --------------------------

      res = re*re
      wts = 0.0
      favg = 0.0

      icount = 0

c      write (*,62) iibeg,iiend,jjbeg,jjend
c 62   format (1x,'in barnes, iibeg= ',i3,' iiend= ',i3,' jjbeg= ',i3
c     &       ,' jjend= ',i3)

      do j=jjbeg,jjend
        do i=iibeg,iiend

          icount = icount + 1

          ! Call a different version of calcdist in order to
          ! get the distance in km....

          call calcdistkm(flon,flat,rlon(i),rlat(j),dist)

c          write (*,65) i,j,flon,flat,rlon(i),rlat(j),dist

          if (dist .gt. ri) cycle

c          write (*,66) i,j,flon,flat,rlon(i),rlat(j),dist,fxy(i,j)
c          print *,'+++ dist= ',dist,' fxy= ',fxy(i,j)

          wt   = exp(-1.0*dist*dist/res)
          wts  = wts + wt
          favg = favg + wt*fxy(i,j)

        enddo
      enddo

c 65   format (1x,'i= ',i3,' j= ',i3,' flon= ',f6.2,' flat= ',f5.2
c     &       ,' rlon(i)= ',f6.2,' rlat(j)= ',f5.2,' dist= ',f7.2)
c 66   format (1x,'i= ',i3,' j= ',i3,' flon= ',f6.2,' flat= ',f5.2
c     &       ,' rlon(i)= ',f6.2,' rlat(j)= ',f5.2,' dist= ',f7.2
c     &       ,' fxy= ',f5.2)
cJ.Peng_04-22-2013  F4.2 to f5.2 the above line
c      write (*,75) 'in barnes before averaging, wts = ',wts,'  favg= '
c     &             ,favg

      if (wts .gt. 1.0E-5) then
         favg = favg/wts
      else
         favg = 0.0
      endif
      iret = 0

c      write (*,75) 'in barnes after  averaging, wts = ',wts,'  favg= '
c     &             ,favg

c 75   format (1x,a34,f9.4,a8,f10.6)

c
      return
      end

c
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
      subroutine calcdistkm(rlonb,rlatb,rlonc,rlatc,xdist)
c
c     ABSTRACT: This subroutine computes the distance between two
c               lat/lon points by using spherical coordinates to
c               calculate the great circle distance between the points.
c                       Figure out the angle (a) between pt.B and pt.C,
c             N. Pole   then figure out how much of a % of a great
c               x       circle distance that angle represents.
c              / \
c            b/   \     cos(a) = (cos b)(cos c) + (sin b)(sin c)(cos A)
c            /     \
c        pt./<--A-->\c     NOTE: The latitude arguments passed to the
c        B /         \           subr are the actual lat vals, but in
c                     \          the calculation we use 90-lat.
c               a      \
c                       \pt.  NOTE: You may get strange results if you:
c                         C    (1) use positive values for SH lats AND
c                              you try computing distances across the
c                              equator, or (2) use lon values of 0 to
c                              -180 for WH lons AND you try computing
c                              distances across the 180E meridian.
c
c     NOTE: In the diagram above, (a) is the angle between pt. B and
c     pt. C (with pt. x as the vertex), and (A) is the difference in
c     longitude (in degrees, absolute value) between pt. B and pt. C.
c
c     !!! NOTE !!! -- THE PARAMETER ecircum IS DEFINED (AS OF THE
c     ORIGINAL WRITING OF THIS SYSTEM) IN KM, NOT M, SO BE AWARE THAT
c     THE DISTANCE RETURNED FROM THIS SUBROUTINE IS ALSO IN KM.
c
      USE trig_vals
c
      if (rlatb < 0.0 .or. rlatc < 0.0) then
        pole = -90.
      else
        pole = 90.
      endif
c
      distlatb = (pole - rlatb) * dtr
      distlatc = (pole - rlatc) * dtr
      difflon  = abs( (rlonb - rlonc)*dtr )
c
      cosanga = ( cos(distlatb) * cos(distlatc) +
     &            sin(distlatb) * sin(distlatc) * cos(difflon))

c     This next check of cosanga is needed since I have had ACOS crash
c     when calculating the distance between 2 identical points (should
c     = 0), but the input for ACOS was just slightly over 1
c     (e.g., 1.00000000007), due to (I'm guessing) rounding errors.

      if (cosanga > 1.0) then
        cosanga = 1.0
      endif

      degrees    = acos(cosanga) / dtr
      circ_fract = degrees / 360.
      xdist      = circ_fract * ecircum

c     NOTE: whether this subroutine returns the value of the distance
c           in km or m depends on the scale of the parameter ecircum.
c           At the original writing of this subroutine (7/97), ecircum
c           was given in km.

      return
      end

c
c-----------------------------------------------------------------------
c
cJ.Peng-04-15-2013  --------------------------------------------
      subroutine read_atcfunix (traklat,traklon,tcvmx,tcpce,
     &  saverec,stormarr,maxmem,fcst_interval,perts)

c
c     ABSTRACT: This subroutine reads in the atcfunix records from
c     the concatenated file that contains all these records for all 
c     storms, all forecast hours and all ensemble members.

      USE define_atcfunix_rec; USE maxparms; USE ensinfo
      USE grid_bounds
c
      type (atcfunix_card) inprec,saverec

      real    traklat(maxstorms,maxmem,maxtimes)
      real    traklon(maxstorms,maxmem,maxtimes)

cJ.Peng-04-15-2013
      real tcvmx(maxstorms,maxmem,maxtimes)
      real tcpce(maxstorms,maxmem,maxtimes)

      integer maxmem,fcst_interval
      character*4   perts(maxmem)
      character*4   atcfid,stormarr(maxstorms)
      character*1   foundstorm,foundmember,foundmodel
c
      stormarr = 'XXXX'

      ict = 0
      istot = 0
      readloop: do while (.true.)

        read (11,85,end=870) inprec

        if (ict == 0) then
          saverec = inprec
          ict = ict + 1
        endif

        atcfid(1:2) = inprec%atx_basin
        atcfid(3:4) = inprec%atx_storm_num

c       Find the array index for the storm that is being read in.  As 
c       we encounter a new storm in the atcfunix file, we add that 
c       storm to the array and increment the index (istix) by 1.

        foundstorm = 'n'
        stormloop: do istorm = 1,maxstorms
          if (atcfid == stormarr(istorm)) then
            foundstorm = 'y'
            istix = istorm
            exit stormloop
          endif
        enddo stormloop
        if (foundstorm == 'n') then
          istot = istot + 1
          istix = istot
          stormarr(istix) = atcfid
        endif

c       Similarly, find the lat/lon array index for the ensemble member
c       that is being read in.  The members are specified in a module.
c       If the program cannot find the member, chances are the ensemble
c       suite has been updated and new members have been added.

        foundmember = 'n'
        memberloop: do imem = 1,maxmem
          if (inprec%atx_model == perts(imem)) then
            foundmodel = 'y'
            imix = imem
            exit memberloop
          endif
        enddo memberloop
        if (foundmodel == 'n') then
          print *,' '
          print *,'!!! ERROR: IN PROGRAM TRAKAVE, A MODEL/MEMBER NAME'
          print *,'!!!        IS NOT RECOGNIZED.'
          print *,'!!!       '
          print *,'!!!        MODEL/MEMBER NAME = --->'
     &                              ,inprec%atx_model,'<----'
          print *,'!!!       '
          print *,'!!! STOPPING PROGRAM....      '
          print *,'!!!       '
          stop 990
        endif 

c       Calculate the index (ifix) for the forecast hour in lat and
c       lon arrays.

        ifix = (inprec%atx_fcsthr / fcst_interval) + 1

c       Check the lat/lon arrays to see if we've already read in the
c       lat/lon positions for this storm/member/forecast_time (traklat
c       and traklon are initialized to -9999.0).  Yes, it is possible
c       to have more than 1 record for a given storm/member/time, 
c       since there are up to 3 different records written out for each
c       forecast time, each containing the 34-, 50- and 64-kt wind 
c       radii.  Note that if the values read in are zeroes, that means
c       that the tracker could not find the storm; consider these 
c       missing values and ignore them. 

        if (traklat(istix,imix,ifix) < -9998.0 .and. 
     &      traklon(istix,imix,ifix) < -9998.0) then

          if (inprec%atx_lat > 0 .and. inprec%atx_lon > 0) then

            ! Float the integer latitudes and convert any positive 
            ! southern latitudes to be negative.

            if (inprec%atx_latns == 'N') then
              traklat(istix,imix,ifix) = float(inprec%atx_lat)/10.0
            else
              traklat(istix,imix,ifix) = float(inprec%atx_lat)/(-10.0)
            endif

            ! Float the integer longitudes and convert any "western" 
            ! longitudes to be "eastern", 0-360.

            if (inprec%atx_lonew == 'E') then
              traklon(istix,imix,ifix) = float(inprec%atx_lon)/10.0
            else
              traklon(istix,imix,ifix) = float(3600 - inprec%atx_lon)
     &                                        / 10.0
            endif

cJ.Peng-04-15-2013
	    tcvmx(istix,imix,ifix) = float(inprec%atx_vmax)
	    tcpce(istix,imix,ifix) = float(inprec%atx_pcen)

          endif
  
        endif

      enddo readloop

C      close (11)

   85 format (a2,2x,a2,2x,a10,2x,a2,2x,a4,2x,i3,2x,i3,a1,2x,i4,a1,2x
     &       ,i3,2x,i4)
c     &       ,i3,2x,i4,2x,a2,2x,i2,2x,a3,4(2x,i3))

  870 goto 1090

 1010 print *,' '
      print *,'!!! ERROR opening input file.  ERROR = ',ios
      print *,' '

 1090 continue

      return
      end

c
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
      subroutine calcdist(rlonb,rlatb,rlonc,rlatc,dtr,xdist)
c
c     ABSTRACT: This subroutine computes the distance between two
c               lat/lon points by using spherical coordinates to
c               calculate the great circle distance between the points.
c                       Figure out the angle (a) between pt.B and pt.C,
c             N. Pole   then figure out how much of a % of a great
c               x       circle distance that angle represents.
c              / \
c            b/   \     cos(a) = (cos b)(cos c) + (sin b)(sin c)(cos A)
c            /     \
c        pt./<--A-->\c     NOTE: The latitude arguments passed to the
c        B /         \           subr are the actual lat vals, but in
c                     \          the calculation we use 90-lat.
c               a      \
c                       \pt.  NOTE: You may get strange results if you:
c                         C    (1) use positive values for SH lats AND
c                              you try computing distances across the
c                              equator, or (2) use lon values of 0 to
c                              -180 for WH lons AND you try computing
c                              distances across the 180E meridian.
c
c     NOTE: In the diagram above, (a) is the angle between pt. B and
c     pt. C (with pt. x as the vertex), and (A) is the difference in
c     longitude (in degrees, absolute value) between pt. B and pt. C.
c
c     !!! NOTE !!! -- THE PARAMETER ecircum IS DEFINED (AS OF THE
c     WRITING OF THIS SYSTEM) IN NM, NOT M or KM, SO BE AWARE THAT
c     THE DISTANCE RETURNED FROM THIS SUBROUTINE IS ALSO IN NM.
c
      real :: dtk = 60.005       ! Dist (nm) over 1 deg lat
                                 ! using erad=6371.0e+3 km
      real :: erad = 6371.0e+3   ! Earth's radius (m)
      real :: ecircum = 21601.8  ! Earth's circumference
                                 ! (nm) using erad=6371.e3
      real    dtr
c
      if (rlatb < 0.0 .or. rlatc < 0.0) then
        pole = -90.
      else
        pole = 90.
      endif
c
      distlatb = (pole - rlatb) * dtr
      distlatc = (pole - rlatc) * dtr
      difflon  = abs( (rlonb - rlonc)*dtr )
c
      cosanga = ( cos(distlatb) * cos(distlatc) +
     &            sin(distlatb) * sin(distlatc) * cos(difflon))

c     This next check of cosanga is needed since I have had ACOS crash
c     when calculating the distance between 2 identical points (should
c     = 0), but the input for ACOS was just slightly over 1
c     (e.g., 1.00000000007), due to (I'm guessing) rounding errors.

      if (cosanga > 1.0) then
        cosanga = 1.0
      endif

      degrees    = acos(cosanga) / dtr
      circ_fract = degrees / 360.
      xdist      = circ_fract * ecircum

c     NOTE: whether this subroutine returns the value of the distance in
c           km, m or nm depends on the scale of the parameter ecircum.
c           At the original writing of this subroutine (2/01), ecircum
c           was given in nm.

      return
      end

c
c---------------------------------------------------------------------
c
c---------------------------------------------------------------------
      subroutine output_atcfunix (stormid,cymdh,outlat,outlon,ifh
     &                           ,ctype,xprob,dthresh,catcf,fcsthrs)
c
c     ABSTRACT: This subroutine  outputs a 1-line message for a given
c     storm at an input forecast hour in the new ATCF UNIX format.
c     Unlike the old atcf DOS format in which you waited until the
c     whole tracking was over to write the  output for all forecast
c     hours, with this atcfunix format, each time we are calling this
c     subroutine, it is to only write out 1 record, which will be the
c     fix info for a particular storm at a given time.  
c
c     While this new atcfunix format contains much more information than
c     the old 1-line atcf dos message, for the purposes of the ensemble
c     mean track, we will only use the slots for storm location and NOT
c     worry about the intensity parameters.  An example is shown below:
c
c     AL, 13, 2000092500, 03, AEMN, 036, 243N, 675W
c
c     NOTE: The longitudes that are passed into this subroutine are
c     given in 0 - 360, increasing eastward.  The format for the
c     atcfunix system requires that the  output be 0-180E or
c     0-180W, so we must adjust the values, if needed.  Also, the
c     values for southern latitudes must be positive (use 'N' and
c     'S' to distinguish Northern/Southern Hemispheres).
c
c     INPUT:
c     stormid   Storm ID (EP02, AL13, etc...)
c     cymdh     initial date/time of forecast (yyyymmddhh)
c     outlat    latitude fix position for this storm at this time
c               which is to be written out to the  output file
c     outlon    longitude fix position for this storm at this time
c               which is to be written out to the  output file
c     ifh       the index for the current forecast hour being output
c     ctype     'mean' for ens mean, 'mode' for ens prob mode
c     xprob     prob mode value for ctype = 'mode'. Not used for mean
c     dthresh   distance threshold used for prob mode value, to be 
c               used in model output name for 'mode' case.
c
c     LOCAL:
c     intlon    integer that holds the value of outlon*10
c     intlat    integer that holds the value of outlat*10
c
      USE ensinfo; USE maxparms

      real    outlon,outlat,xprob
      integer fcsthrs(maxtimes)
      integer intlon,intlat
      character*1  clatns,clonew
      character*2  basinid,stormnum
      character*4  stormid,ctype,cmode,catcf
      character*10 cymdh
      character*37 outatxfile

c     First convert all of the lat/lon values from reals into integers.
c     These integer values must be 10x their real value (eg. 125.4 will
c     be written out as 1254).  Convert the lon values so that they go
c     from 0-180E or 0-180W, and convert the lat values so that they are
c     positive and use 'N' or 'S' to differentiate hemispheres.

C      outatxfile = 'aemn.'//cymdh//'.cyclone.trackatcfunix'
C      open (unit=53,file=outatxfile,status='unknown',position='append')

      if (outlon < -998.0 .or. outlat < -998.0) then
        intlon = 0
        intlat = 0
        clonew = ' '
        clatns = ' '
      else
        if (outlon >= 180.0) then
c          intlon = 3600 - int(outlon * 10. + 0.5)
c          clonew = 'W'
cJ.Peng-----03/27/2015------------------------------------
          if (outlon < 360.0) then
            intlon = 3600 - int(outlon * 10. + 0.5)
            clonew = 'W'
          else
            intlon =int(outlon * 10. + 0.5) - 3600
            clonew = 'E'
          endif
cJ.Peng-----03/27/2015------------------------------------

        else
          intlon = int(outlon * 10. + 0.5)
          clonew = 'E'
        endif
        intlat = int(abs(outlat) * 10. + 0.5)
        if (outlat < 0.0) then
          clatns = 'S'
        else
          clatns = 'N'
        endif
      endif

      basinid = stormid(1:2)
      stormnum = stormid(3:4)

      if (ctype == 'mean') then
        write (53,81) basinid,stormnum,cymdh,catcf,fcsthrs(ifh)
     &        ,intlat,clatns,intlon,clonew
      else
        itemp = int(dthresh)
        cmode(1:1) = catcf(1:1)
        write (cmode(2:4),'(i3.3)') itemp
c        print *,'atcfunix, cmode= ',cmode,' itemp= ',itemp
c        write (*,79) cmode,itemp
   79   format (1x,' cmode= ....',a4,'....  itemp= ',i3.3)
        write (57,83) basinid,stormnum,cymdh,cmode,fcsthrs(ifh)
     &        ,intlat,clatns,intlon,clonew,int(xprob*100.0)
      endif

   81 format (a2,', ',a2,', ',a10,', 03, ',a4,', ',i3.3,', ',i3,a1
     &       ,', ',i4,a1)
   83 format (a2,', ',a2,', ',a10,', 03, ',a4,', ',i3.3,', ',i3,a1
     &       ,', ',i4,a1,', ',i3)
c
C      close (53)
      return
      end

c
c-----------------------------------------------------------------------
      subroutine output_atcf (stormid,xmeanlon,xmeanlat,cymdh,ist,catcf)
c
c     ABSTRACT: This subroutine  outputs a 1-line message for each storm
c     in ATCF format.  This message contains the model identifier, the
c     forecast initial date, and the positions for 12, 24, 36, 48 and
c     72 hours.  This message also contains the intensity
c     estimates (in knots) for those same hours.  The  conversion for
c     m/s to knots is to multiply m/s by 1.9427 (3.281 ft/m,
c     1 naut mile/6080 ft, 3600s/h).
c
c     NOTE: Intensity estimates are NOT included for the ensemble 
c     mean tracks.
c
c     NOTE: The longitudes that are passed into this subroutine are
c     given in 0 - 360, increasing eastward.  The  output of this
c     subroutine is used by the atcf system at TPC for plotting
c     purposes, and the atcf plotting routines need the longitudes in
c     0 - 360, increasing westward.  Thus, a necessary adjustment is
c     made.

      USE maxparms

      character*4  stormid,catcf
      character*10 cymdh
      character*8  yymmddhh
      character*33 outatcffile

      real    xmeanlat(maxstorms,maxtimes)
      real    xmeanlon(maxstorms,maxtimes)
      integer intlon(maxtimes),intlat(maxtimes)
c
C      outatcffile = 'aemn.'//cymdh//'.cyclone.trackatcf'
C      open (unit=54,file=outatcffile,status='unknown',position='append')

      yymmddhh = cymdh(3:10)

      ifhloop: do ifh = 1,maxtimes

        if (xmeanlon(ist,ifh) < -998.0 .or. xmeanlat(ist,ifh) < -998.0)
     &  then

          intlon(ifh) = 0
          intlat(ifh) = 0

        else
          intlon(ifh) = 3600 - int(xmeanlon(ist,ifh) * 10. + 0.5)
          intlat(ifh) = int(abs(xmeanlat(ist,ifh)) * 10. + 0.5)
          if (xmeanlat(ist,ifh) < 0.0) then
            intlat(ifh) = intlat(ifh) * (-1)
          endif

        endif

      enddo ifhloop

      write (54,82) '99',catcf,yymmddhh,intlat(3),intlon(3)
     &           ,intlat(5),intlon(5),intlat(7),intlon(7)
     &           ,intlat(9),intlon(9),intlat(13),intlon(13)
     &           ,0,0,0,0,0,stormid,cymdh(9:10)

   82 format (a2,a4,a8,10i4,5i3,1x,a4,a2)
c

C      close (54)
      return
      end

c
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
      subroutine output_all (stormid,xmeanlon,xmeanlat,cymdh,ist,catcf)
c
c     ABSTRACT: This subroutine  outputs a 1-line message for each
c     storm.  This message contains the model identifier, the forecast
c     initial date, and the positions for 0, 12, 24, 36, 48, 60 and 72
c     hours.  In the case of the regional models (NGM, Eta), which
c     only go out to 48h, zeroes are included for forecast hours
c     60 and 72.
c
c     NOTE: The longitudes that are passed into this subroutine are
c     given in 0 - 360, increasing eastward.  The  output of this
c     subroutine is used by Steve Lord for plotting purposes, and his
c     plotting routines need the longitudes in 0 - 360, increasing
c     westward.  Thus, a necessary adjustment is made.
c
      USE maxparms

      character*4  stormid,catcf
      character*10 cymdh
      character*8  yymmddhh
      character*2  stormnum
      character*3  basinid
      character*32 outallfile

      real    xmeanlat(maxstorms,maxtimes)
      real    xmeanlon(maxstorms,maxtimes)
      integer intlon(maxtimes),intlat(maxtimes)
c
C      outallfile = 'aemn.'//cymdh//'.cyclone.trackall'
C      open (unit=55,file=outallfile,status='unknown',position='append')

      yymmddhh = cymdh(3:10)
      stormnum = stormid(3:4)
c
      ifhloop: do ifh = 1,maxtimes

        if (xmeanlon(ist,ifh) < -998.0 .or. xmeanlat(ist,ifh) < -998.0)
     &  then
          intlon(ifh) = 0
          intlat(ifh) = 0
        else
          intlon(ifh) = 3600 - int(xmeanlon(ist,ifh) * 10. + 0.5)
          intlat(ifh) = int(abs(xmeanlat(ist,ifh)) * 10. + 0.5)
          if (xmeanlat(ist,ifh) < 0.0) then
            intlat(ifh) = intlat(ifh) * (-1)
          endif
        endif

      enddo ifhloop

      select case (stormid(1:2))
        case ('AL');  basinid(3:3) = 'L'
        case ('CP');  basinid(3:3) = 'C'
        case ('EP');  basinid(3:3) = 'E'
        case ('WP');  basinid(3:3) = 'W'
        case ('SC');  basinid(3:3) = 'O'
        case ('EC');  basinid(3:3) = 'T'
        case ('AU');  basinid(3:3) = 'U'
        case ('SP');  basinid(3:3) = 'P'
        case ('SI');  basinid(3:3) = 'S'
        case ('BB');  basinid(3:3) = 'B'
        case ('AA');  basinid(3:3) = 'A'
        case default; basinid(3:3) = '*'
      end select

      basinid(1:2) = stormid(3:4)

      write (55,83) '99',catcf,yymmddhh,intlat(1),intlon(1)
     &           ,intlat(3),intlon(3),intlat(5),intlon(5)
     &           ,intlat(7),intlon(7),intlat(9),intlon(9)
     &           ,intlat(11),intlon(11),intlat(13),intlon(13)
     &           ,basinid

   83 format (a2,a4,a8,14i4,1x,a3)
c
C      close (55)
      return
      end

c---------------------------------------------------------------------
c
c---------------------------------------------------------------------
      subroutine read_nlists (dthresh,cmodel)
c
c     ABSTRACT: This subroutine simply reads in the namelists that are
c     created in the shell script.  Namelist datain currently contains
c     just the value of the dthresh parameter.
c
      real      dthresh
      character*5 :: cmodel
c
      namelist/datain/dthresh,cmodel
c
      read (5,NML=datain,END=801)
  801 continue

      print *,' '
      print *,'After datain namelist read in trakave.f, '
      print *,'namelist parms follow: ' 
      print *,' '
      print *,'Distance threshold (dthresh) = ',dthresh
      print *,'cmodel = ',cmodel
c
      return
      end
c---------------------------------------------------------------------
cJ.Peng-04-15-2013--------------------------
      subroutine output_atcfunix1 (stormid,cymdh,outlat,outlon
     &                           ,stdvmx,stdpce,yspread
     &                           ,outvmx,outpce,ifh
     &                           ,ctype,xprob,dthresh,catcf,fcsthrs)
c
c     ABSTRACT: This subroutine  outputs a 1-line message for a given
c     storm at an input forecast hour in the new ATCF UNIX format.
c     Unlike the old atcf DOS format in which you waited until the
c     whole tracking was over to write the  output for all forecast
c     hours, with this atcfunix format, each time we are calling this
c     subroutine, it is to only write out 1 record, which will be the
c     fix info for a particular storm at a given time.  
c
c     While this new atcfunix format contains much more information than
c     the old 1-line atcf dos message, for the purposes of the ensemble
c     mean track, we will only use the slots for storm location and NOT
c     worry about the intensity parameters.  An example is shown below:
c
c     AL, 13, 2000092500, 03, AEMN, 036, 243N, 675W
c
c     NOTE: The longitudes that are passed into this subroutine are
c     given in 0 - 360, increasing eastward.  The format for the
c     atcfunix system requires that the  output be 0-180E or
c     0-180W, so we must adjust the values, if needed.  Also, the
c     values for southern latitudes must be positive (use 'N' and
c     'S' to distinguish Northern/Southern Hemispheres).
c
c     INPUT:
c     stormid   Storm ID (EP02, AL13, etc...)
c     cymdh     initial date/time of forecast (yyyymmddhh)
c     outlat    latitude fix position for this storm at this time
c               which is to be written out to the  output file
c     outlon    longitude fix position for this storm at this time
c               which is to be written out to the  output file
c     ifh       the index for the current forecast hour being output
c     ctype     'mean' for ens mean, 'mode' for ens prob mode
c     xprob     prob mode value for ctype = 'mode'. Not used for mean
c     dthresh   distance threshold used for prob mode value, to be 
c               used in model output name for 'mode' case.
c
c     LOCAL:
c     intlon    integer that holds the value of outlon*10
c     intlat    integer that holds the value of outlat*10
c
      USE ensinfo; USE maxparms

      real    outlon,outlat,xprob
      real    outvmx,outpce,stdvmx,stdpce,yspread 
      integer fcsthrs(maxtimes)
      integer intlon,intlat
      character*1  clatns,clonew
      character*2  basinid,stormnum
      character*4  stormid,ctype,cmode,catcf
      character*10 cymdh
      character*37 outatxfile

c     First convert all of the lat/lon values from reals into integers.
c     These integer values must be 10x their real value (eg. 125.4 will
c     be written out as 1254).  Convert the lon values so that they go
c     from 0-180E or 0-180W, and convert the lat values so that they are
c     positive and use 'N' or 'S' to differentiate hemispheres.

C      outatxfile = 'aemn.'//cymdh//'.cyclone.trackatcfunix'
C      open (unit=53,file=outatxfile,status='unknown',position='append')

      if (outlon < -998.0 .or. outlat < -998.0) then
        intlon = 0
        intlat = 0
        clonew = ' '
        clatns = ' '
      else
        if (outlon >= 180.0) then
c          intlon = 3600 - int(outlon * 10. + 0.5)
c          clonew = 'W'
cJ.Peng-----03/27/2015------------------------------------
          if (outlon < 360.0) then
            intlon = 3600 - int(outlon * 10. + 0.5)
            clonew = 'W'
          else
            intlon =int(outlon * 10. + 0.5) - 3600
            clonew = 'E'
          endif
cJ.Peng-----03/27/2015------------------------------------

        else
          intlon = int(outlon * 10. + 0.5)
          clonew = 'E'
        endif
        intlat = int(abs(outlat) * 10. + 0.5)
        if (outlat < 0.0) then
          clatns = 'S'
        else
          clatns = 'N'
        endif
      endif

      basinid = stormid(1:2)
      stormnum = stormid(3:4)

      if (ctype == 'mean') then
c        write (53,81) basinid,stormnum,cymdh,catcf,fcsthrs(ifh)
c     &        ,intlat,clatns,intlon,clonew
c     &        ,int(outvmx+0.5),int(outpce+0.5),stdvmx,stdpce
c     &        ,yspread

        write (53,81) basinid,stormnum,cymdh,catcf,fcsthrs(ifh)
     &        ,intlat,clatns,intlon,clonew
     &        ,int(outvmx+0.5),int(outpce+0.5)
     &        ,int(yspread+0.5),int(stdpce+0.5)
     &        ,int(stdvmx+0.5)
      else
        itemp = int(dthresh)
        cmode(1:1) = catcf(1:1)
        write (cmode(2:4),'(i3.3)') itemp
c        print *,'atcfunix, cmode= ',cmode,' itemp= ',itemp
c        write (*,79) cmode,itemp
   79   format (1x,' cmode= ....',a4,'....  itemp= ',i3.3)
        write (57,83) basinid,stormnum,cymdh,cmode,fcsthrs(ifh)
     &        ,intlat,clatns,intlon,clonew,int(xprob*100.0)
      endif

c   81 format (a2,', ',a2,', ',a10,', 03, ',a4,', ',i3.3,', ',i3,a1
c     &       ,', ',i4,a1,', ',i3,', ',i4,', ',f5.1,', ',f5.1
c     &       ,', ',f7.2)

   81 format (a2,', ',a2,', ',a10,', 03, ',a4,', ',i3.3,', ',i3,a1
     &       ,', ',i4,a1,', ',i3,', ',i4,
     &       ', XX,  34, NEQ, 0000, 0000, 0000, 0000',
     &       ', ',i4,', ',i4,', ',i3)

   83 format (a2,', ',a2,', ',a10,', 03, ',a4,', ',i3.3,', ',i3,a1
     &       ,', ',i4,a1,', ',i3)
c
C      close (53)
      return
      end

c
c-----------------------------------------------------------------------
