program mix
  use MixReader

!  Run: gmake; then run: driver

  real, allocatable :: gridX(:,:),gridY(:,:)
  real, allocatable :: potN(:,:),pots(:,:)
  real, allocatable :: facN(:,:),facs(:,:)
  real, allocatable :: pedN(:,:),peds(:,:)
  real, allocatable :: hallN(:,:),halls(:,:)
  real, allocatable :: engN(:,:),engs(:,:)
  real, allocatable :: fluxN(:,:),fluxs(:,:)
  real, allocatable :: nwN(:,:),nws(:,:)
!  integer, parameter :: ndata=365*12*24, iyr=2006
!  integer, parameter :: ndata=365*12*24, iyr=2007
  integer, parameter :: ndata=366*12*24, iyr=2008
  character*33 mixdir07
  character*40 filtim07
  character*73 filenam07
  character*43 mixdir08,filtim08
  character*86 filenam08
  character*53 mixdir06 	! for TIEGCM 1.94.2 versin to 0 UT day 350
! character*47 mixdir06		! for TIEGCM 1.94.1 version to ~5 UT day 349
  character*38 filtim06
! character*85 filenam06	! for 1.94.1 version
  character*91 filenam06	! for 1.94.2 version
  character*20 datehrmn
  real :: radius(27,181),theta(181)
  integer nji(2),njlat,nimlt,j,i,nmph,dmin,n23
  real :: mlatdeg(27),colatrad(27),mltn(181),mlts(181),mlt(181),colatdeg(27)
  real :: pii,dtr,htr
  real :: latnx,latnn,latsx,latsn,latnamx,latnpmn,latsamx,latspmn
  real :: mltnx,mltnn,mltsx,mltsn,mltnamx,mltnpmn,mltsamx,mltspmn
      real :: theta0(2),offa(2),dskofa(2),phid(2),phin(2),rrad(2),offc(2),dskofc(2),weictpoten(2),byimf
  real :: date(ndata),bx(ndata),by(ndata),bz(ndata),swvel(ndata),swden(ndata)
  real :: epotSH(181,27),epotNH(181,27),efluxSH(181,27),efluxNH(181,27)
  real :: offcen(2),dskofcen(2),radcen(2)
  integer :: inn,inx,jnn,jnx,isn,isx,jsn,jsx,ifbad(2),nbad(3)
  real :: hem,mltday,mltnit
! 09/12 bae: eflux=flux*fac_p2e*alfa   fac_p2e=2*1.602e-9 for alfa, so half for eng=2*alfa
      real, parameter:: fac_p2e = 1.602e-9 ! convert from particle to energy flux
      integer :: ifarad
      real :: Re,areatot,dmlon,hpi_sh_mix,hpi_nh_mix
      real :: da(181),mlat2pol(27)
      real :: mlat_sh_max,mlat_nh_max,hpi_sh_max,hpi_nh_max

  pii = 4.*atan(1.)
  dtr = pii/180.
  rtd = 1./dtr
  htr = pii/12.
 
! ifarad = 1,0   if do,not calculate the auroral radius (arad) from efxS,N
        ifarad = 1 

! Open OMNI IMF file  (fails because of diff dims in NC_var_shape in 2 libs!)
!  call imf15m (iyr,ndata,date,bx,by,bz,swvel,swden)

! nh0 = 0 usually and nm1 = 1 usually except for ifwr=1 for 1 time only
  nh0 = 0
  nm1 = 1
! nmph = number of minutes per hour (30x2min WHI, 12x5min AGU Dec06)
! dmin = number of minutes between CMIT files 
! ifwr = 1,0  if do,not write out epot for plotting with plottxt.m
  ifwr = 0

  if (iyr .eq. 2006) then
   open (11,file='outcmit_dec2006_testc')
   open (12,file='ifbad_dec2006c')
!  mixdir06=     "/hao/aim4/schmitt/data/LTR-2_1_4/Agu2006/CMITs/"
   mixdir06=     "/hao/aim4/schmitt/data/LTR-2_1_5b/CMIT_Agu2006/CMITs/"
   dmin = 5
   nd0 = 348
   nd1 = 348
!   nd1 = 349
   nd2 = 349
!   nd2 = 348
  endif
