subroutine kameleon_swmf (istep)

c  Ja Soon Shim:  This is an example program mainly to test the time interpolator
c  04/13 bae:  Revised to read Dec 13-15, 2006 SWMF 5 min files of epot,eflux,eave

c Variables to be used for interpolation and data extraction
c  ep = electric potential in kV
c  eave = electron average energy in keV
c  eflux = electron energy flux in W/m2
c  ex,y,z = electric field in x(phi??), y(theta??), z(up??) directions in mV/m

      use input_module,only:  ! from user input
     |  ctpoten, ! cross-cap potential (kV)      (e.g., 45.)
     |  power    ! hemispheric power   (GW)      (e.g., 16.)
      use input_module,only: kameleon_path,facmw,lbmlat
      use params_module,only: nlonp1,nmlonp1,nlat,nmlon,nmlat,
     |    nlatp1,nlon
      use cons_module,only: pi,dtr,rtd,ylonm,ylatm
      use dynamo_module,only: nmlat0,phihm
      use magfield_module,only: sunlons 
      use magfield_module,only: im,jm,dim,djm
      use init_module,only: iyear,iday,uthr
      use hist_module,only: modeltime
      use aurora_module,only: gekev,gmwpm2	! (nlonp1,0:nlatp1) in geog coords
      use dynamo_module,only: mag2geo
      use aurora_module,only: theta0,offa,dskofa,phid,phin,rrad,offc,
     |                        dskofc
      use addfld_module,only: addfld
      implicit none

!
! Args:
      integer,intent(in) :: istep
!
! Local:
      integer :: iyr,nda,imon,ida,ihh,imm
      integer :: iyear1,iyear2,imon1,imon2,iday1,iday2,ihh1,ihh2,
     |           imm1,imm2
      integer :: j,i,iv,ih,ip,ipp,j1,j2
      real :: rmlt(nmlonp1),mlat,vnx(2,2),hemctpoten(2),hemhp(2)
      real :: xshmwpm2,xnhmwpm2,pole
! Define auroral output fields in TIEGCM distorted magnetic grid (same as phihm)
      real ::  tieefx(nmlonp1,nmlat),tieekv(nmlonp1,nmlat)
      real :: Re,areatot,dmlat,dlonm,da(nmlat)
      real :: rhon,rhox,thetandeg,thetaxdeg,xn,yn,xx,yx,xcen,ycen
      real :: dskofcen(2),offcen(2),radcen(2)

      character*300 cdf_file_path1
      character*300 cdf_file_path2
      character*24 time
      character*100 variable
       
      real :: time_epoch
      integer :: tid
      integer :: iskp2,jnx(2,2),inx(2,2)

! External functions 
      external f_timeinterp_create
      external f_timeinterp_addtimestep
      external f_timeinterp_managemem,f_timeinterp_interpolate
      external f_timeinterp_timetoepoch
      external f_timeinterp_timestrtoepoch
      external f_timeinterp_closeandcleanupmemory

! Find da for HPI estimation
      if (istep<=1) then
        do i=1,nmlonp1  
! sunlons(nlat): sun's longitude in dipole coordinates (see sub sunloc) in rad
          rmlt(i) = (ylonm(i)-sunlons(1)) * rtd / 15. + 12.
        enddo
! Re=radius of the earth in m
      Re = 6371.e+3
! ylatm (97 values) and ylonm (80 values 4.5 deg apart) are arrays of latitudes 
!   and longitudes of the distorted tiegcm magnetic grid in radians - from consdyn.h
!  nmlat0 = 49 for mag equator, nmlat = 97, ylatm from -/-90 to +/-0.94, with dellat>3deg 35-48mlat
!   #77-97 or 21 values for 41.9mlat to 90mlat;  Use kameleon only for mlat>=lbmlat
      dlonm = 2.*pi/float(nmlon)
      areatot = 0.
      write (6,"(1x,'nmlon,dlonm =',i4,f7.2)") nmlon,dlonm/dtr
