! read L1 data    by I. Tamagawa 2008 Dec.
!   mod. to output temp. & vapor density   2009. Nov.  by I. Tamagawa

module  meteor_h

  implicit none
  integer, parameter :: maxn=100000, maxch=64, maxfile=100
  type datetime
   integer :: year, yday, hour, min, ymin
   integer :: month, day
   real    :: r_yday      ! r_yday = yday + (hour + min/60)/24  
  end type datetime

end


program readL1
  use meteor_h

  type(datetime) :: time(maxn), req_time
  real           :: data(maxn, maxch)

  integer :: argc
  character( LEN=80 ) :: argv, file(maxfile)

  integer :: i, j, num=0
  integer, external :: dayofyear 

  argc = iargc()

 if(argc==0) then
   write(0, *) 'read L1 meteor data program for TKC data'
   write(0,*) './program file1 file2 file3 ....'
   stop
 end if


  write(0,*) 'required year month day hour minute'
  read(*,*) req_time%year, req_time%month, req_time%day, &
  				req_time%hour,  req_time%min
  write(0,*) req_time%year, req_time%month,  req_time%day, &
  				req_time%hour,  req_time%min
  req_time%yday = dayofyear( req_time%year, req_time%month, req_time%day)

 if(argc < maxfile) then
  do i=1, argc
  call getarg( i, argv )
  file(i)=argv
  end do
 else
  write(0,*) 'too many file specified', argc
  stop
 end if

 do i=1, argc
   call readdata( file(i), time,  data, num )
 end do


 do i=1, num
 if( req_time%year == time(i)%year    &
     .AND.  req_time%month == time(i)%month  &
     .AND.   req_time%day  == time(i)%day    &
     .AND.   req_time%hour == time(i)%hour   &
     .AND.   req_time%min == time(i)%min )  then
  write(0,*)  time(i)%year, time(i)%r_yday, time(i)%month, time(i)%day, &
   time(i)%hour, time(i)%min
   
  call vapordensity(data(i,20), data(i,30), density)
  write(*,*) data(i, 20), density

  call vapordensity(data(i,21), data(i,31), density)
  write(*,*) data(i, 21), density

!     ika do-bun  konna kanji
  end if

 end do
 

end


subroutine readdata( file, time, data, num )
 use meteor_h
 character( LEN=80 ) :: file
 integer :: num
 real ::    data(maxn, maxch)
 character( LEN=10 ) :: dum
 type(datetime) :: time(maxn)

 character( LEN=10240 ) :: dumline
 integer :: i

 i=num
 open( unit=10, file=file, err=900 ) 

  read(10, '(A)') dumline
 do while( 1==1 )
 num=num+1
 if( num > maxn ) then
   write(0,*) 'Too many data lines!.'
   go to 900
 end if
 read(10,*, end=900)   time(num)%year, time(num)%yday, dum, time(num)%ymin, &
 		(data(num, j), j=1, maxch)
 read(dum, '(I2, 1X, I2)') time(num)%hour, time(num)%min
! write(0,*) 'Check hh:mm', dum, '=>', time(num)%hour, time(num)%min


 call doy2monthday( time(num)%year, time(num)%yday, &
 			time(num)%month, time(num)%day)
 time(num)%r_yday = time(num)%yday + (time(num)%hour+time(num)%min/60.0)/24.0

 if(num==i+1) then
 	write(0,*) file, time(num)%year, time(num)%month, &
		time(num)%day, time(num)%hour, time(num)%min
 end if

 end do
  
 900 close(10)
 num = num -1


end subroutine


subroutine satvap2(T, e)
! calc. saturation vapor pressure (hPa) by Tetens
 implicit none
 real :: T   !  temperature in C deg.
 real :: e   !  saturate vapor pressure in hPa

! for water
 e = 6.1078*10**(7.5*T/(237.3+T))
! for ice
! e = 6.1078*10**(9.5*T/(265.3+T))

end




subroutine vapordensity(T, RH, density)
! T ( C. deg)  RH (0-1)  ==> density (kg/m3)
implicit none

real :: T, RH, density
real :: e, e_Pa
real, parameter :: T0 = 273.15, Rv = 461.0

  call satvap2(T, e)
