!-----------------------------------------------------------------------
!      subroutine calccloc
      subroutine calccloc(potS,potN,efxS,efxN,ifarad,rtd,mltarr,mlatdeg,model_ctpoten,theta0,offa,dskofa,phid,phin,rrad,offc,dskofc,ifbad,offcen,dskofcen,radcen)

!
!  NCAR Nov 02:  Calculate the convection center location (dskofc,offc),
!       radius (theta0), dayside entry for cusp location (phid, poten=0),
!       and nightside exit (phin, poten=0), and the associated auroral
!       radius (arad).  (Leave the aurora center as is, dskofa=0, offa=1.)
! 01/11 bae:  calccloc has a problem with Bz>0, |Bz/By|>1 conditions where
!   multiple cells are possible, so for these conditions, set defaults of:
!   theta0 = 10 deg, offc = 4.2 deg, and dskofc = 0 deg (offc and dskofc from 2005 model)
! Dec 2011 bae: Revise for CMIT tests using MIX potS(27,181) potN(27,181)
!       27 co-lats in radians or from 90 (pole), 89.9,.. to 51.662, 50.199 variable deg
!      181 theta from 0 to 2pi (2deg) where 0=noon MLT for NH and 0=midnight MLT for SH
! Aug 2012 bae: Revise for ifbad (0 if good; 1, 2, or 3 for various failures) each hem
!    and put in opposite hem (flipped for By) if only 1 bad, or default if both bad.
!    NOTE:  These revisions for ifbad.ne.0 can be revised outside this routine as desired
!      where one example of this is CMIT which makes the convection center always at the pole
!      (offc=offa=dskofc=0, dskofa=-2.5deg) outside this routine whether ifbad is 0 or not.
! 09/12 bae:  Add efxS,efxN,ifarad to find auroral radius to replace arad=crad+2 if ifarad=1

! ih=1,2 for SH and NH
!
! Input phihm: High lat model potential in magnetic coordinates (single level).
!
! (dimension 2 is for south, north hemispheres)
!  Calculate crad, offc, dskofc, and phid and phin if possible
!   Use Fig 8 of Heelis et al. [JGR, 85, 3315-3324, 1980]
!     This shows:  arad = 18.7 deg, crad = 16.7 deg (so arad = crad + 2 deg)
!           offa = offc = 3 deg (so offa = offc)
!           dskofc = 2 deg, dskofa = -0.5 deg  (so dskofa = dskofc - 2.5 deg)
!   Parameterization defaults for phid (phid(MLT)=9.39 +/- 0.21By - 12)
!                             and phin (phin(MLT)=23.50 +/- 0.15By - 12)
!   (In aurora_cons, phid=0., phin=180.*rtd)
!     (For zero By, should be phid=21.39MLT*15*rtd, phin=11.5*15*rtd)
! These are the dimensions and descriptions (corrected phid,n) from aurora.F:
!    |  theta0(2), ! convection reversal boundary in radians
!    |  offa(2),   ! offset of oval towards 0 MLT relative to magnetic pole (rad)
!    |  dskofa(2), ! offset of oval in radians towards 18 MLT (f(By))
!    |  phid(2),   ! dayside convection entrance in MLT-12 converted to radians (f(By))
!    |  phin(2),   ! night convection entrance in MLT-12 converted to radians (f(By))
!    |  rrad(2),   ! radius of auroral circle in radians
!    |  offc(2),   ! offset of convection towards 0 MLT relative to mag pole (rad)
!    |  dskofc(2)  ! offset of convection in radians towards 18 MLT (f(By))
! sunlons(nlat): sun's longitude in dipole coordinates (see sub sunloc)
!
! Dec 2011: Changes for CMIT MIX potS,N
      implicit none
     integer, parameter :: nmlat=27,nmlonp1=181
!      real :: potS(nmlat,nmlonp1),potN(nmlat,nmlonp1) !from MIX, but reversed for same mlt,mlatdeg
      real :: efxS(nmlonp1,nmlat),efxN(nmlonp1,nmlat)
      real :: potS(nmlonp1,nmlat),potN(nmlonp1,nmlat),mlatdeg(nmlat),rtd
      real :: theta0(2),offa(2),dskofa(2),phid(2),phin(2),rrad(2),offc(2),dskofc(2),phihm(nmlonp1,nmlat),mltarr(nmlonp1),model_ctpoten(2)
      integer :: ifbad(2),ifarad,ihb,ihg
      real :: dskofcen(2),offcen(2),radcen(2)
! Additional auroral parameters (see sub aurora_cons):
!      use aurora_module,only: theta0,offa,dskofa,phid,phin,rrad,offc,
!     |                        dskofc
!     use magfield_module,only: sunlons 
!      use input_module,only:  ! from user input
!     |  byimf,   ! By component of IMF (nT)      (e.g., 0.)
!     |  bzimf    ! Bz component of IMF (nT)      (e.g., 0.)
!       use cons_module,only: rtd,
!     |  ylonm,ylatm  ! magnetic grid lons, co-lats in radians
!     implicit none
!
! Args:
!     integer,intent(in) :: nmlat0,nmlat,nmlonp1
!     real,dimension(nmlonp1,nmlat),intent(in) :: phihm ! potential in magnetic
!
! Local:
      integer :: i,i1,i2,ih,j,j1,j2,k,ifwr
      real :: vnx(2,2),hem,mltdn,mltdx,mltnn,mltnx,mltd,mltn
      integer :: jnx(2,2),inx(2,2),jinx(nmlonp1,2)
      real :: vinx(nmlonp1,2),latinx(nmlonp1,2),mltinx(nmlonp1,2)
      integer :: iflgnn,iflgdx,inm3,inp3,ixm3,ixp3,i06,i18
      real :: offcn,offcx,offcdeg,dskof,ofdc,crad,arad,crad0,craduse
      real :: cosl06,sinl06,colat06,cosm06,cosl18,sinl18,colat18,cosm18
      real :: asind,dmlthalf,latnm3,latnp3,latxm3,latxp3
      real :: rhon,rhox,thetandeg,thetaxdeg,xn,yn,xx,yx,xcen,ycen
      real :: x06,y06,rho06,x18,y18,rho18,cradcoord
      real :: nxmlt(2,2)
