PRO budget_lev_full

  resp=''
  myinicolors
  if ( !D.NAME EQ 'X' ) then begin
     device,decomposed=0
     csize = 1.5
  endif

  units_sclfac = 1.0e9
  yax_tit = '10!U-9!N psu sec!u-1!N'

  sec_per_day = 3600.0*24.0

  exp_name = 'g.e12.G.T62_t12.001'
  cyear = '0004'
  ;; monthday_str = [ '01-05_0004-01-30', '02-04_0004-03-01', '03-06_0004-03-31', $
  ;;                  '04-05_0004-04-30', '05-05_0004-05-30', '06-04_0004-06-29', $
  ;;                  '07-04_0004-07-29', '08-03_0004-08-28', '09-02_0004-09-27', $
  ;;                  '10-02_0004-10-27', '11-01_0004-11-26', '12-01_0004-12-26']

  monthday_str = $
     ['04-20', '04-25', '04-30']
  ;; monthday_str = $
  ;;  [ '01-05', '01-10', '01-15', '01-20', '01-25', '01-30',          $
  ;;    '02-04', '02-09', '02-14', '02-19', '02-24',                   $
  ;;    '03-01', '03-06', '03-11', '03-16', '03-21', '03-26', '03-31', $
  ;;    '04-05', '04-10', '04-15', '04-20', '04-25', '04-30',          $
  ;;    '05-05', '05-10', '05-15', '05-20', '05-25', '05-30',          $
  ;;    '06-04', '06-09', '06-14', '06-19', '06-24', '06-29',          $
  ;;    '07-04', '07-09', '07-14', '07-19', '07-24', '07-29',          $
  ;;    '08-03', '08-08', '08-13', '08-18', '08-23', '08-28',          $
  ;;    '09-02', '09-07', '09-12', '09-17', '09-22', '09-27',          $
  ;;    '10-02', '10-07', '10-12', '10-17', '10-22', '10-27',          $
  ;;    '11-01', '11-06', '11-11', '11-16', '11-21', '11-26',          $
  ;;    '12-01', '12-06', '12-11', '12-16', '12-21', '12-26', '12-31'  ]

  nsamp = n_elements(monthday_str)

  path_root = '/glade/scratch/bryan/'
  path_work = '/glade/p/work/bryan/'
  grid_type = 'global_0.1_tripole'
  fprefix = ''
  fpostfix = '.pop.h.'
  file_grid = 'grid.tx01_62l.2013-07-13.nc'

  lon_box = [315.0,330.0]
  lat_box=  [19,29]
  zmax = 300.0e2

