! monthly_report_tmp.f90
!   Mar. 2010,  by I. Tamagawa

! 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 = 30
 type awsdata
  integer :: minute
  character(len=80) :: date
  real :: dat( maxch )
 end type

 type(type_aveparam) 
  real:: prev
  real:: foll
 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,parameter :: file1 = '00_All_AWS1_10Min.csv' 
integer,parameter :: ch1 = 9, 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 ::	u       = 10+shift1 ,  v         = 11+shift1
integer,parameter ::	sphum   = 12+shift1 ,  vappress  = 13+shift1
integer,parameter ::	prec24_1= 14+shift1 ,  prec24_2  = 15+shift1
integer,parameter ::	prec2_1 = 16+shift1 ,  prec2_2   = 17+shift1
integer,parameter ::	prec1_1 = 18+shift1 ,  prec1_2   = 19+shift1


character,parameter :: file2 = '00_All_AWS2_10Min.csv'

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

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, dt
type(awsdata), allocatable  :: data(:)
type(stat), allocatable :: sdata(:,:)


 extent_min = 0 ! 24*60  ! extent of data duration for moving average etc...
 dt = 10             ! 10 minute interval data
 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) )

 call initdata( data, total_num )

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

 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)

 call mkstatistics( data, total_num, sdata, days, duration_min, sduration )

! output

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 runing_sum( data, total_num, aveparam, prec_1, prec1_1 )
  call runing_sum( data, total_num, aveparam, prec_2, prec1_2 )

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

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

end



subroutine runing_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_noda(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

 write(0,*) 'Open ', 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), *) data(i)%dat(j+shift1) 
      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)
    write(0,*) file1, 'closed'

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