! 09/12 bae:
      real :: efluxm(nmlonp1,nmlat),mlatefx(nmlonp1),maxefx
      integer :: narad
! 08/12 bae:
      real :: mltdby,mltnby,mltn24
      real :: byloc,byimf

! Dec 2011 for CMIT tests:
      byimf = 0.
      byloc = byimf
      ifwr = 0

!
!  Look at both hemispheres
      do ih=1,2
! Dec 2011 for CMIT tests (j1 and j2 are used often, so must redefine them for each ih)
      j1 = 1
      j2 = nmlat
! nmlat0 = (nmlat+1)/2 for end of SH in TIEGCM!     
	if (ih .eq. 1) then
	  hem = -1.
          do i=1,nmlonp1
           do j=1,nmlat
            phihm(i,j) = potS(i,j)
            efluxm(i,j) = efxS(i,j)
           enddo
          enddo
	else
	  hem = 1.
          do i=1,nmlonp1
           do j=1,nmlat
            phihm(i,j) = potN(i,j)
            efluxm(i,j) = efxN(i,j)
           enddo
          enddo
        endif

! Print out un-revised values:
!       write (6,"(1x,'Original convection/oval params (hem,By,off,dsk',
!    |    ',rad,phid,n=',10f9.4)") hem,byimf,offc(ih)*rtd,offa(ih)*rtd,
!    |    dskofc(ih)*rtd,dskofa(ih)*rtd,theta0(ih)*rtd,rrad(ih)*rtd,
!    |    phid(ih)*rtd/15.+12.,phin(ih)*rtd/15.+12.

!  08/12 bae:  Find parameterization defaults for phid (phid(MLT)=9.39 +/- 0.21By - 12)
!                             and phin (phin(MLT)=23.50 +/- 0.15By - 12)
	  mltdby = 9.39 - hem*0.21*byloc
	  mltnby = 23.50 - hem*0.15*byloc
!  The cusp is put at mltd (phid in radians) at theta0
!  Moving the cusp around a lot in MLT is not good.
!  Thus, restrict both phid and phin to be within 2 hours of mltdby and mltnby, 
!   even though By is often 0 or AMIE or CMIT runs.

!  09/12 bae: find the auroral radius from efx if ifarad>0
      if (ifarad .gt. 0) then
!  Calculate the locations of the maximum efx and get arad from mlatdeg 12h apart ~6-18 MLT
       do i=1,nmlonp1 
         maxefx = -1.
         do j=j1,j2
           if (efluxm(i,j) .gt. maxefx) then
	     maxefx = efluxm(i,j)
             mlatefx(i) = mlatdeg(j)
           endif
         enddo
! TEMP
!      if (mltarr(i).ge. 2. .and. mltarr(i) .le. 6.) write (6,"(1x,'i mlt mlat maxefx =',i4,2f8.2,e12.4)") i,mltarr(i),mlatefx(i),maxefx
!       if (mltarr(i).ge. 14. .and. mltarr(i) .le. 18.) write (6,"(1x,'i mlt mlat maxefx =',i4,2f8.2,e12.4)") i,mltarr(i),mlatefx(i),maxefx
       enddo
!  Calculate for 2-6 MLT only
	 arad = 0.
         narad = 0
         do i=1,nmlonp1-1
	   if (mltarr(i) .ge. 2. .and. mltarr(i) .le. 6.) then
             j = i+(nmlonp1-1)/2
	     if (j .gt. nmlonp1) j = j - nmlonp1 + 1
	     narad = narad + 1
	     arad = arad + 0.5*(180.-abs(mlatefx(i))-abs(mlatefx(j)))
! TEMP
!	   write (6,"(1x,'i mlt mlat j mlt mlat arad =',2(i4,2f8.2),f8.2)") i,mltarr(i),mlatefx(i),j,mltarr(j),mlatefx(j),0.5*(180.-abs(mlatefx(i))-abs(mlatefx(j)))
           endif
         enddo
       if (narad .gt. 0) arad = arad / narad
       rrad(ih) = arad / rtd
       write (6,"(1x,'ih narad arad = ',i2,i4,f8.2)") ih,narad,arad  
      endif		! if ifarad>0

!  Find min/max
	vnx(ih,1) = 0.
	vnx(ih,2) = 0.
	do j=j1,j2
	  do i=1,nmlonp1-1
	    if (phihm(i,j) .gt. vnx(ih,2)) then
	      vnx(ih,2) = phihm(i,j)
	      jnx(ih,2) = j
	      inx(ih,2) = i
              nxmlt(ih,2) = mltarr(i)
              if (abs(mlatdeg(j)) .gt. 89.99) nxmlt(ih,1) = 6.
	    endif
	    if (phihm(i,j) .lt. vnx(ih,1)) then
	      vnx(ih,1) = phihm(i,j)
	      jnx(ih,1) = j
	      inx(ih,1) = i
              nxmlt(ih,1) = mltarr(i)
              if (abs(mlatdeg(j)) .gt. 89.99) nxmlt(ih,1) = 18.
	    endif
	  enddo  !  i=1,nmlonp1-1
	enddo  !  j=j1,j2
          
! 02/10: Calculate model_ctpoten in kV from min/max in V
	model_ctpoten(ih) = 0.001 * (vnx(ih,2) - vnx(ih,1))
!       write (6,"(1x,'ih min/max pot,lat,mlt=',i2,3e12.4,2x,3e12.4))")
!    |    ih,(vnx(ih,i),ylatm(jnx(ih,i)),(ylonm(inx(ih,i))-sunlons(1))
!    |    * rtd/15.+12.,i=1,2)
       if (ifwr .eq. 1) write (6,"(1x,'ih inx jnx min/max pot,lat,mlt=',i2,2i4,e12.4,2f7.2,2x,2i4,e12.4,2f7.2)") ih,(inx(ih,i),jnx(ih,i),vnx(ih,i),mlatdeg(jnx(ih,i)),(nxmlt(ih,i)),i=1,2)

