Index: advance.F
===================================================================
--- advance.F	(revision 489)
+++ advance.F	(working copy)
@@ -16,7 +16,7 @@
 !
       use fields_module
       use hist_module,only: nstep,modeltime,nsource
-      use heelis_module,only: heelis,colath
+      use heelis_module,only: heelis
       use weimer_module,only: weimer01  ! Weimer 2001 model
       use wei05sc,only: weimer05,       ! Weimer 2005 model
      |  wei05sc_fac  ! field-aligned current diagnostic (nmlonp1,nmlat)
Index: colath.F
===================================================================
--- colath.F	(revision 0)
+++ colath.F	(revision 0)
@@ -0,0 +1,89 @@
+!-----------------------------------------------------------------------
+      subroutine colath
+!
+! Calculate pfrac fractional presence of dynamo equation using critical
+!  convection colatitudes crit(2).  (crit is in cons module)  All in radians.
+!
+! 01/11: Move colath from heelis.F to colath.F.  Calculate crit in colath from theta0 
+!       found from ctpoten in aurora_cons, where offc=1.1 deg and dskofc = -0.08 deg 
+!       (converted to radians).  In Weimer, theta0, offc, and dskofc are recalculated.
+!       (For Weimer 2005, offc=4.2 deg, dskofc = 0 deg, and theta0 = crad = bndyfitr/2)
+!       AMIE and CISM could: (1) use calccloc in wei01gcm.F to calculate theta0, offc, dskofc,
+!        from their potential patterns (where if Bz>0 and Bz>|By|, set theta0=10), 
+!        (2) let aurora_cons calculate theta0 from their ctpoten,
+!	 (3) use the old default crit(1,2) = 15,30 deg (in radians) from cons.F by not
+!        setting ctpoten, but using the NAMELIST CTPOTEN = 13. (or less) and BYIMF = 0 
+!        (or do no set BYIMF because if missing, it is set to 0 in input.F) with
+!        POTENTIAL_MODEL = 'HEELIS' or 'NONE' (The low CTPOTEN ensures theta0<10 deg, while
+!	 byimf=0 ensures dskofa in both hemispheres is identical at -0.08 deg.),
+!        (4) set theta0=10 default in aurora_cons as CISM does
+! TEST: NONE, CTPOTEN=13., BYIMF=0 or no BYIMF - crit1,2 = 15,30, theta0=9.917, offc=1.1, dskofc=-0.08
+! TEST: HEELIS, CTPOTEN=13., no BYIMF - crit1,2 = 15,30, theta0=9.917, offc=1.1, dskofc=-0.08
+! TEST bz=+4.4, By=2.0, Vsw=484, 02080: W01 theta0=10 (if Bz>0,Bz>|By|), W05 theta0=13.36 deg
+!
+      use dynamo_module,only: nmlat0,pfrac
+      use params_module,only: nmlonp1
+      use aurora_module,only: theta0, dskofc, offc ! see aurora.F
+      use magfield_module,only: sunlons 
+      use cons_module,only: rtd,
+     |  crit,        ! critical colatitudes crit(2)
+     |  ylonm,ylatm  ! magnetic grid lons, lats
+      implicit none
+!
+! Local:
+      integer :: i,j
+      real :: sinlat,coslat,aslonc,ofdc,cosofc,sinofc,crit1deg
+      real,dimension(nmlonp1,nmlat0) :: colatc
+
+! 01/11 bae:  Revise crit in rad so crit(1)=theta0(=crad in rad)+5deg, crit(2)=crit(1)+15deg
+        crit1deg = max(15.,0.5*(theta0(1)+theta0(2))*rtd + 5.)
+	crit1deg = min(30.,crit1deg) 
+! To make the same as in cons.F, comment out next line
+ 	crit(1) = crit1deg/rtd
+	crit(2) = crit(1) + 15./rtd
+!
+! offc(2), dskofc(2) are for northern hemisphere aurora (see aurora.F)
+! 01/11 bae:  Revised so that use average of both hemispheres instead of NH only
+!     ofdc = sqrt(offc(2)**2+dskofc(2)**2)
+      ofdc = sqrt((0.5*(offc(1)+offc(2)))**2 +
+     |  (0.5*(dskofc(1)+dskofc(2)))**2)
+      cosofc = cos(ofdc)
+      sinofc = sin(ofdc)
+      aslonc = asin(0.5*(dskofc(1)+dskofc(2))/ofdc)
+!     aslonc = asin(dskofc(2)/ofdc)
+! TEMP
+!      write (6,"(1x,'COLATH: crit1,2 theta0 offc dskofc =',8e12.4)")
+!    |   crit(1)*rtd,crit(2)*rtd,theta0(1)*rtd,theta0(2)*rtd,
+!    |   offc(1)*rtd,offc(2)*rtd,dskofc(1)*rtd,dskofc(2)*rtd
+!
+! Define colatc with northern convection circle coordinates
+! sunlons(nlat): sun's longitude in dipole coordinates (see sub sunloc)
+!
+      do j=1,nmlat0
+        sinlat = sin(abs(ylatm(j+nmlat0-1)))
+        coslat = cos(    ylatm(j+nmlat0-1))
+        do i=1,nmlonp1
+          colatc(i,j) = cos(ylonm(i)-sunlons(1)+aslonc)
+          colatc(i,j) = acos(cosofc*sinlat-sinofc*coslat*colatc(i,j))
+        enddo ! i=1,nmlonp1
+
+!       write(6,"('colath: j=',i3,' colatc(:,j)=',/,(6e12.4))")
+!    |    j,colatc(:,j)*rtd
+
+!
+! Calculate fractional presence of dynamo equation at each northern
+! hemisphere geomagnetic grid point. Output in pfrac(nmlonp1,nmlat0)
+!
+	do i=1,nmlonp1
+          pfrac(i,j) = (colatc(i,j)-crit(1))/(crit(2)-crit(1))
+          if (pfrac(i,j) < 0.) pfrac(i,j) = 0.
+          if (pfrac(i,j) >= 1.) pfrac(i,j) = 1.
+	enddo ! i=1,nmlonp1
+
+!       write(6,"('colath: j=',i3,' pfrac(:,j)=',/,(6e12.4))")
+!    |    j,pfrac(:,j)
+
+      enddo ! j=1,nmlat0
+
+      end subroutine colath
+!-----------------------------------------------------------------------
Index: aurora.F
===================================================================
--- aurora.F	(revision 489)
+++ aurora.F	(working copy)
@@ -309,7 +309,11 @@
       rh = (h2 - h1) / (h1 + h2)
 !
       theta0(isouth) = (-3.80+8.48*(ctpoten**0.1875))*dtr
