program wei01_test
	use weimer_module
        implicit none

! 01/10: wei01_test to test out NCAR TIEGCM Weimer 2001 from wei01gcm.F
!  Compile on linux Sun: 
!   pgf90 -r8 -o wei01.exe wei01_nc_subset.F wei01_test.F
!  Compile on unix Sun: f90 -o wei01.exe wei01_bigc_subset.F wei01_test.F 
!  Run as: wei01.exe > tie_wei01_out
!      Comment out #, modules, substitute rmlat,rlonm,etc, delete ?
!   replace calccloc with wei01loc
!   output poten in dmlat and hmlt can be plotted later in matlab

	real,parameter :: re=6371.e5    !  earth radius in cm
	real :: calchp,                 !  calculated HP to cf with power input
     |          ds,                     !  surface area of earth for integration
     |          efluxds                 !  eflux over a surface area ds
	real :: guvihpkp(10),rkp(10)    !  GUVI HP(GW)=f(Kp) from all observations (X Luan); real Kp
	integer :: n,j,i
	real :: emery_hp, plevel
	integer :: kp, year, secs
!       real :: fkp, emery_hp, plevel, f107
!	integer :: kp, year, iday, secs
        real :: bzwkp(10)		!  Effective Bz with Kp (< 'real')
        real :: fac = 2.0        

	pi = 4.*atan(1.)
	dtr = pi/180.
        rtd = 180./pi


!   from Xiaoli Luan (remove dayside noise) and Rachel Miller, HAO/NCAR, Jun 2009
	rkp(1) = 0.19
	rkp(2) = 1.00
	rkp(3) = 1.99
	rkp(4) = 2.99
	rkp(5) = 3.95
	rkp(6) = 4.92
	rkp(7) = 5.92
	rkp(8) = 6.98
	rkp(9) = 7.96
	rkp(10) = 8.76
	guvihpkp(1) = 3.6
	guvihpkp(2) = 8.8
	guvihpkp(3) = 17.9
	guvihpkp(4) = 29.1
	guvihpkp(5) = 45.8
	guvihpkp(6) = 64.6
	guvihpkp(7) = 90.2
	guvihpkp(8) = 130.4
	guvihpkp(9) = 229.8
	guvihpkp(10)= 387.6

! Calculate the effective Bz for Weimer runs from Kp (see aurora TIEGCM notes)
	do kp=1,10
	 bzwkp(kp) = 0.3610 - 0.8353*rkp(kp) + 0.1119*rkp(kp)*rkp(kp)
     |     - 0.0319*rkp(kp)**3
	enddo

!  Use 1 deg mlat from 40-89, 15 min MLT (0-24) grid

	do n=1,nlat
	 dmlat(n) = 49.5 + n*deltamlat
	enddo
	do n=1,nlon
	 hmlt(n) = (n-1)*deltahmlt
	enddo

! Set IMF etc for Weimer 
	byimf = 0.
        swvel = 400.
        swden = 4.
        al = -100.
        aluse = .false.

C Set sample date and UT secs
	year = 2002
	iday = 172
        iday = 355
	secs = 0
	iyear = year
	uthr = secs/3600.

C Set f107 (100 for all GUVI)
	f107 = 100.

        write(6,"(' NH year iday f107 by,zimf= ',2i4,3f7.2)") 
     |   year,iday,f107,byimf,bzimf

	do kp=0,9
         bzimf = bzwkp(kp+1)
	 fkp = rkp(kp+1)
! Formula for potential:
! modified by LQIAN, 2007
! formula given by Wenbin based on data fitting
	ctpoten = 15.+15.*fkp + 0.8*fkp**2
!
! Formula for hemispheric power:
! Emery et al [JGR, 2008]
	emery_hp = 2.5103+9.9207*fkp-2.1825*fkp*fkp+0.3033*fkp*fkp*fkp
! Revise so power=2*Emery (fac=2 in util.F hp_from_bz_swvel)
	power = emery_hp*fac

! modified by LQIAN, 2007
! for fkp<=7: formula given by Zhang Yongliang based on TIMED/GUVI
! for fkp>7: power is 153.13 when fkp=7 from Zhang's formula, 
! assume power is 300.(based on NOAA satellites) when fkp=9 
! do linear interporation in between
            if (fkp <=7.) power = 16.82*exp(0.32*fkp)-4.86
            if (fkp > 7.) power = 153.13 + 
     |                            (fkp-7.)/(9.-7.)*(300.-153.13)
! 07/09 Miller/Luan: Cubic from all 2002-2007 GUVI data
!	plevel = 2.271 + 2.499*fkp - 0.36*fkp*fkp + 0.02372*fkp*fkp*fkp
!	power = exp(plevel/2.09)

! Revise so power=guvihpkp
	power = guvihpkp(kp+1)

! Revise so power=f(Bz,Vsw) from util.F function hp_from_bz_swvel(bz,swvel)
      if (bzimf < 0.) then
        power = 6.0 + 3.3*abs(bzimf) + (0.05 + 0.003*abs(bzimf))*
     |    (min(swvel,700.)-300.)
      else
        power = 5.0 + 0.05 * (min(swvel,700.)-300.)
      endif
      power = max(2.5,power)*fac
      
	call weimer01

! Write out for input to matlab plot
        write(6,"('Kp,power,ctpoten= ',2f7.2,f11.2)")
     |   fkp,power,ctpoten
	write (6,"('hmlt=',/,(6f12.2))") hmlt
	do n=1,nlat
	 write(6,"('mlat=',f6.1,' epot =',/,(6e12.4))") dmlat(n),
     |     (phihm(i,n),i=3,nlon+2)
 	enddo

	enddo ! End of Kp loop

	end program wei01_test