!  Feb 2012 bae:  Find cartesian coordinates for min/max, and assume center of circle fit is the midpoint
!   theta = atan2(y,x), rho = sqrt(x*x+y*y);  x=rho*cos(theta), y=rho*sin(theta)
!   theta = (MLT-6)*360/24, rho=colat
     rhon = 90.-abs(mlatdeg(jnx(ih,1)))
     rhox = 90.-abs(mlatdeg(jnx(ih,2)))
     thetandeg = (mltarr(inx(ih,1)) - 6.)*360./24.
     thetaxdeg = (mltarr(inx(ih,2)) - 6.)*360./24.
     xn = rhon*cos(thetandeg/rtd)
     yn = rhon*sin(thetandeg/rtd)
     xx = rhox*cos(thetaxdeg/rtd)
     yx = rhox*sin(thetaxdeg/rtd)
     dskofcen(ih) = -(xx - 0.5*(xx-xn))
     offcen(ih) = -0.5*(yx+yn)
!  Center in x,y coordinates is y=-offcen(ih) and x=-dskofcen(ih)
     xcen = -dskofcen(ih)
     ycen = -offcen(ih)
     radcen(ih) = sqrt( (xx-xcen)*(xx-xcen) + (yx-ycen)*(yx-ycen) )
       if (ifwr .eq. 1) write (6,"(1x,'rhon,x thetan,x x,yn x,y,x dskof,offc,radc =',11f6.1)") rhon,rhox,thetandeg,thetaxdeg,xn,yn,xx,yx,dskofcen(ih),offcen(ih),radcen(ih)
 
!  Feb 2012 bae:  Calculate ifbad for test of MLT min>12. and MLT max<12. (bad min<12., max>12.)
        ifbad(ih) = 0
        if (nxmlt(ih,1) .lt. 12. .or. nxmlt(ih,2) .gt. 12.) then
          ifbad(ih) = 1
!  Skip the rest of the calculation but put defaults in at the end for this or other ifbad cases
        else
! Feb and Aug 2012:  Need to set ifbad(ih) > 1 for other problems with undef dskof or offcdeg w then bad crad when use the -999 instead of something good!

! 01/11 bae:  Check to see if Bz positive and |Bz/By|>1:  if so, use defaults
!	for crad (theta0), offcdeg and dskof to set variables for colath.
!   	Do not redefine phin and phid, but use defaults from aurora_cons
!	if (bzimf>0 .and. abs(bzimf)>abs(byimf)) then
!	 crad0 = 0.
!	 crad = 10.
!	 offcdeg = 4.2
!         dskof = 0.
!        else
!  Set default values
      do k=1,2
	do i=1,nmlonp1
	  jinx(i,k) = -999.
	  vinx(i,k) = -999.
	  mltinx(i,k) = -999.
	  latinx(i,k) = -999.
	enddo  ! i=1,nmlonp1
      enddo   !  k=1,2
!  MLT is 0.3 MLT apart (24/80=0.3) so find half this for testing near edge
    dmlthalf = 0.5 * 24./(nmlonp1-1)

!  Find min/max +/-8 hrs (nmlonp1/3) from peaks and +/-4 lats away
!  Feb 2012:  If have small cell (less than 6 MLT wide), then latinx remains -999 
!    if go to the sign of the other cell.
!    Therefore, need to find 3 hours or less away from min/max to where latinx is not -999
! sunlons(nlat): sun's longitude in dipole coordinates (see sub sunloc) in rad
! Min:
      k = 1
! NOTE: mltdn and/or mltnn can remain -999
      mltdn = -999.
      mltnn = -999.
      iflgnn = 0
	j1 = jnx(ih,k) - 4
	if (j1 .lt. 1) j1 = 1
	j2 = jnx(ih,k) + 4
	if (j2 .gt. nmlat) j2 = nmlat
	i1 = inx(ih,k) - nmlonp1/3
	if (i1 .lt. 1) i1=1
	i2 = inx(ih,k) + nmlonp1/3
	if (i2 .gt. nmlonp1) i2=nmlonp1
! Look at mid-point part
!       if (ifwr .eq. 1) write (6,"(1x,'k j1,2 i1,2=',i2,2i3,2i4)") k,j1,j2,i1,i2
	do i=i1,i2
	vinx(i,k) = 0.
	mltinx(i,k) = mltarr(i)
!!	mltinx(i,k) = (ylonm(i)-sunlons(1)) * rtd / 15. + 12.
	if (mltinx(i,k) .gt. 24.) mltinx(i,k) = mltinx(i,k) - 24.
	if (mltinx(i,k) .lt.  0.) mltinx(i,k) = mltinx(i,k) + 24.
!  MLT is dmlthalf*2 apart in hours - find edge
	if (abs(mltinx(i,k)-18.) .lt. dmlthalf) i18=i
	  do j=j1,j2
	    if (phihm(i,j) .lt. vinx(i,k)) then
	      vinx(i,k) = phihm(i,j)
	      jinx(i,k) = j
!              latinx(i,k) = colatrad(j)*rtd
              latinx(i,k) = mlatdeg(j)
!!	      latinx(i,k) = ylatm(j) * rtd
	    endif
	  enddo   !  j=j1,j2
	  if (vinx(i,k) .ge. 0.) then
!  Look at vinx=0 for low values of i (decreasing time - phid)
	    if (mltinx(i,k) .gt. 4.5 .and. mltinx(i,k) .lt. 16.5) then
		 mltdn = mltinx(i,k)
	    else
!  Look at vinx=0 for high values of i (increasing time - phin)
	      if (iflgnn .eq. 0) then
		mltnn = mltinx(i,k)
		iflgnn = 1
	      endif
	    endif
	  endif
	enddo   !  i=i1,i2