! Set da for NH first from eq to pole (so set j-1 to abs for SH at j=nmlat0)
      do j=nmlat0,nmlat
       if (j==nmlat) then
        dmlat = 90.*dtr - ylatm(j) + 0.5*(ylatm(j)-ylatm(j-1))
       else
        dmlat = 0.5*(ylatm(j+1)-abs(ylatm(j-1)))
       endif
       da(j) = Re*cos(ylatm(j))*dlonm*(Re*dmlat)
       write (6,"(1x,'j dmlat mlat da =',i3,3e12.4)") j,dmlat/dtr,
     |  ylatm(j)/dtr,da(j)
       areatot = areatot + da(j)*nmlon
      enddo
! Set SH from NH da
      do j=1,nmlat0-1
       da(j) = da(nmlat+1-j)
      enddo
       write (6,"(1x,'area of hem = ',e12.4)") areatot
      endif	! if (istep<=1) then

! Convert from day-of-year to month,day to get names for SWMF 5-min files
!  sub cvt2md in wei05sc.F
!
      iyr   = iyear
      nda = iday
      call cvt2md(6,iyr,nda,imon,ida)
      ihh = int(uthr)
      imm = int((uthr - ihh)*60.)
      write(time,'(i4.4,2("-",i2.2),"T",2(i2.2,":"),"00.000Z")')   
     &   iyear,imon,ida,ihh,imm 
      if (istep<=1) write (6,"(1x,'time = ',a)") time
C Find times to nearest 5 minutes
       imm1=imm/5*5
       imm2=imm1+5
       if(imm2.eq.60) then
          imm2=0         
          ihh1=ihh
          ihh2=ihh1+1
        else
          ihh1=ihh
          ihh2=ihh1 
       endif 
       if(ihh2.eq.24) then
          ihh2=0         
          iday1=ida
          iday2=iday1+1
        else
          iday1=ida
          iday2=iday1
        endif 
         
         imon1=imon
         imon2=imon1  
         iyear1=iyear
         iyear2=iyear1

C Have no 0005UT 16 Dec, so have special case
       iskp2 = 0
       if (iday2.eq.16 .and.ihh.eq.0 .and.imm.lt.5) iskp2 = 1

! Write 2 cdf file paths
      write(cdf_file_path1,'(a,i4.4,2i2.2,"-",2i2.2,"00-000.cdf")')   
     &   trim(kameleon_path),iyear1,imon1,iday1,ihh1,imm1 
      write(cdf_file_path2,'(a,i4.4,2i2.2,"-",2i2.2,"00-000.cdf")')   
     &   trim(kameleon_path),iyear2,imon2,iday2,ihh2,imm2 
      if (istep<=1) write (6,"(1x,'path1,2 = ',a,1x,a)") 
     |   trim(cdf_file_path1),trim(cdf_file_path2)

c     tid will be set when calling f_timeinterp_create. Just initialize here.
      tid=0
      time_epoch = 0.

      call f_timeinterp_create(tid)
CCC
      call f_timeinterp_timestrtoepoch(time, time_epoch)
      call f_timeinterp_addtimestep(tid, cdf_file_path1)
      if (iskp2==0) call f_timeinterp_addtimestep(tid, cdf_file_path2)
      do iv=1,3
       if (iv==1) variable = 'ep'     
       if (iv==2) variable = 'eave'
       if (iv==3) variable = 'eflux'

!      print *, 'variable =',variable
      call f_timeinterp_managemem(tid, time_epoch, variable);
      if (istep<=1) write (6,"(1x,'tid time_epoch =',i5,e12.4)")
     |  tid,time_epoch

! Get SH (j=1,nmlat0) and NH (j=nmlat0+1,nmlat)
!  Pole values are trash for SWMF, so skip j=1 and nmlat and compute later
!     do j=1,nmlat
      do j=2,nmlat-1
       ih = 1
       if (j>nmlat0) ih = 2
!  Find min/max
!	vnx(ih,1) = 0.
!	vnx(ih,2) = 0.
        do i=1,nmlonp1  

          mlat = ylatm(j)*rtd
        if (abs(mlat)<lbmlat) then
          phihm(i,j) = 0.
          tieekv(i,j) = 1.
          tieefx(i,j) = 0.
        else
