subroutine calccloc(model_ctpoten,ifbad)
! ... Change log .......................................................

! 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.)
! Jan 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 11: 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
! Jan 12: bae: Move calccloc from wei01gcm.F to util.F; divide dskofc by 2; 
!       check if 'nogood' and use defaults if: MLT(max)>12, MLT(min)<12, dsckofc<-10 or >+10,
!       offc<-5 or >+10;  set defaults from 2005 Weimer model, but defaults not based on IMF
! Aug 12: 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.
! Sep 12: bae:  Add efxS,efxN,ifarad to find auroral radius and to replace arad=crad+2 if ifarad>0
!       ifarad = 2 (do 2 things:) calc Ra (arad) and use Rc=Ra-4 deg instead of 
!       the default 12 deg for Rc (crad) when both ifbad>0
!       ih=1,2 for SH and NH
! Nov 12 Dec 12 Jan 13: Major cleanup.  See SVN commit log for
!       details. Most of the above comments are now irellevant.
!
! ... Description ......................................................
!
! Calculate crad, offc, dskofc 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)
!
! ... Use Association ..................................................
      use aurora_module,only:  ! dimension (2) is for south, north hemispheres
     |     theta0, ! theta0(2), ! convection reversal boundary in radians
     |     offc,   ! offc(2),   ! offset of convection towards 0 MLT relative to mag pole (rad)
     |     dskofc  ! dskofc(2)  ! offset of convection in radians towards 18 MLT (f(By))
      use magfield_module,only: sunlons  ! sunlons(nlat): sun's longitude in dipole coordinates (see sub sunloc)
      use cons_module,only: 
     |     ylonm,ylatm, ! magnetic grid longitudes (nmlonp1) and latitudes (nmlat) in radians
     |     rtd, dtr
      use params_module,only: 
     |     nmlat,nmlonp1 ! phihm dimensions
      use dynamo_module,only: phihm ! Input phihm(potS,N): High lat model potential in magnetic coordinates (single level).

      use cism_coupling_module,only: validate_potential_parameters
      implicit none
! ... Arguments ........................................................      
      real,intent(out) :: model_ctpoten(2)
      integer,intent(out) :: ifbad(2)
! ... Constants ........................................................      
      integer,parameter :: ifwr=0 ! ifwr=1 is a print flag for extra output
      integer,parameter :: nmlat_subset=nmlat/4-1
! ... Local variables ..................................................      
      real,dimension(nmlonp1) ::
     |     mltsh, ! CMIT/MIX: hours from noon clockwise PM to AM to noon, 
                  ! AMIE: hours from 0 to 24 MLT
     |     mltnh  ! CMIT/MIX: hours from midnight clockwise AM to PM to midnight 
                  ! AMIE: hours from 0 to 24 MLT

      real,dimension(nmlat_subset) :: 
     |     ylatm_deg_nh, ! Magnetic latitude in Northern Hemisphere (degrees)
     |     ylatm_deg_sh  ! Magnetic latitude in Southern Hemisphere (degrees)

      integer :: i,i1,i2,ih,j,j1,j2,k
      real :: phihm_min(2), phihm_max(2)
      integer j_min(2), i_min(2), kmlt_min(2)
      integer j_max(2), i_max(2), kmlt_max(2)      

      integer :: inx(2,2),jinx(nmlonp1,2)
      real :: vinx(nmlonp1,2),latinx(nmlonp1,2),mltinx(nmlonp1,2)
      integer :: inm3,inp3,ixm3,ixp3,i06,i18
      real :: offcn,offcx,offcdeg,dskof,ofdc,crad,crad0,craduse
      real :: asind,dmlthalf,latnm3,latnp3,latxm3,latxp3
      real :: mlatdeg(nmlat_subset),mltarr(nmlonp1)		
      real :: dskofcen(2),offcen(2),radcen(2)
      real :: rhon,rhox,thetandeg,thetaxdeg,xn,yn,xx,yx,xcen,ycen
      real :: x06,y06,rho06, cosl06,sinl06,colat06,cosm06, ! 06 => dawn
     |     x18,y18,rho18, cosl18,sinl18,colat18,cosm18, ! 18 => dusk
     |     cradcoord
      integer :: ihb,ihg
! ... Begin ............................................................

! Calculate magnetic local time.
! Find MLT from sunlons (sunlons(1-nlat) are all the same value so use the first value)
      do i=1,nmlonp1
       mltSh(i) = (ylonm(i)-sunlons(1)) * rtd / 15. + 12.
       if (mltSh(i) .gt. 24.) mltSh(i) = mltSh(i) - 24.
       if (mltSh(i) .lt. 0.) mltSh(i) = mltSh(i) + 24.
       mltNh(i) = mltSh(i)
      enddo

