PRO tgcm_sf
; read 5 minutely sampled TGCM netCDF output 
;dtr = !pi/180.
;Missing value = 1.e36
  missing = 1.e36
; p0 is not a parameter, as it must be changed for mtgcm (see getflds)
; (for tgcm, p0=5.e-4, for mtgcm p0=1.2e-3)
;
    p0 = 5.e-4  ;microbar

    grav = 870.                  ; cm/s2
;k0 is Boltzmann constant
    k0 = 1.38e-16 
; gas constant 
    gask = 8.314e7 ;R
;Ri is Earth's radius at  0 km altitude in cm
    Ri = 6370*1.e5
    dlat = 5.*!pi/180
    dlon = 5.*!pi/180

;file 11111111111111111111111111-------------------------------------------------------------------
;files = file_search('$LINUX/SolarFlare/sfism_all_yr200331*.nc')
files = file_search('$LINUX/SolarFlare/nc2003301/sfism_preflare_2003301.nc')
nfiles = n_elements(files)
print, 'number of files = ',nfiles

mtime = intarr(3,288*nfiles)
mut = fltarr(288*nfiles)
mhpower = fltarr(288*nfiles)
msigma_ped = fltarr(72,36,29,288*nfiles);S/m
meden = fltarr(72, 36, 29, 288*nfiles)  ;cm-3
mden = fltarr(72,36,29,288*nfiles);ilev
mqjoule = fltarr(72,36,29,288*nfiles)   ; (lon,lat,lev,time)
mqjoule_integ = fltarr(72,36,288*nfiles); (lon,lat,time)
msheating = fltarr(72,36,29,288*nfiles)
;mo1 = fltarr(72,36,29,288*nfiles)
;mo2 = fltarr(72,36,29,288*nfiles)
mtn = fltarr(72,36,29,288*nfiles)
;mscht = fltarr(72,36,29,288*nfiles)
mzg = fltarr(72,36,29,288*nfiles);ilev
mz = fltarr(72,36,29,288*nfiles);ilev
mtec = fltarr(72,36,288*nfiles)
mqop2pa = fltarr(72,36,29,288*nfiles)   ; production rate (cm-3s-1)
mqop2da = fltarr(72,36,29,288*nfiles)   ; (lon,lat,lev,time)
mqo2p = fltarr(72,36,29,288*nfiles)   ; (lon,lat,lev,time)
mqn2p = fltarr(72,36,29,288*nfiles)   ; (lon,lat,lev,time)
mqnp = fltarr(72,36,29,288*nfiles)   ; (lon,lat,lev,time)
mqopa = fltarr(72,36,29,288*nfiles)   ; (lon,lat,lev,time)
mqnop = fltarr(72,36,29,288*nfiles)   ; (lon,lat,lev,time)


for i = 0L,nfiles-1 do begin
print, 'i=',i
filein = files(i)
;; open and read netCDF file
    print, 'Reading ',filein
    id=ncdf_open(filein)
    ncdf_varget,id,'mtime',var
    mtime(*,i*288:(i*288+287)) = var
    ncdf_varget,id,'lat',lat   ;(36),[-87.5,87.5]
    ncdf_varget,id,'lon',lon   ;(72),[-180., 175.]
    ncdf_varget,id,'lev',lev   ;
    ncdf_varget,id,'ilev',ilev  ;Z,ZG,DEN
;    ncdf_varget,id,ncdf_varid(id,'CO2_COOL'),vco2 ;erg/g/s
;mco2cool(*,*,*,i*288:(i*288+287)) = vco2
;    ncdf_varget,id,ncdf_varid(id,'NO_COOL'),vno   ;erg/g/s
;mnocool(*,*,*,i*288:(i*288+287)) = vno
    ncdf_varget,id,ncdf_varid(id,'ut'),vut ;universal time (hours)
    mut(i*288:(i*288+287)) = vut
    ncdf_varget,id,ncdf_varid(id,'hpower'),vhpower ;hemispheric power, GW
    mhpower(i*288:(i*288+287)) = vhpower
    ncdf_varget,id,ncdf_varid(id,'SIGMA_PED'),vped ;pederson conductivity
    msigma_ped(*,*,*,i*288:(i*288+287)) = vped
    ncdf_varget,id,ncdf_varid(id,'NE'),veden ;electron density, /cm3
    meden(*,*,*,i*288:(i*288+287)) = veden
    ncdf_varget,id,ncdf_varid(id,'DEN'),vden ;total density, g/cm3
    mden(*,*,*,i*288:(i*288+287)) = vden
    ncdf_varget,id,ncdf_varid(id,'QJOULE'),vqj
    mqjoule(*,*,*,i*288:(i*288+287)) = vqj
    ncdf_varget,id,ncdf_varid(id,'QJOULE_INTEG'),vqj_int ;erg/cm2/s = mW/m2
    mqjoule_integ(*,*,i*288:(i*288+287)) = vqj_int
    ncdf_varget,id,ncdf_varid(id,'SHEATING'),vheat ;solar heating,erg/g/s
    msheating(*,*,*,i*288:(i*288+287)) = vheat
    ;ncdf_varget,id,ncdf_varid(id,'QTOTAL'),qtotal
;    ncdf_varget,id,ncdf_varid(id,'O1'),vo1
;mo1(*,*,*,i*288:(i*288+287)) = vo1
;    ncdf_varget,id,ncdf_varid(id,'O2'),vo2
;mo2(*,*,*,i*288:(i*288+287)) = vo2
    ncdf_varget,id,ncdf_varid(id,'TN'),vtn
    mtn(*,*,*,i*288:(i*288+287)) = vtn
    ncdf_varget,id,ncdf_varid(id,'grav'),grav ;cm/s
;    ncdf_varget,id,ncdf_varid(id,'SCHT'),vsh  ;scale height, km
;mscht(*,*,*,i*288:(i*288+287)) = vsh
    ncdf_varget,id,ncdf_varid(id,'ZG'),vzg
    mzg(*,*,*,i*288:(i*288+287)) = vzg
    ncdf_varget,id,ncdf_varid(id,'Z'),vz
    mz(*,*,*,i*288:(i*288+287)) = vz
    ncdf_varget,id,ncdf_varid(id,'TEC'),vtec
    mtec(*,*,i*288:(i*288+287)) = vtec ;1/cm2
    ncdf_varget,id,ncdf_varid(id,'QOP2Pa'),vqop2pa
    mqop2pa(*,*,*,i*288:(i*288+287)) = vqop2pa
    ncdf_varget,id,ncdf_varid(id,'QOP2Da'),vqop2da
    mqop2da(*,*,*,i*288:(i*288+287)) = vqop2da
    ncdf_varget,id,ncdf_varid(id,'QO2P'),vqo2p
    mqo2p(*,*,*,i*288:(i*288+287)) = vqo2p
    ncdf_varget,id,ncdf_varid(id,'QN2P'),vqn2p
    mqn2p(*,*,*,i*288:(i*288+287)) = vqn2p
    ncdf_varget,id,ncdf_varid(id,'QNP'),vqnp
    mqnp(*,*,*,i*288:(i*288+287)) = vqnp
    ncdf_varget,id,ncdf_varid(id,'QOPa'),vqopa
    mqopa(*,*,*,i*288:(i*288+287)) = vqopa
    ncdf_varget,id,ncdf_varid(id,'QNOP'),vqnop
    mqnop(*,*,*,i*288:(i*288+287)) = vqnop

    print,'Reading Finished~~~~~~~~~~~~~~~'