!	if (ifwr .eq. 1) write (6,"(1x,'latinx(i,k),i=i1,i2)=',10f7.2)") (latinx(i,k),i=i1,i2)
!         write (6,"(1x,'knx i j v mlt lat =',3i3,3e12.4)") (k,i,
!    |     jinx(i,k),vinx(i,k),mltinx(i,k),latinx(i,k),i=i1,i2)
!          write (6,"(1x,'knx1 i j v mlt lat =',3i4,3e12.4)") (k,i,jinx(i,k),vinx(i,k),mltinx(i,k),latinx(i,k),i=i1,i2)
!  Now look at i<1 for dusk side:
	if (inx(ih,k) - nmlonp1/3 .lt. 1) then
	  i1 = inx(ih,k) - nmlonp1/3 + nmlonp1 - 1
	  i2 = nmlonp1
!       if (ifwr .eq. 1) write (6,"(1x,'i<1 dusk: i1,2=',2i4)") i1,i2
	  do i=i1,i2
	  vinx(i,k) = 0.
	mltinx(i,k) = mltarr(i)
!!	  mltinx(i,k) = (ylonm(i)-sunlons(1)) * rtd / 15. + 12.
	  if (mltinx(i,k) .gt. 24.) mltinx(i,k) = mltinx(i,k) - 24.
	  if (mltinx(i,k) .lt.  0.) mltinx(i,k) = mltinx(i,k) + 24.
!  MLT is dmlthalf*2 apart in hours - find edge
	  if (abs(mltinx(i,k)-18.) .lt. dmlthalf) i18=i
	    do j=j1,j2
	      if (phihm(i,j) .lt. vinx(i,k)) then
		vinx(i,k) = phihm(i,j)
		jinx(i,k) = j
!             latinx(i,k) = colatrad(j)*rtd
              latinx(i,k) = mlatdeg(j)
!!		latinx(i,k) = ylatm(j) * rtd
	      endif
	    enddo   !  j=j1,j2
	    if (vinx(i,k) .ge. 0.) then
!  Look at vinx=0 for low values of i (decreasing time - phid)
!	      if (mltinx(i,k) .gt. 4.5 .and. mltinx(i,k) .lt. 16.5)
!     |             mltdn = mltinx(i,k)
	      if (mltinx(i,k) .gt. 4.5 .and. mltinx(i,k) .lt. 16.5) mltdn = mltinx(i,k)
	    endif
	  enddo   !  i=i1,i2
!	if (ifwr .eq. 1) write (6,"(1x,'latinx(i,k),i=i1,i2)=',10f7.2)") (latinx(i,k),i=i1,i2)
!         write (6,"(1x,'knx i j v mlt lat =',3i3,3e12.4)") (k,i,
!    |     jinx(i,k),vinx(i,k),mltinx(i,k),latinx(i,k),i=i1,i2)
!          write (6,"(1x,'knx2 i j v mlt lat =',3i4,3e12.4)") (k,i,jinx(i,k),vinx(i,k),mltinx(i,k),latinx(i,k),i=i1,i2)
	endif
!  Now look at i>nmlonp1 for dusk side:
	if (inx(ih,k) + nmlonp1/3 .gt. nmlonp1) then
	  i2 = inx(ih,k) + nmlonp1/3 - nmlonp1 + 1
	  i1 = 1
!      if (ifwr .eq. 1) write (6,"(1x,'i>nmlonp1 dusk: i1,2=',2i4)") i1,i2
	  do i=i1,i2
	  vinx(i,k) = 0.
	mltinx(i,k) = mltarr(i)
!!	  mltinx(i,k) = (ylonm(i)-sunlons(1)) * rtd / 15. + 12.
	  if (mltinx(i,k) .gt. 24.) mltinx(i,k) = mltinx(i,k) - 24.
	  if (mltinx(i,k) .lt.  0.) mltinx(i,k) = mltinx(i,k) + 24.
!  MLT is dmlthalf*2 apart in hours - find edge
	  if (abs(mltinx(i,k)-18.) .lt. dmlthalf) i18=i
	    do j=j1,j2
	      if (phihm(i,j) .lt. vinx(i,k)) then
		vinx(i,k) = phihm(i,j)
		jinx(i,k) = j
!              latinx(i,k) = colatrad(j)*rtd
              latinx(i,k) = mlatdeg(j)
!!		latinx(i,k) = ylatm(j) * rtd
	      endif
	    enddo   !  j=j1,j2
	  if (vinx(i,k) .ge. 0.) then
!  Look at vinx=0 for high values of i (increasing time - phin)
	    if (mltinx(i,k) .le. 4.5 .or. mltinx(i,k) .ge. 16.5) then
	      if (iflgnn .eq. 0) then
		mltnn = mltinx(i,k)
		iflgnn = 1
	      endif
	    endif
	  endif
	  enddo   !  i=i1,i2
!	if (ifwr .eq. 1) write (6,"(1x,'latinx(i,k),i=i1,i2)=',10f7.2)") (latinx(i,k),i=i1,i2)
!         write (6,"(1x,'knx i j v mlt lat =',3i3,3e12.4)") (k,i,
!    |     jinx(i,k),vinx(i,k),mltinx(i,k),latinx(i,k),i=i1,i2)
!          write (6,"(1x,'knx3 i j v mlt lat =',3i4,3e12.4)") (k,i,jinx(i,k),vinx(i,k),mltinx(i,k),latinx(i,k),i=i1,i2)
	endif	! if (inx(ih,k) + nmlonp1/3 .gt. nmlonp1)
! Max:
      k = 2
! NOTE:  mltnx and/or mltdx can remain -999
      mltnx = -999.
      mltdx = -999.
      iflgdx = 0
	j1 = jnx(ih,k) - 4
	if (j1 .lt. 1) j1 = 1
	j2 = jnx(ih,k) + 4
	if (j2 .gt. nmlat) j2 = nmlat
	i1 = inx(ih,k) - nmlonp1/3
	if (i1 .lt. 1) i1=1
	i2 = inx(ih,k) + nmlonp1/3
	if (i2 .gt. nmlonp1) i2=nmlonp1