! !  Get points from poles to about +/-35 mlat
! nmlat_subset=23 from -90 about every 2 deg to -70, every 3 to -30:
!  23 magS = -0.9000E+02 -0.8812E+02 -0.8624E+02 -0.8433E+02 -0.8240E+02 
!            -0.8043E+02 -0.7841E+02 -0.7633E+02 -0.7419E+02 -0.7197E+02
!            -0.6967E+02 -0.6728E+02 -0.6480E+02 -0.6222E+02 -0.5955E+02 
!            -0.5678E+02 -0.5393E+02 -0.5100E+02 -0.4801E+02 -0.4498E+02
!            -0.4192E+02 -0.3885E+02 -0.3582E+02

      do j=1,nmlat_subset
       ylatm_deg_sh(j) = ylatm(j)*rtd
       ylatm_deg_nh(j) = ylatm(nmlat-j+1) * rtd
      enddo
    
!
!  Look at both hemispheres (ih=1 SH, ih=2 NH)
      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_subset
	if (ih .eq. 1) then
          do j=1,nmlat_subset
	   mlatdeg(j) = abs(ylatm_deg_sh(j))
          enddo
          do i=1,nmlonp1
           mltarr(i) = mltsh(i)
          enddo
	else
          do j=1,nmlat_subset
	   mlatdeg(j) = ylatm_deg_nh(j)
          enddo
          do i=1,nmlonp1
           mltarr(i) = mltnh(i)
          enddo
        endif
!
! Print out un-revised values:
        if (ifwr .eq. 1) write (6,"(1x,
     | 'Original convection/oval params (By,loc,off,dsk',
     |    'n=',11f9.4)") dskofc(ih)*rtd,theta0(ih)*rtd


!  Find min/max values & location of potential
! ih = hemisphere selection
	phihm_min(ih) =  99999999.0
	phihm_max(ih) = -99999999.0
	do j=j1,j2
	  do i=1,nmlonp1-1
	    if (phihm(i,j) .gt. phihm_max(ih)) then
	      phihm_max(ih) = phihm(i,j)
	      j_max(ih) = j
	      i_max(ih) = i
              kmlt_max(ih) = mltarr(i)
              if (abs(mlatdeg(j)) .gt. 89.99) kmlt_min(ih) = 6. !FIXME: Magic numbers
	    endif
	    if (phihm(i,j) .lt. phihm_min(ih)) then
	      phihm_min(ih) = phihm(i,j)
	      j_min(ih) = j
	      i_min(ih) = i
              kmlt_min(ih) = mltarr(i)
              if (abs(mlatdeg(j)) .gt. 89.99) kmlt_min(ih) = 18. !FIXME: Magic numbers
	    endif
	  enddo  !  i=1,nmlonp1-1
	enddo  !  j=j1,j2

! 02/10: Calculate the model ctpoten in kV from model (max - min) in V
	model_ctpoten(ih) = 0.001 * (phihm_max(ih) - phihm_min(ih))

