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

implicit none

!integer, parameter :: parameter_ind=1
!character (len=50), parameter :: parameter_name="surface_air_pressure (hPa)"
!integer, parameter :: parameter_ind=2
!character (len=50), parameter :: parameter_name="air_temperature (deg C)"
!integer, parameter :: parameter_ind=3
!character (len=50), parameter :: parameter_name="dew_point_temperature (deg C)"
!integer, parameter :: parameter_ind=4
!character (len=50), parameter :: parameter_name="relative_humidity (%)"
integer, parameter :: parameter_ind=5
character (len=50), parameter :: parameter_name="specific_humidity (g/kg)"
!!!!integer, parameter :: parameter_ind=6
!!!!character (len=50), parameter :: parameter_name="wind_speed (m/s)"
!!!!integer, parameter :: parameter_ind=7
!!!!character (len=50), parameter :: parameter_name="wind_from_direction (deg)"
!!!!integer, parameter :: parameter_ind=8
!!!!character (len=50), parameter :: parameter_name="u"
!!!!integer, parameter :: parameter_ind=9
!!!!character (len=50), parameter :: parameter_name="v"
!integer, parameter :: parameter_ind=10
!character (len=50), parameter :: parameter_name="precipitation_amount (kg/m2)"
!integer, parameter :: parameter_ind=11
!character (len=50), parameter :: parameter_name="surface_snow_thickness (m)"
!integer, parameter :: parameter_ind=12
!character (len=50), parameter :: parameter_name="surface_downwelling_shortwave_flux_in_air (Watt/m2)"
!integer, parameter :: parameter_ind=13
!character (len=50), parameter :: parameter_name="surface_upwelling_shortwave_flux_in_air (Watt/m2)"
!integer, parameter :: parameter_ind=14
!character (len=50), parameter :: parameter_name="surface_downwelling_longwave_flux_in_air (Watt/m2)"
!integer, parameter :: parameter_ind=15
!character (len=50), parameter :: parameter_name="surface_upwelling_longwave_flux_in_air (Watt/m2)"
!!!!integer, parameter :: parameter_ind=16
!!!!character (len=50), parameter :: parameter_name="surface_net_downward_radiative_flux"
!integer, parameter :: parameter_ind=17
!character (len=50), parameter :: parameter_name="surface_temperature (deg C)"
!!!!integer, parameter :: parameter_ind=18
!!!!character (len=50), parameter :: parameter_name="surface_downwelling_photosynthetic_photon_flux_in_air"
!!!!integer, parameter :: parameter_ind=19
!!!!character (len=50), parameter :: parameter_name="surface_upwelling_photosynthetic_photon_flux_in_air"

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

integer, parameter :: ini_base=107

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=305) :: rline
real :: lat, lon, alt, value, datetime
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_ind,'   ',parameter_name
 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, '(a305)', 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

   if (parameter_ind <= 11) then

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

     read (rline(ia_ini:ia_fin-2),'(f8.2)') value

   else

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

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

   endif

   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)

 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)

stop
end