! CR2060 14 Aug 2007 (1UT!) to 10 Sep 2007 (2358UT) (one last at 11 Sep 0000UT)
  if (iyr .eq. 2007) then
   open (11,file='outcmit_2007cr2060_testc')
   open (12,file='ifbad_2007cr2060c')
   mixdir07=     "/hao/aim3/emery/cmit/CR2060_CMIT/"
   dmin = 2
! Aug07
   nd0 = 226
   nd1 = 226
   nd2 = 243
   open (10,file='ls_CR2060_hdf')
! Sep07
!   nd0 = 244
!   nd1 = 244
!   nd2 = 253
!   open (10,file='ls_CR2060_sep07')
   if (ifwr .eq. 1) then
! Do test of 07238 2156UT
    nd1 = 238
    nd2 = 238
! Do test of 07239 0100UT 
!    nd1 = 239
!    nd2 = 239
   endif
  endif
  if (iyr .eq. 2008) then
   open (11,file='outcmit_whi2008_testc')
   open (12,file='ifbad_whi2008c')
   mixdir08=     "/hao/aim2/wiltbemj/mhd_data/WHI-CMIT-final/"
   dmin = 2
! Open up list of files for WHI since too many are labelled 1sec early (ie, 59Z.hdf!)
!  Mar 20-31, 2008
!  nd0 = 80
!  nd1 = 80
!  nd2 = 91
! open (10,file='ls_whi_mar08_aim2')
!  Apr 1-15, 2008
   nd0 = 92
   nd1 = 92
   nd2 = 106
   open (10,file='ls_whi_apr08_aim2')
  endif
  nmph = 60/dmin
! Change values of nh0 and nm1 to get 1 time only for ifwr=1 starting at the begin time
! TEMP to get epot in ascii file!
   if (ifwr .eq. 1) then
    open (11,file='outcmit_single_time.txt')
    open (12,file='single_time_epot.txt')
! Test 07328 2156 UT
  nh0 = 21
  nm1 = 56/2 + 1
! Do test of 07239 0100UT 
!  nh0 = 1
!  nm1 = 1
! Test 06348 0630 UT
   nh0 = 6
   nm1 = 30/5 + 1
   endif
  write (11,"('yr daynm mo da hr mn SthenN:cp x,lat,mlt n,lat,mlt; HPS HPN ifarad')")
  write (11,"('second line: SthenN:cp,theta0,arad,offc,offa,dskofc,dskofa,phid,phin,offcen,dskofcen,radcen,ifbad')") 
  write (12,"('Bad convection locations: =1 if max<12MLT or min>12MLT; =2 if latinx06 or 18 -999; =3 if latinx+/-3h -999')")
  write (12,"('yr daynm mo da hr mn SthenN:cp x,lat,mlt n,lat,mlt; HPS HPN ifarad')")
  write (12,"('second line: SthenN:cp,theta0,arad,offc,offa,dskofc,dskofa,phid,phin,offcen,dskofcen,radcen,ifbad')") 
  nbad(1) = 0
  nbad(2) = 0
  nbad(3) = 0
  
! WHI 80-106 in 2008, -03-20(80) to -03-31(91) and -04-01(92) to -04-15(106)
nd = 0;
 do ndnum=nd1,nd2
  nd = nd + 1
if (iyr .eq. 2006) then
 ida0 = 14
 imon = 12
endif
if (iyr .eq. 2007) then
 if (ndnum .lt. 244) then
   ida0 = 14
   imon = 8
 else
   ida0 = 1
   imon = 9
 endif
endif
if (iyr .eq. 2008) then
 if (ndnum .lt. 92) then
   ida0 = 20
   imon = 3
 else
   ida0 = 1
   imon = 4
 endif
endif
  ndmo = ida0 + nd1 - nd0 + nd - 1
! Would get segmentation fault if nhr=5 and nm=3,12 in 0530, but find nm=3,5 or nm=6,12??
 n23 = 23
! if (ndnum .eq. 349) n23 = 4		! for 1.94.1 version
 if (ifwr .eq. 0) then
  nh0 = 0
  if (ndnum .eq. 226) nh0 = 1
 endif
! for ifwr=1, change end time (n23 and nmph) to be same as begin time (nh0 and nm1)
 if (ifwr .eq. 1) then
  n23 = nh0
  nmph = nm1
 endif
 do nhr=nh0,n23
