subroutine kameleon_tiegcm (istep,hemhp)

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,offc,dskofc
      use addfld_module,only: addfld
      implicit none		

!
! Args:
      integer,intent(in) :: istep
      real,intent(out) :: hemhp(2)
!
! 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)
      real :: xshmwpm2,xnhmwpm2,polekv,polekev,polemwpm2
! Define auroral output fields in TIEGCM distorted magnetic grid (same as phihm)
      real ::  tieefx(nmlonp1,nmlat),tieekv(nmlonp1,nmlat)
      real :: Re,areatot,dmlat,dlonm
      real,save :: 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
       
      real :: time_epoch
      integer :: tid
      integer :: iskp2

! 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
! 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,'istep<=1 area of hem = ',e12.4)") areatot
      endif	! if (istep<=1) then

! Redefine RMLT with new sunlons(1-nmlonp1) all the same but change each time step
        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

! 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)
      call f_timeinterp_timestrtoepoch(time, time_epoch)
      if (istep<3) write (6,"(1x,'tid time_epoch =',i5,f18.1)")
     |  tid,time_epoch

      call f_timeinterp_addtimestep(tid, cdf_file_path1)
      if (iskp2==0) call f_timeinterp_addtimestep(tid, cdf_file_path2)

C Find 'ep'
C  NOTE in TimeInterpolator.cpp, could put 3 variables in as:
C   manageMemory(tid, time_epoch, 'ep, 'eave', 'eflux', 3)
C  But fortran wrapper (managemem) has only 1 variable passed in at a time
      call f_timeinterp_managemem(tid, time_epoch, 'ep')
! 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
       mlat = ylatm(j)*rtd
       ih = 1
       if (j>nmlat0) ih = 2
        do i=1,nmlonp1  
         if (abs(mlat)<lbmlat) then
          phihm(i,j) = 0.
         else
! Obtain electric potential in kV and convert to V
          call f_timeinterp_interpolate(tid, time_epoch, 'ep', 0.0,
     |     mlat, rmlt(i), phihm(i,j))
          phihm(i,j) = phihm(i,j)*1000.
         endif	!  if (abs(mlat)<lbmlat) then,else
        enddo ! i=1,nmlonp1
      enddo ! j=2,nmlat-1

C Find 'eave'
      call f_timeinterp_managemem(tid, time_epoch, 'eave')
      do j=2,nmlat-1
       mlat = ylatm(j)*rtd
        do i=1,nmlonp1  
         if (abs(mlat)<lbmlat) then
          tieekv(i,j) = 0.5
         else
! Obtain electron mean energy in keV
         call f_timeinterp_interpolate(tid, time_epoch, 'eave', 0.0,
     |    mlat, rmlt(i), tieekv(i,j)) 
         endif	!  if (abs(mlat)<lbmlat) then,else
        enddo ! i=1,nmlonp1
      enddo ! j=2,nmlat-1

C Find 'eflux' and convert to mW/m2 if necessary
      call f_timeinterp_managemem(tid, time_epoch, 'eflux')
      do j=2,nmlat-1
       mlat = ylatm(j)*rtd
        do i=1,nmlonp1  
         if (abs(mlat)<lbmlat) then
          tieefx(i,j) = 0.
         else
! Obtain electron energy flux in w/m2 (SWMF) or in mw/m2 (ASTRA_AMIE)    
          call f_timeinterp_interpolate(tid, time_epoch, 'eflux', 0.0,
     |     mlat, rmlt(i), tieefx(i,j))
          tieefx(i,j) = tieefx(i,j)*facmw 
         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
         polekv = 0. 
         polekev = 0. 
         polemwpm2 = 0. 
         do i=1,nmlon
           polekv = polekv + phihm(i,ipp)
           polekev = polekev + tieekv(i,ipp)
           polemwpm2 = polemwpm2 + tieefx(i,ipp)
         enddo
         polekv = polekv/nmlon
         polekev = polekev/nmlon
         polemwpm2 = polemwpm2/nmlon
         do i=1,nmlonp1
           phihm(i,ip) = polekv
           tieekv(i,ip) = polekev           
           tieefx(i,ip) = polemwpm2
         enddo
         ip = nmlat
         ipp = nmlat-1
         polekv = 0. 
         polekev = 0. 
         polemwpm2 = 0. 
         do i=1,nmlon
           polekv = polekv + phihm(i,ipp)
           polekev = polekev + tieekv(i,ipp)
           polemwpm2 = polemwpm2 + tieefx(i,ipp)
         enddo
         polekv = polekv/nmlon
         polekev = polekev/nmlon
         polemwpm2 = polemwpm2/nmlon
         do i=1,nmlonp1
           phihm(i,ip) = polekv
           tieekv(i,ip) = polekev           
           tieefx(i,ip) = polemwpm2
         enddo
! TEMP
!	if (istep<3) write (6,"(1x,'polekv,kev,mwpm2 =',3e12.4)") 
!     |   polekv,polekev,polemwpm2

!  Calculate hem HP in GW in distorted magnetic grid
        hemhp(1) = 0.
        hemhp(2) = 0.
        xnhmwpm2 = 0.
        xshmwpm2 = 0.
        do j=1,nmlat0
         do i=1,nmlon
          tieefx(i,j) = tieefx(i,j) 
          tieefx(i,nmlat-j) = tieefx(i,nmlat-j) 
          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

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 gekev and gmwpm2  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 
! 8/13 bae:   and calculate theta0 in advance.F

      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)

      end subroutine