-      theta0(inorth) = (-3.80+8.48*(ctpoten**0.1875))*dtr
+! Set theta0 = 10 deg so crit(1,2)=15,30 (old values in cons.F) in colath.F for CISM
+#if defined(INTERCOMM) || defined(CISMAH)
+      theta0(isouth) = 10.*dtr
+#endif
+      theta0(inorth) = theta0(isouth)
       offa(isouth) = 1.0*dtr
       offa(inorth) = 1.0*dtr
       dskofa(isouth) = 0.
Index: heelis.F
===================================================================
--- heelis.F	(revision 489)
+++ heelis.F	(working copy)
@@ -8,7 +8,7 @@
 ! Module used to calculcate the Heelis model potential in both hemispheres
 ! Byimf, Ctpoten and Power at a minimum using paramaters from aurora_cons
 !
-      use dynamo_module,only: nmlat0,pfrac,phihm
+      use dynamo_module,only: phihm
       implicit none
       contains 
 !-----------------------------------------------------------------------
@@ -37,65 +37,6 @@
 !
       end subroutine heelis
 !-----------------------------------------------------------------------
-      subroutine colath
-!
-! Calculate pfrac fractional presence of dynamo equation using critical
-!  convection colatitudes crit(2).  (crit is in cons module)
-!
-      use params_module,only: nmlonp1
-      use aurora_module,only: dskofc, offc ! see aurora.F
-      use magfield_module,only: sunlons 
-      use cons_module,only: rtd,
-     |  crit,        ! critical colatitudes crit(2)
-     |  ylonm,ylatm  ! magnetic grid lons, lats
-      implicit none
-!
-! Local:
-      integer :: i,j
-      real :: sinlat,coslat,aslonc,ofdc,cosofc,sinofc
-      real,dimension(nmlonp1,nmlat0) :: colatc
-!
-! offc(2), dskofc(2) are for northern hemisphere aurora (see aurora.F)
-      ofdc = sqrt(offc(2)**2+dskofc(2)**2)
-      cosofc = cos(ofdc)
-      sinofc = sin(ofdc)
-      aslonc = asin(dskofc(2)/ofdc)
-! TEMP
-!       write (6,"(1x,'COLATH: crit1,2 dskofc offc deg=',6e12.4)")
-!    |   crit(1)*rtd,crit(2)*rtd,dskofc(1)*rtd,offc(1)*rtd,
-!    |   dskofc(2)*rtd,offc(2)*rtd
-!
-! Define colatc with northern convection circle coordinates
-! sunlons(nlat): sun's longitude in dipole coordinates (see sub sunloc)
-!
-      do j=1,nmlat0
-        sinlat = sin(abs(ylatm(j+nmlat0-1)))
-        coslat = cos(    ylatm(j+nmlat0-1))
-        do i=1,nmlonp1
-          colatc(i,j) = cos(ylonm(i)-sunlons(1)+aslonc)
-          colatc(i,j) = acos(cosofc*sinlat-sinofc*coslat*colatc(i,j))
-        enddo ! i=1,nmlonp1
-
-!       write(6,"('colath: j=',i3,' colatc(:,j)=',/,(6e12.4))")
-!    |    j,colatc(:,j)*rtd
-
-!
-! Calculate fractional presence of dynamo equation at each northern
-! hemisphere geomagnetic grid point. Output in pfrac(nmlonp1,nmlat0)
-!
-	do i=1,nmlonp1
-          pfrac(i,j) = (colatc(i,j)-crit(1))/(crit(2)-crit(1))
-          if (pfrac(i,j) < 0.) pfrac(i,j) = 0.
-          if (pfrac(i,j) >= 1.) pfrac(i,j) = 1.
-	enddo ! i=1,nmlonp1
-
-!       write(6,"('colath: j=',i3,' pfrac(:,j)=',/,(6e12.4))")
-!    |    j,pfrac(:,j)
-
-      enddo ! j=1,nmlat0
-
-      end subroutine colath
-!-----------------------------------------------------------------------
       subroutine potm
       use params_module,only: nmlat,nmlon,nmlonp1
       use magfield_module,only: sunlons 