!  do nhr=3,3
! Dec 2006 AGU storm is every 5 min
! WHI 2008 period is every 2 min from 03-20 (80) to 04-15(106)(04-16T00), 86-88 active
 do nm=nm1,nmph
! do nm=1,1
  nmin = (nm-1)*dmin
  isbad = 0
 if (iyr .eq. 2006) then
  write (filtim06,"('AGU_CMITs_mix_2006-12-',i2.2,'T',i2.2,'-',i2.2,'-00Z.hdf')") ndmo,nhr,nmin
  write (datehrmn,"('2006 ',i3.3,' 12 ',i2.2,1x,i2.2,1x,i2.2)") ndnum,ndmo,nhr,nmin
!  write (filenam06,"(a47,a38)") mixdir06,filtim06	! for 1.94.1 version
  write (filenam06,"(a53,a38)") mixdir06,filtim06
!  write (*,*) "nd filenam06 =",nd,filenam06
  call openMixFile(filenam06)
 endif
 if (iyr .eq. 2007) then
  write (6,"(1x,i4.4,1x,i3.3,1x,i2.2,1x,i2.2,1x,i2.2,1x,i2.2)") iyr,ndnum,imon,ndmo,nhr,nmin
  write (datehrmn,"('2007 ',i3.3,1x,i2.2,1x,i2.2,1x,i2.2,1x,i2.2)") ndnum,imon,ndmo,nhr,nmin
!  Some .hdf files are corrupted like 2007-08-18T13-24
  if (ndnum .eq. 230 .and. nhr .eq. 13 .and. nmin .eq. 24) isbad = 1
  if (ifwr .eq. 1) then
   write (filtim07,"('CR2060-CMIT_mix_2007-',i2.2,'-',i2.2,'T',i2.2,'-',i2.2,'-00Z.hdf')") imon,ndmo,nhr,nmin
  else
  read (10,"(a40)") filtim07
  endif
  write (filenam07,"(a33,a40)") mixdir07,filtim07
  write (6,"(a73)") filenam07
  if (isbad .eq. 0) call openMixFile(filenam07)
 endif
 if (iyr .eq. 2008) then
  write (datehrmn,"('2008 ',i3.3,1x,i2.2,1x,i2.2,1x,i2.2,1x,i2.2)") ndnum,imon,ndmo,nhr,nmin
  if (ifwr .eq. 1) then
  write (filtim08,"('WHI-CMIT-final_mix_2008-',i2.2,'-',i2.2,'T',i2.2,'-',i2.2,'-00Z.hdf')") imon,ndmo,nhr,nmin
! Exceptions:  03-21T-00-22-00Z is 03-21T-00-21-59Z  TOO MANY!  Read from ls_whi_mar08
  write (filenam08,"(a43,a43)") mixdir08,filtim08
  else
   read (10,"(a86)") filenam08
  endif
  call openMixFile(filenam08)
 endif

 if (isbad .eq. 0) then
  call readVariable("Grid X", gridX)
!  write(*,*) "Grid X dims", shape(gridX)
  call readVariable("Grid Y", gridY)
!  write(*,*) "Grid Y dims", shape(gridY)
  call readVariable("Potential North [V]",potN)
!  write(*,*) "Potential North [V] dims", shape(potN)
  call readVariable("Potential South [V]",potS)
!  write(*,*) "Potential South [V] dims", shape(potS)
!  call readVariable("FAC North [A/m^2]",facN)
!  write(*,*) "FAC North [A/m^2] dims", shape(facN)
!  call readVariable("FAC South [A/m^2]",facS)
!  write(*,*) "FAC South [A/m^2] dims", shape(facS)
!  call readVariable("Pedersen conductance North [S]",pedN)
!  write(*,*) "Pedersen conductance North [S] dims", shape(pedN)
!  call readVariable("Pedersen conductance South [S]",pedS)
!  write(*,*) "Pedersen conductance South [S] dims", shape(pedS)
!  call readVariable("Hall conductance North [S]",hallN)
!  write(*,*) "Hall conductance North [S] dims", shape(hallN)
!  call readVariable("Hall conductance South [S]",hallS)
!  write(*,*) "Hall conductance South [S] dims", shape(hallS)
  call readVariable("Average energy North [keV]",engN)
