!-----------------------------------------------------------------------
      subroutine wei05loc (ih)
! ih=1,2 for SH,NH called from weimer05

!  Jan 2012 bae:  added byloc to limit size of byimf (probably incorrect assymmetry in hems)
!  NCAR July 08 bae:  Write wei05loc to use values from Weimer 2005
!       setboundary for convection offc=4.2 deg, dskofc=0, theta0=bndyfitr
!       where deg should be in radians.  wei05loc replaces calccloc05.
!
! (dimension 2 is for south, north hemispheres)
!  Calculate offa, dskofa, rrad, phid, and phin from Weimer 2005 offc, dskofc, theta0
!   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)
!  05/08:  But formulae for ra-rc using IMF CP between 5-7 deg (not 2) so use
!          difference of ra(max IMF CP or HP)-rc(IMF CP) as in aurora.F
! 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))
!                      phid is the MLT-12 location of the cusp on the dayside
!    |  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)
!
      use pdynamo_module,only: nmlat0,phihm
! 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: &
        byimf,   & ! By component of IMF (nT)      (e.g., 0.)
        bzimf,   & ! Bz component of IMF (nT)      (e.g., 0.)
        ctpoten, & ! cross-cap potential (kV)      (e.g., 45.)
        power      ! hemispheric power   (GW)      (e.g., 16.)
      use cons_module,only: &
        rtd,ylonm,ylatm  ! magnetic grid lons, lats in radians
      implicit none
!
! Args:
      integer,intent(in) :: ih
!
! Local:
! 05/08:  Added rccp etc
      real(r8) :: rccp,racp,rahp,ramx,diffrac,plevel,tmltmin,tmltmax
      real(r8) :: diffrac2
      real(r8) :: offcdegp(2),dskofcp(2),radegp(2)
      integer :: i,i1,i2,j,j1,j2
      real(r8) :: vnx(2,2),hem,mltd,mltn
      integer :: jnx(2,2),inx(2,2)
      real(r8) :: offcdeg,dskof,arad,crad
      real(r8) :: cosl06,sinl06,colat06,cosm06,cosl18,sinl18,colat18,cosm18
      real(r8) :: byloc

! Limit size of byimf in phin and phid calculations (as in aurora.F) 
!  NOTE:  This byloc is assymetric in hemisphere, which is probably not correct
      byloc = byimf
      if (byloc .gt. 7._r8) byloc = 7._r8
      if (byloc .lt. -11._r8) byloc = -11._r8
!
!  ih=1 is SH, ih=2 is NH
	if (ih .eq. 1) then
	  j1 = 1
	  j2 = nmlat0
	  hem = -1._r8
	else
	  j1 = nmlat0 + 1
	  j2 = nmlat
	  hem = 1._r8
	endif
! Print out un-revised values:
!       write (iulog,"(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._r8+12._r8,phin(ih)*rtd/15._r8+12._r8
!  Find min/max
	vnx(ih,1) = 0._r8
	vnx(ih,2) = 0._r8
	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
	    endif
	    if (phihm(i,j) .lt. vnx(ih,1)) then
	      vnx(ih,1) = phihm(i,j)
	      jnx(ih,1) = j
	      inx(ih,1) = i
	    endif
	  enddo  !  i=1,nmlonp1-1
	enddo  !  j=j1,j2
! 05/08: Calculate weictpoten in kV from Weimer model min/max in V
	weictpoten(ih) = 0.001 * (vnx(ih,2) - vnx(ih,1))
	tmltmin = (ylonm(inx(ih,1))-sunlons(1)) * rtd/15._r8 + 12._r8
	if (tmltmin .gt. 24.) tmltmin = tmltmin - 24._r8
	tmltmax = (ylonm(inx(ih,2))-sunlons(1)) * rtd/15._r8 + 12._r8
	if (tmltmax .gt. 24.) tmltmax = tmltmax - 24._r8
!       write(iulog,"('ih Bz By Hp ctpoten,wei min/max potV,lat,mlt=',i2, &
!         5f8.2,2x,e12.4,2f8.2,2x,e12.4,2f8.2))") ih,bzimf,byimf,power, &
!         ctpoten,weictpoten(ih), &
!         vnx(ih,1),ylatm(jnx(ih,1))*rtd,tmltmin, &
!         vnx(ih,2),ylatm(jnx(ih,2))*rtd,tmltmax
! 05/08: From aurora_cons, calculate convection and aurora radii using IMF convection
!   and power (plevel);  racp (DMSP/NOAA) - rccp (AMIE) = 5.32 (Bz>0) to 6.62 (Bz<0) deg
!  Heelis et al [1980, JGR, 85, pp 3315-3324] Fig 8: ra=rc+2deg, and is 2.5 deg to dusk
	rccp = -3.80_r8+8.48_r8*(weictpoten(ih)**0.1875_r8)
	racp = -0.43_r8+9.69_r8*(weictpoten(ih)**0.1875_r8)
	plevel = 0._r8
	if (power >=1.00_r8) plevel = 2.09_r8*alog(power)
	rahp = 14.20_r8 + 0.96_r8*plevel
	ramx = max(racp,rahp)
	diffrac = ramx - rccp

