pro msis00, date, sec, glat, glon, alt, t_alt, t_exo, nd, md, $
              ft7p=ft7p_usr, ft7a=ft7a_usr, apd=ap_usr, $
              SI=si_yesno, Ap_3hr=ap3_yesno, HELP=h_flag
;
; IDL interface to Al Hedin's MSIS-90 model atmosphere.
;
; B. Knapp, 1993-04-15
; 
; This IDL procedure provides a simplified interface to Al Hedin's
; MSIS-90 2-dimensional model atmosphere.  MSIS-90 requires the F10.7,
; F10.7 81-day smoothed, and Ap solar activity indices as input, and
; this procedure retrieves the appropriate ones automatically before
; calling MSIS itself.  There are several MSIS-90 entry points, but
; only one of them, GTD7, is supported by this procedure.  Also, support
; for varying the MSIS-90 model calculation parameters (via its common
; block SW variables) is not provided.  For full details on the MSIS-90
; model see the Fortran source code MSIS00.F, and the two papers by
; Al Hedin, "MSIS-86 Thermosphere Model" JGR 92, pp 4649-4662 (1987),
; and "Extension of the MSIS Thermosphere Model into the Middle and
; Lower Atmosphere" JGR 96, pp. 1159-1172 (1991).
;
; Required input (positional) arguments:
;
;    date    -- scalar date of the form yyyyddd or yyddd
;
;    sec     -- scalar time of day (UT) in seconds
;
;    glat    -- scalar geographic latitude (deg)
;
;    glon    -- scalar geographic longitude (deg)
;
;    alt     -- scalar or array of requested altitude(s) (km)
;
; Optional input (keyword) arguments:
;
;    ft7p    -- User-supplied value for F10.7 flux value for previous day
;               (date-1)
;
;    ft7a    -- User-supplied value for 81-day mean of F10.7 flux centered
;               on this day (date)
;
;    apd     -- User-supplied value (or array) for Ap index (or 3-hourly
;               Ap indices; see msis00.f source code for format) for this
;               day (date)
;
;    Note: For any of the above indices not supplied by the user, the
;       correct values for the given date are automatically supplied.
;
;    /SI     -- to request results in SI units (m^-3); if not supplied,
;               results are in cgs units (cm^-3)
;
;    /Ap_3hr -- to request 3-hourly, rather than daily, Ap indices
;               be used; if not supplied, daily Ap value is used.
;
;    /Help   -- to request full documentation to be displayed; if this
;               keyword is set, only the documentation is displayed,
;               and no MSIS data are returned.
;
; Output arguments:
;
;    t_alt   -- scalar or array of same type as alt, with the temperature
;               (K) at each of the input altitudes
;
;    t_exo   -- scalar, exospheric temperature
;
;               Note: if only temperatures are desired, omit the following
;                     two output parameters, so that time is not wasted
;                     computing the number densities.
;
;    nd      -- number density, an array of the same length as the input
;               argument alt, each element of which is a structure with
;               seven scalar  fields, one for each of eight atmospheric
;               species:
;                  nd.h  -- atomic hydrogen
;                  nd.he -- helium
;                  nd.n  -- atomic nitrogen
;                  nd.o  -- atomic oxygen
;                  nd.n2 -- molecular nitrogen
;                  nd.o2 -- molecular oxygen
;                  nd.ar -- argon
;                  nd.ao -- anomalous oxygen
;
;    md      -- scalar or array of same type as alt, with the total mass
;               density of all species at each of the requested altitudes.
;
; Change Log:
;
;    1994-02-03, BGK, to make the ap input parameter to the
;       msis00 fortran code a vector of length 7 (even though
;       only the first element is used), so that when gtd7 calls
;       vtst to see if any inputs have changed--and it looks at
;       all 7 elements of the ap array--this array is fully defined.
;
;    1995-08-17, B. Knapp--improved on-line help and added
;       /HELP keyword
;
;    1995-10-04, B. Knapp--add /Ap_3hr keyword to support
;       3-hourly Ap inputs to GTD7
;
;    1998-06-10, B. Knapp, IDL v. 5 compliance
;
;    2000-01-26, B. Knapp, Y2K compliance
;
;    2000-02-24, B. Knapp, Bug fix: YYDDD wrong type--made Long
;
;    2000-02-25, B. Knapp, Change 91-day smooth to 81-day, and
;       add ft7p, ft7a, apd keyword parameters.
;
;    2001-04-30, B. Knapp, to choose 32-bit or 64-bit dll
;       automatically (IDL 5.4).
;
;    2007-05-02, B. Knapp, Linux port (use getenv(), path_sep())
;
  display = 'more '
  sep = path_sep()
  root = getenv('bgkroot')
  src_path = root+sep+'idl'+sep+'msis'+sep