!  write(*,*) "Average energy North [keV] dims", shape(engN)
  call readVariable("Average energy South [keV]",engS)
!  write(*,*) "Average energy South [keV] dims", shape(engS)
  call readVariable("Number flux North [1/cm^2 s]",fluxN)
!  write(*,*) "Number flux North [1/cm^2 s] dims", shape(fluxN)
  call readVariable("Number flux South [1/cm^2 s]",fluxS)
!  write(*,*) "Number flux South [1/cm^2 s] dims", shape(fluxS)
!  call readVariable("Neutral wind speed North [m/s]",nwN)
!  write(*,*) "Neutral wind speed North [m/s] dims", shape(nwN)
!  call readVariable("Neutral wind speed South [m/s]",nwS)
!  write(*,*) "Neutral wind speed South [m/s] dims", shape(nwS)
  call closeMixFile()

  !!! Do something with the data here!!!
  !!! Do something with the data here!!!
  !!! Do something with the data here!!!

  nji = shape(gridX)
  njlat = nji(1)
  nimlt = nji(2)
!  write(*,*) "njlat nimlt =",njlat,nimlt


! Calculate rho and theta polar coordinates
! Calculate radius(27,181) 181 columns all the same, rows 1-27 from 0-.6947(0 to 39.8deg co-lat)
  do j=1,njlat
   do i=1,nimlt
    radius(j,i) = sqrt(gridX(j,i)**2 + gridY(j,i)**2)
   enddo
  enddo
!  write(*,*) "radius1,27 =",radius(1,1),radius(njlat,1)
  do j=1,njlat
    colatrad(j) = radius(j,1)
    mlatdeg(j) = 90. - colatrad(j)/dtr
! colatdeg only used for printout where colat changed from 0-40 (radius) to 40-0deg
! mlat2pol is used with efluxS(N)H and epot 
    colatdeg(j) = radius(28-j,1)/dtr
    mlat2pol(j) = 90. - radius(28-j,1)/dtr
  enddo
    if (nhr==0 .and. nm==1)  write(*,*) "mlat2pol =",mlat2pol   

! theta is 0 to +2pi, where 0=noon counterclockwise to later times for NH, 
!   but is clockwise to earlier times for SH
!  Avoid atan2 for pole, since is a singularity - should be same for all j but j=1
  j=njlat
   do i=1,nimlt
    theta(i) = atan2(gridY(j,i),gridX(j,i))
    if (theta(i) < 0.) theta(i) = theta(i) + 2.*pii
   enddo
!  write(*,*) theta(1),theta(nimlt)
  do i=1,nimlt
   mltn(i) = 12.+theta(i)/htr
   if (mltn(i)>=24.) mltn(i) = mltn(i)-24.
! 09/12 corrected mlts time going backwards from noon
  mlts(i) = 12.-theta(i)/htr;
   if (mlts(i)<0.) mlts(i) = mlts(i)+24.
! Want mlt=0-24 (theta always >=0 so goes 0-24)
   mlt(i) = theta(i)/htr
  enddo
  dmlt = (nimlt-1)/24.
  deltamlt = 1./dmlt
!    if (nhr==0 .and. nm==1)   write(*,*) "nimlt dmlt deltamlt MLT =", nimlt,dmlt,deltamlt,mlt

!  Change epotn,s from 0-40 colat to 40-0 colat;  add 180 deg to epotn so starts 0 MLT;
!    change epotn so goes from 0-24 not noon-0-noon backwards in time
! Switch SH to later times is i(noon_later)=182-i(noon_earlier)
!  Add 12 hours is +90 or -90
!  Also calculate efluxSH,NH from engS,N and fluxS,N using fac_p2e=1.602e-9
  do i=1,181
  do j=1,27
    if (i .lt. 91) then
     epotSH(i,j) = potS(28-j,92-i)
     efluxSH(i,j) = engS(28-j,92-i)*fac_p2e*fluxS(28-j,92-i)
     epotNH(i,j) = potN(28-j,i+90)
     efluxNH(i,j) = engN(28-j,i+90)*fac_p2e*fluxN(28-j,i+90)
    else
     epotSH(i,j) = potS(28-j,272-i)
     efluxSH(i,j) = engS(28-j,272-i)*fac_p2e*fluxS(28-j,272-i)
     epotNH(i,j) = potN(28-j,i-90)
     efluxNH(i,j) = engN(28-j,i-90)*fac_p2e*fluxN(28-j,i-90)
    endif
  enddo
  enddo	!  for i,j