!  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(j_min(ih)))
      rhox = 90.-abs(mlatdeg(j_max(ih)))
      thetandeg = (mltarr(i_min(ih)) - 6.)*360./24.
      thetaxdeg = (mltarr(i_max(ih)) - 6.)*360./24.
      xn = rhon*cos(thetandeg*dtr)
      yn = rhon*sin(thetandeg*dtr)
      xx = rhox*cos(thetaxdeg*dtr)
      yx = rhox*sin(thetaxdeg*dtr)
      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 (kmlt_min(ih) .lt. 12. .or. kmlt_max(ih) .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
	  cycle
        endif
! 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!

!  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 - NOT USED!
! Min:
      k = 1
      i06 = -999
      i18 = -999
	j1 = j_min(ih) - 4
	if (j1 .lt. 1) j1 = 1
	j2 = j_min(ih) + 4
	if (j2 .gt. nmlat_subset) j2 = nmlat_subset
	i1 = j_min(ih) - nmlonp1/3
	if (i1 .lt. 1) i1=1
	i2 = j_min(ih) + 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)
	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) = mlatdeg(j)
	    endif
	  enddo   !  j=j1,j2
	enddo   !  i=i1,i2
      if (ifwr .eq. 2) write (6,"(1x,'latinx(i,k),i=i1,i2)=',10f7.2)") 
     |  (latinx(i,k),i=i1,i2)
       if (ifwr .eq. 2)
     |   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 (i_min(ih) - nmlonp1/3 .lt. 1) then
	  i1 = i_min(ih) - 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)
	  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) = mlatdeg(j)
	      endif
	    enddo   !  j=j1,j2
	  enddo   !  i=i1,i2
      if (ifwr .eq. 2) write (6,"(1x,'latinx(i,k),i=i1,i2)=',10f7.2)") 
     |   (latinx(i,k),i=i1,i2)
        if (ifwr .eq. 2)
     |   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 (i_min(ih) + nmlonp1/3 .gt. nmlonp1) then
	  i2 = i_min(ih) + 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)
	  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) = mlatdeg(j)
	      endif
	    enddo   !  j=j1,j2
	  enddo   !  i=i1,i2
      if (ifwr .eq. 2) write (6,"(1x,'latinx(i,k),i=i1,i2)=',10f7.2)") 
     |   (latinx(i,k),i=i1,i2)
        if (ifwr .eq. 2)
     |   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
	j1 = j_max(ih) - 4
	if (j1 .lt. 1) j1 = 1
	j2 = j_max(ih) + 4
	if (j2 .gt. nmlat_subset) j2 = nmlat_subset
	i1 = i_max(ih) - nmlonp1/3
	if (i1 .lt. 1) i1=1
	i2 = i_max(ih) + 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=',5i3)")k,j1,j2,i1,i2
	do i=i1,i2
	vinx(i,k) = 0.
	mltinx(i,k) = mltarr(i)
	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) = mlatdeg(j)
	    endif
	  enddo   !  j=j1,j2
	enddo   !  i=i1,i2
       if (ifwr .eq. 2)
     |   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 (i_max(ih) - nmlonp1/3 .lt. 1) then
	  i1 = i_max(ih) - nmlonp1/3 + nmlonp1 - 1
	  i2 = nmlonp1
	  do i=i1,i2
	  vinx(i,k) = 0.
	mltinx(i,k) = mltarr(i)
	  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) = mlatdeg(j)
	      endif
	    enddo   !  j=j1,j2
!  Look at vinx=0 for low values of i (decreasing time - phin)
	  enddo   !  i=i1,i2
        if (ifwr .eq. 2)
     |   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 (i_max(ih) + nmlonp1/3 .gt. nmlonp1) then
	  i2 = i_max(ih) + nmlonp1/3 - nmlonp1 + 1
	  i1 = 1
	  do i=i1,i2
	  vinx(i,k) = 0.
	mltinx(i,k) = mltarr(i)
	  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) = mlatdeg(j)
	      endif
	    enddo   !  j=j1,j2
	  enddo   !  i=i1,i2
       if (ifwr .eq. 2)
     |   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 (i_max(ih) + nmlonp1/3 .gt. nmlonp1) 
       if (i06 .lt. 1 .or. i18 .lt. 1) then
	  ifbad(ih) = 2
	  cycle
       endif
       if (ifwr .eq. 1) write (6,"(1x,'lat_i18(1),_i06(2) =',2f6.1)")
     |    latinx(i18,1),latinx(i06,2)
        if (latinx(i06,2) .lt. -990. .or. latinx(i18,1) .lt. -990.) then
          ifbad(ih) = 2
          cycle
       endif
!  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 (ifwr .eq. 1) write (6,"(1x,'dskof =',f6.1)") 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
       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)
         ifbad(ih) = 3
	 cycle
       endif
! 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(latnm3) - abs(latnp3)
	offcx = abs(latxp3) - abs(latxm3)
        offcdeg = 0.5*(offcn+offcx)*0.5/0.7071

       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) !FIXME: magic number. Define/use epsilon from cons module?
	sinl18 = sin(abs(latinx(i18,1))*dtr)
	cosl18 = cos(abs(latinx(i18,1))*dtr)
	cosm18 = cos(mltinx(i18,1)*15.*dtr+asind)
	colat18 = cos(ofdc*dtr)*sinl18-sin(ofdc*dtr)*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))*dtr)
	cosl06 = cos(abs(latinx(i06,2))*dtr)
	cosm06 = cos(mltinx(i06,2)*15.*dtr+asind)
	colat06 = cos(ofdc*dtr)*sinl06-sin(ofdc*dtr)*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)

!  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
	theta0(ih) = craduse / rtd
       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
!
	offc(ih) = offcdeg / rtd
       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
       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)

      if (ifwr .eq. 1) write (6,
     |  "('Revised convection/oval params hem,By,off,dsk,n=',
     |  i2,9e12.4)")ih,offc(ih)*rtd,dskofc(ih)*rtd,
     |  theta0(ih)*rtd
      enddo  ! ih=1,2 (and cycle point for 'cycle' of the do loop to continue ih=2 or end)


      call validate_potential_parameters(offc, dskofc, theta0, ifbad)

      return
      end subroutine calccloc