endfor

      mtime = transpose(mtime)
      ntimes = n_elements(mtime(*,0))  ;# of hours
      smtime = fltarr(ntimes)
      ;for i = 0, ntimes-1 do begin
          smtime(*) = mtime(*,0)+mtime(*,1)/24.+mtime(*,2)/1440.
          ;smtime(i) = mtime(i*8,0)    ;daily average
      ;endfor
      
      nlat = n_elements(lat)
      nlon = n_elements(lon)
      nlev = n_elements(lev)

;fp = where(mco2cool ge missing)
;mco2cool(fp) = !Values.F_NAN ;NaN
;fp = where(mnocool ge missing)
;mnocool(fp) = !Values.F_NAN ;NaN
fp = where(mhpower ge missing)
mhpower(fp) = !Values.F_NAN ;NaN
fp = where(msigma_ped ge missing)
msigma_ped(fp) = !Values.F_NAN ;NaN
fp = where(meden ge missing)
meden(fp) = !Values.F_NAN ;NaN
fp = where(mden ge missing)
mden(fp) = !Values.F_NAN ;NaN
fp = where(mqjoule ge missing)
mqjoule(fp) = !Values.F_NAN ;NaN
fp = where(mqjoule_integ ge missing)
mqjoule_integ(fp) = !Values.F_NAN ;NaN
fp = where(msheating ge missing)
msheating(fp) = !Values.F_NAN ;NaN       
;fp = where(mo1 ge missing)
;mo1(fp) = !Values.F_NAN ;NaN
;fp = where(mo2 ge missing)
;mo2(fp) = !Values.F_NAN ;NaN
fp = where(mtn ge missing)
mtn(fp) = !Values.F_NAN ;NaN
;fp = where(mscht ge missing)
;mscht(fp) = !Values.F_NAN ;NaN
fp = where(mzg ge missing)
mzg(fp) = !Values.F_NAN ;NaN
fp = where(mz ge missing)
mz(fp) = !Values.F_NAN ;NaN
fp = where(mtec ge missing)
mtec(fp) = !Values.F_NAN ;NaN
fp = where(mqop2pa ge missing)
mqop2pa(fp) = !Values.F_NAN ;NaN
fp = where(mqop2da ge missing)
mqop2da(fp) = !Values.F_NAN ;NaN
fp = where(mqo2p ge missing)
mqo2p(fp) = !Values.F_NAN ;NaN
fp = where(mqn2p ge missing)
mqn2p(fp) = !Values.F_NAN ;NaN
fp = where(mqnp ge missing)
mqnp(fp) = !Values.F_NAN ;NaN
fp = where(mqopa ge missing)
mqopa(fp) = !Values.F_NAN 
fp = where(mqnop ge missing)
mqnop(fp) = !Values.F_NAN ;NaN
tpr = mqop2pa+mqop2da+mqo2p+mqn2p+mqnp+mqopa+mqnop ;total production rate

gzg_400 = fltarr(ntimes)
gden_400 = fltarr(ntimes)
geden_400 = fltarr(ntimes)
gtn_400 = fltarr(ntimes)
gsheating_400=fltarr(ntimes)
gjheating_400=fltarr(ntimes)
tn_400 = fltarr(nlon,nlat,ntimes)
zg_400 = fltarr(nlon,nlat,ntimes)
den_400 = fltarr(nlon,nlat,ntimes)
eden_400 = fltarr(nlon,nlat,ntimes)
sheating_400=fltarr(nlon,nlat,ntimes)
jheating_400=fltarr(nlon,nlat,ntimes)
tpr_400 = fltarr(nlon,nlat,ntimes)
gtpr_400 = fltarr(ntimes)
sigma_ped_400 = fltarr(nlon,nlat,ntimes)
gsigma_ped_400 = fltarr(ntimes)

gtpr = fltarr(ntimes,nlev)
gtn = fltarr(ntimes,nlev)  ;global mean
gden = fltarr(ntimes,nlev) ;global mean
geden = fltarr(ntimes,nlev)
gsheating = fltarr(ntimes,nlev)
gjheating = fltarr(ntimes,nlev)
gsigma_ped = fltarr(ntimes, nlev)
gzg = fltarr(ntimes, nlev)
gcm_sh = fltarr(nlon,nlat)
gcm_jh = fltarr(nlon,nlat)

for n=0, ntimes-1 do begin
   for k=0, nlev-1 do begin
      gtpr(n,k) = mean(tpr(*,*,k,n),/nan)
      gtn(n,k) = mean(mtn(*,*,k,n),/nan)            ;K
      gden(n,k) = mean(mden(*,*,k,n),/nan)
      geden(n,k) = mean(meden(*,*,k,n),/nan)
      gcm_sh = msheating(*,*,k,n)*mden(*,*,k,n)
      gsheating(n,k) = mean(gcm_sh,/nan)*100.       ;mW/m3
      ;gsheating(n,k) = mean(msheating(*,*,k,n),/nan)
      gcm_jh = mqjoule(*,*,k,n)*mden(*,*,k,n)
      gjheating(n,k) = mean(gcm_jh,/nan)*100.       ;mW/m3 
      ;gjheating(n,k) = mean(mqjoule(*,*,k,n),/nan)
      gsigma_ped(n,k) = mean(msigma_ped(*,*,k,n),/nan)
      gzg(n,k) = mean(mzg(*,*,k,n),/nan)*1.e-5 ;Km
   endfor ;k:lev
   for i = 0, nlon-1 do begin
      for j=0, nlat-1 do begin
         mind = min(abs(mzg(i,j,*,n)-4.e7))
         ffp = where(abs(mzg(i,j,*,n)-4.e7) eq mind)
         zg_400(i,j,n) = mzg(i,j,ffp,n)
         den_400(i,j,n) = mden(i,j,ffp,n)
         eden_400(i,j,n) = meden(i,j,ffp,n)
         tn_400(i,j,n) = mtn(i,j,ffp,n)
         tpr_400(i,j,n) = tpr(i,j,ffp,n)
         sigma_ped_400(i,j,n) = msigma_ped(i,j,ffp,n)
         sheating_400(i,j,n) = msheating(i,j,ffp,n)*den_400(i,j,n)*100. ;mW/m3
         jheating_400(i,j,n) = mqjoule(i,j,ffp,n)*den_400(i,j,n)*100.   ;mW/m3
      endfor                    ;j:lat
   endfor                       ;i:lon
   gzg_400(n) = mean(zg_400(*,*,n),/nan)
   gden_400(n) = mean(den_400(*,*,n),/nan)
   geden_400(n) = mean(eden_400(*,*,n),/nan)
   gtpr_400(n) = mean(tpr_400(*,*,n),/nan)
   gsigma_ped_400(n) = mean(sigma_ped_400(*,*,n),/nan)
   gtn_400(n) = mean(tn_400(*,*,n),/nan)
   gden_400(n) = mean(den_400(*,*,n),/nan)
   gsheating_400(n) = mean(sheating_400(*,*,n),/nan)   ;mW/m3
   gjheating_400(n) = mean(jheating_400(*,*,n),/nan)   ;mW/m3