if (ifwr .eq. 1) then
  write (12,"(a,' colatdeg, mlthr epotSH(181,27),NH')") datehrmn 
  write (12,"(10f7.2)") colatdeg
  write (12,"(10f7.2)") mlt
! Convert from V to kV for printout
 do j=1,njlat
  write (12,"(10f7.1)") (epotSH(i,j)*0.001,i=1,nimlt)
 enddo
 do j=1,njlat
  write (12,"(10f7.1)") (epotNH(i,j)*0.001,i=1,nimlt)
 enddo
endif

! Calculate HP from efluxS,NH
!  Re = radius of the Earth in m at the level of the ionosphere (use 6370e+3 in MIX)
      Re = 6371.e+3
      areatot = 0.
!  The MIX grid is 2 deg in maglon (irreg in maglat) 
      dmlon = 360./180.
      
      do j=1,26
       da(j) = Re*cos(mlat2pol(j)*dtr)*dmlon*dtr*(Re*dtr*(mlat2pol(j+1)-mlat2pol(j)))
       areatot = areatot + da(j)*180.
      enddo
      hpi_sh_mix = 0.
      hpi_nh_mix = 0.
      do i=1,180
      hpi_sh_max = 0.
      hpi_nh_max = 0.
       do j=1,26
        hpi_sh_mix = hpi_sh_mix + da(j)*efluxSH(i,j)
        hpi_nh_mix = hpi_nh_mix + da(j)*efluxNH(i,j)
        if (efluxSH(i,j) > hpi_sh_max) then
	   hpi_sh_max = efluxSH(i,j)
	   mlat_sh_max = mlat2pol(j)
        endif
        if (efluxNH(i,j) > hpi_nh_max) then
	   hpi_nh_max = efluxNH(i,j)
	   mlat_nh_max = mlat2pol(j)
        endif
       enddo
! TEMP
!    write (6,"(1x,'i mlt S,N_maxmlat S,N_eflux=',i3,3f7.2,2e12.4)") i,mlt(i),mlat_sh_max,mlat_nh_max,hpi_sh_max,hpi_nh_max
      enddo
! Convert from mW/m2 * m2(da) * 1.e-12GW/mW to GW
      hpi_sh_mix = hpi_sh_mix*1.e-12  
      hpi_nh_mix = hpi_nh_mix*1.e-12  
  write (*,*) "HP SH,NH (GW) =",hpi_sh_mix,hpi_nh_mix

  i02mlt = (nimlt-1)/2 + 2*dmlt
  i08mlt = (nimlt-1)/2 + 8*dmlt
  i14mlt = (14-12)*dmlt
  i20mlt = (21-12)*dmlt
  k02mlt = 2*dmlt + 1
  k08mlt = 8*dmlt + 1
  k14mlt = 14*dmlt + 1
  k20mlt = 21*dmlt + 1
  