;  path_data = path_root + '/' + exp_name + '/ocn/proc/tavg/'
  path_data = path_root + '/' + exp_name + '/ocn//hist/'
  path_grid = path_work + '/' + grid_type + '/aux_data/grid/'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  ncid = ncdf_open(path_grid + file_grid)

  id = ncdf_varid(ncid,'ULAT')
  ncdf_varget,ncid,id,ulat
  id = ncdf_varid(ncid,'ULONG')
  ncdf_varget,ncid,id,ulon
  id = ncdf_varid(ncid,'TAREA')
  ncdf_varget,ncid,id,tarea

  id = ncdf_varid(ncid,'TLAT')
  ncdf_varget,ncid,id,tlat
  id = ncdf_varid(ncid,'TLONG')
  ncdf_varget,ncid,id,tlon
  id = ncdf_varid(ncid,'KMT')
  ncdf_varget,ncid,id,KMT
  id = ncdf_varid(ncid,'DZBC')
  ncdf_varget,ncid,id,DZBC
  DZBC = double(DZBC)

  id = ncdf_varid(ncid,'z_t')
  ncdf_varget,ncid,id,z_t
  id = ncdf_varid(ncid,'z_w')
  ncdf_varget,ncid,id,z_w
  id = ncdf_varid(ncid,'dz')
  ncdf_varget,ncid,id,dz

  ncdf_close,ncid

  imt = n_elements(ULON(*,0))
  jmt = n_elements(ULON(0,*))
  km = n_elements(z_t)

  get_bounds,ulon,ulat,lon_box,lat_box,i_box,j_box
  ibeg = i_box(0)
  iend = i_box(1)
  jbeg = j_box(0)
  jend = j_box(1)
  ;; ibeg = 650
  ;; iend = 655
  ;; jbeg = 1400
  ;; jend = 1405
  ilen = iend - ibeg + 1
  jlen = jend - jbeg + 1
  print, ' ibeg,iend = ',ibeg,iend,'  jbeg/jend=',jbeg,jend

  indx = where(z_t GT zmax)
  kmax = indx(0)
  klen=kmax+1
  print,'kmax=',kmax,' zt(kmax) = ',z_t(kmax),' z_w(klen)',z_w(klen)

  DZT = dblarr(imt,jmt,km)
  CELL_VOL = dblarr(imt,jmt,km)
  OMASK3D = lonarr(imt,jmt,km)

  DWORK2D = dblarr(imt,jmt)
  IWORK2D = intarr(imt,jmt)

  ;;; Create 3D cell thickness and cell volume arrays
  for k=0,km-1 do begin
     IWORK2D = IWORK2D*0
     DWORK2D(*,*) = dz(k)
     indx = where(k EQ KMT-1)
     if ( max(indx) GE 0 ) then begin
        DWORK2D(indx)=DZBC(indx)
     endif
     DZT(*,*,k) = DWORK2D

     CELL_VOL(*,*,k) = DWORK2D*TAREA

     indx = where(k LT KMT)
     IWORK2D(indx) = 1
     OMASK3D(*,*,k) = IWORK2D

     ;; print,'k=',k,' z_t=',z_t(k),$
     ;;       ' DZT min/max=',min(DZT(*,*,k)),max(DZT(*,*,k))
  endfor
  lmask_indx = where(OMASK3D LT 1)
  n_land = n_elements(lmask_indx)
  n_total = long(imt)*long(jmt)*long(km)
  n_ocean_kmt = total(KMT,/integer)
  n_ocean_omask = total(OMASK3D,/integer)

  print,' total # points = ',n_total
  print,' total # ocean points from KMT = ',n_ocean_kmt
  print,' # ocean points from OMASK = ',n_ocean_omask
  print,' # land points from OMASK = ',n_land