! Look at mid-point part
!       write (6,"(1x,'k j1,2 i1,2=',5i3)") k,j1,j2,i1,i2
	do i=i1,i2
	vinx(i,k) = 0.
	mltinx(i,k) = mltarr(i)
!!	mltinx(i,k) = (ylonm(i)-sunlons(1)) * rtd / 15. + 12.
	if (mltinx(i,k) .gt. 24.) mltinx(i,k) = mltinx(i,k) - 24.
	if (mltinx(i,k) .lt.  0.) mltinx(i,k) = mltinx(i,k) + 24.
!  MLT is dmlthalf*2 apart in hours - find edge
	if (abs(mltinx(i,k)-6.) .lt. dmlthalf) i06=i
	  do j=j1,j2
	    if (phihm(i,j) .gt. vinx(i,k)) then
	      vinx(i,k) = phihm(i,j)
	      jinx(i,k) = j
!              latinx(i,k) = colatrad(j)*rtd
              latinx(i,k) = mlatdeg(j)
!!	      latinx(i,k) = ylatm(j) * rtd
	    endif
	  enddo   !  j=j1,j2
	  if (vinx(i,k) .le. 0.) then
!  Look at vinx=0 for low values of i (decreasing time - phin)
	    if (mltinx(i,k) .le. 4.5 .or. mltinx(i,k) .ge. 16.5) then
		 mltnx = mltinx(i,k)
	    else
!  Look at vinx=0 for high values of i (increasing time - phid)
	      if (iflgdx .eq. 0) then
		mltdx = mltinx(i,k)
		iflgdx = 1
	      endif
	    endif
	  endif
	enddo   !  i=i1,i2
!         write (6,"(1x,'knx i j v mlt lat =',3i3,3e12.4)") (k,i,
!    |     jinx(i,k),vinx(i,k),mltinx(i,k),latinx(i,k),i=i1,i2)
!          write (6,"(1x,'knx4 i j v mlt lat =',3i4,3e12.4)") (k,i,jinx(i,k),vinx(i,k),mltinx(i,k),latinx(i,k),i=i1,i2)
!  Now look at i<1 for dawn side:
	if (inx(ih,k) - nmlonp1/3 .lt. 1) then
	  i1 = inx(ih,k) - nmlonp1/3 + nmlonp1 - 1
	  i2 = nmlonp1
	  do i=i1,i2
	  vinx(i,k) = 0.
	mltinx(i,k) = mltarr(i)
!!	  mltinx(i,k) = (ylonm(i)-sunlons(1)) * rtd / 15. + 12.
	  if (mltinx(i,k) .gt. 24.) mltinx(i,k) = mltinx(i,k) - 24.
	  if (mltinx(i,k) .lt.  0.) mltinx(i,k) = mltinx(i,k) + 24.
!  MLT is dmlthalf*2 apart in hours - find edge
	  if (abs(mltinx(i,k)-6.) .lt. dmlthalf) i06=i
	    do j=j1,j2
	      if (phihm(i,j) .gt. vinx(i,k)) then
		vinx(i,k) = phihm(i,j)
		jinx(i,k) = j
!              latinx(i,k) = colatrad(j)*rtd
              latinx(i,k) = mlatdeg(j)
!!		latinx(i,k) = ylatm(j) * rtd
	      endif
	    enddo   !  j=j1,j2
!  Look at vinx=0 for low values of i (decreasing time - phin)
	  if (vinx(i,k) .le. 0.) then
!	    if (mltinx(i,k) .le. 4.5 .or. mltinx(i,k) .ge. 16.5)
!     |           mltnx = mltinx(i,k)
	    if (mltinx(i,k) .le. 4.5 .or. mltinx(i,k) .ge. 16.5) mltnx = mltinx(i,k)
	  endif
	  enddo   !  i=i1,i2
!         write (6,"(1x,'knx i j v mlt lat =',3i3,3e12.4)") (k,i,
!    |     jinx(i,k),vinx(i,k),mltinx(i,k),latinx(i,k),i=i1,i2)
!          write (6,"(1x,'knx5 i j v mlt lat =',3i4,3e12.4)") (k,i,jinx(i,k),vinx(i,k),mltinx(i,k),latinx(i,k),i=i1,i2)
	endif
!  Now look at i>nmlonp1 for dawn side:
	if (inx(ih,k) + nmlonp1/3 .gt. nmlonp1) then
	  i2 = inx(ih,k) + nmlonp1/3 - nmlonp1 + 1
	  i1 = 1
	  do i=i1,i2
	  vinx(i,k) = 0.
	mltinx(i,k) = mltarr(i)
!!	  mltinx(i,k) = (ylonm(i)-sunlons(1)) * rtd / 15. + 12.
	  if (mltinx(i,k) .gt. 24.) mltinx(i,k) = mltinx(i,k) - 24.
	  if (mltinx(i,k) .lt.  0.) mltinx(i,k) = mltinx(i,k) + 24.
!  MLT is dmlthalf*2 apart in hours - find edge
	  if (abs(mltinx(i,k)-6.) .lt. dmlthalf) i06=i
	    do j=j1,j2
	      if (phihm(i,j) .gt. vinx(i,k)) then
		vinx(i,k) = phihm(i,j)
		jinx(i,k) = j
!              latinx(i,k) = colatrad(j)*rtd
              latinx(i,k) = mlatdeg(j)
!!		latinx(i,k) = ylatm(j) * rtd
	      endif
	    enddo   !  j=j1,j2
	  if (vinx(i,k) .le. 0.) then
!  Look at vinx=0 for high values of i (increasing time - phid)
	    if (mltinx(i,k) .gt. 4.5 .and. mltinx(i,k) .lt. 16.5) then
	      if (iflgdx .eq. 0) then
		mltdx = mltinx(i,k)
		iflgdx = 1
	      endif
	    endif
	  endif
	  enddo   !  i=i1,i2