!  e_Pa = (e*100)*(RH/100.0)
  e_Pa = (e*100.0)*RH

  density = e_Pa / (Rv * (T+T0))
 
end subroutine






! program diff_date_sample
! !  Sample program for calculating date in 'year month day' form
! !    2008. 5. 29   by I. Tamagawa	, (mod 12.12)
! !
! implicit none
! 
!   integer :: year1, month1, day1, year2, month2, day2, diffday, addday
!   integer :: doy
! 
!   write(*, *) 'yyyy mm dd yyyy mm dd'
!   read(*,*) year1, month1, day1, year2, month2, day2
!   call datedif( year1, month1, day1, year2, month2, day2, diffday )
!   write(*,*) diffday, ' from ', year2, month2, day2, 'to', year1, month1, day1
! 
! 
!   write(*,*) 'year day_of_year'
!   read(*,*)   year1, doy
!   call  doy2monthday( year1, doy, month1, day1)
!   write(*,*) year1, month1, day1
! 
! 
!   write(*, *) 'yyyy mm dd addday'
!   read(*,*) year1, month1, day1, addday
!   call dateadd( year1, month1, day1, addday, year2, month2, day2)
!   write(*,*) year2, month2, day2
! 
! 
! 
! end
! 

subroutine dateadd( year1, month1, day1, addday, year2, month2, day2)
! add  addday to year1, month1, day1
! results in    year2, month2, day2
  implicit none
  integer :: year1, month1, day1
  integer :: year2, month2, day2
  integer :: addday
  integer :: doy1, doy2, ddd
! functions
  integer :: dayofyear, days_in_year
  logical :: leepyear

  doy1 = dayofyear( year1, month1, day1 )

  doy2 = doy1 + addday
  if( doy2 > 0 )  then
    year2 = year1

    do while( .TRUE. )
      ddd = days_in_year(year2)
      if( doy2 > ddd) then
    	year2 = year2 + 1
	doy2  = doy2 - ddd
      else
        exit
      end if
    end do
   else    ! doy2 <= 0
    year2 = year1 
    do while( doy2 <= 0 )
      year2 = year2 - 1
      ddd = days_in_year(year2)
      doy2 = doy2 + ddd
    end do
   end if

   call doy2monthday( year2, doy2, month2, day2 )
   return

end


subroutine doy2monthday( year, doy, month, day)
! convert year and doy (day of year) into month and day
 implicit none
 integer ::  year, doy, month, day
 integer ::  tmp
! functions
 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

 tmp = 0; month=0
 do while( tmp < doy )
  month = month + 1
  tmp = tmp + m(month)
 end do

 tmp = tmp - m(month)
 day = doy - tmp
end






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

  if( leepyear( year ) ) then
  	days_in_year = 366
  else
  	days_in_year = 365
  end if
end




subroutine datedif( year1a, month1a, day1a, year2a, month2a, day2a, diffday )
!     difference between year1/month1/day1 and year2/month2/day2
!        
  implicit none
  integer :: year1a, month1a, day1a
  integer :: year2a, month2a, day2a
  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

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


  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
 
 if( year2 < year1 ) then
  do tmp=year2, year1-1
   if( leepyear(tmp) ) then
    diffday = diffday + 366
   else
    diffday = diffday + 365
   end if
  end do
 end if
  


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

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

   diffday = fact*(diffday + d1 - d2)
   return
end



logical function leepyear( year )
 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

   
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


 
real function hourdif( year1, month1, day1, hour1, min1, sec1, &
		year2, month2, day2, hour2, min2, sec2 )
!
!  hour1:min1:sec1  -  hour2:min2:sec2      output unit : hour
 implicit none
 integer :: year1, month1, day1, year2, month2, day2
 integer :: hour1, min1, sec1, hour2, min2, sec2
 integer :: diff_day
 real :: h1, h2

 h1 = hour1 + min1/60.0 + sec1/3600.0
 h2 = hour2 + min2/60.0 + sec2/3600.0

 hourdif = h1 - h2

 call datedif( year1, month1, day1, year2, month2, day2, diff_day )

 hourdif = diff_day*24.0+hourdif
 return
end