!  Set default values
!  Use parameterization defaults for phid (phid(MLT)=9.39 +/- 0.21By - 12)
!                             and phin (phin(MLT)=23.50 +/- 0.15By - 12)
	mltd = 9.39_r8 - hem*0.21_r8*byloc
	mltn = 23.50_r8 - hem*0.15_r8*byloc
	phid(ih) = (mltd-12._r8) * 15._r8 / rtd
	phin(ih) = (mltn-12._r8) * 15._r8 / rtd
! 05/18/08:  Note that phid,phin are only for Heelis and are irrelevant for Weimer
!       write (iulog,"(1x,'mltd mltn phid,n =',4f8.2)") &
!        mltd,mltn,phid(ih)*rtd/15._r8,phin(ih)*rtd/15._r8
!  Use default constant value of offcdegp from setboundary in Weimer 2005
	offcdeg = 4.2_r8
	offcdegp(ih) = offcdeg
	offc(ih) = offcdegp(ih) / rtd
	offa(ih) = offcdegp(ih) / rtd
!       write (iulog,"(1x,'offcdeg,rad =',2e12.4)") offcdeg,offc(ih)
	dskof = 0._r8
	dskofcp(ih) = dskof
	dskofc(ih) = dskof / rtd
!  oval offset is 2.5 deg towards dawn (more neg dskof)
	dskofa(ih) = (dskof-2.5_r8) / rtd
!       write(iulog,"(1x,'dskof,c,a=',3f8.2)") &
!         dskof,dskofc(ih)*rtd,dskofa(ih)*rtd
! Set crad from bndyfitr/2 of setboundary of Weimer 2005
	crad = bndyfitr/2._r8
!      write(iulog,"(1x,'wei05loc: ih,bz,y,crad =',i2,3f8.2)") &
!         ih,bzimf,byimf,crad
!  Fig 8 Heelis et al [1980]: ra=rc+2deg, and shifted 2.5 deg to dusk
	arad = crad + 2._r8
! 05/08:  Make ra=rc+diffrac(=ramx-rccp) - same difference as in aurora.F
! Choose to have arad=crad(Weimer) + diffrac(same diff as in aurora.F)
	arad = crad + diffrac
! 08/08: OR make ra=ramx=max(racp,rahp) so diffrac=arad-crad
	diffrac2 = ramx - crad
! Choose to have arad=ramx (same as in aurora.F as determined by P/CP)
!       arad = ramx
	theta0(ih) = crad / rtd
	radegp(ih) = arad
	rrad(ih) = arad / rtd
!       write (iulog,"(1x,'radius: crad,rccp,racp,rahp diffa-c', &
!        '(aurF,ramx-Weic) ramx,Weic+d,arad deg=',9f8.2)") crad,rccp, &
!        racp,rahp,diffrac,diffrac2,ramx,crad+diffrac,arad

! Print out revised values (revised 05/08):
!       write (iulog,"(1x,'Revised convection/oval params (off,dsk,', &
!         'rad,phid,n=',8f9.4)")offc(ih)*rtd,offa(ih)*rtd, &
!         dskofc(ih)*rtd,dskofa(ih)*rtd,theta0(ih)*rtd,rrad(ih)*rtd, &
!         phid(ih)*rtd/15._r8+12._r8,phin(ih)*rtd/15._r8+12._r8

      end subroutine wei05loc
!-----------------------------------------------------------------------