endfor                          ;ntimes

;openw,13,'gden_400_time_preflare.dat',width=160
;printf,13,'smtime(day) gsheating(mW/m3) gjheating(mW/m3) gden_400(g/cm3) gtn_400(K) gped_400(S/m) gtpr_400(cm-3s-1) geden_400(cm-3)'
;for n = 0, n_elements(smtime)-1 do begin
;printf,13,smtime(n),gsheating_400(n),gjheating_400(n),gden_400(n),gtn_400(n),gsigma_ped_400(n),gtpr_400(n),geden_400(n)
;endfor
;close,13
;stop

int_jh = fltarr(ntimes) & int_jh(*)=0.
int_sh = fltarr(ntimes) & int_sh(*)=0.
height_sh = fltarr(nlon, nlat,ntimes) & height_sh(*,*,*) =0. ;height integrated solar heating
for n=0,ntimes-1 do begin
   for i = 0, nlon-1 do begin
      for j=0, nlat-1 do begin
         for k = 1, nlev-2 do begin
            dmsh = msheating(i,j,k,n)*mden(i,j,k,n)*(mZ(i,j,k+1,n)-mZ(i,j,k-1,n))/2.
            height_sh(i,j,n) = height_sh(i,j,n)+dmsh
            int_sh(n) = int_sh(n) + dmsh*cos(lat(j)*!pi/180.)             ;mW/m2
         endfor                 ;k:lev
         int_jh(n) = int_jh(n)+mqjoule_integ(i,j,n)*cos(lat(j)*!pi/180.)  ;mW/m2
      endfor                    ;j:lat
   endfor                       ;i:lon
endfor                          ;n:ntimes

; Z was in km, and convert from ergs/s*1e5 to GW by a factor of 1.e-11
area = (Ri+mZ(0,0,nlev-1,ntimes-1))^2*dlat*dlon
int_sh = int_sh*area*1.e-16     ;GW
int_jh = int_jh*area*1.e-16     ;GW
     


;file 22222222222222222222222222222222222--------------------------------------------------------------
files2 = file_search('$LINUX/SolarFlare/nc2003301/sfism_0_195nm_2003301.nc')
fn = '2003301_0_195nm'
nfiles2 = n_elements(files2)
print, 'number of files = ',nfiles2

mtime2 = intarr(3,288*nfiles2)
mut2 = fltarr(288*nfiles2)
mhpower2 = fltarr(288*nfiles2)
msigma_ped2 = fltarr(72,36,29,288*nfiles2);S/m
meden2 = fltarr(72,36,29,288*nfiles2);ilev
mden2 = fltarr(72,36,29,288*nfiles2);ilev
mqjoule2 = fltarr(72,36,29,288*nfiles2)   ; (lon.lat,lev,time)
mqjoule_integ2 = fltarr(72,36,288*nfiles2); (lon,lat,time)
msheating2 = fltarr(72,36,29,288*nfiles2)
mtn2 = fltarr(72,36,29,288*nfiles2)
mzg2 = fltarr(72,36,29,288*nfiles2);ilev
mz2 = fltarr(72,36,29,288*nfiles2);ilev
mtec2 = fltarr(72,36,288*nfiles2)
mqop2pa2 = fltarr(72,36,29,288*nfiles2)   ; production rate (cm-3s-1)
mqop2da2 = fltarr(72,36,29,288*nfiles2)   ; (lon,lat,lev,time)
mqo2p2 = fltarr(72,36,29,288*nfiles2)   ; (lon,lat,lev,time)
mqn2p2 = fltarr(72,36,29,288*nfiles2)   ; (lon,lat,lev,time)
mqnp2 = fltarr(72,36,29,288*nfiles2)   ; (lon,lat,lev,time)
mqopa2 = fltarr(72,36,29,288*nfiles2)   ; (lon,lat,lev,time)
mqnop2 = fltarr(72,36,29,288*nfiles2)   ; (lon,lat,lev,time)