! Find max/min Epot
! Find max Epot and position of max for 2-8 MLT
! Find min Epot and position of min for 14-20 MLT
  potnx = 0.
  potnn = 0.
  potnamx = 0.
  potnpmn = 0.
  potsx = 0.
  potsn = 0.
  potsamx = 0.
  potspmn = 0.
  do j=1,njlat
   do i=1,nimlt
     if (potN(j,i)>potnx*1.e+3) then
      potnx = potN(j,i)*1.e-3
      latnx = mlatdeg(j)
      mltnx = mltn(i)
      if (abs(latnx)>89.99) mltnx = 6.
      inx = i
      jnx = j
     endif
     if (potN(j,i)<potnn*1.e+3) then
      potnn = potN(j,i)*1.e-3
      latnn = mlatdeg(j)
      mltnn = mltn(i)
      if (abs(latnn)>89.99) mltnn = 18.
      inn = i
      jnn = j
     endif
     if (potS(j,i)>potsx*1.e+3) then
      potsx = potS(j,i)*1.e-3
      latsx = mlatdeg(j)
      mltsx = mlts(i)
      if (abs(latsx)>89.99) mltsx = 6.
      isx = i
      jsx = j
     endif
     if (potS(j,i)<potsn*1.e+3) then
      potsn = potS(j,i)*1.e-3
      latsn = mlatdeg(j)
      mltsn = mlts(i)
      if (abs(latsn)>89.99) mltsn = 18.
      isn = i
      jsn = j
     endif
     if (potN(j,i)>potnamx*1.e+3 .and. i>i02mlt .and. i<i08mlt ) then
      potnamx = potN(j,i)*1.e-3
      latnamx = mlatdeg(j)
      mltnamx = mltn(i)
     endif
     if (potN(j,i)<potnpmn*1.e+3 .and. i>i14mlt .and. i<i20mlt ) then
      potnpmn = potN(j,i)*1.e-3
      latnpmn = mlatdeg(j)
      mltnpmn = mltn(i)
     endif
     if (potS(j,i)>potsamx*1.e+3 .and. i>k02mlt .and. i<k08mlt ) then
      potsamx = potS(j,i)*1.e-3
      latsamx = mlatdeg(j)
      mltsamx = mlts(i)
     endif
     if (potS(j,i)<potspmn*1.e+3 .and. i>k14mlt .and. i<k20mlt ) then
      potspmn = potS(j,i)*1.e-3
      latspmn = mlatdeg(j)
      mltspmn = mlts(i)
     endif
   enddo
  enddo
  cpn = potnx - potnn
  cps = potsx - potsn
  cpampmn = potnamx - potnpmn
  cpampms = potsamx - potspmn
!  write(*,*) "CP(S,N) =", cps,cpn,cpampms,cpampmn
  write (*,*) "potx(S,N) latx(S,N) mltx(S,N) i,jsx i,jnx =",potsx,potnx,latsx,latnx,mltsx,mltnx,isx,jsx,inx,jnx
  write(*,*) "potn(S,N) latn(S,N) mltn(S,N) = i,jsn i,jnn",potsn,potnn,latsn,latnn,mltsn,mltnn,isn,jsn,inn,jnn
 
!   write(11,"(7f7.2,' SH CP X,lat,MLT N,lat,MLT')") cps,potsx,latsx,mltsx,potsn,latsn,mltsn
!   write(11,"(7f7.2,' NH CP X,lat,MLT N,lat,MLT')")cpn,potnx,latnx,mltnx,potnn,latnn,mltnn
!  write (11,"(a20,2x,2f7.1,2f7.2,f7.1,2f7.2,2x,2f7.1,2f7.2,f7.1,2f7.2)") datehrmn, cps,potsx,latsx,mltsx,potsn,latsn,mltsn, cpn,potnx,latnx,mltnx,potnn,latnn,mltnn


!  Call calccloc  (CALCulate the Convection LOCation)
    call calccloc (epotSH,epotNH,efluxSH,efluxNH,ifarad,rtd,mlt,mlat2pol,weictpoten,theta0,offa,dskofa,phid,phin,rrad,offc,dskofc,ifbad,offcen,dskofcen,radcen)
!!    call calccloc (potS,potN,rtd,mltn,mlts,mlatdeg,weictpoten,theta0,offa,dskofa,phid,phin,rrad,offc,dskofc,ifbad,offcen,dskofcen,radcen)
!    call calccloc (cpmin,latmin,latmax,potS,potN,rtd,mltn,mlts,mlatdeg,weictpoten,theta0,offa,dskofa,phid,phin,rrad,offc,dskofc)
 write (*,*) datehrmn,weictpoten,hpi_sh_mix,hpi_nh_mix,' SH,NH CP(kV) HP(GW)'
  write (11,"(a20,2x,2f7.1,2f7.2,f7.1,2f7.2,2x,2f7.1,2f7.2,f7.1,4f7.2,i2)") datehrmn, cps,potsx,latsx,mltsx,potsn,latsn,mltsn, cpn,potnx,latnx,mltnx,potnn,latnn,mltnn, hpi_sh_mix,hpi_nh_mix,ifarad
  write (11,"(22x,2f7.1,22f7.2,2i2)")  weictpoten,theta0*rtd,rrad*rtd,offc*rtd,offa*rtd,dskofc*rtd,dskofa*rtd,phid*rtd/15.+12.,phin*rtd/15.+12.,offcen,dskofcen,radcen,ifbad
