! monthly_report_tmp.f90
!   Mar. 2010,  by I. Tamagawa
!     - Apr.2010
!   Jan. 2011      AWS1 added Snow depth and 2 more data column
!   Dec. 2011      AWS1 add Tg WC Qsu

! make monthly report of TAKAYAMA AWS
!  http://sateco-archive.green.gifu-u.ac.jp/AWS1/

module typedef
 type date_time
   integer :: year
   integer :: month
   integer :: day
   integer :: hour
   integer :: min
   integer :: sec
 end type

 type stat
  real :: max
  real :: min
  real :: ave
  double precision  :: sum
  integer :: num
 end type

 integer, parameter :: maxch = 50
 type awsdata
  integer :: minute
  character(len=80) :: date
  real :: dat( maxch )
 end type

 type type_aveparam
  integer:: prev
  integer:: foll
 end type


 integer, parameter :: hist_nmax = 100
 type histdata
  real    :: dx, xmin, xmax
  integer :: hist_n
  integer :: hist(hist_nmax)
 end type

 integer, parameter :: windrose_class=16
 type type_windrose
  integer :: n = windrose_class
  real    :: sum_ws(windrose_class)
 end type

end module

module mod_nodata
implicit none
 integer, parameter :: NODATA = -9999
 logical, external  :: is_nodata
end module

logical function is_nodata( data )
implicit none
 integer, parameter :: NODATA = -9999
  real :: data
  real, parameter :: eps = 0.5

  if( abs( data - NODATA ) < eps ) then
      is_nodata = .TRUE.
  else
      is_nodata = .FALSE.
  end if
end function


module fileinfo
implicit none

character(len=1024) :: file1 = '00_All_AWS1_10Min.csv' 
integer,parameter :: ch1 = 12, shift1 = 0
integer,parameter ::   ws      =  2+shift1 ,  wdir      = 3+shift1 
integer,parameter ::  ta      =  4+shift1 ,  rh        = 5+shift1
integer,parameter ::  qs      =  6+shift1 
integer,parameter ::  prec_1  =  7+shift1 ,  press     = 8+shift1
integer,parameter ::  prec_2  =  9+shift1
integer,parameter ::    snow    = 10 + shift1, qsnow = 12+shift1
integer,parameter :: tg10  = 13+shift1, tg20  =14+shift1
integer,parameter :: qsu  = 15+shift1
integer,parameter :: sm10  = 17+shift1, sm20  = 18+shift1


integer,parameter ::  u       = 19+shift1 ,  v         = 20+shift1
integer,parameter ::  sphum   = 21+shift1 ,  vappress  = 22+shift1
integer,parameter ::  prec24_1= 23+shift1 ,  prec24_2  = 24+shift1
integer,parameter ::  prec2_1 = 25+shift1 ,  prec2_2   = 26+shift1
integer,parameter ::  prec1_1 = 27+shift1 ,  prec1_2   = 28+shift1
integer,parameter ::    ws_v    = 29+shift1 , wdir_v     = 30+shift1
integer,parameter :: albedo  = 31+shift1


character(len=1024) :: file2 = '00_All_AWS2_10Min.csv'

integer,parameter :: ch2 = 5, shift2 = 35
integer,parameter :: qs2     = 2+shift2  , ir         = 4+shift2


integer,parameter    :: dt=10  ! 10minutes
end module


program main
use typedef
use mod_nodata
use fileinfo
implicit none

integer :: command_argument_count
integer :: argc
character(len=80) :: argv
integer :: year, month                   ! Target Year and Month
integer :: duration_min, extent_min, days
type(date_time) :: start

integer :: total_ch, total_num,  sduration
type(awsdata), allocatable  :: data(:)
type(stat), allocatable :: sdata(:,:), msdata(:,:)

integer :: i, j

type(histdata) :: hist(maxch)
type(type_windrose) :: windrose

 extent_min = 24*60  ! extent of data duration for moving average etc...
 argc = command_argument_count()
 if( argc >= 2 ) then
  call get_command_argument(1, argv)
  read(argv, *) year
  call  get_command_argument(2, argv)
  read(argv, *) month
 else
  write(0,*) 'Please add year and month like  ./a.out 2010 3'
  stop
 end if
 write(0,*) 'year/momth =', year, '/', month
 call setparam( year, month, start, duration_min, days )


 total_num = ( duration_min + 2*extent_min )/dt
 allocate( data(1:total_num) )
 allocate( sdata(1:days, 1:maxch) )
 allocate( msdata(1, 1:maxch) )

 call initdata( data, total_num )

 call read_copy( file1, ch1, shift1,  &
   start, duration_min, extent_min, dt, total_num, data )

do i=1, total_num
 call checksnow(  data(i)%dat(snow), data(i)%dat(qsnow) )
end do

 call read_copy( file2, ch2, shift2,  &
   start, duration_min, extent_min, dt, total_num, data )

! call testout( data, total_num, ch1, shift1, ch2, shift2 )


 call mkcalcdata( data, total_num )  

! make statistics (daily)
 write(0,*) 'make statistics'
 sduration = 24*60
 call mkstat( data, total_num, sdata, days, duration_min, sduration )
! make daily albedo
 call daily_albedo( sdata, days )



 sduration = 24*60*days
 call mkstat( data, total_num, msdata, 1, duration_min, sduration )
 call daily_albedo( msdata, 1 )


! make histgrams

 write(0,*) 'make histgrams'
 call mkhist(  data, total_num, 0, duration_min, hist )

 write(0,*) 'make windrose'
 call mkwindrose( data, total_num, 0, duration_min, windrose)

! output

call fileout( year, month, days, sdata, msdata, hist, windrose )

end


subroutine fileout( year, month, days, sdata, msdata, hist, windrose )
use fileinfo
use typedef
use mod_nodata
implicit none

integer :: year, month, days
type(stat)  :: sdata(days,maxch), msdata(1,maxch)
type(histdata) :: hist(maxch)
type(type_windrose) :: windrose
character(len=1024) :: out, direc, command, out1, outm, outd, out2
integer :: iday, i



write(direc,'(I0.4A1I0.2)') year,'/',month
command = 'mkdir -p '