for i = 0L,nfiles2-1 do begin
print, 'i=',i
filein2 = files2(i)
;; open and read netCDF file
    print, 'Reading ',filein2
    id=ncdf_open(filein2)
    ncdf_varget,id,'mtime',var
    mtime2(*,i*288:(i*288+287)) = var
    ncdf_varget,id,ncdf_varid(id,'ut'),vut ;universal time (hours)
    mut2(i*288:(i*288+287)) = vut
    ncdf_varget,id,ncdf_varid(id,'hpower'),vhpower ;hemispheric power, GW
    mhpower2(i*288:(i*288+287)) = vhpower
    ncdf_varget,id,ncdf_varid(id,'SIGMA_PED'),vped ;pederson conductivity
    msigma_ped2(*,*,*,i*288:(i*288+287)) = vped
    ncdf_varget,id,ncdf_varid(id,'NE'),veden ;total density, g/cm3
    meden2(*,*,*,i*288:(i*288+287)) = veden
    ncdf_varget,id,ncdf_varid(id,'DEN'),vden ;total density, g/cm3
    mden2(*,*,*,i*288:(i*288+287)) = vden
    ncdf_varget,id,ncdf_varid(id,'QJOULE'),vqj
    mqjoule2(*,*,*,i*288:(i*288+287)) = vqj
    ncdf_varget,id,ncdf_varid(id,'QJOULE_INTEG'),vqj_int    ;erg/cm2/s = mW/m2
    mqjoule_integ2(*,*,i*288:(i*288+287)) = vqj_int
    ncdf_varget,id,ncdf_varid(id,'SHEATING'),vheat  ;solar heating,erg/g/s
    msheating2(*,*,*,i*288:(i*288+287)) = vheat
    ncdf_varget,id,ncdf_varid(id,'TN'),vtn
    mtn2(*,*,*,i*288:(i*288+287)) = vtn
    ncdf_varget,id,ncdf_varid(id,'ZG'),vzg
    mzg2(*,*,*,i*288:(i*288+287)) = vzg
    ncdf_varget,id,ncdf_varid(id,'Z'),vz
    mz2(*,*,*,i*288:(i*288+287)) = vz
    ncdf_varget,id,ncdf_varid(id,'TEC'),vtec
    mtec2(*,*,i*288:(i*288+287)) = vtec ;1/cm2
    ncdf_varget,id,ncdf_varid(id,'QOP2Pa'),vqop2pa
    mqop2pa2(*,*,*,i*288:(i*288+287)) = vqop2pa
    ncdf_varget,id,ncdf_varid(id,'QOP2Da'),vqop2da
    mqop2da2(*,*,*,i*288:(i*288+287)) = vqop2da
    ncdf_varget,id,ncdf_varid(id,'QO2P'),vqo2p
    mqo2p2(*,*,*,i*288:(i*288+287)) = vqo2p
    ncdf_varget,id,ncdf_varid(id,'QN2P'),vqn2p
    mqn2p2(*,*,*,i*288:(i*288+287)) = vqn2p
    ncdf_varget,id,ncdf_varid(id,'QNP'),vqnp
    mqnp2(*,*,*,i*288:(i*288+287)) = vqnp
    ncdf_varget,id,ncdf_varid(id,'QOPa'),vqopa
    mqopa2(*,*,*,i*288:(i*288+287)) = vqopa
    ncdf_varget,id,ncdf_varid(id,'QNOP'),vqnop
    mqnop2(*,*,*,i*288:(i*288+287)) = vqnop

    print,'Reading Finished~~~~~~~~~~~~~~~'
endfor

      mtime2 = transpose(mtime2)
      ntimes2 = n_elements(mtime2(*,0))  ;# of hours
      smtime2 = fltarr(ntimes2)
      smtime2(*) = mtime2(*,0)+mtime2(*,1)/24.+mtime2(*,2)/1440.
      
fp = where(mhpower2 ge missing)
mhpower2(fp) = !Values.F_NAN ;NaN
fp = where(msigma_ped2 ge missing)
msigma_ped2(fp) = !Values.F_NAN ;NaN
fp = where(meden2 ge missing)
meden2(fp) = !Values.F_NAN ;NaN
fp = where(mden2 ge missing)
mden2(fp) = !Values.F_NAN ;NaN
fp = where(mqjoule2 ge missing)
mqjoule2(fp) = !Values.F_NAN ;NaN
fp = where(mqjoule_integ2 ge missing)
mqjoule_integ2(fp) = !Values.F_NAN ;NaN
fp = where(msheating2 ge missing)
msheating2(fp) = !Values.F_NAN ;NaN       
fp = where(mtn2 ge missing)
mtn2(fp) = !Values.F_NAN ;NaN
fp = where(mzg2 ge missing)
mzg2(fp) = !Values.F_NAN ;NaN
fp = where(mz2 ge missing)
mz2(fp) = !Values.F_NAN ;NaN
fp = where(mtec2 ge missing)
mtec2(fp) = !Values.F_NAN ;NaN
fp = where(mqop2pa2 ge missing)
mqop2pa2(fp) = !Values.F_NAN ;NaN
fp = where(mqop2da2 ge missing)
mqop2da2(fp) = !Values.F_NAN ;NaN
fp = where(mqo2p2 ge missing)
mqo2p2(fp) = !Values.F_NAN ;NaN
fp = where(mqn2p2 ge missing)
mqn2p2(fp) = !Values.F_NAN ;NaN
fp = where(mqnp2 ge missing)
mqnp2(fp) = !Values.F_NAN ;NaN
fp = where(mqopa2 ge missing)
mqopa2(fp) = !Values.F_NAN 
fp = where(mqnop2 ge missing)
mqnop2(fp) = !Values.F_NAN ;NaN
tpr2 = mqop2pa2+mqop2da2+mqo2p2+mqn2p2+mqnp2+mqopa2+mqnop2 ;total production rate

gzg2_400 = fltarr(ntimes2)
gden2_400 = fltarr(ntimes2)
geden2_400 = fltarr(ntimes2)
zg2_400 = fltarr(nlon,nlat,ntimes2)
den2_400 = fltarr(nlon,nlat,ntimes2)
eden2_400 = fltarr(nlon,nlat,ntimes2)
gtn2_400 = fltarr(ntimes2)
gsheating2_400=fltarr(ntimes2)
gjheating2_400=fltarr(ntimes2)
tn2_400 = fltarr(nlon,nlat,ntimes2)
sheating2_400=fltarr(nlon,nlat,ntimes2)
jheating2_400=fltarr(nlon,nlat,ntimes2)
tpr2_400 = fltarr(nlon,nlat,ntimes2)
gtpr2_400 = fltarr(ntimes2)
sigma_ped2_400 = fltarr(nlon,nlat,ntimes2)
gsigma_ped2_400 = fltarr(ntimes2)

gtpr2 = fltarr(ntimes2,nlev)
gtn2 = fltarr(ntimes2,nlev)  ;global mean
gden2 = fltarr(ntimes2,nlev) ;global mean
geden2 = fltarr(ntimes2,nlev) ;global mean
gsheating2 = fltarr(ntimes2,nlev)
gjheating2 = fltarr(ntimes2,nlev)
gsigma_ped2 = fltarr(ntimes2, nlev)
gzg2 = fltarr(ntimes2, nlev)
gcm_sh2 = fltarr(nlon,nlat)
gcm_jh2 = fltarr(nlon,nlat)