;
; Identify shareable object library
; dll = shareable_object('msis00_')
;
; 11/17/08 btf:
; Sys var !msis_so is set by startup.pro, using env var UTPROC_MSO,
; which is set by the utproc exec script.
;
  dll = !msis_so
  ext = '_pass_'
  
; Print usage?
  if n_params() eq 0 and (not keyword_set( h_flag )) then begin
     print," "
     print," MSIS00 is a simplified IDL interface to the Fortran code of"
     print," Al Hedin's two-dimensional thermospheric model MSIS-90, which"
     print," provides estimates of temperature and of number density of seven"
     print," prevalent species (H, He, N, O, N2, O2, and Ar).  It uses the"
     print," 10.7 cm radio flux (for the previous day, and a 81-day smooth"
     print," centered on the requested day) and the Ap geomagnetic indices"
     print," as input, and these data are obtained automatically.  Complete"
     print," documentation is avaialable with the /HELP keyword."
     print," "
     print," Usage:"
     print," "
     print,"     msis00, date, sec, glat, glon, alt, t_alt, t_exo, nd, md, $"
     print,"             ft7p=<ft7p_value>, ft7a=<ft7a_value>, "+$
                                                       "apd=<ap_value[s]>, $"
     print,"             /SI, /Ap_3hr, /HELP"
     print," "
     return
  endif
;
; Print full documentation?
  if keyword_set( h_flag ) then begin
      spawn, display+src_path+'msis00.readme'
      return
  endif
  
; Check inputs
  sz = size( date )
  nz = n_elements(sz)
  if ( sz[0] ne 0 ) or ( sz[nz-2] eq 0 ) then begin
     message,'argument  date  must be a scalar'
     return
  endif
;
  sz = size( sec )
  nz = n_elements(sz)
  if ( sz[0] ne 0 ) or ( sz[nz-2] eq 0 ) then begin
     message,'argument  sec  must be a scalar'
     return
  endif
;
  sz = size( glat )
  nz = n_elements(sz)
  if ( sz[0] ne 0 ) or ( sz[nz-2] eq 0 ) then begin
     message,'argument  glat  must be a scalar'
     return
  endif
;
  sz = size( glon )
  nz = n_elements(sz)
  if ( sz[0] ne 0 ) or ( sz[nz-2] eq 0 ) then begin
     message,'argument  glon  must be a scalar'
     return
  endif
;
  szalt = size( alt )
  nz = n_elements(szalt)
  if ( szalt[0] gt 1 ) or ( szalt[nz-2] eq 0 ) then begin
     message,'argument  alt  must be a scalar or vector'
     return
  endif
;
; Is the number density structure type defined?
  common msis00_sav, nd_type
  if n_elements( nd_type ) eq 0 then nd_type = $
     { _msis00_, h:0., he:0., n:0., o:0., n2:0., o2:0., ar:0., ao:0. }
;
; Are we computing only temperatures, or also number densities?
  req_code = default_int((n_params() lt 8) ? 0 : 48)
;
; Convert inputs to correct types and ranges
  ldate = y2toy4( date )
  yyddd = default_int( ldate mod 100000L )

; print,'msis00: yyddd=',yyddd

  fsec = default_real( sec )
  flat = default_real( glat )
  flon = ( (default_real( glon )+180.) mod 360. ) - 180.
  if flon lt -180. then $
     flon = flon + 360. $
  else if flon ge  180. then $
     flon = flon - 360.
;      
  if szalt[0] eq 0 then $
     falt = [ default_real( alt ) ] $
  else $
     falt = default_real( alt )
;
; Get the current GTD7 switch values
  sw = default_real_array(25)
  entry = 'tretrv'+ext

; print,'msis00: dll=',dll,' entry=',entry
; print,'msis00: sw=',sw

  dummy = call_external( dll, entry, sw )

