@ $TGCMROOT/readhist/readfile.pro
;-----------------------------------------------------------------------
function mklevels, f,interval ; interval is output arg
;
; Return a reasonable contour levels array, given field f. 
; Also return interval as output arg.
;
  nlev = 12 ; default to start
  fmin = min(f) & fmax = max(f)
  ifmin = round(fmin-0.5) & ifmax = round(fmax+0.5)
  clevel0 = float(ifmin)
  interval = round((fmax-fmin)/float(nlev))
  nlev = 0
  while clevel0+float(nlev*interval) le fmax do begin
    nlev = nlev+1
  endwhile
  nlev = nlev+1
  clevels = fltarr(nlev)
  for i=0,nlev-1 do begin
    clevels[i] = clevel0+float(i*interval)
  endfor
  return,clevels
end
;-----------------------------------------------------------------------
function slt2glon, slt,ut,mag=mag
;
; Return geo longitude, given slt and ut.
;
  glon = (slt-ut) * 15.
  if keyword_set(mag) then glon = glon+71.
  if (glon ge 180.) then glon = glon-360.
  if (glon lt -180.) then glon = glon+360.
  return,glon
end
;-----------------------------------------------------------------------
function ixfind, z,val,del
;
; Find nearest index to val in z, with increment del
; If not found, return -1
; If del==0, look for exact val
;
nz = n_elements(z)
if del gt 0 then begin
  for i=0,nz-1 do begin
    if (val ge z(i)-0.5*del and val le z(i)+0.5*del) then return,i
  endfor
endif else begin
  for i=0,nz-1 do begin
    if (val eq z(i)) then return,i
  endfor
endelse
return,-1
end
;-----------------------------------------------------------------------
function ifind_mtime, mtime,ntime,mtime_find
;
; itime_beg = ifind_mtime(file.mtime,file.ntime,mtime_beg)
;
for i=0,ntime-1 do begin
  if (mtime[0,i] eq mtime_find[0] and $
      mtime[1,i] eq mtime_find[1] and $
      mtime[1,2] eq mtime_find[2]) then return,i
endfor
print,format="('>>> Could not find mtime ',3i4')",mtime_find
print,format="('ntime=',i4,' mtime=',/,(3i4)",ntime,mtime
return,-1
end
;-----------------------------------------------------------------------
pro meso,ps=ps
;
; User options:
;   path                 ; History file (or directory) to be read and processed
;   mtime_beg,mtime_end  ; Beginning and ending nth histories to be proc'd
;   nthmean              ; history-interval in which to average zm tn,ht
;
; Days of month, for convenience:
 months = ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec']
