pro storage_lc
!p.multi=[0,2,2]
print,'Enter start date in yyyymmdd format:'
startdate=0L
read,startdate
if startdate lt 19900000L then begin
   print,'Sorry. Please enter start date in 4-digit year format: yyyymmdd.'
   goto,fini
endif

print,'Enter number of days to calculate storage over:'
print,'You must make sure that data files exist in the period you select'
ndays=0
read,ndays

print,'Do you wish to save plots: y/n?'
plt=''
read,plt
if plt eq 'y' then begin
   plt=1
endif else begin
   plt=0
endelse

print,'Do you wish to save  results to storage files: y/n?'
sav=''
read,sav
if sav eq 'y' then begin
   sav=1
endif else begin
   sav=0
endelse

IF sav EQ 1 THEN BEGIN

print,'please select way of saving,'
print,' 0--simple adding in the end,1--update one to one 2--directly write'

way=0
read,way

endif


date = startdate

top:

;Get strings and integers related to the date chosen
yyyymmdd = strtrim(string(date),2) ;string of startdate
yrss = strmid(yyyymmdd,2,2) ;string of 2-digit year
yrls = strmid(yyyymmdd,0,4) ;string of 4-digit year
yrl=fix(yrls)               ;integer of 4-digit year
jday = jdate(date)          ;integer jday
sndays=strtrim(string(ndays),2); string of ndays


print,'calculated Date:',' yyyymmdd=',yyyymmdd
print,'jday=',jday
print,'sndays=',sndays


;initialize the halt program flag
hltprgm=0


;   cs - `slow' co2  from Vaisala (ppm)
;   r - a vector of indicies telling where all the good data are
;   tm - the time vector (in seconds from the start time

; Load the data

storload_lc,jday,yrss,ndays,cs,hltprgm





;------------------------------------------------------------------------------
; This lets user choose a new date - see bottom of this prgm

if hltprgm eq 1 then begin
    print,'warning: there is data full of badval'
;   print,'DATA PROBLEMS FOR DATES CHOSEN'
;   goto, fini
endif
;------------------------------------------------------------------------------
badval=-999.0
c=fltarr(6)
total=480L*ndays
cint=fltarr(total)
dcint=fltarr(total)
nave=total/10
average=fltarr(nave)+badval
cint_30=fltarr(nave)+badval

uair=28.9


h=[0,1.2,1.8,3.0,5.5,10.2];level height (m)


print,'total=',total

;print,cs


for ntot=0L,total-1 do begin
   for i=1,5 do begin
    c(i)=cs(1+i,ntot)
   endfor
    c(0)=c(1)

;integration
; count badval
   badflag=0
    for i=0,5 do begin
     if(c(i) eq badval)then begin
      badflag=badflag+1
     endif
    endfor

     if(badflag ge 3)then begin
      sum=badval
      goto, intsum
     endif else begin
;;;;;;;         for j=0,7 do begin
          sum1=0.0
          nsum1=0.0
          for j=0,5 do begin
          if(c(j) ne badval) then begin
           sum1=sum1+c(j)
           nsum1=nsum1+1
          endif
         endfor

          sum1=sum1/nsum1
         for j=0,5 do begin
          if(c(j) eq badval)then begin
           c(j)=sum1
          endif
         endfor
     endelse      

   sum=0.0
    for i=1,5 do begin
     sum=sum+(c(i-1)+c(i))*(h(i)-h(i-1))/2.0
    endfor

intsum: cint(ntot)=sum
endfor

print,cint


;calculate 30-min-averaged co2 concentration

  for i=0L,nave-1 do begin

    i1=i*10
    i2=i1+10-1
    cint_30(i)=meanbadval(cint(i1:i2))
  endfor
    

; cint is the integration of every 3 mins.

; calculate the d(cint)/dt

;;  for i=1L,total-1 do begin
;;   if(cint(i) ne badval and cint(i-1) ne badval) then begin
;;   dcint(i)=(cint(i)-cint(i-1))/(3.0*60.0)
;;   endif else begin
;;   dcint(i)=badval
;;   endelse
;;   endfor

; set dcint(0) equal to dcint(1)
;;   dcint(0)=dcint(1)

; calculate 30-min average

;;   for i=1L,nave do begin
 ;;         sum2=0.0
 ;;         nsum2=0.0
;;        for j=(i-1)*10+1,(i)*10 do begin
;;         kk=j-1
;;;         sum2=0.0
;;;         nsum2=0.0
;;          if(dcint(kk) ne badval)then begin
;;           nsum2=nsum2+1
;;           sum2=sum2+dcint(kk)
 ;;         endif
;;           endfor

;;          if(nsum2 ge 1)then begin
;;            average(i-1)=sum2/nsum2
;;          endif else begin
;;            average(i-1)=badval
;;          endelse
;;   endfor


       for i=1L,nave-1 do begin

         if(cint_30(i) ne badval and cint_30(i-1) ne badval) then begin
           average(i)=(cint_30(i)-cint_30(i-1))/(30.0*60.0)
         endif
       endfor




;        average(0)=average(1)
;00, 30 , 00, 30, 

; convert unit,

;load air temp, air pressure, mixr....

print,sndays
help,yrss
help,jday
help,sndays

;hrmet_wc,yrss,jday,sndays,airp,mixr
; changed to hrmet_lc_psu,  dec17 2000 by wwg
hrmet_lc_psu,yrss,jday,sndays,airp,mixr

;Mar 22, 2002, use licor mixratio instead of bad mairatio of water
;from qmix file


;Mar 25, 2002, use lic temp replace bad temperature
loadlicmixr_lc,yrls,jday,ndays,mixr_lic,temp_lic



;print,mixr
;compute the airtemp every half-hour

storairtemp_lc,jday,yrss,ndays,airt

plot,mixr,min_val=-100

print,'airt',airt(0:47)

   for ind=0,nave-1 do begin

     airpress=airp(ind)
     mixrat=mixr(ind)
     IF(mixrat EQ -999.) THEN mixrat=mixr_lic(ind)
     airtemp=airt(ind)    
      IF(airtemp EQ -999.) THEN airtemp=temp_lic(ind)
print,airpress,mixrat,mixr_lic(ind),airtemp,ind
; Compute dry air density (DRHO - kg/m^3) averaged over the time of
; the flux calculation
   ddensity,airpress,mixrat,airtemp,drho,badval         
print,drho,airtemp,mixrat,airpress
;convert ppm---> micromole/m2/s
   if(average(ind) ne badval and drho ne badval) then begin
    average(ind)=average(ind)*drho/uair*10^3
   endif else begin
    average(ind)=badval
   endelse

   endfor





; calculate the year, julian days and time to save

sa_yr=intarr(nave)
sa_day=intarr(nave)
sa_time=fltarr(nave)
atime=0.0-0.5
aday=0
   for i=0,nave-1 do begin
    sa_yr(i)=yrl
     atime=atime+0.5
     if(atime eq 24) then begin
      aday=aday+1
      atime=0
     endif
    sa_day(i)=aday+jday
    sa_time(i)=atime
print,sa_yr(i),sa_day(i),sa_time(i),average(i)

    endfor

    
           
if (sav eq 1) then begin
  storsave_wc,yrls,sa_yr,sa_day,sa_time,average,way
endif

;***************************************************************
;Plot fits in .ps file.
if plt eq 1 then begin

storplot_lc,plt,yrss,jday,ndays,average

endif



fini: return
end