Index: wei05sc.F
===================================================================
--- wei05sc.F	(revision 489)
+++ wei05sc.F	(working copy)
@@ -80,9 +80,8 @@
 ! 07/08:  Added ctpoten to store average weictpoten from wei05loc
       use input_module,only: ctpoten
       use cons_module,only: rtd,ylonm,ylatm
-      use dynamo_module,only: nmlat0,pfrac,phihm
+      use dynamo_module,only: nmlat0,phihm
       use magfield_module,only: sunlons 
-      use heelis_module,only: colath
       use init_module,only: iyear,iday,uthr
       implicit none
 !
Index: input.F
===================================================================
--- input.F	(revision 489)
+++ input.F	(working copy)
@@ -2167,6 +2167,14 @@
         endif ! imf run
       endif ! weimer potential model
 !
+! 12/3/08 btf: User may set BY to test BY effect on Heelis (see aurora.F)
+! 
+      if (potential_model(1:4)=='NONE' .and. byimf == spval) then
+        write(6,"('Note input: Setting BY to 0 with no (NONE) ',
+     |      'potential model.')")
+        byimf = 0.
+      endif
+!
 ! Calendar must be advanced if using gpi or imf data:
       if (len_trim(gpi_ncfile) > 0) then
         if (calendar_advance <= 0) then