;                    J  F  M   A   M   J   J   A   S   O   N   D
 days_per_month  = [31,28,31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
 first_day_month = [ 1,32,60, 91,121,152,182,213,244,274,305,335]
 last_day_month  = [31,59,90,120,151,181,212,243,273,304,334,365]
;
; User set: path to file or directory (if dir, pickfile dialogue will pop-up)
;
 path='/hao/tgcm/timegcm1.2_yr04/py4n/TNZ_yr04.nc'                 ; full 2004 yr (daily)
;path='/hao/tgcm/timegcm1.2_yr04/py4n/ROBLE.timegcm1.2.py4n001.nc'
;path='/hao/tgcm/timegcm1.2_yr04/sy4n/ROBLE.timegcm1.2.sy4n001.nc'
;path='/e/foster/tiegcm1.8/tiegcm1.8-linux/FOSTER.tiegcm1.8.stest001.nc'
;
; User set: beginning and ending histories to be proc'd
;
;mtime_begin = intarr(3) & mtime_begin = [60,0,0] ; Mar
;mtime_endin = intarr(3) & mtime_endin = [90,0,0] ; Mar
;mtime_begin = intarr(3) & mtime_begin = [152,0,0] ; Jun
;mtime_endin = intarr(3) & mtime_endin = [181,0,0] ; Jun
;mtime_begin = intarr(3) & mtime_begin = [244,0,0] ; Sep
;mtime_endin = intarr(3) & mtime_endin = [273,0,0] ; Sep
 mtime_begin = intarr(3) & mtime_begin = [335,0,0] ; Dec
 mtime_endin = intarr(3) & mtime_endin = [365,0,0] ; Dec
;
; Full year:
;mtime_begin = intarr(3) & mtime_begin = [1,0,0]
;mtime_endin = intarr(3) & mtime_endin = [365,0,0]
;
; User set: nthmean = history interval in which to take means of zm tn,ht
; If calendar monthly means are desired, init nthmean to days_per_month[0],
;   and set nthmean_flag = 'monthly' (otherwise nthmean_flag is ignored)
; Set nthmean=0 to turn this off.
;
 nthmean = 0  ; average zmthmin and zmht_tnmin every nthmean histories 
;nthmean = days_per_month[0]
;nthmean_flag = 'monthly'
;
; User set: zp0,zp1 = bottom,top levels of column to search for min tn.
; Mesopause (min tn) zp is typically -7 to -8. (~80-95 km)
; (can set -999.,999. to get full column)
;
 zp0 = -12.
 zp1 = -5.
;zp0 = -999.
;zp1 =  999.
;
;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
;
; Read history file:
;
file = readfile(path=path)
;
;ncdf_attget,file.ncid,/global,'missing_value',missing_value
;print,'missing_value=',missing_value
;
; Read TN and Z from file:
;
print,'Reading TN...'
showvar,file.ncid,'TN'
TN = getvar(file.ncid,'TN',tn_info) ; (lon,lat,lev,time)
print,'Reading Z...'
showvar,file.ncid,'Z'
Z  = getvar(file.ncid,'Z' ,z_info)  ; (lon,lat,lev,time)
;
ut = getvar(file.ncid,'ut',ut_info)
time_mins = getvar(file.ncid,'time',time_info)
;
nlev = n_elements(file.lev)
nlat = n_elements(file.lat)
TN[*,*,nlev-1,*] = TN[*,*,nlev-2,*] ; set top boundary tn to nlev-2
dz = file.lev[1]-file.lev[0]
dlon = file.lon[1]-file.lon[0]
;
if (keyword_set(ps)) then begin
  set_plot,'PS'
  device,/color
endif else begin
  set_plot,'X'
endelse
;
tnlt   = fltarr(nlat,nlev)    ; tn slice at current local time
htlt   = fltarr(nlat,nlev)    ; height slice at current local time
tnmin       = fltarr(nlat,24) ; min of tn at hourly local times
ht_tnmin    = fltarr(nlat,24) ; ht at tnmin at hourly local times
av_tnmin    = fltarr(nlat,24) ; history-averaged tnmin (init to zero)
av_ht_tnmin = fltarr(nlat,24) ; history-averaged ht_tnmin (init to zero)
;
tnzm    = fltarr(nlat,nlev)    ; zonal mean tn slice
htzm    = fltarr(nlat,nlev)    ; zonal mean height slice
av_tnzm = fltarr(nlat,nlev)    ; time-averaged zm tn slice
av_htzm = fltarr(nlat,nlev)    ; time-averaged zm ht slice
;
; Beginning and ending histories to process:
;
mtime_beg = long(mtime_begin) & mtime_end = long(mtime_endin)
itime_beg = ifind_mtime(file.mtime,file.ntime,mtime_beg)
if (itime_beg lt 0) then return
itime_end = ifind_mtime(file.mtime,file.ntime,mtime_end)
if (itime_end lt 0) then return
;
if (itime_beg gt file.ntime-1 or itime_end gt file.ntime-1 or $
    itime_beg gt itime_end) then begin
  print,format="(/,'>>> meso: bad itime_beg=',i4,' or bad itime_end=',i4)",$
    itime_beg,itime_end
  print,format="('There are ',i4,' histories on this file.',/)",file.ntime
  return
endif
ntime = itime_end-itime_beg+1           ; number of histories to process
if (ntime lt 2) then begin
  print,format="(/,'>>> meso: must process at least 2 histories: ntime=',i4,/)",$
    ntime
  return
endif
;
; nthmean = history interval in which to take means
;
if (nthmean gt 0) then begin
  if (ntime mod nthmean ne 0 and nthmean_flag ne 'monthly') then begin
    print,format=$
      "(/,'>>> meso: ntime mod nthmean /= 0: ntime=',i4,' nthmean=',i4)",$
      ntime,nthmean
    print,format="('itime_beg=',i4,' itime_end=',i4)",itime_beg,itime_end
    return
  endif
  if (nthmean_flag ne 'monthly') then $
    nmean = ntime/nthmean $ ; number of averaged time (hist) points (y-axis)
  else nmean = 12           ; average monthly
endif
;
; Validate search column zp0,zp1 (user set above):
;
if zp0 ge zp1 then begin
  print,format="(/,'>>> zp0 must be < zp1: zp0=',f8.2,' zp1=',f8.2)",zp0,zp1
  return
endif
if zp0 lt file.lev[0] then zp0 = file.lev[0]
if zp1 gt file.lev[nlev-1] then zp1 = file.lev[nlev-1]
izp0 = ixfind(file.lev,zp0,dz)
izp1 = ixfind(file.lev,zp1,dz)
if (izp0 lt 0) then begin
  print,format="(/,'>>> bad izp0=',i3,' zp0=',f8.2)",izp0,zp0
  return
endif
if (izp1 lt 0) then begin
  print,format="(/,'>>> bad izp1=',i3,' zp1=',f8.2)",izp1,zp1
  return
endif
print,format="('Search column: zp0,zp1=',f8.2,',',f8.2,' izp0,izp1=',2i4)",$
  zp0,zp1,izp0,izp1
;
zmtnmin       = fltarr(nlat)       ; min of zonal mean tn
zmht_tnmin    = fltarr(nlat)       ; ht at min of zonal mean tn
if (nthmean gt 0) then begin
  av_zmtnmin    = fltarr(nlat,nmean) ; min of zonal mean tn (ave every nthmean) 
  av_zmht_tnmin = fltarr(nlat,nmean) ; ht at min of zm tn (ave every nthmean)
  dday          = fltarr(nmean)      ; decimal days for y-axis of av zm tn,ht
endif
;
nhist = 0                          ; init nth history
imean = 0                          ; init index to av_zmtnmin, av_zmht_tnmin
imonth = 0
;
; History loop:
;
for itime=itime_beg,itime_end do begin  ; history loop
  nhist = nhist+1                       ; nth history being processed (1-based)
  decday = float(file.mtime[0,itime])+float(file.mtime[1,itime])+$
           float(file.mtime[2,itime])/60.
  print,format="('History ',i4,' (',i4,' of ',i4,'): mtime=',3i4,' decday=',f8.3)",$
    itime,nhist,ntime,file.mtime[*,itime],decday
;
; Local time loop:
;
  for ltime=0,23 do begin
    glon = slt2glon(ltime,ut[itime])    ; longitude at current local time
    ixlon = ixfind(file.lon,glon,dlon)  ; index to glon
    for k=0,nlev-1 do begin
      for j=0,nlat-1 do begin
        tnlt[j,k] = TN[ixlon,j,k,itime]       ; tn at local time ltime
        htlt[j,k] = Z [ixlon,j,k,itime]*1.e-5 ; ht at local time ltime (km)
      endfor
    endfor ; k
;
; Find height at min tn:
;
    avzp = 0.
    for j=0,nlat-1 do begin
      tnmin[j,ltime] = min(tnlt[j,izp0:izp1])            ; tn minimum of column 
      kmin = ixfind(tnlt[j,izp0:izp1],tnmin[j,ltime],0.) ; k index at tnmin
      kmin = kmin+izp0
      ht_tnmin[j,ltime] = htlt[j,kmin] ; ht at tn min
      avzp = avzp+file.lev[kmin]/float(nlat)

; Plot profile at current lat and lt:
;     title=string(format=$
;       "('ltime=',i2,' j=',i2,' zp0,1=',f6.2,',',f6.2,' zpmin=',f6.2,' tnmin=',f8.2,' htmin=',f6.2)",$
;       ltime,j,zp0,zp1,file.lev[kmin],tnmin[j,ltime],ht_tnmin[j,ltime])
;     plot,tnlt[j,*],file.lev,title=title,xtitle='TN',ytitle='Zp',xstyle=1,ystyle=1
;     if (not keyword_set(ps)) then cursor,x,y,/up

    endfor ; j
;   print,format="('ltime=',i3,' ht_tnmin (lat ave)=',f8.2,' avzp=',f8.2)",$
;     ltime,mean(ht_tnmin[*,ltime]),avzp
  endfor ; ltime
;
; Update time-average over histories:
;
  av_tnmin    = av_tnmin    + tnmin/float(ntime)
  av_ht_tnmin = av_ht_tnmin + ht_tnmin/float(ntime)
;
; Zonal mean tn and ht:
;
  if (nthmean eq 0) then goto,endloop
  for k=0,nlev-1 do begin
    for j=0,nlat-1 do begin
      tnzm[j,k] = mean(TN[*,j,k,itime])        ; zonal mean tn 
      htzm[j,k] = mean(Z [*,j,k,itime])*1.e-5  ; zonal mean ht (km)
;
; Update zm hist means:
;
      av_tnzm[j,k] = av_tnzm[j,k] + tnzm[j,k]/float(nthmean)
      av_htzm[j,k] = av_htzm[j,k] + htzm[j,k]/float(nthmean)
    endfor ; j
  endfor ; k
;
; Re-init ave zm every nthmean:
;
  if (nthmean_flag ne 'monthly' and nhist mod nthmean eq 0) then begin
    dday[imean] = decday
    if (imean lt nmean) then begin
      print,'Re-init means: nhist=',nhist,' nthmean=',nthmean,' imean=',imean,$
        ' dday=',dday[imean]
;
; Find zm height at min zm tn:
;
      for j=0,nlat-1 do begin
        av_zmtnmin[j,imean] = min(av_tnzm[j,izp0:izp1])
        kmin = ixfind(av_tnzm[j,izp0:izp1],av_zmtnmin[j,imean],0.)
        kmin = kmin+izp0
        av_zmht_tnmin[j,imean] = av_htzm[j,kmin]
      endfor
    endif
    av_tnzm[*,*] = 0.
    av_htzm[*,*] = 0.
    imean = imean+1                          ; index for next set
;
; Doing monthly means:
;
  endif else if (nthmean_flag eq 'monthly') then begin
;
; Re-init at end of each month:
;
    if (imean+1 eq days_per_month[imonth]) then begin
      dday[imonth] = float(imonth+1) ; dday is 1=based, imonth is 0-based
      format="('Re-init means: nhist=',i4,' nthmean=',i4,' imonth=',i2,' (',a,')')"
      print,format=format,nhist,nthmean,imonth,months[imonth]

;
; Find zm height at min zm tn:
;
      for j=0,nlat-1 do begin
        av_zmtnmin[j,imonth] = min(av_tnzm[j,izp0:izp1])
        kmin = ixfind(av_tnzm[j,izp0:izp1],av_zmtnmin[j,imonth],0.)
        kmin = kmin+izp0
        av_zmht_tnmin[j,imonth] = av_htzm[j,kmin]
      endfor
      av_tnzm[*,*] = 0.
      av_htzm[*,*] = 0.
      imean = 0
      imonth = imonth+1
      if (imonth lt 12) then nthmean = days_per_month[imonth]
    endif else begin    ; end of month
      imean = imean+1
    endelse
  endif ; doing monthly means
  endloop:
endfor ; itime=itime_beg,itime_end
;
; Load color table:
;
loadct,get_names=ctnames
for i=0,n_elements(ctnames)-1 do begin
  if ctnames[i] eq 'Rainbow + white' then begin
    print,'Loading color table ',ctnames[i]
    loadct,i
  endif 
endfor
;
; Plot position:
;
xur = .90 & yur = .93
xll = .15 & yll = .20
position  = [xll,yll,xur,yur]
xmid = .5*(xur+xll)
chsize = 1.2
c_chsize = 1.2
;
; mtime title:
;
dtime = 0. ; delta time (mins)
if (ntime gt 1) then $
  dtime = time_mins[itime_beg+1]-time_mins[itime_beg] ; mins (assume uniform dtime)
mtime_title = 'Averaged over '+string(format="(i4)",ntime)+' model times: '+$
  string(format="(i4,2i3)",file.mtime[*,itime_beg])+' to '+$
  string(format="(i4,2i3)",file.mtime[*,itime_end])+' by '+$
  string(format="(i5)" ,dtime)+' mins'
mtime_title_pos = yll-.11  & if (keyword_set(ps)) then mtime_title_pos  = yll-.14
minmax_title_pos = yll-.16 & if (keyword_set(ps)) then minmax_title_pos = yll-.19
;
; Contour history-averaged tnmin at lt:
;
title = 'MESOPAUSE TEMPERATURE (DEG K)'
levels_tn = mklevels(av_tnmin,interval_tn)
;print,'levels for av_tnmin:' & print,levels_tn
minmax_title = 'Min,max = '+string(format="(f8.2)",min(av_tnmin))+', '+$
                            string(format="(f8.2)",max(av_tnmin))+' cint = '+$
                            string(format="(i3)",interval_tn)
localtime = findgen(24)+1.
; Monochrome:
;contour,av_tnmin,file.lat,localtime,/follow,xstyle=1,ystyle=1,$
;  xtitle='LATITUDE (DEG NORTH)',ytitle='LOCAL TIME (HRS)',$
;  nlevels=nclev,position=position,/norm,title=title,charsize=chsize
;
; Color fill:
 contour,av_tnmin,file.lat,localtime,xstyle=5,ystyle=5,/follow,/fill,$
   position=position,/norm,levels=levels_tn
 contour,av_tnmin,file.lat,localtime,/noerase,/follow,xstyle=1,ystyle=1,$
   xtitle='LATITUDE (DEG NORTH)',ytitle='LOCAL TIME (HRS)',c_charsize=c_chsize,$
   position=position,/norm,title=title,charsize=chsize,levels=levels_tn
;
; Add bottom labels:
xyouts,xmid,mtime_title_pos,mtime_title,align=0.5,/norm,charsize=chsize
xyouts,xmid,minmax_title_pos,minmax_title,align=0.5,/norm,charsize=chsize
;
; Contour history-averaged ht_tnmin at lt:
;
if (not keyword_set(ps)) then cursor,x,y,/up
title = 'MESOPAUSE HEIGHT (KM)'
levels_ht = mklevels(av_ht_tnmin,interval_ht)
;print,'levels for av_ht_tnmin:' & print,levels_ht
minmax_title = 'Min,max = '+string(format="(f8.2)",min(av_ht_tnmin))+', '+$
                            string(format="(f8.2)",max(av_ht_tnmin))+' cint = '+$
                            string(format="(i3)",interval_ht)
;
; Monochrome:
;contour,av_ht_tnmin,file.lat,localtime,/follow,xstyle=1,ystyle=1,$
;  xtitle='LATITUDE (DEG NORTH)',ytitle='LOCAL TIME (HRS)',$
;  nlevels=nclev,position=position,/norm,title=title
;
; Color fill:
 contour,av_ht_tnmin,file.lat,localtime,xstyle=5,ystyle=5,/follow,/fill,$
   levels=levels_ht,position=position,/norm
 contour,av_ht_tnmin,file.lat,localtime,/noerase,/follow,xstyle=1,ystyle=1,$
   xtitle='LATITUDE (DEG NORTH)',ytitle='LOCAL TIME (HRS)',c_charsize=c_chsize,$
   levels=levels_ht,position=position,/norm,title=title,charsize=chsize
;
xyouts,xmid,mtime_title_pos,mtime_title,align=0.5,/norm,charsize=chsize
xyouts,xmid,minmax_title_pos,minmax_title,align=0.5,/norm,charsize=chsize
;
; Contour history-averaged zm tnmin:
;
if (nthmean eq 0) then goto,done
if (not keyword_set(ps)) then cursor,x,y,/up
title = 'MESOPAUSE TEMPERATURE (DEG K)'
levels_tn = mklevels(av_zmtnmin,interval_tn)
;print,'levels for av_zmtnmin:' & print,levels_tn
minmax_title = 'Min,max = '+string(format="(f8.2)",min(av_zmtnmin))+', '+$
                            string(format="(f8.2)",max(av_zmtnmin))+' cint = '+$
                            string(format="(i3)",interval_tn)
if (nthmean_flag eq 'monthly') then begin
  ytitle = 'DAY OF YEAR'
  ytitle = 'MONTH'
  mtime_title = 'Monthly average'
endif else begin
  mtime_title = 'Averaged every '+string(format="(i4)",nthmean)+' model times'
  mtime_title = mtime_title+' ('+string(format="(i3)",nmean)+' points on y-axis)'
endelse
;
; Monochrome:
;contour,av_zmtnmin,file.lat,dday,/follow,xstyle=1,ystyle=1,$
;  xtitle='LATITUDE (DEG NORTH)',ytitle=ytitle,$
;  nlevels=nclev,position=position,/norm,title=title,charsize=chsize
;
; Color fill:
 contour,av_zmtnmin,file.lat,dday,/follow,/fill,xstyle=5,ystyle=5,$
   levels=levels_tn,position=position,/norm
 contour,av_zmtnmin,file.lat,dday,/follow,/noerase,xstyle=1,ystyle=1,$
    xtitle='LATITUDE (DEG NORTH)',ytitle=ytitle,c_charsize=c_chsize,$
    levels=levels_tn,position=position,/norm,title=title,charsize=chsize
;
xyouts,xmid,mtime_title_pos,mtime_title,align=0.5,/norm,charsize=chsize
xyouts,xmid,minmax_title_pos,minmax_title,align=0.5,/norm,charsize=chsize
;
; Contour history-averaged zm ht at zm tnmin:
;
if (not keyword_set(ps)) then cursor,x,y,/up
title = 'MESOPAUSE HEIGHT (KM)'
minmax_title = 'Min,max = '+string(format="(f8.2)",min(av_zmht_tnmin))+', '+$
                            string(format="(f8.2)",max(av_zmht_tnmin))+' cint = '+$
                            string(format="(i3)",interval_ht)
levels_ht = mklevels(av_zmht_tnmin,interval_ht)
;print,'levels for av_zmht_tnmin:' & print,levels_ht
;
; Monochrome:
;contour,av_zmht_tnmin,file.lat,dday,/follow,xstyle=1,ystyle=1,$
;  xtitle='LATITUDE (DEG NORTH)',ytitle=ytitle,$
;  nlevels=nclev,position=position,/norm,title=title,charsize=chsize
;
; Color:
 contour,av_zmht_tnmin,file.lat,dday,/follow,/fill,xstyle=5,ystyle=5,$
   levels=levels_ht,position=position,/norm
 contour,av_zmht_tnmin,file.lat,dday,/follow,/noerase,xstyle=1,ystyle=1,$
   xtitle='LATITUDE (DEG NORTH)',ytitle=ytitle,c_charsize=c_chsize,$
   levels=levels_ht,position=position,/norm,title=title,charsize=chsize
;
xyouts,xmid,mtime_title_pos,mtime_title,align=0.5,/norm,charsize=chsize
xyouts,xmid,minmax_title_pos,minmax_title,align=0.5,/norm,charsize=chsize
done:
;
; Clean up:
;
free_lun,file.ncid
if (keyword_set(ps)) then device,/close
end ; pro meso