!
! Obtain electric potential in kV
        if (iv==1) then
         call f_timeinterp_interpolate(tid, time_epoch, variable, 0.0,
     |    mlat, rmlt(i), phihm(i,j))
!  This code fails to work here (vnx remains zero), but works outside in kameleon/fortran_test but not later????
!	    if (phihm(i,j) .lt. vnx(ih,1)) vnx(ih,1) = phihm(i,j)
!	    if (phihm(i,j) .gt. vnx(ih,2)) vnx(ih,2) = phihm(i,j)
!        phihm(i,j) = phihm(i,j)*1000.
        endif  
        if (iv==2) then 
! Obtain electron mean energy in keV
         call f_timeinterp_interpolate(tid, time_epoch, variable, 0.0,
     |    mlat, rmlt(i), tieekv(i,j)) 
       endif
! Obtain electron energy flux in w/m2 (SWMF) or in mw/m2 (ASTRA_AMIE)
       if (iv==3) then     
        call f_timeinterp_interpolate(tid, time_epoch, variable, 0.0,
     |   mlat, rmlt(i), tieefx(i,j))
!        tieefx(i,j) = tieefx(i,j)*facmw 
       endif

       endif	!  if (abs(mlat)<lbmlat) then,else
        enddo ! i=1,nmlonp1
      enddo ! j=2,nmlat-1

! Find pole value as average of closest latitude for SWMF (and for AMIE also)
         ip = 1
         ipp = 2
         pole = 0. 
         do i=1,nmlon
           if (iv==1) pole = pole + phihm(i,ipp)
           if (iv==2) pole = pole + tieekv(i,ipp)
           if (iv==3) pole = pole + tieefx(i,ipp)
         enddo
         pole = pole/nmlon
         do i=1,nmlonp1
           if (iv==1) phihm(i,ip) = pole
           if (iv==2) tieekv(i,ip) = pole           
           if (iv==3) tieefx(i,ip) = pole
         enddo
         ip = nmlat
         ipp = nmlat-1
         pole = 0. 
         do i=1,nmlon
           if (iv==1) pole = pole + phihm(i,ipp)
           if (iv==2) pole = pole + tieekv(i,ipp)
           if (iv==3) pole = pole + tieefx(i,ipp)
         enddo
         pole = pole/nmlon
         do i=1,nmlonp1
           if (iv==1) phihm(i,ip) = pole
           if (iv==2) tieekv(i,ip) = pole           
           if (iv==3) tieefx(i,ip) = pole
         enddo

      enddo  ! iv variable

!  Calculate hem HP in GW in distorted magnetic grid and convert to mw/m2 if necessary
        hemhp(1) = 0.
        hemhp(2) = 0.
        xnhmwpm2 = 0.
        xshmwpm2 = 0.
        do j=1,nmlat0
         do i=1,nmlon
          tieefx(i,j) = tieefx(i,j)*facmw 
          tieefx(i,nmlat-j) = tieefx(i,nmlat-j)*facmw 
          hemhp(1) = hemhp(1) + da(j)*tieefx(i,j)
          hemhp(2) = hemhp(2) + da(j)*tieefx(i,nmlat-j)
          xshmwpm2 = max(xshmwpm2,tieefx(i,j))
          xnhmwpm2 = max(xnhmwpm2,tieefx(i,nmlat-j))
         enddo
!  Set nmlonp1=1
         tieefx(nmlonp1,j) = tieefx(1,j)
         tieefx(nmlonp1,nmlat-j) = tieefx(1,nmlat-j)
	enddo
! Convert from mW/m2 * m2(da) * 1.e-12GW/mW
        hemhp(1) = hemhp(1)*1.e-12
        hemhp(2) = hemhp(2)*1.e-12
! Return the average of the hemhp from the SH and NH in power
        power = 0.5*(hemhp(1) + hemhp(2))


C  Convert aurora from TGCM distorted magnetic grid to geographic one for aurora.F
       call mag2geo(tieekv(1,1),gekev(1,0),im(1,0),jm(1,0),
     |    dim(1,0),djm(1,0),nlonp1,nmlonp1,nlon,nlat+2,nmlon,nmlat)
       call mag2geo(tieefx(1,1),gmwpm2(1,0),im(1,0),jm(1,0),
     |    dim(1,0),djm(1,0),nlonp1,nmlonp1,nlon,nlat+2,nmlon,nmlat)