Index: cons.F
===================================================================
--- cons.F	(revision 489)
+++ cons.F	(working copy)
@@ -125,9 +125,8 @@
      |  ylatg(0:nlatp1), ! geographic latitudes (radians)
      |  ylong(nlonp1)    ! geographic longitudes (radians)
 !
-! Critical colatitude limits for use of Heelis potential in dynamo:
-      real,parameter :: 
-     |  crit(2) = (/0.261799387, 0.523598775/)
+! Critical colatitude limits (15,30 deg) for use of Heelis potential in dynamo:
+      real :: crit(2) = (/0.261799387, 0.523598775/)
 !
 ! Kut is used in filtering longitudinal waves (see filter.F):
 !
Index: wei01gcm.F
===================================================================
--- wei01gcm.F	(revision 489)
+++ wei01gcm.F	(working copy)
@@ -9,7 +9,7 @@
 !  hemispheres from Bz, By, Sws, Swd, and AL.
 !
       use params_module,only: nmlat,nmlon,nmlonp1
-      use dynamo_module,only: nmlat0,pfrac,phihm
+      use dynamo_module,only: nmlat0,phihm
       implicit none
 !
 ! Coefficients for Weimer 2001 from reading coef file:
@@ -42,13 +42,15 @@
 
       contains 
 !-----------------------------------------------------------------------
-!     subroutine calccloc(phihm)
       subroutine calccloc
 !
 !  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)
 !
 ! Input phihm: High lat model potential in magnetic coordinates (single level).
 !
@@ -78,7 +80,8 @@
      |                        dskofc
       use magfield_module,only: sunlons 
       use input_module,only:  ! from user input
-     |  byimf    ! By component of IMF (nT)      (e.g., 0.)
+     |  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, lats in radians
       implicit none
@@ -136,6 +139,15 @@
 !    |    ih,(vnx(ih,i),ylatm(jnx(ih,i)),(ylonm(inx(ih,i))-sunlons(1))
 !    |    * rtd/15.+12.,i=1,2)
 
+! 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
@@ -382,12 +394,6 @@
 
 !  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))
-	dskofc(ih) = dskof / rtd
-!  oval offset is 2.5 deg towards dawn (more neg dskof)
-	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)
 
 !  Estimate offc from lat of peak +/-3 hrs from each maximum
 !   (In colat, is nighside-dayside, but in lat is dayside-nightside)
@@ -404,11 +410,6 @@
 !  07/03:  Correction to make sure offset is at least 0.5 deg towards 0 MLT
 !       offcdeg = 0.5*(offcn+offcx)
 	offcdeg = max(0.5,0.5*(offcn+offcx))
-	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)
 
 !  Estimate theta0 from 6-18 MLT line first
 	crad0 = 90. - 0.5*abs(latinx(i18,1)+latinx(i06,2))
@@ -430,6 +431,7 @@
 !       write (6,"(1x,'06 sinl,cosl,cosm,colat asin=',5e12.4)")
 !    |   sinl06,cosl06,cosm06,colat06,asin(dskof/ofdc)
 	crad = 0.5*(colat18+colat06)
+	endif		! of using defaults for Bz>0, |Bz|>|By|
 !  Make sure crad is largest of crad and crad0
 	craduse = max(crad,crad0)
 	if (craduse .gt. crad) write (6,"(1x,'Used crad0 from 6-18',
@@ -440,6 +442,19 @@
 !       write (6,"(1x,'radius: 0,18,06,c,a deg,rad=',7e12.4)")crad0,
 !    |    colat18,colat06,crad,arad,theta0(ih),rrad(ih)
 !
+	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),dskofc(ih) = dskof / rtd
+
+!  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)
+
 ! 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,
@@ -450,12 +465,10 @@
       end subroutine calccloc
 !-----------------------------------------------------------------------
 !
-!     subroutine weimer01(pfrac,phihm,nmlat0,nmlat,nmlon,nmlonp1)
       subroutine weimer01
 
 ! 01/10:  Added ctpoten to store average weictpoten from calccloc
       use input_module,only: ctpoten