for n=0, ntimes2-1 do begin
   for k=0, nlev-1 do begin
      gtpr2(n,k) = mean(tpr2(*,*,k,n),/nan)           
      gtn2(n,k) = mean(mtn2(*,*,k,n),/nan)            ;K
      gden2(n,k) = mean(mden2(*,*,k,n),/nan)
      geden2(n,k) = mean(meden2(*,*,k,n),/nan)
      gcm_sh2 = msheating2(*,*,k,n)*mden2(*,*,k,n)
      gsheating2(n,k) = mean(gcm_sh2,/nan)*100.       ;mW/m3
      ;gsheating2(n,k) = mean(msheating2(*,*,k,n),/nan)
      gcm_jh2 = mqjoule2(*,*,k,n)*mden2(*,*,k,n)
      gjheating2(n,k) = mean(gcm_jh2,/nan)*100.       ;mW/m3 
      ;gjheating(n,k) = mean(mqjoule2(*,*,k,n),/nan)
      gsigma_ped2(n,k) = mean(msigma_ped2(*,*,k,n),/nan)
      gzg2(n,k) = mean(mzg2(*,*,k,n),/nan)*1.e-5 ;Km
   endfor                       ;k:lev
   for i = 0, nlon-1 do begin
      for j=0, nlat-1 do begin
         mind = min(abs(mzg2(i,j,*,n)-4.e7))
         ffp = where(abs(mzg2(i,j,*,n)-4.e7) eq mind)
         zg2_400(i,j,n) = mzg2(i,j,ffp,n)
         den2_400(i,j,n) = mden2(i,j,ffp,n)
         eden2_400(i,j,n) = meden2(i,j,ffp,n)
         tn2_400(i,j,n) = mtn2(i,j,ffp,n)
         tpr2_400(i,j,n) = tpr2(i,j,ffp,n)
         sigma_ped2_400(i,j,n) = msigma_ped2(i,j,ffp,n)
         sheating2_400(i,j,n) = msheating2(i,j,ffp,n)*den2_400(i,j,n)*100. ;mW/m3
         jheating2_400(i,j,n) = mqjoule2(i,j,ffp,n)*den2_400(i,j,n)*100.   ;mW/m3
;if ((i eq 0) and (j eq 1) and (n eq 0)) then print, 'mind =',mind, 'ffp=',ffp,zg2_400(i,j,n) 
      endfor                    ;j:lat
   endfor                       ;i:lon
   gzg2_400(n) = mean(zg2_400(*,*,n),/nan)
   gden2_400(n) = mean(den2_400(*,*,n),/nan)
   geden2_400(n) = mean(eden2_400(*,*,n),/nan)
   gtpr2_400(n) = mean(tpr2_400(*,*,n),/nan)
   gsigma_ped2_400(n) = mean(sigma_ped2_400(*,*,n),/nan)
   gtn2_400(n) = mean(tn2_400(*,*,n),/nan)
   gsheating2_400(n) = mean(sheating2_400(*,*,n),/nan)   ;mW/m3
   gjheating2_400(n) = mean(jheating2_400(*,*,n),/nan)   ;mW/m3
endfor ;ntimes

gazg2 = fltarr(nlev)
for k = 0, nlev-1 do begin
gazg2(k) = mean(mzg2(*,*,k,*),/nan)
endfor

int_jh2 = fltarr(ntimes2) & int_jh2(*)=0.
int_sh2 = fltarr(ntimes2) & int_sh2(*)=0.
height_sh2 = fltarr(nlon, nlat,ntimes2) & height_sh2(*,*,*) =0. ;height integrated solar heating
for n=0,ntimes2-1 do begin
   for i = 0, nlon-1 do begin
      for j=0, nlat-1 do begin
         for k = 1, nlev-2 do begin
            dmsh = msheating2(i,j,k,n)*mden2(i,j,k,n)*(mZ2(i,j,k+1,n)-mZ2(i,j,k-1,n))/2.
            height_sh2(i,j,n) = height_sh2(i,j,n)+dmsh
            int_sh2(n) = int_sh2(n) + dmsh*cos(lat(j)*!pi/180.)
         endfor                 ;k:lev
         int_jh2(n) = int_jh2(n)+mqjoule_integ2(i,j,n)*cos(lat(j)*!pi/180.)
      endfor                    ;j:lat
   endfor                       ;i:lon
endfor                          ;n:ntimes

; Z was in km, and convert from ergs/s*1e5 to GW by a factor of 1.e-11
area = (Ri+mZ2(0,0,nlev-1,ntimes-1))^2*dlat*dlon
int_sh2 = int_sh2*area*1.e-16     ;GW
int_jh2 = int_jh2*area*1.e-16

;difference--------------------------------------------------------
dgjheating=gjheating2-gjheating
dgsheating=gsheating2-gsheating
dgeden=geden2-geden
dgden=gden2-gden
dgtn=gtn2-gtn
dgped=gsigma_ped2-gsigma_ped
dgtpr = gtpr2-gtpr

dgtpr_400 = gtpr2_400-gtpr_400
dgped_400 = gsigma_ped2_400-gsigma_ped_400
dgden_400 = gden2_400-gden_400
dgeden_400 = geden2_400-geden_400
dgtn_400 = gtn2_400-gtn_400
dgsheating_400 = gsheating2_400-gsheating_400
dgjheating_400 = gjheating2_400-gjheating_400

;openw,13,'dgeden_400_time'+fn+'.dat',width=160
;printf,13,'smtime(day) dgsheating_400(mW/m3) dgjheating_400(mW/m3) dgden_400(g/cm3) dgtn_400(K) dgped_400(S/m) dgtpr_400(cm-3s-1) dgeden_400(cm-3)'
;for n = 0, n_elements(smtime)-1 do begin
;printf,13,smtime(n),dgsheating_400(n),dgjheating_400(n),dgden_400(n),dgtn_400(n),dgped_400(n),dgtpr_400(n),dgeden_400(n)
;endfor
;close,13

;stop
;maximum perturbation in electron density: dgeden(288, 29)--------------------------------------
max_eden = max(dgeden,i,/nan)
max_eden_smtime = i mod 288
max_eden_lev = i/288
max_ped = max(dgped,i,/nan)
max_ped_smtime = i mod 288
max_jh = max(dgjheating,i,/nan)
max_jh_smtime = i mod 288
max_sh = max(dgsheating,i,/nan)
max_sh_smtime = i mod 288
print, smtime(max_ped_smtime),smtime(max_ped_smtime),smtime(max_jh_smtime),smtime(max_sh_smtime)
!p.multi = [0,1,2]
plot, dgeden(max_eden_smtime,*),lev
plot, dgped(max_ped_smtime,*),lev
!p.multi =0

openw,33,'dgeden_lev'+fn+'.dat',width = 160
printf,33,'lev zg(km) dgeden(cm-3)  dgped  dgjheating  dgsheating  dgtpr(cm-3s-1)'
for n = 0, 28 do begin
printf,33,lev(n),gazg2(n),dgeden(max_sh_smtime,n),dgped(max_sh_smtime,n),dgjheating(max_sh_smtime,n),dgsheating(max_sh_smtime,n),dgtpr(max_sh_smtime, n)
endfor
close,33

