! Read and write dataset of meteorological and soil data
! of Baltex project in Lindenberg (Germany) 2002-2009
!
! Parameters at levels! Not surface!
!
! Read and write 1 parameter only! Launch this program some times for
! different parameters.
!

implicit none

!! for BALTEX_LINDENBERG_20020101_20091231.flx
!
!integer, parameter :: parameter_ind=1
!character (len=50), parameter :: parameter_name="surface_upward_sensible_heat_flux W/m2"
!
!!integer, parameter :: parameter_ind=2
!!character (len=50), parameter :: parameter_name="surface_upward_latent_heat_flux W/m2"
!
!real, parameter :: parameter_level=2.40 ! FALKENBERG
!!real, parameter :: parameter_level=30.50 ! FOREST

! for BALTEX_LINDENBERG_20021001_20091231.stm

!integer, parameter :: parameter_ind=1
!character (len=50), parameter :: parameter_name="soil_temperature (deg C)"

integer, parameter :: parameter_ind=2
character (len=50), parameter :: parameter_name="volume_fraction_of_water_in_soil percent"

!real, parameter :: parameter_level=-1.50 ! FALKENBERG t
!real, parameter :: parameter_level=-1.20 ! FALKENBERG t
!real, parameter :: parameter_level=-1.00 ! FALKENBERG t
!real, parameter :: parameter_level=-0.90 ! FALKENBERG t q
!real, parameter :: parameter_level=-0.60 ! FALKENBERG t q
!real, parameter :: parameter_level=-0.50 ! FALKENBERG t
!real, parameter :: parameter_level=-0.45 ! FALKENBERG t q
real, parameter :: parameter_level=-0.30 ! FALKENBERG t q
!real, parameter :: parameter_level=-0.20 ! FALKENBERG t
!real, parameter :: parameter_level=-0.15 ! FALKENBERG t q
!real, parameter :: parameter_level=-0.10 ! FALKENBERG t
!real, parameter :: parameter_level=-0.08 ! FALKENBERG q
!real, parameter :: parameter_level=-0.05 ! FALKENBERG t

!real, parameter :: parameter_level=-1.50 ! FOREST t q
!real, parameter :: parameter_level=-0.90 ! FOREST t q
!real, parameter :: parameter_level=-0.60 ! FOREST t q
!real, parameter :: parameter_level=-0.30 ! FOREST t q
!real, parameter :: parameter_level=-0.20 ! FOREST t q
!real, parameter :: parameter_level=-0.10 ! FOREST t q
!real, parameter :: parameter_level=-0.05 ! FOREST t

integer, parameter :: ini_base=115

character (len=30), parameter :: point_name="FALKENBERG"
!character (len=30), parameter :: point_name="FOREST"

integer :: unit_inp=21, unit_out=31
character (len=30) :: file_name_inp="file_inp.txt", file_name_out="file_out.txt", name_curr
character (len=159) :: rline
real :: lat, lon, alt, level, datetime, value, value_min=1.e8, value_max=-1.e8
integer, dimension(2) :: year, month, day, hour, minute
integer :: val_miss=999, ierr_open=0, ierr_read=0, ini_flag=0, ia_ini, ia_fin, flag, imon, jday, dd_feb, ndyear
character (len=1) :: ind_quality, quality1="G", quality2="I", a1
integer, dimension(12) :: dd_mon=(/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/)

 print *
 print *,parameter_name,'  ',parameter_level
 print *
 print *,'   ',point_name
 print *

 open (unit_inp, file=file_name_inp, form="formatted", status="old", iostat=ierr_open)

 if (ierr_open /= 0) then
   print *
   print *,"Requested file ",file_name_inp," not found, stop"
   stop
 endif

 open (unit_out, file=file_name_out, form="formatted", status="unknown")

 ini_flag=0

 do while (.true.)

   read (unit_inp, '(a159)', iostat=ierr_read) rline
   if (ierr_read /= 0) exit

   read (rline(62:76),*) name_curr
   if (trim(name_curr) /= trim(point_name)) cycle

   if (ini_flag == 0) then
! Definition of point geografical coordinates 

     read (rline(77:107),'(2(1x,f10.5),1x,f7.2)') lat, lon, alt

     write (unit_out,'(a20,2(1x,f10.5),1x,f7.2,3x,a50)') point_name, lat, lon, alt, parameter_name

     ini_flag=val_miss

   endif

   read (rline(109:115),'(f7.2)') level

   if (level /= parameter_level) cycle

   ia_ini=ini_base+(parameter_ind-1)*11+1
   ia_fin=ia_ini+10

   read (rline(ia_ini:ia_fin-2),'(f9.2)') value
   read (rline(ia_fin:ia_fin),'(a1)') ind_quality

   if (int(value) == val_miss) cycle
   flag=-1
   if (ind_quality == quality1) flag=1
   if (ind_quality == quality2) flag=1
   if (flag /= 1) cycle

! Data value is good: decoding of data and time, writting data

   read (rline( 1:16), '(i4,4(a1,i2))' ) year(1),a1,month(1),a1,day(1),a1,hour(1),a1,minute(1)
   read (rline(18:33), '(i4,4(a1,i2))' ) year(2),a1,month(2),a1,day(2),a1,hour(2),a1,minute(2)

   dd_feb=28
   ndyear=365
   if (mod(year(1), 4) == 0) then
     dd_feb=29
     ndyear=366
   endif
   dd_mon(2)=dd_feb

   jday=0
   do imon=2,month(1)
     jday=jday+dd_mon(imon-1)
   enddo
   jday=jday+day(1)
   jday=jday-1

   datetime=(float(jday)*1440.+float(hour(1))*60.+float(minute(1)))/(float(ndyear)*1440.)

   write (unit_out, 101) year(1),datetime,value,year(1),month(1),day(1),hour(1),':',minute(1)

   if (value > value_max) value_max=value
   if (value < value_min) value_min=value

 enddo

 close (unit_inp)
 close (unit_out)

 101 format (i4.4,f9.6,1x,f8.2,5x,i4.4,2i2.2,1x,i2.2,a1,i2.2)

 print *
 print *,' Max and min values: ',value_max,value_min

stop
end