C     ****
C     ****     INSERT PERIODIC POINTS
C     ****
C 4/13 bae: use i=1,nlon of ekvg and efxg in aurora.F and below, BUT nlonp1 is needed for alfa,eflux!
      DO j = 1,nlat
        gekev(nlonp1,j) = gekev(1,j)
        gmwpm2(nlonp1,j) = gmwpm2(1,j)
      ENDDO
! Could calculate theta0 (cusplat and colath) and set MLT for cusp also from calccloc
!  or could use defaults in aurora.F or could revise as in amie.F
! 4/13 bae:  Use defaults in aurora.F for now for SWMF

!  Find min/max in electric potentail and convert from kV to V
      do j=1,nmlat
        do i=1,nmlonp1
          phihm(i,j) = phihm(i,j)*1000.
        enddo
      enddo
! Get SH (j=1,nmlat0) and NH (j=nmlat0+1,nmlat)
      do ih=1,2
	vnx(ih,1) = 0.
	vnx(ih,2) = 0.  
        if (ih==1) then
         j1 = 1
         j2 = nmlat0
        else
         j1 = nmlat0+1
         j2 = nmlat
        endif     
      do j=j1,j2
        do i=1,nmlonp1
	    if (phihm(i,j) .lt. vnx(ih,1)) then
              vnx(ih,1) = phihm(i,j)
	      jnx(ih,1) = j
	      inx(ih,1) = i
            endif
	    if (phihm(i,j) .gt. vnx(ih,2)) then
              vnx(ih,2) = phihm(i,j)
	      jnx(ih,2) = j
	      inx(ih,2) = i
            endif
        enddo
      enddo

!  Calculate hemctpoten in kV from kameleon model min/max in V
	hemctpoten(ih) = (vnx(ih,2) - vnx(ih,1))*0.001
      enddo
!
! Return the average of the hemctpoten from the SH and NH in ctpoten
      ctpoten = 0.5*(hemctpoten(1)+hemctpoten(2))
!      write (6,"(1x,'kameleon: istep,mtime, CPS,N ctpoten=',i8,3i4,
!     |   3f8.1)") istep,modeltime(1:3),hemctpoten,ctpoten
!
! Calculate pfrac from theta0 in colath in advance.F

! Calculate center from min/max of epot using code from calccloc.F
!  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     
      do ih=1,2
      rhon = 90.-abs(ylatm(jnx(ih,1))*rtd)
      rhox = 90.-abs(ylatm(jnx(ih,2))*rtd)
      thetandeg = (rmlt(inx(ih,1)) - 6.)*360./24.
      thetaxdeg = (rmlt(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) )
      theta0(ih) = radcen(ih)*dtr
      enddo

      if (istep<=1) write (6,"(1x,'s,nefx HP CP pot1,p1 rad=',12e12.4)")
     |  xshmwpm2,xnhmwpm2,power,hemhp(1),hemhp(2),ctpoten,hemctpoten(1)
     |  ,hemctpoten(2),phihm(1,nmlat-2),phihm(nmlonp1,nmlat-2),radcen

      call f_timeinterp_closeandcleanupmemory(tid)

        call addfld('PHIHM',' ',' ',phihm,'mlon',1,nmlonp1,
     |    'mlat',1,nmlat,0)
        call addfld('TIEEKV',' ',' ',tieekv,'mlon',1,nmlonp1,
     |    'mlat',1,nmlat,0)
        call addfld('TIEEFX',' ',' ',tieefx,'mlon',1,nmlonp1,
     |    'mlat',1,nmlat,0)
        call addfld('GEKEV',' ',' ',gekev,'lon',1,nlonp1,
     |    'lat',0,nlatp1,0)
        call addfld('GMWPM2',' ',' ',gmwpm2,'lon',1,nlonp1,
     |    'lat',0,nlatp1,0)

!      print *, 'finished'
      end subroutine