stop
;--------------------------------------
dsh = int_sh2-int_sh
djh = int_jh2-int_jh
fp = finite(dsh,/nan)  ; return1 for NaN
dsh(where(fp eq 1)) = 0.
fp = finite(djh,/nan)  ; return1 for NaN
djh(where(fp eq 1)) = 0.
dt = 5./60./24. ; 5-minute resolution
total_sh = 0.
total_jh = 0.
for i = 0, n_elements(dsh)-2 do begin
total_sh = total_sh+0.5*dt*(dsh(i)+dsh(i+1))
;print, total_sh
total_jh = total_jh+0.5*dt*(djh(i)+djh(i+1))
endfor
total_sh = total_sh*1.e9 ;J
total_jh = total_jh*1.e9 ;J


;openw,3,'Energyinputs_'+fn+'.dat'
;;printf,3,'time   dsh(GW)  djh(GW)'
;for i = 0, n_elements(smtime)-1 do begin
;printf, 3, smtime(i), dsh(i), djh(i)
;endfor
;close, 3

stop
!p.multi = [0,1,2]
set_plot,'ps'
device,filename=fn+'_3.eps', /color;/ENCAPSUL, BITS_PER_PIXEL=8;, /COLOR
xrange = [301,302]
conmake,dgsheating,smtime,lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title = 'GSH(mW/m3)'+fn;, overp=dgsheating 
plot, smtime, dgsheating_400,xran=xrange,ytitle='GSH at 400 km(mW/m3)',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
conmake,dgjheating,smtime,lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title = 'GJH(mw/m3)'+fn;, overp=dgjheating
plot, smtime, dgjheating_400,xran=xrange,ytitle='GJH at 400 km(mW/m3)',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
conmake,dgden,smtime,lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title= 'Gdensity (g/cm3)'+fn;,minlev=-2.0e-14,maxlev =2.0e-14;, overp=dgden
plot, smtime, dgden_400,xran=xrange,ytitle='Gdensity at 400 km(g/cm3)',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
conmake,dgtn,smtime,lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title = 'Gtemperature(K)'+fn;, overp=dgtn
plot, smtime, dgtn_400,xran=xrange,ytitle='Gtemperature at 400 km(K)',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
conmake,dgped, smtime, lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title = 'GPederson(S/m)'+fn;, overp=dgtn
plot, smtime, dgped_400,xran=xrange,ytitle='GPederson 400 km(S/m)',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
conmake,dgtpr, smtime, lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title = 'GProduction rate(cm-3s-1)'+fn
plot, smtime, dgtpr_400,xran=xrange,ytitle='GProduction rate at 400 km(cm-3s-1)',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
conmake,dgeden, smtime, lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title = 'GElectron density(cm-3)'+fn
plot, smtime, dgeden_400,xran=xrange,ytitle='GElectron density at 400 km(cm-3)',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
device,/close
set_plot,'x'

!p.multi = [0,1,2]
set_plot,'ps'
device,filename=fn+'_4.eps', /color;/ENCAPSUL, BITS_PER_PIXEL=8;, /COLOR
xrange = [301,302]
conmake,dgsheating/gsheating*100.,smtime,lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title = 'GSH %increase'+fn;, overp=dgsheating 
plot, smtime, dgsheating_400/gsheating_400*100.,xran=xrange,ytitle='GSH %increase at 400 km',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
conmake,dgjheating/gjheating*100.,smtime,lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title = 'GJH %increase'+fn,maxlev =30.;, overp=dgjheating
plot, smtime, dgjheating_400/gjheating_400*100.,xran=xrange,ytitle='GJH %increase at 400 km',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
conmake,dgden/gden*100.,smtime,lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title= 'Gdensity %increase'+fn;,minlev=-2.0e-14,maxlev =2.0e-14;, overp=dgden
plot, smtime, dgden_400/gden_400*100.,xran=xrange,ytitle='Gdensity %increase at 400 km',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
conmake,dgtn/gtn*100.,smtime,lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title = 'Gtemperature %increase'+fn;, overp=dgtn
plot, smtime, dgtn_400/gtn_400*100.,xran=xrange,ytitle='Gtemperature %increase at 400 km',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
conmake,dgped/gsigma_ped*100., smtime, lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title = 'GPederson %increase'+fn,maxlev=30.;, overp=dgtn
plot, smtime, dgped_400/gsigma_ped_400*100.,xran=xrange,ytitle='GPederson %increase 400 km',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
conmake,dgtpr/gtpr*100., smtime, lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title = 'GProduction rate %increase'+fn
plot, smtime, dgtpr_400/gtpr_400*100.,xran=xrange,ytitle='GProduction rate %increase at 400 km',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
conmake,dgeden/geden*100., smtime, lev,256,/cbar,xlabel=5,ylabel=2,xran=xrange,yran=[-6.25, 6.25],title = 'GNe %increase'+fn
plot, smtime, dgeden_400/geden_400*100.,xran=xrange,ytitle='GNe %increase at 400 km',position=[[0.14,0.08],[0.82,0.42]]
oplot, [250,311],[0,0],linestyle=2 
device,/close
set_plot,'x'

stop
;total integration------------------------------------
;xrange = [301.,302.]
!p.multi=[0,1,2]
set_plot,'ps'
device,filename='Energyinputs_'+fn+'_3.eps',/portrait
;plot, smtime2,int_sh2,ytitle='solar flux (GW)',title = fn
;oplot, smtime,int_sh,linestyle=2
plot, smtime, dsh,ytitle='dSH (GW)',title=fn,xran=xrange
xyouts, 301.7, max(dsh,/nan)/2.,string(total_sh)+'J'
plot, smtime, djh,ytitle='dJH (GW)',xran=xrange
xyouts, 301.7, max(djh,/nan)/2., string(total_jh)+'J'
device,/close
set_plot,'x'

stop
;polar contour----------------------------------------------------------
!p.multi= 0
lon2 = fltarr(nlon+1) & lon2(0:71) = lon & lon2(72) = 180.
mqjoule_integ_t = fltarr(nlon+1,nlat,ntimes) & mqjoule_integ_t(0:71,*,*)= mqjoule_integ & mqjoule_integ_t(72,*,*)=mqjoule_integ(0,*,*)

set_plot,'ps'
device,filename='polar.eps',/color
nl =256
LOADCT, 33,ncolors = nl
DEVICE, DECOMPOSED = 0
map_set,90,-100,0,/orthographic,/isotropic,/grid,/continents,/horizon,e_horizon={fill:1},limit=[50.,-180,90.,180]
contour,mqjoule_integ_t(*,18:35,150),lon2,lat(18:35),$
        nlevels=nl,C_COLORS = INDGEN(nl),xstyle=1,ystyle=1, $
        /overplot,/fill,/isotropic,POSITION=[ 0.14, 0.08, 0.82, 0.92 ]