;;   DZT = fltarr(imt,jmt,km)
;;   CELL_VOL = fltarr(imt,jmt,km)
;;   dummy = fltarr(imt,jmt)
;;   for k=0,km-1 do begin
;;      dummy(*,*) = dz(k)
;;      indx = where(k EQ KMT-1)
;;      if ( max(indx) GE 0 ) then begin
;;         dummy(indx)=DZBC(indx)
;;      endif
;;      dzt(*,*,k) = dummy
;;      work = TAREA
;;      work(where(k GT KMT-1))=0.
;;      CELL_VOL(*,*,k) = DZT(*,*,k)*work
;; ;    print,'k=',k,' z_t=',z_t(k),' dz=',dzt(ibeg,jbeg,k),CELL_VOL(ibeg,jbeg,k)
;;   endfor

  KMTbox = KMT(ibeg:iend,jbeg:jend)
  TAREAbox = TAREA(ibeg:iend,jbeg:jend)
  DZTbox = DZT(ibeg:iend,jbeg:jend,*)

  loc_str = ' @ ' + $
            string(lon_box(0),format="(f5.1,'E')") + $
            string(lat_box(0),format="(f5.1,'N')") + $
            ' to ' + $
            string(lon_box(1),format="(f5.1,'E')") + $
            string(lat_box(1),format="(f5.1,'N')") + $
            string(long(zmax/100.),format="(' / 0-',i3,'m')")
  print,loc_str
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  MLD_area_avg = fltarr(nsamp)
  SFWF_area_avg =  fltarr(nsamp)
  FLUX_TOP_area_avg =  fltarr(nsamp)
  FLUX_TOP_vol_avg =  fltarr(nsamp)
  FLUX_BOT_vol_avg =  fltarr(nsamp)

  S_vol_avg =  fltarr(nsamp)

  H_adv_vol_avg = fltarr(nsamp)
  V_adv_vol_avg = fltarr(nsamp)
  T_adv_vol_avg = fltarr(nsamp)

  H_diff_vol_avg = fltarr(nsamp)
  V_diff_vol_avg = fltarr(nsamp)
  T_diff_vol_avg = fltarr(nsamp)
  KPP_SRC_vol_avg = fltarr(nsamp)

  RHS_vol_avg = fltarr(nsamp)
  ERR_vol_avg = fltarr(nsamp)
  TOT_TEND_vol_avg = fltarr(nsamp)

  Area_avg =  fltarr(nsamp)
  Vol_avg =  fltarr(nsamp)
  time = fltarr(nsamp)

  for n=0,nsamp-1 do begin

     date_str = cyear + '-' + monthday_str(n)

     file_name = path_data + fprefix + exp_name + fpostfix + date_str + '.nc'
     print,'file=',file_name
     ncid = ncdf_open(file_name)

     id = ncdf_varid(ncid,'time')
     ncdf_varget,ncid,id,indat
     time(n) = indat
     print,'time = ',time(n)


     ;;;; Get some scale factors
     if ( n EQ 0 ) then begin

        id = ncdf_varid(ncid,'SALT')
        result = ncdf_attinq(ncid,id,'scale_factor')
        if ( result.length NE 0 ) then begin
           ncdf_attget,ncid,id,'scale_factor',salt_scale_factor
        endif else begin
           salt_scale_factor = 1.0
        endelse

        id = ncdf_varid(ncid,'rho_fw')
        ncdf_varget,ncid,id,rho_fw
        rho_fw = rho_fw*1.0e3   ; kg/m^3

        id = ncdf_varid(ncid,'days_in_norm_year')
        ncdf_varget,ncid,id,days_per_year
        sec_per_year = sec_per_day*days_per_year

        id = ncdf_varid(ncid,'salinity_factor')
        ncdf_varget,ncid,id,salinity_factor

        print,'salinity factor before adjust=',salinity_factor
        salinity_factor = salinity_factor*1.0e3 ; convert from Model sal. untis to PSU
        print,'salinity factor after adjust=',salinity_factor

        fw_flux_units_factor = sec_per_year/rho_fw
        print,'fw_flux_units_factor = ',fw_flux_units_factor

     endif

     var = 'SALT'
     id = ncdf_varid(ncid,var)
     ncdf_varget,ncid,id,s,count=[ilen,jlen,klen,1],offset=[ibeg,jbeg,0,0]
     s = s*salt_scale_factor
     print,' s = ',min(s),' to ',max(s)

     var = 'SFWF'
     id = ncdf_varid(ncid,var)
     ncdf_varget,ncid,id,SFWF,count=[ilen,jlen,1],offset=[ibeg,jbeg,0]
     print,' SFWF = ',min(SFWF),' to ',max(SFWF)

     srf_salt_flux = sfwf*salinity_factor
     print,' SRF FLUX (salt units)= ',min(srf_salt_flux),' to ',max(srf_salt_flux)

     sfwf = sfwf*fw_flux_units_factor
     print,' SFWF (m/yr)= ',min(sfwf),' to ',max(sfwf)

     var = 'TEND_SALT'
     id = ncdf_varid(ncid,var)
     ncdf_varget,ncid,id,TOT_TEND,count=[ilen,jlen,klen,1],offset=[ibeg,jbeg,0,0]
     TOT_TEND = TOT_TEND*salt_scale_factor
     print,' TOT_TEND = ',min(TOT_TEND),' to ',max(TOT_TEND)

     id_ues = ncdf_varid(ncid,'UES')
     ncdf_varget,ncid,id_ues,RWORK3D
     RWORK3D(lmask_indx) = 0.
     print,' UES min/max=',min(RWORK3D),max(RWORK3D)
     RWORK3D = RWORK3D*CELL_VOL
     RWORK3D = -(RWORK3D - shift(RWORK3D,1,0,0))/CELL_VOL
     Z_ADV = RWORK3D(ibeg:iend,jbeg:jend,*)
     Z_ADV = Z_ADV*salt_scale_factor
     print,' Z_ADV min/max=',min(Z_ADV),max(Z_ADV)

     id_vns = ncdf_varid(ncid,'VNS')
     ncdf_varget,ncid,id_vns,RWORK3D
     RWORK3D(lmask_indx) = 0.
     print,' VNS min/max=',min(RWORK3D),max(RWORK3D)
     RWORK3D = RWORK3D*CELL_VOL
     RWORK3D = -(RWORK3D - shift(RWORK3D,0,1,0))/CELL_VOL
     M_ADV = RWORK3D(ibeg:iend,jbeg:jend,*)
     M_ADV = M_ADV*salt_scale_factor
     print,' M_ADV min/max=',min(M_ADV),max(M_ADV)

     ;; z_adv = fltarr(ilen,jlen,klen)
     ;; var = 'UES'
     ;; id = ncdf_varid(ncid,var)
     ;; for k=0,klen-1 do begin
     ;;    ncdf_varget,ncid,id,RWORK2D,count=[imt,jmt,1,1],offset=[0,0,k,0]
     ;;    RWORK2D = RWORK2D*CELL_VOL(*,*,k)
     ;;    RWORK2D = -(RWORK2D - shift(RWORK2D,1,0))/CELL_VOL(*,*,k)
     ;;    z_adv(*,*,k) = RWORK2D(ibeg:iend,jbeg:jend)
     ;; endfor
     ;; z_adv = z_adv*salt_scale_factor

     ;; m_adv = fltarr(ilen,jlen,klen)
     ;; var = 'VNS'
     ;; id = ncdf_varid(ncid,var)
     ;; for k=0,klen-1 do begin
     ;;    ncdf_varget,ncid,id,RWORK2D,count=[imt,jmt,1,1],offset=[0,0,k,0]
     ;;    RWORK2D = RWORK2D*CELL_VOL(*,*,k)
     ;;    RWORK2D = -(RWORK2D - shift(RWORK2D,0,1))/CELL_VOL(*,*,k)
     ;;    m_adv(*,*,k) = RWORK2D(ibeg:iend,jbeg:jend)
     ;; endfor
     ;; m_adv = m_adv*salt_scale_factor

     var = 'WTS'
     id = ncdf_varid(ncid,var)
     ncdf_varget,ncid,id,WTS,count=[ilen,jlen,klen+1,1],offset=[ibeg,jbeg,0,0]
     for k=0,klen do WTS(*,*,k) = WTS(*,*,k)*dz(k)
     v_adv = fltarr(ilen,jlen,klen)
     for k=0,kmax do v_adv(*,*,k) = -(WTS(*,*,k)-WTS(*,*,k+1))/DZT(ibeg:iend,jbeg:jend,k)
     v_adv = v_adv*salt_scale_factor

     z_hdif = fltarr(ilen,jlen,klen)
     var = 'HDIFE_SALT'
     id = ncdf_varid(ncid,var)
     for k=0,klen-1 do begin
        ncdf_varget,ncid,id,RWORK2D,count=[imt,jmt,1,1],offset=[0,0,k,0]
        RWORK2D = RWORK2D*CELL_VOL(*,*,k)
        RWORK2D = (RWORK2D - shift(RWORK2D,1,0))/CELL_VOL(*,*,k)
        z_hdif(*,*,k) = RWORK2D(ibeg:iend,jbeg:jend)
     endfor
     z_hdif = z_hdif*salt_scale_factor

     m_hdif = fltarr(ilen,jlen,klen)
     var = 'HDIFN_SALT'
     id = ncdf_varid(ncid,var)
     for k=0,klen-1 do begin
        ncdf_varget,ncid,id,RWORK2D,count=[imt,jmt,1,1],offset=[0,0,k,0]
        RWORK2D = RWORK2D*CELL_VOL(*,*,k)
        RWORK2D = (RWORK2D - shift(RWORK2D,0,1))/CELL_VOL(*,*,k)
        m_hdif(*,*,k) = RWORK2D(ibeg:iend,jbeg:jend)
     endfor
     m_hdif = m_hdif*salt_scale_factor

     var = 'DIA_IMPVF_SALT'
     id = ncdf_varid(ncid,var)
     ncdf_varget,ncid,id,VMIX_FLUX,count=[ilen,jlen,klen+1,1],offset=[ibeg,jbeg,0,0]
     print,' VMIX_FLUX = ',min(VMIX_FLUX),' to ',max(VMIX_FLUX)
     VMIX_FLUX = VMIX_FLUX*salt_scale_factor

     v_diff = fltarr(ilen,jlen,klen)
     k = 0
     v_diff(*,*,k) = (srf_salt_flux - VMIX_FLUX(*,*,k))/DZT(*,*,k)
     for k=1,klen-1 do v_diff(*,*,k) = (VMIX_FLUX(*,*,k-1) - VMIX_FLUX(*,*,k))/DZT(*,*,k)

     var = 'KPP_SRC_SALT'
     id = ncdf_varid(ncid,var)
     ncdf_varget,ncid,id,KPP_SRC,count=[ilen,jlen,klen,1],offset=[ibeg,jbeg,0,0]
     print,' KPP_SRC = ',min(KPP_SRC),' to ',max(KPP_SRC)
     KPP_SRC = KPP_SRC*salt_scale_factor
     
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

     h_adv = z_adv + m_adv
     t_adv = h_adv + v_adv
     h_diff = z_hdif + m_hdif
     t_diff = v_diff + kpp_src + h_diff
     rhs = t_adv + t_diff
     err = tot_tend - rhs

     if ( n EQ 0 ) then begin
        !P.MULTI=[0,1,2]
        plot,z_t(0:klen-1)/100.,z_adv(0,0,*)*units_sclfac,title='Zonal Adv'
        plot,z_t(0:klen-1)/100.,m_adv(0,0,*)*units_sclfac,title='Merid Adv'

        read,prompt='Continue ?',resp

        plot,z_t(0:klen-1)/100.,v_adv(0,0,*)*units_sclfac,title='Vert Adv'
        plot,z_t(0:klen-1)/100.,t_adv(0,0,*)*units_sclfac,title='Total Adv'

        read,prompt='Continue ?',resp

        plot,z_t(0:klen-1)/100.,z_hdif(0,0,*)*units_sclfac,title='Zonal Hdiff'
        plot,z_t(0:klen-1)/100.,m_hdif(0,0,*)*units_sclfac,title='Merid Hdiff'

        read,prompt='Continue ?',resp

        plot,z_t(0:klen-1)/100.,h_diff(0,0,*)*units_sclfac,title='Total Hdiff'
        plot,z_t(0:klen-1)/100.,kpp_src(0,0,*)*units_sclfac,title='KPP SRC'

        read,prompt='Continue ?',resp

        plot,z_t(0:klen-1)/100.,v_diff(0,0,*)*units_sclfac,title='Vert diff'
        plot,z_t(0:klen-1)/100.,t_diff(0,0,*)*units_sclfac,title='Tot diff'

        read,prompt='Continue ?',resp

        plot,z_t(0:klen-1)/100.,tot_tend(0,0,*)*units_sclfac,title='TOT_TEND'
        plot,z_t(0:klen-1)/100.,rhs(0,0,*)*units_sclfac,title='RHS'
        oplot,z_t(0:klen-1)/100.,tot_tend(0,0,*)*units_sclfac,linestyle=2

        read,prompt='Continue ?',resp

        plot,z_t(0:klen-1)/100.,err(0,0,*)*units_sclfac,title='ERR'
        plot,z_t(0:klen-1)/100.,err(0,0,*)/abs(tot_tend(0,0,*)),/ystyle,yrange=[-0.025,0.025],$
             title='Relative error'

        read,prompt='Continue ?',resp
        !P.MULTI=0
     endif

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

     omask = KMTbox
     omask(where(omask GE 1))=1

     work = TAREAbox

     Area_avg(n) = total(work,/DOUBLE)
     SFWF_area_avg(n) = total(SFWF*work,/DOUBLE)/Area_avg(n)
     FLUX_TOP_area_avg(n) = total(srf_salt_flux(*,*,0)*work,/DOUBLE)/Area_avg(n)

     for k=0,kmax do begin
        print,' summing k=',k,' zt=',z_t(k),$
              ' zw(k)=',z_w(k),' z_w(k+1)=',z_w(k+1)

        indx = where(KMTbox LT k+1) 
        if ( max(indx) GE 0 ) then omask(indx)=0.

        indx = where(omask GT 0.)
        if ( max(indx) LT 0 ) then break
        n_in = n_elements(indx)

        print,'k=',k,' # in ML = ',n_in

        work = TAREAbox*DZTbox(*,*,k)*omask

        Vol_avg(n) = Vol_avg(n) + total(work,/DOUBLE)
        S_vol_avg(n) = S_vol_avg(n) + total(s(*,*,k)*work,/DOUBLE)

        TOT_TEND_vol_avg(n) = TOT_TEND_vol_avg(n) + total(TOT_TEND(*,*,k)*work,/DOUBLE)
        RHS_vol_avg(n) = RHS_vol_avg(n) + total(RHS(*,*,k)*work,/DOUBLE)
        ERR_vol_avg(n) = ERR_vol_avg(n) + total(err(*,*,k)*work,/DOUBLE)

        H_adv_vol_avg(n) = H_adv_vol_avg(n) +total(H_adv(*,*,k)*work,/DOUBLE)
        V_adv_vol_avg(n) = V_adv_vol_avg(n) +total(V_adv(*,*,k)*work,/DOUBLE)
        T_adv_vol_avg(n) = T_adv_vol_avg(n) +total(T_adv(*,*,k)*work,/DOUBLE)

        V_diff_vol_avg(n) = V_diff_vol_avg(n) +total(V_diff(*,*,k)*work,/DOUBLE)
        H_diff_vol_avg(n) = H_diff_vol_avg(n) +total(h_diff(*,*,k)*work,/DOUBLE)
        kpp_src_vol_avg(n) = kpp_src_vol_avg(n) +total(kpp_src(*,*,k)*work,/DOUBLE)
        T_diff_vol_avg(n) = T_diff_vol_avg(n) +total(T_diff(*,*,k)*work,/DOUBLE)
     endfor
     S_vol_avg(n) = S_vol_avg(n)/Vol_avg(n)

     H_ADV_vol_avg(n) = H_ADV_vol_avg(n)/Vol_avg(n)
     V_ADV_vol_avg(n) = V_ADV_vol_avg(n)/Vol_avg(n)
     T_ADV_vol_avg(n) = T_ADV_vol_avg(n)/Vol_avg(n)

     H_DIFF_vol_avg(n) = H_DIFF_vol_avg(n)/Vol_avg(n)
     V_DIFF_vol_avg(n) = V_DIFF_vol_avg(n)/Vol_avg(n)
     KPP_SRC_vol_avg(n) = KPP_SRC_vol_avg(n)/Vol_avg(n)
     T_DIFF_vol_avg(n) = T_DIFF_vol_avg(n)/Vol_avg(n)

     RHS_vol_avg(n) = RHS_vol_avg(n)/Vol_avg(n)
     TOT_TEND_vol_avg(n) = TOT_TEND_vol_avg(n)/Vol_avg(n)
     ERR_vol_avg(n) = ERR_vol_avg(n)/Vol_avg(n)

     FLUX_TOP_vol_avg(n) = Area_avg(n)*FLUX_TOP_area_avg(n)/Vol_avg(n)
     FLUX_BOT_vol_avg(n) = FLUX_TOP_vol_avg(n) - V_DIFF_vol_avg(n)

     print,' Area_avg   = ',Area_avg(n)*(1.0e-5)^2,' km^2'
     print,' Vol_avg(n) = ',Vol_avg(n)*(1.0e-5)^3,' km^3'

     print,' SFWF_area_avg = ',SFWF_area_avg(n), 'm/yr'
     print,' FLUX_TOP_area_avg = ',FLUX_TOP_area_avg(n)
     print,' '
     print,' S_vol_avg = ',S_vol_avg(n)
     print,' '
     print,' H_ADV_vol_avg = ',H_ADV_vol_avg(n)
     print,' V_ADV_vol_avg = ',V_aDV_vol_avg(n)
     print,' T_ADV_vol_avg = ',T_ADV_vol_avg(n)
     print,' '
     print,' H_Diff_vol_avg  = ',H_Diff_vol_avg(n)
     print,' V_Diff_vol_avg  = ',V_Diff_vol_avg(n)
     print,' KPP_Src_vol_avg = ',KPP_Src_vol_avg(n)
     print,' T_Diff_vol_avg  = ',T_Diff_vol_avg(n)
     print,' '
     print,' FLUX_TOP_vol_avg = ',FLUX_TOP_vol_avg(n)
     print,' FLUX_BOT_vol_avg = ',FLUX_BOT_vol_avg(n)
     print,' '
     print,' TOT_TEND_vol_avg = ',TOT_TEND_vol_avg(n)
     print,' RHS_vol_avg = ',RHS_vol_avg(n)
     print,' ERR_vol_avg = ',ERR_vol_avg(n)

  endfor

  year0 = long(time(0)/365.0)
  print,'year0=',year0
  time = time - year0*365
  t_range=[min(time),max(time)]
  print,'t range=',t_range

  !P.MULTI=[0,1,2]
  dummy = fltarr(nsamp)
  print,S_vol_avg
  plot,time,S_vol_avg,thick=2,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       ytitle='psu',ythick=2,$
       /ynozero,$
       title = 'Salinity ' + loc_str,$
       charsize=csize,charthick=2

  plot,time,-(SFWF_area_avg-mean(SFWF_area_avg)),thick=2,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       ytitle='m/yr',ythick=2,$
       /ynozero,$
       title = 'E-P (anom)' + loc_str,$
       charsize=csize,charthick=2

  if (!D.NAME EQ 'X' ) then read,prompt='Continue?',resp


  plot,time,TOT_TEND_vol_avg*units_sclfac,thick=1,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       /ynozero,ytitle=yax_tit,ythick=2,$
       title = 'Salinity Tendency' + loc_str,$
       charsize=csize,charthick=2
  oplot,time,RHS_vol_avg*units_sclfac,thick=2,linestyle=4

  plot,time,ERR_vol_avg*units_sclfac,thick=2,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       /ynozero,ytitle=yax_tit,ythick=2,$
       title = 'ERROR' + loc_str,$
       charsize=csize,charthick=2

  if (!D.NAME EQ 'X' ) then read,prompt='Continue?',resp

  plot,time,H_ADV_vol_avg*units_sclfac,thick=2,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       /ynozero,ytitle=yax_tit,ythick=2,$
       title = 'Horiz. Adv.' + loc_str,$
       charsize=csize,charthick=2

  plot,time,V_ADV_vol_avg*units_sclfac,thick=2,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       /ynozero,ytitle=yax_tit,ythick=2,$
       title = 'Vert. Adv.' + loc_str,$
       charsize=csize,charthick=2

  if (!D.NAME EQ 'X' ) then read,prompt='Continue?',resp

  plot,time,H_Diff_vol_avg*units_sclfac,thick=2,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       /ynozero,ytitle=yax_tit,ythick=2,$
       title = 'Horiz. Diff.' + loc_str,$
       charsize=csize,charthick=2

  plot,time,V_Diff_vol_avg*units_sclfac,thick=2,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       /ynozero,ytitle=yax_tit,ythick=2,$
       title = 'Vert. Diff.' + loc_str,$
       charsize=csize,charthick=2

  if (!D.NAME EQ 'X' ) then read,prompt='Continue?',resp

  plot,time,T_ADV_vol_avg*units_sclfac,thick=2,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       /ynozero,ytitle=yax_tit,ythick=2,$
       title = 'Total Adv.' + loc_str,$
       charsize=csize,charthick=2

  plot,time,T_DIFF_vol_avg*units_sclfac,thick=2,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       /ynozero,ytitle=yax_tit,ythick=2,$
       title = 'Total DIFF.' + loc_str,$
       charsize=csize,charthick=2

  if (!D.NAME EQ 'X' ) then read,prompt='Continue?',resp

  plot,time,FLUX_TOP_vol_avg*units_sclfac,thick=2,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       /ynozero,ytitle=yax_tit,ythick=2,$
       title = 'DIFF FLUX TOP' + loc_str,$
       charsize=csize,charthick=2

  plot,time,-FLUX_BOT_vol_avg*units_sclfac,thick=2,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       /ynozero,ytitle=yax_tit,ythick=2,$
       title = 'DIFF FLUX BOT' + loc_str,$
       charsize=csize,charthick=2

  if (!D.NAME EQ 'X' ) then read,prompt='Continue?',resp

  !P.MULTI=0

  dummy = fltarr(nsamp)
  for n=0,nsamp-1 do dummy(n) = max([abs(TOT_TEND_vol_avg(n)),$
                                     abs(FLUX_TOP_vol_avg(n)),$
                                     abs(FLUX_BOT_vol_avg(n)),$
                                     abs(T_ADV_vol_avg(n))])*units_sclfac
  plot,time,TOT_TEND_vol_avg*units_sclfac,/nodata,$
       /xstyle,xrange=t_range,xtitle='Day of Year',xthick=2,$
       /ystyle,yrange=[-max(dummy),max(dummy)],ytitle=yax_tit,ythick=2,$
       title = 'Salinity Tendency' + loc_str,$
       charsize=csize,charthick=2
  oplot,time,TOT_TEND_vol_avg*units_sclfac,thick=2,color=6
  oplot,time,RHS_vol_avg*units_sclfac,thick=2,color=0,linestyle=2
  oplot,time,FLUX_TOP_vol_avg*units_sclfac,thick=2,color=4
  oplot,time,-FLUX_BOT_vol_avg*units_sclfac,thick=2,color=4,linestyle=2
  oplot,time,(T_ADV_vol_avg+H_Diff_vol_avg)*units_sclfac,thick=2,color=5
  oplot,time,T_ADV_vol_avg*units_sclfac,thick=2,color=5
  oplot,time,H_Diff_vol_avg*units_sclfac,thick=2,color=5,linestyle=1


  plots,t_range,[0,0],linestyle=1
  return


end