; Get the solar indices for the requested date
  if n_elements( ft7p_usr ) eq 0 or n_elements( ft7a_usr ) eq 0 or $
     n_elements( ap_usr ) eq 0 then $
     get_t7_ap, ldate, ft7p_auto, ft7a_auto, ap_auto
;
; Use user-supplied indices?
  if n_elements( ft7p_usr ) eq 0 then $
     ft7 = default_real( ft7p_auto ) $
  else $
     ft7 = default_real( ft7p_usr )
;
  if n_elements( ft7a_usr ) eq 0 then $
     ft7a = default_real( ft7a_auto ) $
  else $
     ft7a = default_real( ft7a_usr )

; print,'msis00: date=',date,' ft7=',ft7,' ft7a=',ft7a
; print,format="('msis00: date=',i8,' ft7=',f10.3,' ft7a=',f10.3)",$
;   date,ft7,ft7a
;
  fap=default_real_array(7)
  if n_elements( ap_usr ) eq 0 then $
     fap[0] = default_real( ap_auto ) $
  else $
     fap[0] = default_real( ap_usr )

;
; If 3-hourly Ap values are requested, compute the remaining entries of fap
  if keyword_set( ap3_yesno ) and n_elements( ap_usr ) eq 0 then begin
;
;    Get today's and 3 previous days' 3-hourly Ap values
     ap32 = default_real_array(32)
     ap3_dat = intarr(8)
     for j=-3,0 do begin
        yd = yd2( ldate, j )
        get_ap3, yd, ap3_dat, ap3_status
        if ap3_status ne 0 then goto,quit
        ap32[(j+3)*8] = ap3_dat
     endfor
     quit:
     if ap3_status ne 0 then begin
        print,' *** Ap 3-hourly data not available; using daily average ***'
        sw[8] = 1.
     endif else begin
;
;       Compute the fap input values for GTD7
        k = 24+floor(fsec/10800.)
        fap[1] = ap32[k]
        fap[2] = ap32[k-1]
        fap[3] = ap32[k-2]
        fap[4] = ap32[k-3]
        fap[5] = total(ap32[k-11:k- 4])/8.0
        fap[6] = total(ap32[k-19:k-12])/8.0
;
;       Signal 3-hourly inputs are supplied
        sw[8] = -1.
     endelse
  endif else begin
;
;    Signal only daily Ap value supplied
     sw[8] = 1.
  endelse
;
; Reset the switches
  entry = 'tselec'+ext
  dummy = call_external( dll, entry, sw )
;
; Prepare the outputs
  if szalt[0] eq 0 then nalt=1 else nalt=szalt[1]
  md = fltarr( nalt )
  t_alt = fltarr( nalt )
  t_exo = 0.
  nd = replicate( nd_type, nalt )
;
; Variables for GTD7
  d = default_real_array( 9 )
  t = default_real_array( 2 )
  stl = fsec/3600. + flon/15.  ;local solar time
;
; Select cgs or SI units
  si = default_int(keyword_set( si_yesno ) ? -1 : 0)
  entry = 'meters'+ext
  dummy = call_external( dll, entry, si )
;
; For each requested altitude, call MSIS model
  entry = 'gtd7'+ext
  for j = 0,nalt-1 do begin
;
;    if (yyddd eq 2001 and j eq 0) then begin
;      print,'msis00 call gtd7: j=',j,' yyddd=',yyddd
;      print,'fap=',fap
;      help,yyddd & help,fsec & help,falt[j]
;      help,flat & help,flon & help,stl & help,ft7a & help,ft7
;      help,fap & help,req_code & help,d & help,t
;    endif

     dummy = call_external( dll, entry, $
        yyddd, fsec, falt[j], flat, flon, stl, ft7a, ft7, fap, req_code, $
        d, t )
;
     t_alt[j] = t[1]
     nd[j].h  = d[6]
     nd[j].he = d[0]
     nd[j].n  = d[7]
     nd[j].o  = d[1]
     nd[j].n2 = d[2]
     nd[j].o2 = d[3]
     nd[j].ar = d[4]
     nd[j].ao = d[8]
     md[j]    = d[5]
  endfor
  t_exo = t[0]
  if szalt[0] eq 0 then begin
     t_alt = t_alt[0]
     md = md[0]
  endif
  return
  end