;,position=position,zrange=[0.,2.6],zstyle= 1,c_colors = c_colors
map_grid,/label,latlab = 175,color=130
map_continents,color=130
;color bar
mxx = MAX( mqjoule_integ_t, MIN = mnn ,/nan) * 0.98
delta_lev = ( mxx - mnn ) / FLOAT( nl )
level = FINDGEN( nl ) * delta_lev + mnn
;colorbar array consists of levels previously determined
cbar_arr = FLTARR( 2, N_ELEMENTS( level ) )
cbar_arr[0,*] = level
cbar_arr[1,*] = level
CONTOUR, cbar_arr, [1,2], level, $
         NLEVELS = nl, C_COLORS = INDGEN( nl ), XSTYLE = 5, YSTYLE = 5, $
         /ZSTYLE, /FILL, BACKGROUND = nl - 1, COLOR = 0, $ 
         POSITION = [ 0.82 + 0.03, 0.03, $ 
                     0.82 + 0.06, 0.945],$
         LEVELS = level, /NOERASE, XTICKS = 1, XTICKNAME = [' ', ' ', ' ']
AXIS, /YAXIS, charsize=1., charthick=2, /DATA, COLOR = 0, $
      /ZSTYLE, /YSTYLE
device,/close
set_plot,'x'


lon2 = fltarr(nlon+1) & lon2(0:71) = lon & lon2(72) = 180.
height_sh_t = fltarr(nlon+1,nlat,ntimes) & height_sh_t(0:71,*,*)= height_sh & height_sh_t(72,*,*)=height_sh(0,*,*)
set_plot,'ps'
device,filename='polar_sh.eps',/color
nl =256
LOADCT, 33,ncolors = nl
DEVICE, DECOMPOSED = 0
map_set,90,-100,0,/orthographic,/isotropic,/grid,/continents,/horizon,e_horizon={fill:1},limit=[50.,-180,90.,180]
contour,height_sh_t(*,18:35,0),lon2,lat(18:35),$
        nlevels=nl,C_COLORS = INDGEN(nl),xstyle=1,ystyle=1, $
        /overplot,/fill,/isotropic,POSITION=[ 0.14, 0.08, 0.82, 0.92 ]
;,position=position,zrange=[0.,2.6],zstyle= 1,c_colors = c_colors
map_grid,/label,latlab = 175,color=130
map_continents,color=130
;color bar
mxx = MAX(height_sh, MIN = mnn ,/nan) * 0.98
delta_lev = ( mxx - mnn ) / FLOAT( nl )
level = FINDGEN( nl ) * delta_lev + mnn
;colorbar array consists of levels previously determined
cbar_arr = FLTARR( 2, N_ELEMENTS( level ) )
cbar_arr[0,*] = level
cbar_arr[1,*] = level
CONTOUR, cbar_arr, [1,2], level, $
         NLEVELS = nl, C_COLORS = INDGEN( nl ), XSTYLE = 5, YSTYLE = 5, $
         /ZSTYLE, /FILL, BACKGROUND = nl - 1, COLOR = 0, $ 
         POSITION = [ 0.82 + 0.03, 0.03, $ 
                     0.82 + 0.06, 0.945],$
         LEVELS = level, /NOERASE, XTICKS = 1, XTICKNAME = [' ', ' ', ' ']
AXIS, /YAXIS, charsize=1., charthick=2, /DATA, COLOR = 0, $
      /ZSTYLE, /YSTYLE
device,/close
set_plot,'x'

lon2 = fltarr(nlon+1) & lon2(0:71) = lon & lon2(72) = 180.
tec_t = fltarr(nlon+1,nlat,ntimes) & tec_t(0:71,*,*)= mtec & tec_t(72,*,*)=mtec(0,*,*)
tec_t2 = fltarr(nlon+1,nlat,ntimes2) & tec_t2(0:71,*,*)= mtec2 & tec_t2(72,*,*)=mtec2(0,*,*)
dtec_t = tec_t-tec_t2
set_plot,'ps'
device,filename='polar_tec.eps',/color
nl =256
LOADCT, 33,ncolors = nl
DEVICE, DECOMPOSED = 0
map_set,90,-100,0,/orthographic,/isotropic,/grid,/continents,/horizon,e_horizon={fill:1},limit=[50.,-180,90.,180]
contour,dtec_t(*,*,150),lon2,lat,$
        nlevels=nl,C_COLORS = INDGEN(nl),xstyle=1,ystyle=1, $
        /overplot,/fill,/isotropic,POSITION=[ 0.14, 0.08, 0.82, 0.92 ]
;,position=position,zrange=[0.,2.6],zstyle= 1,c_colors = c_colors
map_grid,/label,latlab = 175,color=130
map_continents,color=130
;color bar
mxx = MAX(dtec_t(*,*,150), MIN = mnn ,/nan) * 0.98
delta_lev = ( mxx - mnn ) / FLOAT( nl )
level = FINDGEN( nl ) * delta_lev + mnn
;colorbar array consists of levels previously determined
cbar_arr = FLTARR( 2, N_ELEMENTS( level ) )
cbar_arr[0,*] = level
cbar_arr[1,*] = level
CONTOUR, cbar_arr, [1,2], level, $
         NLEVELS = nl, C_COLORS = INDGEN( nl ), XSTYLE = 5, YSTYLE = 5, $
         /ZSTYLE, /FILL, BACKGROUND = nl - 1, COLOR = 0, $ 
         POSITION = [ 0.82 + 0.03, 0.03, $ 
                     0.82 + 0.06, 0.945],$
         LEVELS = level, /NOERASE, XTICKS = 1, XTICKNAME = [' ', ' ', ' ']
AXIS, /YAXIS, charsize=1., charthick=2, /DATA, COLOR = 0, $
      /ZSTYLE, /YSTYLE
device,/close
set_plot,'x'


lon2 = fltarr(nlon+1) & lon2(0:71) = lon & lon2(72) = 180.
ped2 = fltarr(nlon+1,nlat) & ped2(0:71,*)= msigma_ped2(*,*,10,134) & ped2(72,*)=msigma_ped2(0,*,10,134)
ped = fltarr(nlon+1,nlat) & ped(0:71,*)= msigma_ped(*,*,10,134) & ped(72,*)=msigma_ped(0,*,10,134)
set_plot,'ps'
device,filename='polar_pedersen.eps',/color
nl =256
LOADCT, 33,ncolors = nl
DEVICE, DECOMPOSED = 0
map_set,90,-100,0,/orthographic,/isotropic,/grid,/continents,/horizon,e_horizon={fill:1},limit=[30.,-180,90.,180]
contour,ped2-ped,lon2,lat,$
        nlevels=nl,C_COLORS = INDGEN(nl),xstyle=1,ystyle=1, $
        /overplot,/fill,/isotropic;,POSITION=[ 0.14, 0.08, 0.82, 0.92 ]