-      use heelis_module,only: colath
       use input_module,only: weimer_ncfile
 ! sunlons(nlat): sun's longitude in dipole coordinates (see sub sunloc)
       use magfield_module,only: sunlons 
@@ -480,7 +493,7 @@
 !
 ! Args:
 ! pfrac:  Fractional presence of dynamo equation in NH given critical
-!           convection colatitudes crit(2) (15,30).
+!           convection colatitudes crit(2) (15,30 in cons.F) re-calc in colath
 ! phihm:  High lat model potential in magnetic coordinates (single level).
 !
 !     integer,intent(in) :: nmlat0,nmlon,nmlat,nmlonp1
@@ -584,8 +597,6 @@
       enddo ! j=1,nmlat
 
 !  Re-calculate the values of theta0,offc,dskofc,arad,phid,phin,offa,dskofa
-!     call calccloc(phihm,nmlat,nmlonp1,nmlat0)
-!     call calccloc(phihm)
       call calccloc
 ! 01/10 bae:  Return the average of the weictpoten from the SH and NH in ctpoten
       ctpoten = 0.5*(weictpoten(1)+weictpoten(2))
@@ -593,9 +604,8 @@
 !    |    ctpoten
 !
 ! Calculate pfrac fractional presence of dynamo equation using critical
-!  convection colatitudes crit(2).  (crit is in cons module)
+!  convection colatitudes crit(2).  (crit is in cons module, re-calc in colath)
 !
-!     call colath(pfrac,nmlat0,nmlonp1)
       call colath
 !
 ! From dynamo.F:
@@ -603,10 +613,6 @@
 !    |   pfrac, ! NH fraction of potential
 !    |   phihm  ! potential in magnetic
 !
-!     do j=1,nmlat0
-!       write(6,"('weimer01: j=',i3,' pfrac=',/,(6e12.4))") j,pfrac(:,j)
-!       write(6,"('weimer01: j=',i3,' phihm=',/,(6e12.4))") j,phihm(:,j)
-!     enddo ! j=1,nmlat0
       end subroutine weimer01
 !-----------------------------------------------------------------------
 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
Index: cism_coupling.F
===================================================================
--- cism_coupling.F	(revision 489)
+++ cism_coupling.F	(working copy)
@@ -53,7 +53,9 @@
       end subroutine initialize
 !-----------------------------------------------------------------------
       subroutine import
-      integer ii,jj
+      use input_module only: ctpoten
+      integer ii,jj,js,jn
+      real cpmaxsh,cpminsh,cpmaxnh,cpminnh
 
       write(6,*)'inside receive'
 #if defined(INTERCOMM)
@@ -67,6 +69,7 @@
 
 !
 ! Add period points
+
 !
       do jj=1,nlat
        do ii=1,nlonp1-1
@@ -87,6 +90,29 @@
        enddo
       enddo
 !
+! 01/11:  Find ctpoten (kV) or min/max average from both hemispheres (gpot in V)
+!         ctpoten is NOT (yet) used to find theta0 in aurora_cons, where theta0
+!         is used in colath for crit(1,2).  Set theta0=10 so crit1,2=15,30 deg (old).
+!
+      cpmaxsh = -1000000.
+      cpmaxnh = -1000000.
+      cpminsh =  1000000.
+      cpminnh =  1000000.
+
+      do js=1,nlat/2
+       jn=nlat/2+js
+       do ii=1,nlonp1+1
+	cpmaxsh = max(cpmaxsh,gpot(ii,js))
+	cpminsh = min(cpminsh,gpot(ii,js))
+	cpmaxnh = max(cpmaxnh,gpot(ii,jn))
+	cpminnh = min(cpminnh,gpot(ii,jn))
+       enddo
+      enddo
+      ctpoten = 0.5*(cpmaxsh-cpminsh+cpmaxnh-cpminnh)*0.001
+!     write (6,"(1x,'cism CP (SH,NH,av) =',3f8.2)") 
+!   |  (cpmaxsh-cpminsh)*0.001,(cpmaxnh-cpminnh)*0.001,ctpoten	
+	
+!
 ! Add pole values
 !
       polev1=0.