! Count ifbad for each hem and type of bad and print-out to 'bad' file
  if (ifbad(1)+ifbad(2)>0) then
   if (ifbad(1)==1 .or. ifbad(2)==1) nbad(1) = nbad(1) + 1
   if (ifbad(1)==2 .or. ifbad(2)==2) nbad(2) = nbad(2) + 1
   if (ifbad(1)==3 .or. ifbad(2)==3) nbad(3) = nbad(3) + 1
   if (ifbad(1)==2 .or. ifbad(2)==2 .or. ifbad(1)==3 .or. ifbad(2)==3) write (12,"('NOTE')")
  write (12,"(a20,2x,2f7.1,2f7.2,f7.1,2f7.2,2x,2f7.1,2f7.2,f7.1,4f7.2,i2)") datehrmn, cps,potsx,latsx,mltsx,potsn,latsn,mltsn, cpn,potnx,latnx,mltnx,potnn,latnn,mltnn, hpi_sh_mix,hpi_nh_mix,ifarad
  write (12,"(22x,2f7.1,22f7.2,2i2)")  weictpoten,theta0*rtd,rrad*rtd,offc*rtd,offa*rtd,dskofc*rtd,dskofa*rtd,phid*rtd/15.+12.,phin*rtd/15.+12.,offcen,dskofcen,radcen,ifbad
  endif


  if (allocated(gridX)) deallocate(gridX)
  if (allocated(gridY)) deallocate(gridY)
  if (allocated(potN)) deallocate(potN)
  if (allocated(potS)) deallocate(potS)
!  if (allocated(facN)) deallocate(facN)
!  if (allocated(facS)) deallocate(facS)
!  if (allocated(pedN)) deallocate(pedN)
!  if (allocated(pedS)) deallocate(pedS)
!  if (allocated(hallN)) deallocate(hallN)
!  if (allocated(hallS)) deallocate(hallS)
  if (allocated(engN)) deallocate(engN)
  if (allocated(engS)) deallocate(engS)
  if (allocated(fluxN)) deallocate(fluxN)
  if (allocated(fluxS)) deallocate(fluxS)
!  if (allocated(nwN)) deallocate(nwN)
!  if (allocated(nwS)) deallocate(nwS)
else
! Have corrupted CR2060 2007-08-18T13-24 .hdf so print out bogus stuff
  write (11,"(a20,'     50.0   25.0 -75.00   6.00  -25.0 -75.00  18.00     50.0   25.0 -75.00   6.00  -25.0 -75.00  18.00   5.00   5.00 1')") datehrmn
  write (11,"(22x,'   50.0   50.0 -15.00 -15.00  20.00  20.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00  11.00  11.00  24.00  24.00   0.00   0.00   0.00   0.00  15.00  15.00 4 4')") 
endif
 
! TEMP
!   if (nm .eq. 1) stop
 enddo		! do nm
 enddo		! do nhr
 enddo		! do ndnum

! AGU 2006 348 times (22 MLT, 4 lat0618 of 24 bad times or 6.9%)
! Aug 2007 12930 times (2223 MLT, 84 lat0618, 45 lat+/-3h (129 NOTEs) of 2318+1(8-18-1324UT) bad times or 17.9%!!! so 34 w bad MLT also bad 2 or 3)
! Sep 2007 7200 times (250 MLT, 10 lat0618, 11 lat+/-3h of 262 bad times or 3.6% where 21 =2 or 3 so 9 with bad MLT also have ifbad=2 or 3)
! Mar 2008 8640 times (744 MLT, 30 lat0618, 20 lat+/-3h (50 NOTEs) of 766 bad or 8.9% where 28 =2 or 3 also have bad MLT)
! Apr 2008 10800 times (1355 MLT, 68 lat0618, 22 lat+/-3h (90 NOTEs) of 1415 bad times or 13.1% with 30 =2 or 3 with bad MLT also)
! TOTAL: 39918 times, 4786 bad or 12% (mostly lower CP values)
 
  write (12,"(3i5)") nbad

end program mix