command = (command(1:len_trim(command))//' ')//direc(1:len_trim(direc))
call system( command )

! output monthly statistics
write(out, '(A7,A1,I0.4A1I0.2A)') direc,'/',year,'_', month, '_montly_stat.csv'
write(outm, '(I0.4A1I0.2A)') year,'_', month, '_montly_stat.csv'
write(0,*) 'output file name=', out(1:len_trim(out))
open(unit=11, file=out(1:len_trim(out)))
write(11,'(a)') '"year","month",&
  &"wind_spped(min)","wind_speed(ave)","wind_speed(max)",&
  &"wind_speed(vec_ave)","wind_direc(vec_ave)",&
  &"Air_temp(min)","Air_temp(ave)","Air_temp(Max)",&
  &"Rel_Hum(min)","Rel_Hum(ave)","Rel_Hum(max)",&
  &"Spec_Hum(min)","Spec_Hum(ave)","Spec_Hum(max)",&
  &"Vap_Pres(min)","Vap_Pres(ave)","Vap_Pres(max)",&
  &"Solor_Rad(ave)","Solor_Rad(max)",&
  &"Precise_Solor_Rad(ave)",&
        &"Precise_Solor_Rad(max)",&
  &"Longwave_Rad(min)","Longwave_Rad(ave)","Longwave_Rad(max)",&
  &"Press(min)","Press(ave)","Press(max)",&
  &"Rainfall(sum)","Rainfall(max/10min.)",&
  &"Precipitation(sum)","Precipitation(max/10min.)",&
  &"Rainfall(max/1hour)","Precipitation(max/1hour)",&
  &"Rainfall(max/2hour)","Precipitation(max/2hour)",&
  &"Rainfall(max/24hour)","Precipitation(max/24hour)",&
  &"Snow depth(min)","Snow depth(ave)","Snow depth(max)",&
  &"Ground_Temp_10cm(min)","Ground_Temp_10cm(ave)","Ground_Temp_10cm(max)",&
  &"Ground_Temp_20cm(min)","Ground_Temp_20cm(ave)","Ground_Temp_20cm(max)",&
  &"WaterContent_10cm(min)","WaterContent_10cm(ave)","WaterContent_10cm(max)",&
  &"WaterContent_20cm(min)","WaterContent_20cm(ave)","WaterContent_20cm(max)",&
  &"Upward SW radiation(ave)","Upward SW radiation(max)",&
  &"Daily albedo"'




write(11,'(a)') ',"(1-12)","(m/s)","(m/s)","(m/s)",&
    &"(m/s)","(deg)",&
    &"(C deg)","(C deg)","(C deg)",&
    &"(%)","(%)","(%)",&
    &"(g/kg)","(g/kg)","(g/kg)",&
    &"(hPa)","(hPa)","(hPa)",&
    &"(W/m2)","(W/m2)",&
    &"(W/m2)","(W/m2)",&
    &"(W/m2)","(W/m2)","(W/m2)",&
    &"(hPa)","(hPa)","(hPa)",&
    &"(mm/month)","(mm/10min)",&
    &"(mm/month)","(mm/10min)",&
    &"(mm/1hour)","(mm/1hour)",&
    &"(mm/2hour)","(mm/2hour)",&
    &"(mm/24hour)","(mm/24hour)"&
    &"(cm)","(cm)","(cm)"&
    &"(C deg)","(C deg)","(C deg)",&
    &"(C deg)","(C deg)","(C deg)",&
    &"(%)","(%)","(%)",&
    &"(%)","(%)","(%)",&
    &"(W/m2)","(W/m2)",&
    &"(-)"'








write(11, '(I4,",",I2,54(",",F10.1),F10.2)') year, month, &
  msdata(1, ws)%min, msdata(1, ws)%ave, msdata(1, ws)%max, &
  msdata(1, ws_v)%ave, msdata(1, wdir_v)%ave,              &
  msdata(1, ta)%min, msdata(1, ta)%ave, msdata(1, ta)%max, &
  msdata(1, rh)%min, msdata(1, rh)%ave, msdata(1, rh)%max, &
  1000*msdata(1, sphum)%min, 1000*msdata(1, sphum)%ave, &
                                  1000*msdata(1, sphum)%max, &
  msdata(1, vappress)%min, msdata(1, vappress)%ave, &
          msdata(1, vappress)%max, &
        msdata(1, qs)%ave, msdata(1,qs)%max, &
        msdata(1, qs2)%ave, msdata(1,qs2)%max, &
  msdata(1,ir)%min, msdata(1,ir)%ave, msdata(1,ir)%max, &
  msdata(1,press)%min, msdata(1,press)%ave, msdata(1,press)%max, &
  msdata(1, prec_1)%sum, msdata(1, prec_1)%max, &
  msdata(1, prec_2)%sum, msdata(1, prec_2)%max, &
  msdata(1, prec1_1)%max,  msdata(1, prec1_2)%max, &
  msdata(1, prec2_1)%max, msdata(1, prec2_2)%max, &
  msdata(1, prec24_1)%max, msdata(1, prec24_2)%max, &
  msdata(1, snow)%min, msdata(1, snow)%ave, msdata(1, snow)%max, &
  msdata(1, tg10)%min, msdata(1, tg10)%ave, msdata(1,tg10)%max, &
  msdata(1, tg20)%min, msdata(1, tg20)%ave, msdata(1,tg20)%max, &
  msdata(1, sm10)%min, msdata(1, sm10)%ave, msdata(1, sm10)%max, &
  msdata(1, sm20)%min, msdata(1, sm20)%ave, msdata(1, sm20)%max, &
  msdata(1, qsu)%ave, msdata(1, qsu)%max, &
  msdata(1, albedo)%ave



close(11)


! output daily statistics
write(out, '(A7,A1,I0.4A1I0.2A)') direc,'/',year,'_', month, '_daily_stat.csv'
write(outd, '(I0.4A1I0.2A)') year,'_', month, '_daily_stat.csv'
write(0,*) 'output file name=', out(1:len_trim(out))

open(unit=11, file=out(1:len_trim(out)))
write(11,'(a)') '"year", "month","day",&
  &"wind_spped(min)","wind_speed(ave)","wind_speed(max)",&
  &"wind_speed(vec_ave)","wind direc(vec_ave)",&
  &"Air_temp(min)","Air_temp(ave)","Air_temp(Max)",&
  &"Rel_Hum(min)","Rel_Hum(ave)","Rel_Hum(max)",&
  &"Spec_Hum(min)","Spec_Hum(ave)","Spec_Hum(max)",&
  &"Vap_Pres(min)","Vap_Pres(ave)","Vap_Pres(max)",&
  &"Solor_Rad(ave)","Solor_Rad(max)",&
  &"Precise_Solor_Rad(ave)",&
        &"Precise_Solor_Rad(max)",&
  &"Longwave_Rad(min)","Longwave_Rad(ave)","Longwave_Rad(max)",&
  &"Press(min)","Press(ave)","Press(max)",&
  &"Rainfall(sum)","Rainfall(max/10min.)",&
  &"Precipitation(sum)","Precipitation(max/10min.)",&
  &"Rainfall(max/1hour)","Precipitation(max/1hour)",&
  &"Rainfall(max/2hour)","Precipitation(max/2hour)",&
  &"Rainfall(max/24hour)","Precipitation(max/24hour)",&
  &"Snow depth(min)","Snow depth(ave)","Snow depth(max)",&
  &"Ground_Temp_10cm(min)","Ground_Temp_10cm(ave)","Ground_Temp_10cm(max)",&
  &"Ground_Temp_20cm(min)","Ground_Temp_20cm(ave)","Ground_Temp_20cm(max)",&
  &"WaterContent_10cm(min)","WaterContent_10cm(ave)","WaterContent_10cm(max)",&
  &"WaterContent_20cm(min)","WaterContent_20cm(ave)","WaterContent_20cm(max)",&
  &"Upward SW radiation(ave)","Upward SW radiation(max)",&
  &"Daily albedo"'




write(11,'(a)') ',"(1-12)",,"(m/s)" ,"(m/s)","(m/s)",&
    &"(m/s)","(deg)",&
    &"(C deg)","(C deg)","(C deg)",&
    &"(%)","(%)","(%)",&
    &"(g/kg)","(g/kg)","(g/kg)",&
    &"(hPa)","(hPa)","(hPa)",&
    &"(W/m2)","(W/m2)",&
    &"(W/m2)","(W/m2)",&
    &"(W/m2)","(W/m2)","(W/m2)",&
    &"(hPa)","(hPa)","(hPa)",&
    &"(mm/day)","(mm/10min)",&
    &"(mm/day)","(mm/10min)",&
    &"(mm/1hour)","(mm/1hour)",&
    &"(mm/2hour)","(mm/2hour)",&
    &"(mm/24hour)","(mm/24hour)"&
    &"(cm)","(cm)","(cm)"&
    &"(C deg)","(C deg)","(C deg)",&
    &"(C deg)","(C deg)","(C deg)",&
    &"(%)","(%)","(%)",&
    &"(%)","(%)","(%)",&
    &"(W/m2)","(W/m2)",&
    &"(-)"'



do iday=1, days


write(11, '(I4,",",I3, ",", I3, 54(",",F10.1),",",F10.2)') year, month, iday, &
  sdata(iday, ws)%min, sdata(iday, ws)%ave, sdata(iday, ws)%max, &
  sdata(iday, ws_v)%ave, sdata(iday, wdir_v)%ave,              &
  sdata(iday, ta)%min, sdata(iday, ta)%ave, sdata(iday, ta)%max, &
  sdata(iday, rh)%min, sdata(iday, rh)%ave, sdata(iday, rh)%max, &
  1000*sdata(iday, sphum)%min, 1000*sdata(iday, sphum)%ave, &
                                  1000*sdata(iday, sphum)%max, &
  sdata(iday, vappress)%min, sdata(iday, vappress)%ave, &
          sdata(iday, vappress)%max, &
        sdata(iday, qs)%ave, sdata(iday,qs)%max, &
        sdata(iday, qs2)%ave, sdata(iday,qs2)%max, &
  sdata(iday,ir)%min, sdata(iday,ir)%ave, sdata(iday,ir)%max, &
  sdata(iday, press)%min, sdata(iday, press)%ave,sdata(iday, press)%max,&
  sdata(iday, prec_1)%sum, sdata(iday, prec_1)%max, &
  sdata(iday, prec_2)%sum, sdata(iday, prec_2)%max, &
  sdata(iday, prec1_1)%max,  sdata(iday, prec1_2)%max, &
  sdata(iday, prec2_1)%max, sdata(iday, prec2_2)%max, &
  sdata(iday, prec24_1)%max, sdata(iday, prec24_2)%max, &
  sdata(iday, snow)%min, sdata(iday, snow)%ave, sdata(iday, snow)%max, &
  sdata(iday, tg10)%min, sdata(iday, tg10)%ave, sdata(iday,tg10)%max, &
  sdata(iday, tg20)%min, sdata(iday, tg20)%ave, sdata(iday,tg20)%max, &
  sdata(iday, sm10)%min, sdata(iday, sm10)%ave, sdata(iday, sm10)%max, &
  sdata(iday, sm20)%min, sdata(iday, sm20)%ave, sdata(iday, sm20)%max, &
  sdata(iday, qsu)%ave, sdata(iday, qsu)%max, &
  sdata(iday, albedo)%ave
end do


close(11)



!  for HTML
write(out, '(A7,A11)') direc,'/index.html'

open(unit=11, file=out(1:len_trim(out)))
write(11,*) '<html>'
write(11,*) '<title>', year, month, 'statistics for TAKAYAMA AWS','</title>'
write(11,*) '<body>'

write(11,*) '<h1> TAKAYAMA AWS', year,'/',month,'</h1>'
write(11,*) 'River Basin Research Center, Gifu University'
write(11,*) '<a href="http://sateco-archive.green.gifu-u.ac.jp/AWS1/"> Takayama AWS WWW Page </a><br>',&
  ' Please check the error inforamtion on this page before using the following data.'
write(11,*) '<p>'
write(11,*) '<h2> Monthly Statistics ', year, month,'</h2>'
write(11,*) '<a href="',outm(1:len_trim(outm)),'"> CSV File </a><br>'
write(11,*) '<table border="1">'
write(11, '(a)') '<tr><td>year</td><td>month</td><td>&
  &wind_spped<br>(min)</td><td>wind_speed<br>(ave)</td><td>wind_speed<br>(max)</td><td>&
  &wind_speed<br>(vec_ave)</td><td>wind_direc<br>(vec_ave)</td><td>&
  &Air_temp<br>(min)</td><td>Air_temp<br>(ave)</td><td>Air_temp<br>(Max)</td><td>&
  &Rel_Hum<br>(min)</td><td>Rel_Hum<br>(ave)</td><td>Rel_Hum<br>(max)</td><td>&
  &Spec_Hum<br>(min)</td><td>Spec_Hum<br>(ave)</td><td>Spec_Hum<br>(max)</td><td>&
  &Vap_Pres<br>(min)</td><td>Vap_Pres<br>(ave)</td><td>Vap_Pres<br>(max)</td><td>&
  &Solor_Rad<br>(ave)</td><td>Solor_Rad<br>(max)</td><td>&
  &Precise_Solor_Rad<br>(ave)</td><td>&
        &Precise_Solor_Rad<br>(max)</td><td>&
  &Longwave_Rad<br>(min)</td><td>Longwave_Rad<br>(ave)</td><td>Longwave_Rad<br>(max)</td><td>&
  &Press<br>(min)</td><td>Press<br>(ave)</td><td>Press<br>(max)</td><td>&
  &Rainfall<br>(sum)</td><td>Rainfall<br>(max/10min.)</td><td>&
  &Precipitation<br>(sum)</td><td>Precipitation<br>(max/10min.)</td><td>&
  &Rainfall<br>(max/1hour)</td><td>Precipitation<br>(max/1hour)</td><td>&
  &Rainfall<br>(max/2hour)</td><td>Precipitation<br>(max/2hour)</td><td>&
  &Rainfall<br>(max/24hour)</td><td>Precipitation<br>(max/24hour)</td><td>&
  &Snow depth<br>(min)</td><td>Snow depth<br>(ave)</td><td>Snow depth<br>(max)</td>&
  &<td>Ground_Temp_10cm<br>(min)</td><td>Ground_Temp_10cm<br>(ave)</td>&
    &<td><td>Ground_Temp_10cm<br>(max)</td>&
  &<td>Ground_Temp_20cm<br>(min)</td><td>Ground_Temp_20cm<br>(ave)</td>&
    &<td><td>Ground_Temp_20cm<br>(max)</td>&
  &<td>WaterContent_10cm<br>(min)</td><td>WaterContent_10cm<br>(ave)</td>&
    &<td>WaterContent_10cm<br>(max)</td>&
  &<td>WaterContent_20cm<br>(min)</td><td>WaterContent_20cm<br>(ave)</td>&
    &<td>WaterContent_20cm<br>(max)</td>&
  &<td>Upward SW radiation<br>(ave)</td><td>Upward SW radiation(max)</td>&
  &<td>Daily albedo</td></tr>'




write(11,'(a)') '<tr><td></td><td>(1-12)</td>&
    &<td>(m/s)</td><td>(m/s)</td><td>(m/s)</td><td>&
    &(m/s)</td><td>(deg)</td><td>&
    &(C deg)</td><td>(C deg)</td><td>(C deg)</td><td>&
    &(%)</td><td>(%)</td><td>(%)</td><td>&
    &(g/kg)</td><td>(g/kg)</td><td>(g/kg)</td><td>&
    &(hPa)</td><td>(hPa)</td><td>(hPa)</td><td>&
    &(W/m2)</td><td>(W/m2)</td><td>&
    &(W/m2)</td><td>(W/m2)</td><td>&
    &(W/m2)</td><td>(W/m2)</td><td>(W/m2)</td><td>&
    &(hPa)</td><td>(hPa)</td><td>(hPa)</td><td>&
    &(mm/month)</td><td>(mm/10min)</td><td>&
    &(mm/month)</td><td>(mm/10min)</td><td>&
    &(mm/1hour)</td><td>(mm/1hour)</td><td>&
    &(mm/2hour)</td><td>(mm/2hour)</td><td>&
    &(mm/24hour)</td><td>(mm/24hour)</td>&
    &(cm)</td><td>(cm)</td><td>(cm)</td>&
    &<td>(C deg)</td><td>(C deg)</td><td>(C deg)</td>&
    &<td>(C deg)</td><td>(C deg)</td><td>(C deg)</td>&
    &<td>(%)</td><td>(%)</td><td>(%)</td>&
    &<td>(%)</td><td>(%)</td><td>(%)</td>&
    &<td>(W/m2)</td><td>(W/m2)</td>&
    &<td>(-)</td></tr>'






write(11, '("<tr><td>",I4,"</td><td>",I2,54("</td><td>",F10.1)),"</td><td>",F10.2,"</td></tr>"')&
  year, month,&
  msdata(1, ws)%min, msdata(1, ws)%ave, msdata(1, ws)%max,&
  msdata(1, ws_v)%ave, msdata(1, wdir_v)%ave,&
  msdata(1, ta)%min, msdata(1, ta)%ave, msdata(1, ta)%max, &
  msdata(1, rh)%min, msdata(1, rh)%ave, msdata(1, rh)%max, &
  1000*msdata(1, sphum)%min, 1000*msdata(1, sphum)%ave, &
                                  1000*msdata(1, sphum)%max, &
  msdata(1, vappress)%min, msdata(1, vappress)%ave, &
          msdata(1, vappress)%max, &
        msdata(1, qs)%ave, msdata(1,qs)%max, &
        msdata(1, qs2)%ave, msdata(1,qs2)%max, &
  msdata(1,ir)%min, msdata(1,ir)%ave, msdata(1,ir)%max, &
  msdata(1,press)%min, msdata(1,press)%ave, msdata(1,press)%max, &
  msdata(1, prec_1)%sum, msdata(1, prec_1)%max, &
  msdata(1, prec_2)%sum, msdata(1, prec_2)%max, &
  msdata(1, prec1_1)%max,  msdata(1, prec1_2)%max, &
  msdata(1, prec2_1)%max, msdata(1, prec2_2)%max, &
  msdata(1, prec24_1)%max, msdata(1, prec24_2)%max, &
  msdata(1, snow)%min, msdata(1, snow)%ave, msdata(1, snow)%max, &
  msdata(1, tg10)%min, msdata(1, tg10)%ave, msdata(1,tg10)%max, &
  msdata(1, tg20)%min, msdata(1, tg20)%ave, msdata(1,tg20)%max, &
  msdata(1, sm10)%min, msdata(1, sm10)%ave, msdata(1, sm10)%max, &
  msdata(1, sm20)%min, msdata(1, sm20)%ave, msdata(1, sm20)%max, &
  msdata(1, qsu)%ave, msdata(1, qsu)%max, &
  msdata(1, albedo)%ave

write(11,*) '</table>'
write(11,*) '</p>'


! output daily statistics
write(11,*) '<h2> Daily Statistics </h2>'
write(11,*) '<p>'
write(11,*) '<a href="',outd(1:len_trim(outm)),'"> CSV File </a><br>'
write(11,*) '<table border="1">'

write(11,'(a)') '<tr><td>year</td><td> month</td><td>day</td><td>&
  &wind_spped<br>(min)</td><td>wind_speed<br>(ave)</td><td>wind_speed<br>(max)</td><td>&
  &wind_speed<br>(vec_ave)</td><td>wind direc<br>(vec_ave)</td><td>&
  &Air_temp<br>(min)</td><td>Air_temp<br>(ave)</td><td>Air_temp<br>(Max)</td><td>&
  &Rel_Hum<br>(min)</td><td>Rel_Hum<br>(ave)</td><td>Rel_Hum<br>(max)</td><td>&
  &Spec_Hum<br>(min)</td><td>Spec_Hum<br>(ave)</td><td>Spec_Hum<br>(max)</td><td>&
  &Vap_Pres<br>(min)</td><td>Vap_Pres<br>(ave)</td><td>Vap_Pres<br>(max)</td><td>&
  &Solor_Rad<br>(ave)</td><td>Solor_Rad<br>(max)</td><td>&
  &Precise_Solor_Rad<br>(ave)</td><td>&
        &Precise_Solor_Rad<br>(max)</td><td>&
  &Longwave_Rad<br>(min)</td><td>Longwave_Rad<br>(ave)</td><td>Longwave_Rad<br>(max)</td><td>&
  &Press<br>(min)</td><td>Press<br>(ave)</td><td>Press<br>(max)</td><td>&
  &Rainfall<br>(sum)</td><td>Rainfall<br>(max/10min.)</td><td>&
  &Precipitation<br>(sum)</td><td>Precipitation<br>(max/10min.)</td><td>&
  &Rainfall<br>(max/1hour)</td><td>Precipitation<br>(max/1hour)</td><td>&
  &Rainfall<br>(max/2hour)</td><td>Precipitation<br>(max/2hour)</td><td>&
  &Rainfall<br>(max/24hour)</td><td>Precipitation<br>(max/24hour)</td><td>&
  &Snow depth<br>(min)</td><td>Snow depth<br>(ave)</td><td>Snow depth<br>(max)</td>&
  &<td>Ground_Temp_10cm<br>(min)</td><td>Ground_Temp_10cm<br>(ave)</td><td><td>Ground_Temp_10cm<br>(max)</td>&
  &<td>Ground_Temp_20cm<br>(min)</td><td>Ground_Temp_20cm<br>(ave)</td><td><td>Ground_Temp_20cm<br>(max)</td>&
  &<td>WaterContent_10cm<br>(min)</td><td>WaterContent_10cm<br>(ave)</td><td>WaterContent_10cm<br>(max)</td>&
  &<td>WaterContent_20cm<br>(min)</td><td>WaterContent_20cm<br>(ave)</td><td>WaterContent_20cm<br>(max)</td>&
  &<td>Upward SW radiation<br>(ave)</td><td>Upward SW radiation(max)</td>&
  &<td>Daily albedo</td></tr>'



write(11,'(a)') '<tr><td></td><td>(1-12)</td><td></td><td>(m/s) </td><td>(m/s)</td><td>(m/s)</td><td>&
    &(m/s)</td><td>(deg)</td><td>&
    &(C deg)</td><td>(C deg)</td><td>(C deg)</td><td>&
    &(%)</td><td>(%)</td><td>(%)</td><td>&
    &(g/kg)</td><td>(g/kg)</td><td>(g/kg)</td><td>&
    &(hPa)</td><td>(hPa)</td><td>(hPa)</td><td>&
    &(W/m2)</td><td>(W/m2)</td><td>&
    &(W/m2)</td><td>(W/m2)</td><td>&
    &(W/m2)</td><td>(W/m2)</td><td>(W/m2)</td><td>&
    &(hPa)</td><td>(hPa)</td><td>(hPa)</td><td>&
    &(mm/day)</td><td>(mm/10min)</td><td>&
    &(mm/day)</td><td>(mm/10min)</td><td>&
    &(mm/1hour)</td><td>(mm/1hour)</td><td>&
    &(mm/2hour)</td><td>(mm/2hour)</td><td>&
    &(mm/24hour)</td><td>(mm/24hour)</td>&
    &(cm)</td><td>(cm)</td><td>(cm)</td>&
    &<td>(C deg)</td><td>(C deg)</td><td>(C deg)</td>&
    &<td>(C deg)</td><td>(C deg)</td><td>(C deg)</td>&
    &<td>(%)</td><td>(%)</td><td>(%)</td>&
    &<td>(%)</td><td>(%)</td><td>(%)</td>&
    &<td>(W/m2)</td><td>(W/m2)</td>&
    &<td>(-)</td></tr>'


do iday=1, days


write(11, '("<tr><td>",I4,2("</td><td>",I2),54("</td><td>",F10.1)),"</td><td>",F10.2,"</td></tr>"')&
  year, month, iday, &
  sdata(iday, ws)%min, sdata(iday, ws)%ave, sdata(iday, ws)%max, &
  sdata(iday, ws_v)%ave, sdata(iday, wdir_v)%ave,              &
  sdata(iday, ta)%min, sdata(iday, ta)%ave, sdata(iday, ta)%max, &
  sdata(iday, rh)%min, sdata(iday, rh)%ave, sdata(iday, rh)%max, &
  1000*sdata(iday, sphum)%min, 1000*sdata(iday, sphum)%ave, &
                                  1000*sdata(iday, sphum)%max, &
  sdata(iday, vappress)%min, sdata(iday, vappress)%ave, &
          sdata(iday, vappress)%max, &
        sdata(iday, qs)%ave, sdata(iday,qs)%max, &
        sdata(iday, qs2)%ave, sdata(iday,qs2)%max, &
  sdata(iday,ir)%min, sdata(iday,ir)%ave, sdata(iday,ir)%max, &
  sdata(iday, press)%min, sdata(iday, press)%ave,sdata(iday, press)%max,&
  sdata(iday, prec_1)%sum, sdata(iday, prec_1)%max, &
  sdata(iday, prec_2)%sum, sdata(iday, prec_2)%max, &
  sdata(iday, prec1_1)%max,  sdata(iday, prec1_2)%max, &
  sdata(iday, prec2_1)%max, sdata(iday, prec2_2)%max, &
  sdata(iday, prec24_1)%max, sdata(iday, prec24_2)%max, &
  sdata(iday, snow)%min, sdata(iday, snow)%ave, sdata(iday, snow)%max, &
  sdata(iday, tg10)%min, sdata(iday, tg10)%ave, sdata(iday,tg10)%max, &
  sdata(iday, tg20)%min, sdata(iday, tg20)%ave, sdata(iday,tg20)%max, &
  sdata(iday, sm10)%min, sdata(iday, sm10)%ave, sdata(iday, sm10)%max, &
  sdata(iday, sm20)%min, sdata(iday, sm20)%ave, sdata(iday, sm20)%max, &
  sdata(iday, qsu)%ave, sdata(iday, qsu)%max, &
  sdata(iday, albedo)%ave

end do

write(11,*) '</table>'
write(11,*) '</p>'

write(11,*) '<h2> daily plots </h2>'
write(out, '(A7,A1,I0.4A1I0.2A)') direc,'/',year,'_', month, '_daily_stat.csv'
call plotdaily( year, month, out )


! histgrams 
write(11,*) '<p>'
write(11,*) '<h2>Histgrams for several variables</h2>'

! output histgrams & windrose
write(out1, '(A7,A1,I0.4A1I0.2)') direc,'/',year,'_', month
write(out2, '(I0.4A1I0.2)') year,'_', month

 call windroseout( windrose, out1 )
 write(11,*) '<br>'



write(11,*) '<table>'
write(11,*) '<tr>'
 call histout( hist(ta), out1, '_Ta_hist.csv', 'Air Temperature')
 call histout( hist(rh), out1, '_RH_hist.csv', 'Relative Humidity')
 call histout( hist(sphum), out1, '_SpHum_hist.csv', 'Specific Humidity')
write(11,*) '</tr><tr>'
 call histout( hist(vappress), out1, '_VapPress_hist.csv', 'Vapor Pressure')
 call histout( hist(qs), out1, '_Qs_hist.csv', 'Solar Radiation(simple)')
 call histout( hist(qs2), out1, '_PreciseQs_hist.csv', 'Solar Radiation(precise)')
write(11,*) '</tr><tr>'
 call histout( hist(ws), out1, '_Ws_hist.csv', 'Wind Speed')
 call histout( hist(ir), out1, '_IR_hist.csv', 'Downward Long Wave Radiation')
 call histout( hist(prec_1), out1, '_Rainfall_hist.csv', 'Rainfall/10min.')
write(11,*) '</tr><tr>'
 call histout( hist(prec1_1), out1, '_Rainfall_1hr_hist.csv', 'Rainfall/1hor')
 call histout( hist(prec2_1), out1, '_Rainfall_2hr_hist.csv', 'Rainfall/2hour')
 call histout( hist(prec24_1), out1, '_Rainfall_24hr_hist.csv','Rainfall/24hour' )
write(11,*) '</tr><tr>'
 call histout( hist(prec_2), out1, '_Precipitation_hist.csv', 'Precipitation/10min.')
 call histout( hist(prec1_2), out1, '_Precipitation_1hr_hist.csv', 'Precipitation/1hour')
 call histout( hist(prec2_2), out1, '_Precipitation_2hr_hist.csv', 'Precipitation/2hour')
write(11,*) '</tr><tr>'
 call histout( hist(prec24_2), out1, '_Precipitation_24hr_hist.csv', 'Precipitation/24hour')
 call histout( hist(press), out1, '_Pressure_hist.csv', 'Pressure')
 call histout( hist(snow), out1, '_SnowDepth_hist.csv', 'Snow Depth')
write(11,*) '</tr><tr>'
 call histout( hist(tg10), out1, '_TG10_hist.csv', 'Ground Temp 10cm')
 call histout( hist(tg20), out1, '_TG20_hist.csv', 'Ground Temp 20cm')
 call histout( hist(sm10), out1, '_SM10_hist.csv', 'Soil Water Content 10cm')
write(11,*) '</tr><tr>'
 call histout( hist(sm20), out1, '_SM20_hist.csv', 'Soil Water Content 20cm')
 call histout( hist(qsu), out1, '_Qsu_hist.csv', 'Upward SW Radiation')
write(11, *) '</tr>'
write(11,*) '</table>'
write(11,*) '</p>'

write(11,*) '<a href="../../"> Up to Monthly Statistics </a>,'
write(11,*) '<a href="http://sateco-archive.green.gifu-u.ac.jp/AWS1/"> Up to TAKAYAMA AWS </a>'
write(11,*) '</body>'
write(11,*) '</html>'

close(11)

end


subroutine plotdaily( year, month, out )
use typedef
use mod_nodata
use fileinfo
implicit none
character(len=*) ::  out
integer :: year, month

character(len=1024) :: direc, body, ext
character(len=1024) :: gfile, jfile, scrfile, command

call resolvefile( out, direc, body, ext,  gfile, jfile, scrfile )

! write(command,*) 'tail +2 ', trim(out), ' | sed s/-9999.0/"-"/g | sed s/,//g > ', &
 write(command,*) 'tail -n +3 ', trim(out), ' | sed s/-9999.0/"-"/g | sed s/,//g > ', &
    trim(direc)//'/tmp1'
call system(command)

scrfile = trim(body)//'-1.gpl'
gfile   = trim(body)//'-1.eps'
jfile   = trim(body)//'-1.jpg'

open(unit=12, file=direc(1:len_trim(direc))//'/'//trim(scrfile))
write(12,*) 'set term post eps noenhanced color "Helvetica" 12'
write(12,*) 'set output "',trim(gfile), '"' 
write(12,*) 'set multiplot title "', year, month, '"'
write(12,*) 'set size 1,0.2'
write(12,*) 'set origin 0, 0.75'
write(12,*) 'set datafile miss "-"'
write(12,*) 'set ylabel "m/s"'
write(12,*) 'plot "tmp1" using 3:5:4:6 title "wind speed(max,ave,min)"  w  yerrorbars,&
  &"tmp1" using 3:7 title "vector ave"  w l'
write(12,*) 'set origin 0, 0.50'
write(12,*) 'set ylabel "deg."'
write(12,*) 'plot "tmp1" using 3:8 title "Wind Dir"'
write(12,*) 'set origin 0,0.25'
write(12,*) 'set ylabel "C.deg"'
write(12,*) 'plot "tmp1" using 3:10:9:11 title "Temp.(max,ave.min)" &
  &w yerrorbars,"tmp1" using 3:45 title "Tg10" with line, "tmp1" using 3:48&
  & title "Tg20" with line'
write(12,*) 'set origin 0,0.0'
write(12,*) 'set ylabel "%"'
write(12,*) 'plot "tmp1" using 3:13:12:14  title "RH(max,ave.min)" w yerrorbars'
write(12,*) 'unset multiplot'
close(12)

write(command, *) '( cd ', direc(1:len_trim(direc)),'; gnuplot ',&
  scrfile(1:len_trim(scrfile)), '; convert -density 120 ', trim(gfile), &
  '  ', trim(jfile), ')'
call system( command )

write(11,*) '<a href="',trim(gfile),'"><img src="', &
    trim(jfile),'"></a>'

scrfile = trim(body)//'-2.gpl'
gfile   = trim(body)//'-2.eps'
jfile   = trim(body)//'-2.jpg'

open(unit=12, file=trim(direc)//'/'//trim(scrfile))
write(12,*) 'set term post eps noenhanced color "Helvetica" 12'
write(12,*) 'set output "',trim(gfile), '"' 
write(12,*) 'set multiplot title "', year, month, '"'
write(12,*) 'set size 1,0.2'
write(12,*) 'set origin 0, 0.75'
write(12,*) 'set datafile miss "-"'
write(12,*) 'set ylabel "g/kg or hPa"'
write(12,*) 'plot "tmp1" using 3:16 title "Spec. Hum" w l,  "tmp1" using 3:19 title "Vapor Pres" w l'
write(12,*) 'set origin 0, 0.50'
write(12,*) 'set ylabel "W/m^2"'
write(12,*) 'plot "tmp1" using 3:21 title "Qs" w l,&
  &"tmp1" using 3:23 title "Qs(precise)" w l, "tmp1" using 3:26 title "Ql" w l,&
  & "tmp1" using 3:56 title "Qs(up)" w l'
write(12,*) 'set origin 0,0.25'
write(12,*) 'set ylabel "hPa"'
write(12,*) 'plot "tmp1" using 3:29 title "P" w l'
write(12,*) 'set origin 0,0.0'
write(12,*) 'set ylabel "mm/day"'
write(12,*) 'plot "tmp1" using 3:31 title "Rainfall" w impulses, "tmp1" using 3:33 title "Preciptation" w impulses'
write(12,*) 'unset multiplot'
close(12)

write(command, *) '( cd ', trim(direc),'; gnuplot ',&
      trim(scrfile), '; convert -density 120 ', trim(gfile), &
          '  ', trim(jfile), ')'
call system( command )

write(11,*) '<a href="',trim(gfile),'"><img src="', &
    trim(jfile),'"></a>'


scrfile = trim(body)//'-3.gpl'
gfile   = trim(body)//'-3.eps'
jfile   = trim(body)//'-3.jpg'

open(unit=12, file=trim(direc)//'/'//trim(scrfile))
write(12,*) 'set term post eps noenhanced color "Helvetica" 12'
write(12,*) 'set output "',trim(gfile), '"' 
write(12,*) 'set multiplot title "', year, month, '"'
write(12,*) 'set size 1,0.2'
write(12,*) 'set origin 0, 0.75'
write(12,*) 'set datafile miss "-"'
write(12,*) 'set ylabel "cm"'
write(12,*) 'plot "tmp1" using 3:42:41:43 title "Snow Depth" w  yerrorbars'



write(12,*) 'set origin 0, 0.5'
write(12,*) 'set ylabel "%"'
write(12,*) 'plot "tmp1" using 3:51 title "WaterContent10" w l, &
  &"tmp1" using 3:53 title "WaterContent20" w l'

write(12,*) 'set origin 0, 0.25'
write(12,*) 'set ylabel ""'
write(12,*) 'plot "tmp1" using 3:58 title "Albedo" w l'

write(12,*) 'unset multiplot'

close(12)

write(command, *) '( cd ', trim(direc),'; gnuplot ',&
      trim(scrfile), '; convert -density 120 ', trim(gfile), &
                '  ', trim(jfile), ')'
call system( command )

write(11,*) '<a href="',trim(gfile),'"><img src="', &
                trim(jfile),'"></a>'


end


subroutine windroseout( windrose, out1 )
use typedef
use mod_nodata
implicit none
type(type_windrose)    :: windrose
character(len=*) :: out1
character(len=1024) :: out, direc, body, ext
character(len=1024) :: gfile, jfile, scrfile, command
integer :: i
real    :: delta

out=trim(out1)//'_windrose.csv'
write(0,*) 'output file name=', trim(out)
open(unit=12, file=trim(out))
delta = 360.0 / windrose%n 

do i=1, windrose%n
write(12,*) delta*(i-1.5),',',delta*(i-1),',', delta*(i-0.5),',',windrose%sum_ws(i)
end do
i=1
write(12,*) delta*(i-1.5),',',delta*(i-1),',', delta*(i-0.5),',',windrose%sum_ws(i)
close(12)


write(11,*) '<p>'
write(11,*) '<h2> Windrose </h2>'

call resolvefile( out, direc, body, ext,  gfile, jfile, scrfile )

! make gnuplot output
!  make tmporal data file : tmp
open(unit=12, file=trim(direc)//'/tmp')
write(12,*) 
do i=1, windrose%n
write(12,*) delta*(i-1),windrose%sum_ws(i)
end do
i=1
write(12,*) delta*(i-1),windrose%sum_ws(i)
close(12)

open(unit=12, file=trim(direc)//'/'//trim(scrfile))
write(12,*) 'set term post eps noenhanced color "Helvetica" 18'
write(12,*) 'set output "',trim(gfile),'"' 
write(12,*) 'set title "', 'Windrose', '(', trim(body),')"'
write(12,*) 'set polar'
write(12,*) 'set grid polar'
write(12,*) 'set size square'
write(12,*) 'set angles degree'
write(12,*) 'plot "tmp" notitle with linespoints'
close(12)

write(command, *) '( cd ', trim(direc),'; gnuplot ',&
  trim(scrfile), '; convert ', trim(gfile), &
  '  ', trim(jfile), ')'
call system( command )


write(11,*) '<a href="',trim(gfile),'"><img src="', &
    trim(jfile),'"></a>, <a href="', &
    trim(body)//trim(ext),&
    '"> data file (csv) </a> <wbr>'





end

subroutine histout( hist, out1, out2, titlename )
use typedef
use mod_nodata
implicit none

type(histdata) :: hist
character(len=*) :: out1, out2, titlename
character(len=1024) :: out, direc, body, ext
character(len=1024) :: gfile, jfile, scrfile, command
integer :: i

out=trim(out1)//trim(out2)
write(0,*) 'output file name=', trim(out)
open(unit=12, file=trim(out))

write(12,*) NODATA, ',',hist%xmin-0.5*hist%dx,',',hist%xmin,',', hist%hist(1)
do i=2, hist%hist_n-1
  write(12,*) hist%xmin+(i-2)*hist%dx,',', hist%xmin+(i-1.5)*hist%dx, &
    ',', hist%xmin+(i-1)*hist%dx, ',',hist%hist(i)
end do
write(12,*)  hist%xmin+(hist%hist_n - 2)*hist%dx,',', &
    hist%xmin+(hist%hist_n -1.5)*hist%dx,',', & 
    NODATA, ',',hist%hist(hist%hist_n)

close(12)

! make graph with gnulot

call resolvefile( out, direc, body, ext,  gfile, jfile, scrfile )

! make gnuplot output
!  make tmporal data file : tmp
open(unit=12, file=trim(direc)//'/tmp')
write(12,*) hist%xmin-0.5*hist%dx,hist%hist(1)
do i=2,  hist%hist_n
write(12, *) hist%xmin+(i-1.5)*hist%dx, hist%hist(i)
end do
close(12)

open(unit=12, file=trim(direc)//'/'//trim(scrfile))
write(12,*) 'set term post eps noenhanced color "Helvetica" 18'
write(12,*) 'set output "',trim(gfile),'"' 
write(12,*) 'set title "', trim(titlename), '(', trim(body),')"'
write(12,*) 'plot "tmp" notitle with boxes'
close(12)

write(command, *) '( cd ', trim(direc),'; gnuplot ',&
  trim(scrfile), '; convert ', trim(gfile), &
  '  ', trim(jfile), ')'
call system( command )


write(11,*) '<td><a href="',trim(gfile),'"><img src="', &
    trim(jfile),'"></a>,<br> <a href="', &
    trim(body)//trim(ext),&
    '"> data file (csv) </a> </td>'



end 

subroutine resolvefile( file, direc, body, ext, gfile, jfile, scrfile )
implicit none
character(len=*) :: file, direc, body, ext, gfile, jfile, scrfile
integer :: i, i1, i2


i1 = 0
do i=1, len_trim(file)
 if( file(i:i)=='/' ) then
   i1 =i
 end if
end do

if( i1 > 1 ) then
 direc=file(1:i1-1)
else if (i1==1) then
 direc='/'
else
 direc=''
end if



body = file(i1+1:(len_trim(file)-4))
ext  = file((len_trim(file)-3):len_trim(file))


gfile =   trim(body)//'.eps'
jfile =   trim(body)//'.jpg'
scrfile = trim(body)//'.gpl'
end




subroutine  mkwindrose( data, total_num, &
  duration_start, duration_end, windrose)
use typedef
use mod_nodata
use fileinfo
implicit none

integer :: total_num
type(awsdata) :: data(total_num)
integer :: duration_start, duration_end
type(type_windrose)    :: windrose


real :: delta
integer :: i, k

delta = 360.0 / windrose%n 
windrose%sum_ws = 0

  do i=1, total_num
   if( data(i)%minute  > duration_start &
        .AND.  data(i)%minute <= duration_end &
  .AND.  (.NOT. is_nodata( data(i)%dat(ws)   ) ) &
  .AND.  (.NOT. is_nodata( data(i)%dat(wdir) ) )) then
  

     k = (data(i)%dat(wdir) - 0.5*delta)/delta + 2
     if( k > windrose_class ) then
        k = 1
     end if

     windrose%sum_ws(k) = windrose%sum_ws(k)+data(i)%dat(ws)*dt/60.0 
                   ! m/s * min / 60
    end if

   end do
end








subroutine mkhist(  data, total_num, duration_start, duration_end, hist )
use typedef
use fileinfo
use mod_nodata
implicit none

integer :: total_num
type(awsdata) :: data(total_num)
type(histdata) :: hist(maxch)
integer :: i, j, class
integer :: duration_start, duration_end

  do i=1, maxch
   hist(i)%dx = 0;  hist(i)%xmin = 0; hist(i)%xmax = 0
   hist(i)%hist_n = 0 ; hist(i)%hist = 0;
  end do

  hist(ta)%dx = 1; hist(ta)%xmin = -20; hist(ta)%xmax = 30
  hist(ws)%dx = 0.25; hist(ws)%xmin = 0; hist(ws)%xmax = 10
  hist(rh)%dx = 5; hist(rh)%xmin = 0; hist(rh)%xmax = 100
  hist(vappress)%dx = 1; hist(vappress)%xmin = 0; hist(vappress)%xmax = 30
  hist(sphum)%dx = 0.001; hist(sphum)%xmin = 0; hist(sphum)%xmax = 0.03
  hist(qs)%dx = 20; hist(qs)%xmin = 0; hist(qs)%xmax = 1000
  hist(qs2)%dx = 20; hist(qs2)%xmin = 0; hist(qs2)%xmax = 1000
  hist(ir)%dx = 20; hist(ir)%xmin = 0; hist(ir)%xmax = 1000
  hist(prec_1)%dx = 0.5; hist(prec_1)%xmin = -0.5; hist(prec_1)%xmax = 10
  hist(prec_2)%dx = 0.5; hist(prec_2)%xmin = -0.5; hist(prec_2)%xmax = 10
  hist(prec1_1)%dx = 0.5; hist(prec1_1)%xmin = -0.5; hist(prec1_1)%xmax = 10
  hist(prec1_2)%dx = 0.5; hist(prec1_2)%xmin = -0.5; hist(prec1_2)%xmax = 10
  hist(prec2_1)%dx = 0.5; hist(prec2_1)%xmin = -0.5; hist(prec2_1)%xmax = 10
  hist(prec2_2)%dx = 0.5; hist(prec2_2)%xmin = -0.5; hist(prec2_2)%xmax = 10
  hist(prec24_1)%dx = 0.5; hist(prec24_1)%xmin = -0.5; hist(prec24_1)%xmax = 10
  hist(prec24_2)%dx = 0.5; hist(prec24_2)%xmin = -0.5; hist(prec24_2)%xmax = 10
  hist(press)%dx = 5; hist(press)%xmin   = 830; hist(press)%xmax = 880

  hist(snow)%dx = 5; hist(snow)%xmin = 0; hist(snow)%xmax = 150


 do j=1, maxch
  if( hist(j)%dx > 0 ) then
  hist(j)%hist_n = (hist(j)%xmax - hist(j)%xmin)/hist(j)%dx + 2

   if( hist_nmax < hist(j)%hist_n ) then
         hist(j)%hist_n = hist_nmax
   write(0,*) 'numbers of histgram data for', j, 'ch is not sufficient.'
   end if


  hist(j)%xmax = hist(j)%xmin + hist(j)%dx * ( hist(j)%hist_n - 2 )

  do i=1, total_num
   if( data(i)%minute  > duration_start &
        .AND.  data(i)%minute <= duration_end &
  .AND.  (.NOT. is_nodata( data(i)%dat(j) ) )) then

      class = (data(i)%dat(j) - hist(j)%xmin)/hist(j)%dx + 2
      if( class < 1 ) then
         class = 1
      else if( class > hist(j)%hist_n ) then
         class = hist(j)%hist_n
      end if

      hist(j)%hist(class) = hist(j)%hist(class) + 1

  end if
 end do
 end if
 end do 

end



subroutine mkstat( data, total_num, sdata, sdata_num, duration_min, sduration )
use typedef
use mod_nodata
use fileinfo
implicit none

integer ::  total_num, sdata_num, duration_min,  sduration
type(awsdata) :: data(total_num)
type(stat) :: sdata(sdata_num,maxch)

integer :: i, j, k, m1, m2
!logical :: is_nodata


 do i=1, sdata_num
  do j=1, maxch
   sdata(i, j)%sum = 0 ; sdata(i, j)%num = 0
  end do
 end do


 do i=1, sdata_num

  m1 = (i-1)*sduration; m2 = i*sduration


    do k=1, total_num
    if( data(k)%minute > m1 .AND. data(k)%minute <= m2 ) then

      do j=1, maxch

     if( .NOT. is_nodata(data(k)%dat(j)) ) then
       if(sdata(i, j)%num == 0) then
         sdata(i, j)%max = data(k)%dat(j)
         sdata(i, j)%min = data(k)%dat(j)
       else
         if(sdata(i, j)%max < data(k)%dat(j)) then
     sdata(i, j)%max = data(k)%dat(j)
   end if
         if(sdata(i, j)%min > data(k)%dat(j)) then
     sdata(i, j)%min = data(k)%dat(j)
   end if
       end if
       sdata(i, j)%num = sdata(i, j)%num +1
       sdata(i, j)%sum = sdata(i, j)%sum + data(k)%dat(j)
       end if
     end do
    end if

    end do

  do j=1, maxch
   if(sdata(i,j)%num == 0 ) then
     sdata(i,j)%sum = NODATA
     sdata(i,j)%max = NODATA
     sdata(i,j)%min = NODATA
     sdata(i,j)%ave = NODATA
   else
    sdata(i,j)%ave = sdata(i,j)%sum/sdata(i,j)%num
  end if
  end do

  call uv2wswdir( sdata(i, u)%ave, sdata(i,v)%ave, &
      sdata(i, ws_v)%ave, sdata(i, wdir_v)%ave )


 end do

end


subroutine daily_albedo( sdata, sdata_num )
use typedef
use mod_nodata
use fileinfo
implicit none

integer ::  sdata_num
type(stat) :: sdata(sdata_num,maxch)

integer :: i
!logical :: is_nodata


 do i=1, sdata_num

  if( (.NOT. is_nodata(sdata(i, qs2)%ave)) &
  .AND. (.NOT. is_nodata(sdata(i, qsu)%ave )) ) then
  sdata(i, albedo)%ave = sdata(i, qsu)%ave/sdata(i, qs2)%ave

   else
    sdata(i,albedo)%ave = NODATA

  end if

 end do

end

subroutine mkcalcdata( data, total_num )
use typedef
use fileinfo
use mod_nodata
implicit none
integer :: total_num
type(awsdata) :: data(total_num)

integer :: i, err
type(type_aveparam)  :: aveparam

! make (u,v)  & humidity related values
do i=1, total_num
 call wswdir2uv( data(i)%dat(ws), data(i)%dat(wdir), &
    data(i)%dat(u), data(i)%dat(v) )
 call vaporcalc( data(i)%dat(ta), data(i)%dat(rh), data(i)%dat(press), &
    data(i)%dat(sphum), data(i)%dat(vappress) )
end do

! make 1hour 2hour 24hour running averaged precipitation

  call make_aveparam( 60, dt, aveparam, err ) 
  call running_sum( data, total_num, aveparam, prec_1, prec1_1 )
  call running_sum( data, total_num, aveparam, prec_2, prec1_2 )

  call make_aveparam( 120, dt, aveparam, err ) 
  call running_sum( data, total_num, aveparam, prec_1, prec2_1 )
  call running_sum( data, total_num, aveparam, prec_2, prec2_2 )

  call make_aveparam( 1440, dt, aveparam, err ) 
  call running_sum( data, total_num, aveparam, prec_1, prec24_1 )
  call running_sum( data, total_num, aveparam, prec_2, prec24_2 )

end



subroutine running_sum( data, total_num, aveparam, prec_1, prec1_1 )
use typedef
use mod_nodata

type(awsdata) :: data(*)
integer :: total_num
type(type_aveparam) :: aveparam
integer :: prec_1, prec1_1

! logical :: is_nodata
integer :: i1, i2

double precision :: sum
integer          :: num


  do i1=(aveparam%prev + 1), total_num - aveparam%foll
   sum = 0; num=0
   do i2 = i1-aveparam%prev, i1+aveparam%foll
    if( .NOT. is_nodata(data(i2)%dat(prec_1)) ) then
      sum = sum + data(i2)%dat(prec_1)
      num = num +1
    end if
   end do

   if( num == (aveparam%prev+aveparam%foll+1) ) then
      data(i1)%dat(prec1_1) = sum
   else
      data(i1)%dat(prec1_1) = NODATA
   end if

  end do


end




subroutine make_aveparam( ave_duration, dt, aveparam, err )
use typedef
implicit none
type(type_aveparam)  :: aveparam
integer :: ave_duration, dt, err


integer :: j

   j = ave_duration/dt
  if( mod( j, 2 ) == 1 ) then
    err = 1
   write(0,*) 'Average duration :', ave_duration, &
    'min is wrong for ', dt, 'min interval data'
  else
   err = 0
  end if

   j = j/2
   aveparam%prev = j-1
   aveparam%foll = j
end


subroutine vaporcalc( ta, rh, p, sphum, vappress )
use mod_nodata
implicit none
  real, intent(in) :: ta, rh, p
  real, intent(out):: sphum, vappress
! logical :: is_nodata
!
!  ta : temperature in C. deg
!  rh : relative humidity (%) 0-100
!  p  : pressure  (hPa)
! sphum : specific humidity  kg/kg
! vappress : vapor pressure (hPa)
!


 if( is_nodata(ta) .OR. is_nodata(rh) .OR. is_nodata(p) ) then
    sphum = NODATA
    vappress = NODATA
 else
  call vaporess_hPa( ta, rh, vappress )
  call specific_humidity( ta, vappress, p, sphum )
 end if

end

subroutine specific_humidity( ta, vappress, p, sphum )
 implicit none
 real :: ta, vappress, p, sphum
!
!  ta : temperature in C. deg
!  rh : relative humidity (%) 0-100
!  p  : pressure  (hPa)
! sphum : specific humidity  kg/kg
! vappress : vapor pressure (hPa)
!

  real, parameter :: Rd=287.0, Rw = 461.0
  real :: r

  r = vappress/(p-vappress) * (Rd/Rw)
 sphum = r/(1+r)
end

   

subroutine vaporess_hPa( T, rh, vappress )
implicit none
 real :: T, rh, vappress

!  ta : temperature in C. deg
!  rh : relative humidity (%) 0-100
! vappress : vapor pressure (hPa)

! USE saturation vapor pressure over water

  vappress = (0.01*rh)*  6.107799961 + T * (4.436518521e-1 &
                 + T * (1.428945805e-2 &
                 + T * (2.650648471e-4 &
                 + T * (3.031240396e-6 &
                 + T * (2.034080948e-8 &
                 + T * (6.136820929e-11))))))
!   from Paul.R.Lowe 1977, J. Appl. Met. 16, 100-103

end




subroutine wswdir2uv( ws, wdir, u, v )
!  ws: wind speed, wdir : wid dir. (0-360: 0 N, 90 E, ...)
! ws, wdir ==> u, v
!
use mod_nodata
real :: ws, wdir, u, v
real :: theta, u1, v1

 if( is_nodata( ws ) .OR. is_nodata( wdir ) ) then
    u = NODATA; v = NODATA
 else
                theta = (270-wdir)*3.14152965/180.0
                u1 = ws * cos(theta)
                v1 = ws * sin(theta)
                u = u1; v=v1
 end if
end subroutine

subroutine uv2wswdir( u, v, ws, wdir )
!  ws: wind speed, wdir : wid dir. (0-360: 0 N, 90 E, ...)
! ws, wdir <== u, v
!
use mod_nodata
real :: ws, wdir, u, v
real :: theta, ws1, wdir1

  if( is_nodata(u) .OR. is_nodata(v) ) then
    ws = NODATA; wdir = NODATA
  else
        ws1 = sqrt( u*u + v*v )
        theta = atan2( v, u )*180/3.14152965
        wdir1 = 270 - theta
        if( wdir1 < 0 ) then
         wdir1 = wdir1 + 360
        end if
        if( wdir1 >= 360) then
         wdir1 = wdir1 - 360
        end if

        ws = ws1; wdir =wdir1
  end if
end subroutine






subroutine testout( data, total_num, ch1, shift1, ch2, shift2 )
use typedef
implicit none
integer :: total_num, ch1, shift1, ch2, shift2
type(awsdata) :: data(total_num)

integer :: i, j

do i=1, total_num
 write(*,*) i, data(i)%date(1:len_trim(data(i)%date)), &
   data(i)%minute, (data(i)%dat(j), j=shift1+1, shift1+ch1), &
    (data(i)%dat(j), j=shift2+1, shift2+ch2)
! write(*,*) i, data(i)%date(1:len_trim(data(i)%date)), &
!   data(i)%minute, (data(i)%dat(j), j=1, maxch)
   
end do

end


subroutine initdata( data, total_num )
use typedef
use mod_nodata
implicit none
integer :: total_num
type(awsdata) :: data(total_num)

integer :: i, j

do i=1, total_num
  data(i)%dat = NODATA
end do


end



subroutine  read_copy( file1, ch1, shift1,  &
    start, duration_min, extent_min, dt, total_num, data )
use typedef
use mod_nodata
implicit none
character(len=*) :: file1
integer :: ch1, shift1
type(date_time) :: start, now
integer :: duration_min, extent_min, dt
integer :: total_num
type(awsdata)   :: data(total_num)

character(len=1024) :: buf, eachcol(maxch)
integer :: i, j, colnum, calc_diffmin, diffmin, ios

 write(0,*) 'Open ', trim(file1)
 open(unit=10, file=file1, status='old', err=999)
! skip 4 lines headder
 do i=1,4
   read(10, '(A)',  end=900, err=900) buf
 end do


 do while( .TRUE.  )
   read(10, '(A)',  end=900, err=900) buf
   call devidecsv( buf(1:len_trim(buf)), eachcol, maxch, colnum )

   read(eachcol(1), "(1x,I4,1x,I2,1x,I2,1x,I2,1x,I2,1x,I2)") &
         now%year, now%month, now%day, now%hour, now%min, now%sec 

   diffmin = calc_diffmin( now, start )

   if(diffmin + extent_min > 0 .AND. diffmin - extent_min <= duration_min) then
   !  COPY
      i = (diffmin+extent_min)/dt

      data(i)%minute = diffmin;
      data(i)%date = eachcol(1)(1:len_trim( eachcol(1) ) )
      do j=1, colnum-1
       read( eachcol(1+j), *, IOSTAT=ios) data(i)%dat(j+shift1) 
       if(ios /= 0) then
             data(i)%dat(j+shift1) = NODATA
       end if
      end do
      if( colnum < ch1 ) then
       do j=colnum+1, ch1
        data(i)%dat(j) = NODATA
       end do
      end if

    else if( diffmin - extent_min > duration_min ) then
     exit
    end if

 end do

900 close(10)

return

999 write(0,*) 'Error! Can not open ', file1
stop

end


integer function calc_diffmin( now, start )
!   now - start     in minute
use typedef
implicit none
type(date_time) :: start, now

 integer :: diffday

  call datedif( now%year, now%month, now%day, &
    start%year, start%month, start%day, diffday )

  calc_diffmin = (diffday*24 + (now%hour-start%hour))*60 &
      +(now%min - start%min)   !  + (now%sec - start%sec)/60
 
end






subroutine devidecsv( buf, eachcol, maxch, num )
implicit none
integer   :: maxch, num
character(len=*) :: buf, eachcol(maxch)
integer :: i, j, prev, now, now1, ll
logical :: quoted


ll = len_trim(buf)
quoted = .FALSE.
prev=0; now=1; i=1
do j=1, ll
 if(buf(j:j) == '"') then
    if( quoted ) then
      quoted = .FALSE.
    else
        quoted = .TRUE.
    end if
 end if

 if(buf(j:j) == ',' .AND. (.NOT.quoted)) then
    now = j
    now1 = now -1
    if( now1 <prev+1 ) then
      eachcol(i) = ''
    else
      eachcol(i) = buf(prev+1:now1)
    end if
    i=i+1
    if( i > maxch ) then
        write(0,*) 'maxch is too small'
        stop
    end if

    prev=now
  end if
end do

if( prev < ll  ) then
 eachcol(i) =  buf(prev+1:ll)
 num = i
else
 num=i-1
end if

end



subroutine setparam( year, month, start, duration_min, days )

use typedef

implicit none
integer :: year, month
type(date_time) :: start
integer :: s_year, s_month, s_day, s_hour, s_min, duration_min
integer :: days
integer, parameter :: min_of_day = 60*24

  start%year = year; start%month = month
  start%day = 1
  start%hour = 0; start%min = 0; start%sec=0

 call days_of_month( year, month, days )
 duration_min = days * min_of_day

end
 



subroutine days_of_month(year, month, days)
! number of days in year/month
implicit none
integer :: year, month, days
  logical, external :: leepyear

 integer :: m(12)
  m(1) = 31; m(2) = 28; m(3) = 31; m(4) = 30
  m(5) = 31; m(6) = 30; m(7) = 31; m(8) = 31
  m(9) = 30; m(10) = 31; m(11) = 30; m(12) = 31

 if(leepyear(year)) then
  m(2) = 29
 end if

 days = m(month)


end


logical function leepyear( year )
implicit none
 integer :: year

 if( mod(year, 4) > 0 ) then
   leepyear  = .FALSE.
 else if( mod(year, 100) > 0 ) then
   leepyear = .TRUE.
 else if( mod(year, 400) > 0 ) then
   leepyear = .FALSE.
 else
   leepyear = .TRUE.
 end if

return
end


subroutine datedif( year1, month1, day1, year2, month2, day2, diffday )
!     difference between year1/month1/day1 and year2/month2/day2
!        
  implicit none
  integer :: year1, month1, day1
  integer :: year2, month2, day2
  integer :: diffday             ! diffday =  (y/m/d)1 - (y/m/d)2
  integer :: dayofyear
  integer :: fact, tmp, tmp2, d1, d2
  logical :: leepyear
  integer :: year1a, month1a, day1a
  integer :: year2a, month2a, day2a

  year1a = year1; month1a = month1; day1a = day1
  year2a = year2; month2a = month2; day2a = day2


  if( year1 > year2 ) then
    fact = 1
  else if( year1 < year2 ) then
    fact = -1
  else 
    if( month1 > month2 ) then
      fact = 1
    else if( month1 < month2 ) then
      fact = -1
    else
      if( day1 > day2 ) then
        fact = 1
      else if( day1 < day2 ) then
        fact = -1
      else
        fact = 0
      end if
    end if
  end if


  if( fact == 0 ) then
    diffday = 0
    return
  else if( fact == -1 ) then

    tmp = year1; year1 = year2; year2 = tmp
    tmp = month1; month1 = month2; month2 = tmp
    tmp = day1; day1 = day2; day2 = tmp

  end if

 diffday = 0
 do tmp=year2, year1-1
  if( leepyear(tmp) ) then
   diffday = diffday + 366
  else
   diffday = diffday + 365
  end if
 end do
  


! count 
! write(*,*) 'year :::', diffday

   d1 = dayofyear( year1, month1, day1 )
   d2 = dayofyear( year2, month2, day2 )

   diffday = fact*(diffday + d1 - d2)
  year1 = year1a; month1 = month1a; day1 = day1a
  year2 = year2a; month2 = month2a; day2 = day2a

   return
end


   
integer function dayofyear(year, month, day)
 integer :: year, month, day
 logical :: leepyear
 integer :: i
 integer :: m(12)
  m(1) = 31; m(2) = 28; m(3) = 31; m(4) = 30
  m(5) = 31; m(6) = 30; m(7) = 31; m(8) = 31
  m(9) = 30; m(10) = 31; m(11) = 30; m(12) = 31

 if( leepyear( year) ) then
  m(2) = 29
 end if

 dayofyear = 0

 do i=1, month-1
  dayofyear = dayofyear + m(i)
 end do

 dayofyear = dayofyear + day
 return

end



subroutine checksnow( dsnow, dqsnow )
use mod_nodata
real :: dsnow, dqsnow

 if( is_nodata( dsnow ) .OR. is_nodata(dqsnow) ) then
    dsnow = NODATA;  dqsnow = NODATA
 else if( dqsnow < 152 .OR. dqsnow > 300 ) then
    dsnow = NODATA
 end if

end subroutine