!         write (6,"(1x,'knx i j v mlt lat =',3i3,3e12.4)") (k,i,
!    |     jinx(i,k),vinx(i,k),mltinx(i,k),latinx(i,k),i=i1,i2)
!          write (6,"(1x,'knx6 i j v mlt lat =',3i4,3e12.4)") (k,i,jinx(i,k),vinx(i,k),mltinx(i,k),latinx(i,k),i=i1,i2)
	endif	! if (inx(ih,k) + nmlonp1/3 .gt. nmlonp1) 
!  Now look at vinx=0 to find phid,n
!  Have mltdx,n from 4.5 to 16.5 MLT (or -999.)
	if (mltdx .ge. 0. .or. mltdn .ge. 0.) then
	  if (mltdx*mltdn .ge. 0.) mltd = 0.5*(mltdx+mltdn)
	  if (mltdx .ge. 0. .and. mltdn .lt. 0.) mltd = mltdx
	  if (mltdn .ge. 0. .and. mltdx .lt. 0.) mltd = mltdn
	else
!  Use parameterization defaults for phid (phid(MLT)=9.39 +/- 0.21By - 12)
!                             and phin (phin(MLT)=23.50 +/- 0.15By - 12)
	  mltd = 9.39 - hem*0.21*byimf
	endif
!  08/12:  Check to make sure mltd is within +/2h of mltdby
!	if (abs(mltd-mltdby) .gt. 2.) then
!	   if (mltd .gt. mltdby) mltd = mltdby + 2.
!	   if (mltd .lt. mltdby) mltd = mltdby - 2.
!	endif
	phid(ih) = (mltd-12.) * 15. / rtd
	if (mltnx .ge. 0. .or. mltnn .ge. 0.) then
!  Make mltnx,n from 16.5 to 28.5 MLT (or -999.)
	  if (mltnx .ge. 0. .and. mltnx .lt. 12.) mltnx = mltnx + 24.
	  if (mltnn .ge. 0. .and. mltnn .lt. 12.) mltnn = mltnn + 24.
	  if (mltnx*mltnn .ge. 0.) mltn = 0.5*(mltnx+mltnn)
	  if (mltnx .ge. 0. .and. mltnn .lt. 0.) mltn = mltnx
	  if (mltnn .ge. 0. .and. mltnx .lt. 0.) mltn = mltnn
	else
!  Use parameterization defaults for phid (phid(MLT)=9.39 +/- 0.21By - 12)
!                             and phin (phin(MLT)=23.50 +/- 0.15By - 12)
	  mltn = 23.50 - hem*0.15*byimf
	endif
!  08/12:  Check to make sure mltn is within +/2h of mltnby
	mltn24 = mltn
	if (mltn24 .lt. 12.) mltn24 = mltn24 + 24.
	if (abs(mltn24-mltnby) .gt. 2.) then
	   if (mltn24 .gt. mltnby) mltn24 = mltnby + 2.
	   if (mltn24 .lt. mltnby) mltn24 = mltnby - 2.
           if (mltn24 .ge. 24.) mltn24 = mltn24 - 24.
!!           mltn = mltn24
	endif
	phin(ih) = (mltn-12.) * 15. / rtd
       if (ifwr .eq. 1) write (6,"(1x,'mltd,n,x  mltn,n,x phid,n =',8e12.4)")  mltd,mltdn,mltdx,mltn,mltnn,mltnx,phid(ih)*rtd,phin(ih)*rtd

!  Estimate dskofc from lat of peak at 6 and 18 MLT (colat(18-6), lat(6-18))
!	dskof = abs(latinx(i06,2)) - abs(latinx(i18,1))
! Dec 2011 should be HALF this difference!
	dskof = (abs(latinx(i06,2)) - abs(latinx(i18,1)))/2.
        if (latinx(i06,2) .lt. -990. .or. latinx(i18,1) .lt. -990.) ifbad(ih) = 2
       if (ifwr .eq. 1) write (6,"(1x,'lat_i18(1),_i06(2),dskof =',3f6.1)") latinx(i18,1),latinx(i06,2),dskof

!  Estimate offc from lat of peak +/-3 hrs (nmlonp1-1)/8 from each maximum
!   (In colat, is nightside-dayside, but in lat is dayside-nightside)
	inm3 = inx(ih,1) - (nmlonp1-1)/8
	inp3 = inx(ih,1) + (nmlonp1-1)/8
	ixm3 = inx(ih,2) - (nmlonp1-1)/8
	ixp3 = inx(ih,2) + (nmlonp1-1)/8
	if (inm3 .lt. 1) inm3 = inm3 + nmlonp1 - 1
	if (inp3 .gt. nmlonp1) inp3 = inp3 - nmlonp1 + 1
	if (ixm3 .lt. 1) ixm3 = ixm3 + nmlonp1 - 1
	if (ixp3 .gt. nmlonp1) ixp3 = ixp3 - nmlonp1 + 1
	latnm3 = latinx(inm3,1)
! Feb 2012:  3 hours is too long a time if the cell is small so find the lesser of 
!   3 hr or when latinx is not -999 for cases when the cell is less than 6 hours wide
	if (latnm3 .lt. -990.) then
          do i=inm3+1,inx(ih,1)
	    if (latinx(i,1) .gt. -990. .and. latinx(i-1,1) .lt. -990.) latnm3 = latinx(i,1)
	  enddo
        endif
        latnp3 = latinx(inp3,1)
	if (latnp3 .lt. -990.) then
          do i=inx(ih,1)+1,inp3
	    if (latinx(i-1,1) .gt. -990. .and. latinx(i,1) .lt. -990.) latnp3 = latinx(i-1,1)
	  enddo
        endif
        latxm3 = latinx(ixm3,2)
	if (latxm3 .lt. -990.) then
          do i=ixm3+1,inx(ih,2)
	    if (latinx(i,2) .gt. -990. .and. latinx(i-1,2) .lt. -990.) latxm3 = latinx(i,2)
	  enddo
        endif
	latxp3 = latinx(ixp3,2)
	if (latxp3 .lt. -990.) then
          do i=inx(ih,2)+1,ixp3
	    if (latinx(i-1,2) .gt. -990. .and. latinx(i,2) .lt. -990.) latxp3 = latinx(i-1,2)
	  enddo
        endif