;,position=position,zrange=[0.,2.6],zstyle= 1,c_colors = c_colors
map_grid,/label,latlab = 175,color=130
map_continents,color=130
;color bar
mxx = MAX(ped2-ped, MIN = mnn ,/nan) * 0.98
delta_lev = ( mxx - mnn ) / FLOAT( nl )
level = FINDGEN( nl ) * delta_lev + mnn
;colorbar array consists of levels previously determined
cbar_arr = FLTARR( 2, N_ELEMENTS( level ) )
cbar_arr[0,*] = level
cbar_arr[1,*] = level
CONTOUR, cbar_arr, [1,2], level, $
         NLEVELS = nl, C_COLORS = INDGEN( nl ), XSTYLE = 5, YSTYLE = 5, $
         /ZSTYLE, /FILL, BACKGROUND = nl - 1, COLOR = 0, $ 
         POSITION = [ 0.82 + 0.03, 0.03, $ 
                     0.82 + 0.06, 0.945],$
         LEVELS = level, /NOERASE, XTICKS = 1, XTICKNAME = [' ', ' ', ' ']
AXIS, /YAXIS, charsize=1., charthick=2, /DATA, COLOR = 0, $
      /ZSTYLE, /YSTYLE
device,/close
set_plot,'x'

;openw,33,'gden_400km_1996_M107_3.dat'
;printf,33,'global density at 400km (ng/m3)'
;for n = 0, ntimes-1 do begin
;printf,33,gden(n)
;endfor
;close,33

stop

;stru = {smtime:0.,density:fltarr(72,36,29),tn:fltarr(72,36,29),qjoule:fltarr(72,36,29),$
;        qjoule_integ:fltarr(72,36),sheating:fltarr(72,36,29),zg:fltarr(72,36,29)}
;openw, 111, '2003300_311fism_all.dat'
;printf, 111, 'smtime(day)  DEN  TN  qjoule  qjoule_integ  sheating  zg'
;for i = 0, ntimes-1 do begin
;stru.smtime = smtime(i)
;stru.density = mden(*,*,*,i)
;stru.tn = mtn(*,*,*,i)
;stru.qjoule = mqjoule(*,*,*,i)
;stru.qjoule_integ = mqjoule_integ(*,*,i)
;stru.sheating = msheating(*,*,*,i)
;stru.zg = mzg(*,*,*,i)
;printf, 111, stru
;endfor
;close, 111


    int_jh = fltarr(ntimes) & int_jh(*)=0.
    int_jh2 = fltarr(ntimes) & int_jh2(*)=0.
    height_jh = fltarr(nlon, nlat,ntimes) & height_jh(*,*,*) =0.     ;height integrated joule heating

  for n=0,ntimes-1 do begin
  print,n
     for i = 0, nlon-1 do begin
       for j=0, nlat-1 do begin
         done = 0
         for k = 1, nlev-2 do begin

          height_jh(i,j,n) = height_jh(i,j,n)+mqjoule(i,j,k,n)*mden(i,j,k,n)*(mZ(i,j,k+1,n)-mZ(i,j,k-1,n))/2.

          int_jh(n) = int_jh(n) + $
            mqjoule(i,j,k,n)*mden(i,j,k,n)*(mZ(i,j,k+1,n)-mZ(i,j,k-1,n))/2.*cos(lat(j)*!pi/180.)
         
      endfor ;k:lev
area2 = (Ri+mZ(i,j,nlev-1,n))^2*dlat*dlon
int_jh2(n) = int_jh2(n)+mqjoule_integ(i,j,n)*cos(lat(j)*!pi/180.)*area2
  endfor ;j:lat

endfor ;i:lon
endfor ;n:ntimes

int_jh2 = int_jh2*1.e-16   ; GW

; Z was in km, and convert from ergs/s*1e5 to GW by a factor of 1.e-11
     area = (Ri+mZ(0,0,nlev-1,ntimes-1))^2*dlat*dlon
;YH 06/13/11
     int_jh = int_jh*area*1.e-16  ;GW


openw, 11, 'hourly_global_JH_1996_M107_300_366.dat'
printf, 11, 'smtime(day)  int_jh2(GW)  int_jh(GW)'
for i = 0, ntimes-1 do begin
printf, 11, smtime(i), int_jh2(i), int_jh(i)
endfor
close, 11

;set_plot,'ps'
;device,filename = 'ijh1996_new.eps'
;plot, smtime_daily,int_jh_daily,ytitle = 'Joule heating(GW)',xtitle='Days',xr=[0,365],xstyle =1
;device,/close
;set_plot,'X'


set_plot,'ps'
device, filename ='height_jh_smax.ps'
;device,filename='../tiegcm_plots/JH_smin_2008_F107_70_hp_20_cpcp_50_F_cut_2.ps'
device,decomposed=0,/color
loadct, 39,ncolors=100,bottom=1
c_colors = indgen(100)+1
position = [0.1,0.1,0.9,0.75]

imax = max(height_jh(*,18:35,ntimes-1))
imin = min(height_jh(*,18:35,ntimes-1))
lon2 = fltarr(nlon+1) & lon2(0:71) = lon & lon2(72) = 180.
height_jh2 = fltarr(nlon+1,nlat,ntimes) & height_jh2(0:71,*,*)= height_jh & height_jh2(72,*,*)=height_jh(0,*,*)

map_set,90,-100,0,/orthographic,/isotropic,/grid,/continents,/horizon,e_horizon={fill:1},limit=[50.,-180,90.,180]
contour, height_jh2(*,18:35,ntimes-1),lon2,lat(18:35),/overplot,/fill,/isotropic,nlevels=100,position=position,zrange=[0.,2.6],zstyle = 1,c_colors = c_colors
map_grid,/label,latlab = 175,color=130
map_continents,color=130

;colorbar,bottom=1,ncolors=50,position=[.2,0.85,0.75,0.9]
;colorbar,bottom=1,divisions=20,ncolors=100,position=[0.2,0.85,0.75,0.9],pscolor=pscolor,title='Joule heating(W/m2)'
;contour, height_jh(*,18:35,ntimes-1),lon,lat(18:35),/overplot,/noerase,nlevels=12

device,/close
set_plot,'X'

stop
;close,/all
end