!      if (ifwr .eq. 1) write (6,"(1x,'inx1,2 inm,p3 ixm,p3 =',6i4)") inx(ih,1),inx(ih,2),inm3,inp3,ixm3,ixp3
! Jan 2012 should be HALF this difference! 
! Feb 2012:  offc=0.5*(abs(lat_12MLT)-abs(lat_24MLT), but have to look away from noon-midnight
!   if min/max at 18/06MLT (is not), then +/-3h is 45 deg or 0.7071 in y 
!    so want 0.5*(offcn+offcx)*0.5/0.7071
!   Better if look at 17-19MLT and 7-5MLT y differences, and NOT lat diffs!
!	offcn = abs(latinx(inm3,1)) - abs(latinx(inp3,1))
!	offcx = abs(latinx(ixp3,2)) - abs(latinx(ixm3,2))
!	offcn = 0.5*(abs(latnm3) - abs(latnp3))
!	offcx = 0.5*(abs(latxp3) - abs(latxm3))
	offcn = abs(latnm3) - abs(latnp3)
	offcx = abs(latxp3) - abs(latxm3)
!  07/03:  Correction to make sure offset is at least 0.5 deg towards 0 MLT - delete Dec11
        offcdeg = 0.5*(offcn+offcx)*0.5/0.7071
!	offcdeg = max(0.5,0.5*(offcn+offcx))
        if (latnm3 .lt. -990. .or. latnp3 .lt. -990. .or. latxp3 .lt. -990 .or. latxm3 .lt. -990.) then
!  Do not revise when ifbad(ih) already is 2 (should not see 1)
         if (ifbad(ih) .eq. 0) ifbad(ih) = 3
       endif
       if (ifwr .eq. 1) write (6,"(1x,'lat,mlt_inm3,inp3,ixp3,ixm3 offcn,x,deg=',11f6.1)") latnm3,mltinx(inm3,1),latnp3,mltinx(inp3,1),latxp3,mltinx(ixp3,2),latxm3,mltinx(ixm3,2),offcn,offcx,offcdeg

!  Estimate theta0 from 6-18 MLT line first
	crad0 = 90. - 0.5*abs(latinx(i18,1)+latinx(i06,2))
!  Estimate theta0 from 6-18 MLT line in 'convection circle coordinates'
! Transform to convection circle coordinates:
	ofdc = sqrt(offcdeg**2+dskof**2)
! Jan 2012: If ofdc=0 (center on magnetic pole), don't divide by 0 or get NaN!
        asind = 0.
        if (abs(ofdc) .gt. 1.e-5) asind = asin(dskof/ofdc)
	sinl18 = sin(abs(latinx(i18,1))/rtd)
	cosl18 = cos(abs(latinx(i18,1))/rtd)
	cosm18 = cos(mltinx(i18,1)*15./rtd+asind)
	colat18 = cos(ofdc/rtd)*sinl18-sin(ofdc/rtd)*cosl18*cosm18
	colat18 = acos(colat18)*rtd
       if (ifwr .eq. 1) write (6,"(1x,'18 sinl,cosl,cosm,colat asin=',5e12.4)")  sinl18,cosl18,cosm18,colat18,asind
!  Feb 2012 bae:  Find rho from ofdc midpoint where rho = sqrt(x*x+y*y)
        y18 = yn + offcdeg
        x18 = xn + dskof
        rho18 = sqrt(x18*x18+y18*y18)
       if (ifwr .eq. 1) write (6,"(1x,'x18 y18 rho18 =',3f6.1)") x18,y18,rho18 
!   theta = atan2(y,x), rho = sqrt(x*x+y*y);  x=rho*cos(theta), y=rho*sin(theta)
!   theta = (MLT-6)*360/24, rho=colat
	sinl06 = sin(abs(latinx(i06,2))/rtd)
	cosl06 = cos(abs(latinx(i06,2))/rtd)
	cosm06 = cos(mltinx(i06,2)*15./rtd+asind)
	colat06 = cos(ofdc/rtd)*sinl06-sin(ofdc/rtd)*cosl06*cosm06
	colat06 = acos(colat06)*rtd
        if (ifwr .eq. 1) write (6,"(1x,'06 sinl,cosl,cosm,colat asin=',5e12.4)")  sinl06,cosl06,cosm06,colat06,asind
!  Feb 2012 bae:  Find rho from ofdc midpoint where rho = sqrt(x*x+y*y)
        y06 = yx + offcdeg
        x06 = xx + dskof
        rho06 = sqrt(x06*x06+y06*y06)
       if (ifwr .eq. 1) write (6,"(1x,'x06 y06 rho06 =',3f6.1)") x06,y06,rho06 
	cradcoord = 0.5*(colat18+colat06)
        crad = 0.5*(rho06+rho18)
!!	endif		! of using defaults for Bz>0, |Bz|>|By|
!  Make sure crad is largest of crad and crad0 (within 0.001 deg in Jan 2012 to avoid printout)
	craduse = max(crad,crad0-0.001)
!	if (craduse .gt. crad+0.001) write (6,"(1x,'Used crad0 from 6-18',
!     |   'instead of calc crad =',2e12.4)") crad0,crad
	if (craduse .gt. crad) write (6,"(1x,'Used crad0 from 6-18 instead of calc crad =',2e12.4)") crad0,crad
	theta0(ih) = craduse / rtd
        if (ifarad .eq. 0) then
	 arad = craduse + 2.
	 rrad(ih) = arad / rtd
        endif
!       write (6,"(1x,'radius: 0,18,06,c,a deg,rad=',7e12.4)")crad0,
!     |    colat18,colat06,crad,arad,theta0(ih),rrad(ih)
       if (ifwr .eq. 1) write (6,"(1x,'radius: 0,18,06,c rho18,06,c,a deg=',8f6.1)")crad0,colat18,colat06,cradcoord,rho18,rho06,crad,arad
!
	offc(ih) = offcdeg / rtd
	offa(ih) = offcdeg / rtd
!       write (6,"(1x,'min/max latd3,n3 offc =',8e12.4)")
!     |    latinx(inm3,1),latinx(inp3,1),offcn,latinx(ixp3,2),
!     |    latinx(ixm3,2),offcx,offcdeg,offc(ih)
       if (ifwr .eq. 1) write (6,"(1x,'min/max latd3,n3 offc =',8e12.4)") latinx(inm3,1),latinx(inp3,1),offcn,latinx(ixp3,2),latinx(ixm3,2),offcx,offcdeg,offc(ih)

!  oval offset is 2.5 deg towards dawn (more neg dskof)
        dskofc(ih) = dskof / rtd
	dskofa(ih) = (dskof-2.5) / rtd
!       write (6,"(1x,'18,06 mlt,lat dskof,c,a=',7e12.4)")
!     |    mltinx(i18,1),latinx(i18,1),mltinx(i06,2),latinx(i06,2),
!     |    dskof,dskofc(ih),dskofa(ih)
       if (ifwr .eq. 1) write (6,"(1x,'18,06 mlt,lat dskof,c,a=',7e12.4)") mltinx(i18,1),latinx(i18,1),mltinx(i06,2),latinx(i06,2),dskof,dskofc(ih),dskofa(ih)

! Print out revised values:
!       write (6,"('Revised convection/oval params hem,By,off,dsk,',
!     |    'rad,phid,n=',i2,9e12.4)")ih,byimf,offc(ih)*rtd,offa(ih)*rtd,
!     |    dskofc(ih)*rtd,dskofa(ih)*rtd,theta0(ih)*rtd,rrad(ih)*rtd,
!     |    phid(ih)*rtd/15.+12.,phin(ih)*rtd/15.+12.
      if (ifwr .eq. 1) write (6,"('Revised convection/oval params hem,By,off,dsk,rad,phid,n=',i2,9e12.4)")ih,byimf,offc(ih)*rtd,offa(ih)*rtd,dskofc(ih)*rtd,dskofa(ih)*rtd,theta0(ih)*rtd,rrad(ih)*rtd,phid(ih)*rtd/15.+12.,phin(ih)*rtd/15.+12.
      endif	! for ifbad(ih)=0 or 1 from MLT min,max test
      enddo  ! ih=1,2

!  Sep 2012:  Check arad if one hem much smaller than other (can have arad=0 if
!    have a large precip at pole like on 06348 2045 UT in SH assoc w big neg cell)
!   Just set both to largest arad (Ra) if difference > 4 deg (crad~arad-4deg)
      if (abs(rrad(1)-rrad(2))*rtd .gt. 4.) then
! Changed 42/576=7.3% of time for 06348-349
       write (6,"(1x,'Ra values more than 4 deg different! Change to max Ra =',2f6.2)") rrad(1)*rtd,rrad(2)*rtd
       if (rrad(1)*rtd .lt. rrad(2)*rtd-4.) then
         rrad(1) = rrad(2)
       endif
       if (rrad(2)*rtd .lt. rrad(1)*rtd-4.) then
         rrad(2) = rrad(1)
       endif
      endif     
!  Aug 2012:  Check if bad.  If both bad, leave the defaults found in aurora_cons
!  SET DEFAULTS using radcen or min crad of 10. used in colath for crit(1,2)
!  These defaults can be revised later if desired  
!  NOTE:  If 1 hemisphere only is bad (ihb) and one is good (ihg), set ihb to ihg with neg By      
      if (ifbad(1) .gt. 0 .and. ifbad(2) .gt. 0) then
       do ih=1,2
!  Can use default of crad=radcen if ifarad=0, or crad=arad-4 (or arad-2) for ifarad=1
	  craduse = max(10.,radcen(ih))
! CANNOT have offc=dskofc=0 or divide by zero in potm in heelis.F
          offcdeg = 0.1
	  dskof = 0.
	  theta0(ih) = craduse / rtd
          if (ifarad .eq. 0) then
	   arad = craduse + 2.
	   rrad(ih) = arad / rtd
          else
	   theta0(ih) = rrad(ih)-4./rtd
          endif
	  offc(ih) = offcdeg / rtd
	  offa(ih) = offcdeg / rtd
!  oval offset is 2.5 deg towards dawn (more neg dskof)
          dskofc(ih) = dskof / rtd
	  dskofa(ih) = (dskof-2.5) / rtd
!  Use parameterization defaults for phid (phid(MLT)=9.39 +/- 0.21By - 12)
!                             and phin (phin(MLT)=23.50 +/- 0.15By - 12)
	  mltd = 9.39 - hem*0.21*byimf
	  phid(ih) = (mltd-12.) * 15. / rtd
	  mltn = 23.50 - hem*0.15*byimf
	  phin(ih) = (mltn-12.) * 15. / rtd
       enddo	! ih=1,2
      else
! Now revise if only 1 hemisphere is bad (and do nothing if both are good)
       if (ifbad(1)+ifbad(2)>0) then
       if (ifbad(1)>0) then
        ihb = 1
        ihg = 2 
        hem = 1.
       else
        ihb = 2
        ihg = 1
        hem = -1.
       endif
       theta0(ihb) = theta0(ihg)
       if (ifarad .eq. 0) rrad(ihb) = rrad(ihg)
       offc(ihb) = offc(ihg)
       offa(ihb) = offa(ihg)
       dskofc(ihb) = -dskofc(ihg)
       dskofa(ihb) = dskofc(ihb) - 2.5/rtd
       phid(ihb) = phid(ihg) + hem*0.21*byimf*2.*15./rtd
       phin(ihb) = phin(ihg) + hem*0.15*byimf*2.*15./rtd
       endif        ! for revisions for 1 bad hemispheres
      endif         ! for revisions for both bad hemispheres or other

      return
      end subroutine calccloc
!-----------------------------------------------------------------------