diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 26b3529ad6..b7ddc8a552 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -22,6 +22,23 @@ individual files. The changes are now listed with the most recent at the top. +**September 18 2023 :: Fluxnet observation converter and obs_def_rttov13_mod.f90 bug-fixes. Tag v10.8.4** + +Fluxnet obs converter: + +- Generates a new observation converter (Fluxnetfull_to_obs) for eddy + covariance flux tower data (carbon, water energy fluxes) +- Documentation changes made to the older, deprecated ameriflux + converter (level4_to_obs) and the broken links have been fixed +- New flux tower observation types added to accomodate the forward + operator approach for time aggregated fluxes (daily through monthly) + +obs_def_rttov13_mod.f90 bug-fixes: + +- Added public get_channel to obs_def_rttov13_mod.f90 to compile WRF + successfully with rttov13. +- Removed cloud_overlap (integer) from the function: get_rttov_option_logical + **August 21 2023 :: CAM-FV shell scripts. Tag v10.8.3** Performance improvements for CAM-FV shell scripts: diff --git a/conf.py b/conf.py index 90264588af..7e5f834bae 100644 --- a/conf.py +++ b/conf.py @@ -21,7 +21,7 @@ author = 'Data Assimilation Research Section' # The full version, including alpha/beta/rc tags -release = '10.8.3' +release = '10.8.4' root_doc = 'index' # -- General configuration --------------------------------------------------- diff --git a/guide/important-capabilities-dart.rst b/guide/important-capabilities-dart.rst index 7abb5044e0..520fa61236 100644 --- a/guide/important-capabilities-dart.rst +++ b/guide/important-capabilities-dart.rst @@ -76,7 +76,9 @@ is as follows: +------------------------------------------------------------------------------------------------------+-------------------+-----------------------------------+ | `Aviso `__: satellite derived sea surface height | Aviso | netCDF | +------------------------------------------------------------------------------------------------------+-------------------+-----------------------------------+ -| Level 4 Flux Tower data from `AmeriFlux `__ | Ameriflux | Comma-separated text | +| Level 4 Flux Tower data | Ameriflux | Comma-separated text | ++------------------------------------------------------------------------------------------------------+-------------------+-----------------------------------+ +| Ameriflux Fullset Flux Tower data from `AmeriFlux `__ | Ameriflux | Comma-separated text | +------------------------------------------------------------------------------------------------------+-------------------+-----------------------------------+ | Level 2 soil moisture from `COSMOS `__ | COSMOS | Fixed-width text | +------------------------------------------------------------------------------------------------------+-------------------+-----------------------------------+ @@ -104,6 +106,8 @@ is as follows: +------------------------------------------------------------------------------------------------------+-------------------+-----------------------------------+ | Sea surface temperature | SST | netCDF | +------------------------------------------------------------------------------------------------------+-------------------+-----------------------------------+ +| Solar-Induced Fluorescence | SIF | netCDF | ++------------------------------------------------------------------------------------------------------+-------------------+-----------------------------------+ | Special Sensor Ultraviolet Spectrographic Imager `(SSUSI) `__ retrievals | SSUSI | netCDF | +------------------------------------------------------------------------------------------------------+-------------------+-----------------------------------+ | World Ocean Database `(WOD) `__ | WOD | World Ocean Database packed ASCII | diff --git a/observations/forward_operators/obs_def_rttov13_mod.f90 b/observations/forward_operators/obs_def_rttov13_mod.f90 index ebcf033d5e..6403995051 100644 --- a/observations/forward_operators/obs_def_rttov13_mod.f90 +++ b/observations/forward_operators/obs_def_rttov13_mod.f90 @@ -397,7 +397,8 @@ module obs_def_rttov_mod write_rttov_metadata, & interactive_rttov_metadata, & get_expected_radiance, & - get_rttov_option_logical + get_rttov_option_logical, & + get_channel ! The rttov_test.f90 program uses these, but no one else should. @@ -4238,8 +4239,6 @@ function get_rttov_option_logical(field_name) result(p) p = USER_CLD_OPT_PARAM case('GRID_BOX_AVG_CLOUD') p = GRID_BOX_AVG_CLOUD - case('CLOUD_OVERLAP') - p = cloud_overlap case('ADDPC') p = ADDPC case('ADDRADREC') diff --git a/observations/obs_converters/Ameriflux/fluxnetfull_to_obs.f90 b/observations/obs_converters/Ameriflux/fluxnetfull_to_obs.f90 new file mode 100644 index 0000000000..655eb990d3 --- /dev/null +++ b/observations/obs_converters/Ameriflux/fluxnetfull_to_obs.f90 @@ -0,0 +1,1078 @@ +! DART software - Copyright UCAR. This open source software is provided +! by UCAR, "as is", without charge, subject to all terms of use at +! http://www.image.ucar.edu/DAReS/DART/DART_download + +program fluxnetfull_to_obs + +!------------------------------------------------------------------------ +! +! fluxnetfull_to_obs - a program that converts Ameriflux/Fluxnet FULLSET eddy +! covariance tower data of NEE, GPP, RE, latent heat and sensible +! heat fluxes into DART obs_seq formatted files. Works on native time +! resolution (HH, HR) that ED tower data or collected or at coarser +! aggregated time resolution files (DD,WW,MM) based on ONEFLUX +! methodology (Pastorello et al., 2020) + + +use types_mod, only : r8, MISSING_R8 + +use utilities_mod, only : initialize_utilities, finalize_utilities, & + error_handler, E_MSG, E_ERR, & + open_file, close_file, do_nml_file, do_nml_term, & + check_namelist_read, find_namelist_in_file, & + nmlfileunit, logfileunit + +use time_manager_mod, only : time_type, set_calendar_type, GREGORIAN, & + set_date, set_time, get_time, print_time, & + print_date, operator(-), operator(+), operator(>), & + operator(<), operator(==), operator(<=), operator(>=), & + operator(/) + +use location_mod, only : VERTISHEIGHT + +use obs_sequence_mod, only : obs_sequence_type, obs_type, & + static_init_obs_sequence, init_obs, write_obs_seq, & + init_obs_sequence, get_num_obs, & + set_copy_meta_data, set_qc_meta_data + +use obs_utilities_mod, only : add_obs_to_seq, create_3d_obs + +use obs_kind_mod, only : TOWER_SENSIBLE_HEAT_FLUX, & + TOWER_NETC_ECO_EXCHANGE, & + TOWER_LATENT_HEAT_FLUX, & + TOWER_ER_FLUX, & + TOWER_GPP_FLUX + +implicit none + +character(len=*), parameter :: source = 'fluxnetfull_to_obs.f90' + +!----------------------------------------------------------------------- +! Namelist with default values +!----------------------------------------------------------------------- + +character(len=256) :: text_input_file = 'textdata.input' +character(len=256) :: obs_out_file = 'obs_seq.out' +real(r8) :: timezoneoffset = -1.0_r8 +real(r8) :: latitude = -1.0_r8 +real(r8) :: longitude = -1.0_r8 +real(r8) :: elevation = -1.0_r8 +real(r8) :: flux_height = -1.0_r8 +! A maxgooqc=3 allows for good=1, medium=2, and poor=3 quality gap-filled data +real(r8) :: maxgoodqc = 3.0_r8 +! Always true except for latent,sensible heat and NEE for hourly time periods +! Can only be false of HH or HR time resolution +logical :: gap_filled = .true. +! Option for energy balance correction for latent and sensible heat +! Recommend to keep false as these values are typically missing +logical :: energy_balance = .false. +character(len=2) :: time_resolution = 'HH' +logical :: verbose = .false. + +namelist /fluxnetfull_to_obs_nml/ text_input_file, obs_out_file, & + timezoneoffset, latitude, longitude, elevation, & + flux_height, maxgoodqc, gap_filled, energy_balance, & + time_resolution, verbose + +!----------------------------------------------------------------------- +! globally-scoped variables +!----------------------------------------------------------------------- + +character(len=5000) :: input_line, bigline +character(len=512) :: string1, string2, string3 +integer :: iline, nlines, nwords +logical :: first_obs +integer :: oday, osec, rcio, iunit +integer :: num_copies, num_qc, max_obs +real(r8) :: oerr, qc +real(r8) :: sig_gppdt, sig_gppnt, sig_recodt, sig_recont +type(obs_sequence_type) :: obs_seq +type(obs_type) :: obs, prev_obs +type(time_type) :: prev_time, offset +real(r8), parameter :: umol_to_gC = (1.0_r8/1000000.0_r8) * 12.0_r8 +logical :: res + +! Initialize with default tower strings, modify later + +type towerdata + type(time_type) :: time_obs + character(len=20) :: startstring = 'TIMESTAMP_START' + character(len=20) :: endstring = 'TIMESTAMP_END' + character(len=20) :: neestring = 'NEE_VUT_REF' + character(len=20) :: neeUNCstring = 'NEE_VUT_REF_JOINTUNC' + character(len=20) :: neeQCstring = 'NEE_VUT_REF_QC' + character(len=20) :: lestring = 'LE_F_MDS' + character(len=20) :: leUNCstring = 'LE_RANDUNC' + character(len=20) :: leQCstring = 'LE_F_MDS_QC' + character(len=20) :: hstring = 'H_F_MDS' + character(len=20) :: hUNCstring = 'H_RANDUNC' + character(len=20) :: hQCstring = 'H_F_MDS_QC' + character(len=20) :: gppDTstring = 'GPP_DT_VUT_REF' + character(len=20) :: gppNTstring = 'GPP_NT_VUT_REF' + character(len=20) :: recoDTstring = 'RECO_DT_VUT_REF' + character(len=20) :: recoNTstring = 'RECO_NT_VUT_REF' + character(len=20) :: gppNTUNC16string = 'GPP_NT_VUT_16' + character(len=20) :: gppNTUNC84string = 'GPP_NT_VUT_84' + character(len=20) :: gppDTUNC16string = 'GPP_DT_VUT_16' + character(len=20) :: gppDTUNC84string = 'GPP_DT_VUT_84' + character(len=20) :: recoNTUNC16string = 'RECO_NT_VUT_16' + character(len=20) :: recoNTUNC84string = 'RECO_NT_VUT_84' + character(len=20) :: recoDTUNC16string = 'RECO_DT_VUT_16' + character(len=20) :: recoDTUNC84string = 'RECO_DT_VUT_84' + integer :: startindex + integer :: endindex + integer :: neeindex + integer :: neeUNCindex + integer :: neeQCindex + integer :: leindex + integer :: leUNCindex + integer :: leQCindex + integer :: hindex + integer :: hUNCindex + integer :: hQCindex + integer :: gppDTindex + integer :: gppNTindex + integer :: recoDTindex + integer :: recoNTindex + integer :: gppDTUNC16index + integer :: gppDTUNC84index + integer :: gppNTUNC16index + integer :: gppNTUNC84index + integer :: recoDTUNC16index + integer :: recoDTUNC84index + integer :: recoNTUNC16index + integer :: recoNTUNC84index + + character(len=12) :: start_time + character(len=12) :: end_time + real(r8) :: nee + real(r8) :: neeUNC + integer :: neeQC + real(r8) :: neeQCfrac + real(r8) :: le + real(r8) :: leUNC + integer :: leQC + real(r8) :: h + real(r8) :: hUNC + integer :: hQC + real(r8) :: gppNT + real(r8) :: gppDT + real(r8) :: gpp + integer :: gppNTQC + integer :: gppDTQC + real(r8) :: gppNTUNC84 + real(r8) :: gppNTUNC16 + real(r8) :: gppDTUNC84 + real(r8) :: gppDTUNC16 + real(r8) :: recoNT + real(r8) :: recoDT + real(r8) :: reco + integer :: recoNTQC + integer :: recoDTQC + real(r8) :: recoNTUNC84 + real(r8) :: recoNTUNC16 + real(r8) :: recoDTUNC84 + real(r8) :: recoDTUNC16 +end type towerdata + +type(towerdata) :: tower + + +!----------------------------------------------------------------------- +! start of executable code +!----------------------------------------------------------------------- + +call initialize_utilities('fluxnetfull_to_obs') + +! Read the namelist entry +call find_namelist_in_file("input.nml", "fluxnetfull_to_obs_nml", iunit) +read(iunit, nml = fluxnetfull_to_obs_nml, iostat = rcio) +call check_namelist_read(iunit, rcio, "fluxnetfull_to_obs_nml") + +! Record the namelist values used for the run +if (do_nml_file()) write(nmlfileunit, nml=fluxnetfull_to_obs_nml) +if (do_nml_term()) write( * , nml=fluxnetfull_to_obs_nml) + +! time setup +call set_calendar_type(GREGORIAN) +offset = set_time(nint(abs(timezoneoffset)*3600.0_r8),0) +prev_time = set_time(0, 0) + +write(string1, *) 'tower located at lat, lon, elev =', latitude, longitude, elevation +write(string2, *) 'flux observations taken at =', flux_height,'m' +if (verbose) call error_handler(E_MSG,'fluxnetfull_to_obs',string1,text2=string2) + +! check the lat/lon values to see if they are ok +if (longitude < 0.0_r8) longitude = longitude + 360.0_r8 + +if (( latitude > 90.0_r8 .or. latitude < -90.0_r8 ) .or. & + (longitude < 0.0_r8 .or. longitude > 360.0_r8 )) then + + write (string2,*)'latitude should be [-90, 90] but is ',latitude + write (string3,*)'longitude should be [ 0,360] but is ',longitude + + string1 ='tower location error in input.nml &fluxnetfull_to_obs_nml' + call error_handler(E_ERR, source, string1, & + text2=string2,text3=string3) +endif + +! Provide gap-filling warning, and provide error if gap-filling +! turned on with DD or coarser time resolution + +if (gap_filled .eqv. .false.) then + + write(string1,*) 'WARNING!: gap filling is turned off. This is only recommended for' + write(string2,*) 'NEE,LE,SH observations 0=measured,gap_filled: 1=good gf,2=medium,3=poor. This approach' + write(string3,*) 'assigns QC values > maxgoodQC thus gap-filled data will not be written to obs_seq' + call error_handler(E_MSG, source, string1, & + text2=string2,text3=string3) + + if (time_resolution == 'DD' .or. time_resolution == 'MM' .or. & + time_resolution == 'WW') then + write(string1,*) 'ERROR!: gap filling can only be turned off for native time' + write(string2,*) 'resolution data (HH, HR)' + write(string3,*) 'Coarser time resolution (DD, WW, MM) must have gap_filled = .true.' + call error_handler(E_ERR, source, string1, & + text2=string2,text3=string3) + endif + +endif + + +! Modify values to towerdata strings based on user input.nml +! 1) time_resolution: DD(daily),MM(monthly),YY(yearly) uses 'TIMESTAMP' header + +if (time_resolution == 'DD' .or. time_resolution == 'MM' .or. & + time_resolution == 'YY') then + + tower%startstring = 'TIMESTAMP' + tower%endstring = 'TIMESTAMP' + res = .false. + + write(string1, *) 'Time resolution is set to =', time_resolution + write(string2, *) 'Using TIMESTAMP to set DART observation time' + if (verbose) call error_handler(E_MSG,'fluxnetfull_to_obs',string1,text2=string2) + +elseif (time_resolution == 'HH' .or. time_resolution == 'HR' .or. & + time_resolution == 'WW') then + + res = .true. + write(string1, *) 'Time resolution is set to =', time_resolution + write(string2, *) 'Using TIMESTAMP_START and TIMESTAMP_END to set: DART observation time' + if (verbose) call error_handler(E_MSG,'fluxnetfull_to_obs',string1,text2=string2) + +else + write(string1,*) 'time_resolution set incorrectly within input.nml' + write(string2,*) 'time_resolution must be HR,HH,DD,WW,MM, or YY' + call error_handler(E_ERR, source, string1,text2=string2) + +endif + +! 2) energy_balance: .true. changes sensible and latent heat strings +! Note: There are no qc values for energy_balance correction +! The qc values are manually set later in code +if (energy_balance .eqv. .true.) then + + tower%lestring = 'LE_CORR' + tower%leUNCstring = 'LE_CORR_JOINTUNC' + tower%hstring = 'H_CORR' + tower%hUNCstring = 'H_CORR_JOINTUNC' + write(string1,*) 'WARNING! Energy balance correction data turned on for LE and H' + write(string2,*) 'Check to make sure LE_CORR and H_CORR data is not missing' + call error_handler(E_MSG, source, string1,text2=string2) +else + + write(string1,*) 'Standard LE (LE_F_MDS) and H (H_F_MDS) data being used' + call error_handler(E_MSG, source, string1) +endif + +! Specify the maximum number of observations in the input file, +! but only the actual number created will be written out. +! Each line has 5 flux observations available. +! Each observation in this series will have a single +! observation value and a quality control flag. +! Initialize two empty observations - one to track location +! in observation sequence - the other is for the new observation. + +iunit = open_file(text_input_file, 'formatted', 'read') +if (verbose) call error_handler(E_MSG,'fluxnetfull_to_obs','opened input file '//trim(text_input_file)) + +nlines = count_file_lines(iunit) +max_obs = 5*nlines +num_copies = 1 +num_qc = 1 +first_obs = .true. + +call static_init_obs_sequence() +call init_obs(obs, num_copies, num_qc) +call init_obs(prev_obs, num_copies, num_qc) +call init_obs_sequence(obs_seq, num_copies, num_qc, max_obs) + +! the first one needs to contain the string 'observation' and the +! second needs the string 'QC'. +call set_copy_meta_data(obs_seq, 1, 'observation') +call set_qc_meta_data( obs_seq, 1, 'Fluxnet QC') + +! Subroutine decode_header reads obs file header and identifies the columns +! where flux variables of interest are located +rewind(iunit) +call decode_header(iunit, nwords) + +obsloop: do iline = 2,nlines + + ! read in entire text line into a buffer + read(iunit,'(A)',iostat=rcio) bigline + if (rcio < 0) exit obsloop + if (rcio > 0) then + write (string1,'(''Cannot read (error '',i3,'') line '',i8,'' in '',A)') & + rcio, iline, trim(text_input_file) + call error_handler(E_ERR,'main', string1, source) + endif + + input_line = adjustl(bigline) + + ! Parse the header line into the tower structure (including the observation time) + call stringparse(input_line, nwords, iline) + + if (iline <= 2) then + write(*,*)'' + call print_date(tower%time_obs, ' first observation date (local time) is') + call print_time(tower%time_obs, ' first observation time (local time) is') + write(*,*)'first observation raw values: (column,string,value) timezone not applied' + write(*,*)tower%hindex , tower%hstring , tower%h + write(*,*)tower%hUNCindex , tower%hUNCstring , tower%hUNC + write(*,*)tower%hQCindex , tower%hQCstring , tower%hQC + write(*,*)tower%leindex , tower%lestring , tower%le + write(*,*)tower%leUNCindex , tower%leUNCstring , tower%leUNC + write(*,*)tower%leQCindex , tower%leQCstring , tower%leQC + write(*,*)tower%neeindex , tower%neestring , tower%nee + write(*,*)tower%neeUNCindex , tower%neeUNCstring , tower%neeUNC + write(*,*)tower%neeQCindex , tower%neeQCstring , tower%neeQC + write(*,*)tower%gppDTindex , tower%gppDTstring , tower%gppDT + write(*,*)tower%gppDTUNC16index , tower%gppDTUNC16string , tower%gppDTUNC16 + write(*,*)tower%gppDTUNC84index , tower%gppDTUNC84string , tower%gppDTUNC84 + write(*,*)tower%gppDTQC + write(*,*)tower%gppNTindex , tower%gppNTstring , tower%gppNT + write(*,*)tower%gppNTUNC16index , tower%gppNTUNC16string , tower%gppNTUNC16 + write(*,*)tower%gppNTUNC84index , tower%gppNTUNC84string , tower%gppNTUNC84 + write(*,*)tower%gppNTQC + write(*,*)tower%recoDTindex , tower%recoDTstring , tower%recoDT + write(*,*)tower%recoDTUNC16index, tower%recoDTUNC16string, tower%recoDTUNC16 + write(*,*)tower%recoDTUNC84index, tower%recoDTUNC84string, tower%recoDTUNC84 + write(*,*)tower%recoDTQC + write(*,*)tower%recoNTindex , tower%recoNTstring , tower%recoNT + write(*,*)tower%recoNTUNC16index, tower%recoNTUNC16string, tower%recoNTUNC16 + write(*,*)tower%recoNTUNC84index, tower%recoNTUNC84string, tower%recoNTUNC84 + write(*,*)tower%recoNTQC + write(*,*)'' + + write(logfileunit,*)'' + call print_date(tower%time_obs, ' first observation date (local time) is',logfileunit) + call print_time(tower%time_obs, ' first observation time (local time) is',logfileunit) + write(logfileunit,*)'first observation raw values: (column,string,value) timezone not applied' + write(logfileunit,*)tower%hindex , tower%hstring , tower%h + write(logfileunit,*)tower%hUNCindex , tower%hUNCstring , tower%hUNC + write(logfileunit,*)tower%hQCindex , tower%hQCstring , tower%hQC + write(logfileunit,*)tower%leindex , tower%lestring , tower%le + write(logfileunit,*)tower%leUNCindex , tower%leUNCstring , tower%leUNC + write(logfileunit,*)tower%leQCindex , tower%leQCstring , tower%leQC + write(logfileunit,*)tower%neeindex , tower%neestring , tower%nee + write(logfileunit,*)tower%neeQCindex , tower%neeQCstring , tower%neeQC + write(logfileunit,*)tower%neeQCindex , tower%neeQCstring , tower%neeQCfrac + write(logfileunit,*)tower%gppDTUNC16index , tower%gppDTUNC16string , tower%gppDTUNC16 + write(logfileunit,*)tower%gppDTUNC84index , tower%gppDTUNC84string , tower%gppDTUNC84 + write(logfileunit,*)tower%gppDTQC + write(logfileunit,*)tower%gppNTindex , tower%gppNTstring , tower%gppNT + write(logfileunit,*)tower%gppNTUNC16index , tower%gppNTUNC16string , tower%gppNTUNC16 + write(logfileunit,*)tower%gppNTUNC84index , tower%gppNTUNC84string , tower%gppNTUNC84 + write(logfileunit,*)tower%gppNTQC + write(logfileunit,*)tower%recoDTindex , tower%recoDTstring , tower%recoDT + write(logfileunit,*)tower%recoDTUNC16index, tower%recoDTUNC16string, tower%recoDTUNC16 + write(logfileunit,*)tower%recoDTUNC84index, tower%recoDTUNC84string, tower%recoDTUNC84 + write(logfileunit,*)tower%recoDTQC + write(logfileunit,*)tower%recoNTindex , tower%recoNTstring , tower%recoNT + write(logfileunit,*)tower%recoNTUNC16index, tower%recoNTUNC16string, tower%recoNTUNC16 + write(logfileunit,*)tower%recoNTUNC84index, tower%recoNTUNC84string, tower%recoNTUNC84 + write(logfileunit,*)tower%recoNTQC + write(logfileunit,*)'' + end if + + call get_time(tower%time_obs, osec, oday) + + if (verbose) then + write(string1,*)'obs time (UTC) is (seconds,days) ',osec, oday,' obs date (UTC) is ' + call print_date(tower%time_obs, trim(string1)) + call print_date(tower%time_obs, trim(string1),logfileunit) + + + write(*,*)'' + write(string1, *) 'Display tower%start_time and tower%end_time (LTC) =', tower%start_time,' ', tower%end_time + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + write(*,*)'' + write(string1, *) 'Display tower%nee tower%neeQC tower%neeQCfrac =', tower%nee, tower%neeQC, tower%neeQCfrac + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + write(*,*)'' + write(string1, *) 'Display tower%le tower%leQC =', tower%le, tower%leQC + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + write(*,*)'' + write(string1, *) 'Display tower%h tower%hQC =', tower%h, tower%hQC + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + write(*,*)'' + write(string1, *) 'Display tower%gppDT tower%gppDTQC =', tower%gppDT, tower%gppDTQC + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + write(*,*)'' + write(string1, *) 'Display tower%gppNT tower%gppNTQC =', tower%gppNT, tower%gppNTQC + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + write(*,*)'' + write(string1, *) 'Display tower%recoDT tower%recoDTQC =', tower%recoDT, tower%recoDTQC + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + write(*,*)'' + write(string1, *) 'Display tower%recoNT tower%recoNTQC =', tower%recoNT, tower%recoNTQC + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + + endif + + ! Create and add observation and uncertainty (1 SD) to obs_seq file + ! Assign the observation the appropriate obs type + if (tower%hQC <= maxgoodqc) then ! Sensible Heat Flux [W m-2] + oerr = tower%hUNC + qc = real(tower%hQC,r8) + ! Check for missing uncertainty value, if needed assign % error + ! based on empirical examination of data + if (oerr <=0) then + select case(time_resolution) + case ('HR', 'HH') + oerr= tower%h*0.2 + case ('DD', 'WW') + oerr= tower%h*0.1 + case ('MM') + oerr= tower%h*0.05 + case default + write(string1, *) 'ERROR, time_resolution must be HH,HR,DD,WW,MM, value is:', time_resolution + call error_handler(E_ERR,'fluxnetfull_to_obs',string1) + end select + endif + call create_3d_obs(latitude, longitude, flux_height, VERTISHEIGHT, tower%h, & + TOWER_SENSIBLE_HEAT_FLUX, oerr, oday, osec, qc, obs) + call add_obs_to_seq(obs_seq, obs, tower%time_obs, prev_obs, prev_time, first_obs) + + if (verbose) then + write(*,*)'' + write(string1, *) 'Display tower%h, tower%hUNC (1 SD) =', tower%h,' ', tower%hUNC + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + endif + + endif + + if (tower%leQC <= maxgoodqc) then ! Latent Heat Flux [W m-2] + oerr = tower%leUNC + qc = real(tower%leQC,r8) + if (oerr <=0) then + select case( time_resolution ) + case ('HR', 'HH') + oerr= tower%le*0.2 + case ('DD', 'WW') + oerr= tower%le*0.1 + case ('MM') + oerr= tower%le*0.05 + case default + write(string1, *) 'ERROR, time_resolution must be HH,HR,DD,WW,MM, value is:', time_resolution + call error_handler(E_ERR,'fluxnetfull_to_obs',string1) + end select + endif + + call create_3d_obs(latitude, longitude, flux_height, VERTISHEIGHT, tower%le, & + TOWER_LATENT_HEAT_FLUX, oerr, oday, osec, qc, obs) + call add_obs_to_seq(obs_seq, obs, tower%time_obs, prev_obs, prev_time, first_obs) + + if (verbose) then + write(*,*)'' + write(string1, *) 'Display tower%le, tower%leUNC (1 SD) =', tower%le,' ', tower%leUNC + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + endif + + endif + + if (tower%neeQC <= maxgoodqc) then ! Net Ecosystem Exchange [umol m-2 s-1] + oerr = tower%neeUNC * umol_to_gC + tower%nee = -tower%nee * umol_to_gC ! Matches units in CLM [gC m-2 s-1] + qc = real(tower%neeQC,r8) + if (oerr <=0) then + select case( time_resolution ) + case ('HR', 'HH') + oerr= abs(tower%nee)*0.2 + case ('DD', 'WW') + oerr= abs(tower%nee)*0.1 + case ('MM') + oerr= abs(tower%nee)*0.05 + case default + write(string1, *) 'ERROR, time_resolution must be HH,HR,DD,WW,MM, value is:', time_resolution + call error_handler(E_ERR,'fluxnetfull_to_obs',string1) + end select + endif + call create_3d_obs(latitude, longitude, flux_height, VERTISHEIGHT, tower%nee, & + TOWER_NETC_ECO_EXCHANGE, oerr, oday, osec, qc, obs) + call add_obs_to_seq(obs_seq, obs, tower%time_obs, prev_obs, prev_time, first_obs) + + if (verbose) then + write(*,*)'' + write(string1, *) 'Display tower%nee, tower%neeUNC (1 SD) =', tower%nee,' ', tower%neeUNC + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + endif + + endif + + + + if (tower%gppDTQC <=maxgoodqc .and. tower%gppNTQC <= maxgoodqc) then ! Gross Primary Production [umol m-2 s-1] + sig_gppdt = (((tower%gppDTUNC84-tower%gppDTUNC16) / 2)**2)**0.5 ! Ustar unc contribution + sig_gppnt = (((tower%gppNTUNC84-tower%gppNTUNC16) / 2)**2)**0.5 + + oerr = (0.25 * (sig_gppdt)**2 + 0.25 * (sig_gppnt)**2)**0.5 ! Combine Ustar and partitioning unc + oerr = oerr * umol_to_gC + ! Take average of night and day partitioning methods + tower%gpp = ((tower%gppDT + tower%gppNT) / 2) * umol_to_gC ! Matches units in CLM [gC m-2 s-1] + qc = maxval((/real(tower%gppDTQC,r8),real(tower%gppNTQC,r8)/)) + if (oerr <=0) then + select case( time_resolution ) + case ('HR', 'HH') + oerr= tower%gpp*0.2 + case ('DD', 'WW') + oerr= tower%gpp*0.1 + case ('MM') + oerr= tower%gpp*0.05 + case default + write(string1, *) 'ERROR, time_resolution must be HH,HR,DD,WW,MM, value is:', time_resolution + call error_handler(E_ERR,'fluxnetfull_to_obs',string1) + end select + endif + call create_3d_obs(latitude, longitude, flux_height, VERTISHEIGHT, tower%gpp, & + TOWER_GPP_FLUX, oerr, oday, osec, qc, obs) + call add_obs_to_seq(obs_seq, obs, tower%time_obs, prev_obs, prev_time, first_obs) + + if (verbose) then + write(*,*)'' + write(string1, *) 'Display tower%gpp, gpp uncertainty (1 SD) =', tower%gpp,' ', oerr + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + endif + + endif + + if (tower%recoDTQC <= maxgoodqc .and. tower%recoNTQC <= maxgoodqc) then ! Gross Primary Production [umol m-2 s-1] + sig_recodt = (((tower%recoDTUNC84-tower%recoDTUNC16) / 2)**2)**0.5 ! Ustar unc contribution + sig_recont = (((tower%recoNTUNC84-tower%recoNTUNC16) / 2)**2)**0.5 + + oerr = (0.25 * (sig_recodt)**2 + 0.25 * (sig_recont)**2)**0.5 ! Combine Ustar and partitioning unc + oerr = oerr * umol_to_gC + ! Take average of night and day partitioning methods + tower%reco = ((tower%recoDT + tower%recoNT) / 2) * umol_to_gC ! Matches units in CLM [gC m-2 s-1] + qc = maxval((/real(tower%recoDTQC,r8),real(tower%recoNTQC,r8)/)) + if (oerr <=0) then + select case( time_resolution ) + case ('HR', 'HH') + oerr= tower%reco*0.2 + case ('DD', 'WW') + oerr= tower%reco*0.1 + case ('MM') + oerr= tower%reco*0.05 + case default + write(string1, *) 'ERROR, time_resolution must be HH,HR,DD,WW,MM, value is:', time_resolution + call error_handler(E_ERR,'fluxnetfull_to_obs',string1) + end select + endif + call create_3d_obs(latitude, longitude, flux_height, VERTISHEIGHT, tower%reco, & + TOWER_ER_FLUX, oerr, oday, osec, qc, obs) + call add_obs_to_seq(obs_seq, obs, tower%time_obs, prev_obs, prev_time, first_obs) + + if (verbose) then + write(*,*)'' + write(string1, *) 'Display tower%reco, reco uncertainty (1 SD) =', tower%reco,' ', oerr + call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + endif + + endif + + + + +end do obsloop + +call close_file(iunit) + +! If obs added to the sequence, write it out to a file now. +if ( get_num_obs(obs_seq) > 0 ) then + write(string1,*)'writing obs_seq, obs_count = ', get_num_obs(obs_seq) + if (verbose) call error_handler(E_MSG,'fluxnetfull_to_obs',string1) + call write_obs_seq(obs_seq, obs_out_file) +endif + +! end of main program +call finalize_utilities() + +contains + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! +! count_file_lines -- +! count the lines in a text file. +! rewinds the unit after counting. +! +! iunit - handle to the already-open text file +! +! created May 2, 2012 Tim Hoar, NCAR/IMAGe +! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +function count_file_lines(iunit) + +integer, intent(in) :: iunit +integer :: count_file_lines + +integer :: i +character(len=128) :: oneline + +integer, parameter :: tenmillion = 10000000 +rewind(iunit) + +count_file_lines = 0 +countloop : do i = 1,tenmillion + + read(iunit,'(A)',iostat=rcio) oneline + + if (rcio < 0) exit countloop ! end of file + if (rcio > 0) then + write (string1,'('' read around line '',i8)')i + call error_handler(E_ERR,'count_file_lines', string1, source) + endif + count_file_lines = count_file_lines + 1 + +enddo countloop +rewind(iunit) + +if (count_file_lines >= tenmillion) then + write (string1,'('' suspiciously large number of lines '',i8)')count_file_lines + call error_handler(E_MSG,'count_file_lines', string1, source) +endif + +end function count_file_lines + + + + +subroutine decode_header(iunit,ncolumns) +! Reads the first line of the obs header and identifies which columns +! the flux variable of interest is located + +integer, intent(in) :: iunit +integer, intent(out) :: ncolumns + +integer, parameter :: maxwordlength = 30 +integer :: i,charcount,columncount,wordlength +character(len=maxwordlength), dimension(:), allocatable :: columns +integer, dimension(23) :: qc = 0 + +! Read the line and strip off any leading whitespace. + + +read(iunit,'(A)',iostat=rcio) bigline +if (rcio /= 0) then + write(string1,*)'Cannot parse header. Begins <',trim(bigline(1:40)),'>' + call error_handler(E_ERR,'decode_header',string1, source) +endif + +input_line = adjustl(bigline) + +! Comma separated file, thus count commas to determine the number of columns + +charcount = CountChar(input_line,',') +ncolumns = charcount + 1 +allocate(columns(ncolumns)) + +columncount = 0 ! track the number of columns +wordlength = 0 ! number of characters in the column descriptor +charcount = 0 ! the position of the (last) comma +do i = 1,len_trim(input_line) + if (input_line(i:i) == ',') then + columncount = columncount + 1 + if (wordlength > maxwordlength) then + write(string1,*)'unexpected long word ... starts <',& + input_line((i-wordlength):(i-1)),'>' + call error_handler(E_ERR,'decode_header',string1, source) + endif + columns(columncount) = input_line((i-wordlength):(i-1)) + write(string1,*) 'word(',columncount,') is ',columns(columncount) + if (verbose) call error_handler(E_MSG,'decode_header',string1) + wordlength = 0 + charcount = i + else + wordlength = wordlength + 1 + endif +enddo + +! There is one more column after the last comma + +if ((columncount+1) /= ncolumns) then + write(string1,*)'parsed wrong number of words ...' + write(string2,*)'expected ',ncolumns,' got ',columncount+1 + call error_handler(E_ERR,'decode_header',string1,source, & + text2=trim(string2), text3=trim(input_line)) +endif + +columns(ncolumns) = input_line((charcount+1):len_trim(input_line)) + +write(string1,*)'word(',ncolumns,') is ',columns(ncolumns) +if (verbose) call error_handler(E_MSG,'decode_header',string1) + +! Finally, identify column index based on string name +tower%startindex = Match(columns, tower%startstring) +tower%endindex = Match(columns, tower%endstring) +tower%neeindex = Match(columns, tower%neestring) +tower%neeUNCindex = Match(columns, tower%neeUNCstring) +tower%neeQCindex = Match(columns, tower%neeQCstring) +tower%neeUNCindex = Match(columns, tower%neeUNCstring) +tower%leindex = Match(columns, tower%lestring) +tower%leQCindex = Match(columns, tower%leQCstring) +tower%leUNCindex = Match(columns, tower%leUNCstring) +tower%hindex = Match(columns, tower%hstring) +tower%hQCindex = Match(columns, tower%hQCstring) +tower%hUNCindex = Match(columns, tower%hUNCstring) +tower%gppDTindex = Match(columns, tower%gppDTstring) +tower%gppNTindex = Match(columns, tower%gppNTstring) +tower%recoDTindex = Match(columns, tower%recoDTstring) +tower%recoNTindex = Match(columns, tower%recoNTstring) +tower%gppDTUNC16index = Match(columns, tower%gppDTUNC16string) +tower%gppDTUNC84index = Match(columns, tower%gppDTUNC84string) +tower%gppNTUNC16index = Match(columns, tower%gppNTUNC16string) +tower%gppNTUNC84index = Match(columns, tower%gppNTUNC84string) +tower%recoDTUNC16index = Match(columns, tower%recoDTUNC16string) +tower%recoDTUNC84index = Match(columns, tower%recoDTUNC84string) +tower%recoNTUNC16index = Match(columns, tower%recoNTUNC16string) +tower%recoNTUNC84index = Match(columns, tower%recoNTUNC84string) + +! Confirm indices were found successfully +qc( 1) = CheckIndex( tower%startindex , tower%startstring) +qc( 2) = CheckIndex( tower%endindex , tower%endstring) +qc( 3) = CheckIndex( tower%neeindex , tower%neestring) +qc( 4) = CheckIndex( tower%neeUNCindex , tower%neeUNCstring) +qc( 5) = CheckIndex( tower%neeQCindex , tower%neeQCstring) +qc( 6) = CheckIndex( tower%leindex , tower%lestring) +qc( 7) = CheckIndex( tower%leQCindex , tower%leQCstring) +qc( 8) = CheckIndex( tower%leUNCindex , tower%leUNCstring) +qc( 9) = CheckIndex( tower%hindex , tower%hstring) +qc(10) = CheckIndex( tower%hQCindex , tower%hQCstring) +qc(11) = CheckIndex( tower%hUNCindex , tower%hUNCstring) +qc(12) = CheckIndex( tower%gppDTindex , tower%gppDTstring) +qc(13) = CheckIndex( tower%gppNTindex , tower%gppNTstring) +qc(14) = CheckIndex( tower%recoDTindex , tower%recoDTstring) +qc(15) = CheckIndex( tower%recoNTindex , tower%recoNTstring) +qc(16) = CheckIndex( tower%gppDTUNC16index , tower%gppDTUNC16string) +qc(17) = CheckIndex( tower%gppDTUNC84index , tower%gppDTUNC84string) +qc(18) = CheckIndex( tower%gppNTUNC16index , tower%gppNTUNC16string) +qc(19) = CheckIndex( tower%gppNTUNC84index , tower%gppNTUNC84string) +qc(20) = CheckIndex( tower%recoDTUNC16index, tower%recoDTUNC16string) +qc(21) = CheckIndex( tower%recoDTUNC84index, tower%recoDTUNC84string) +qc(22) = CheckIndex( tower%recoNTUNC16index, tower%recoNTUNC16string) +qc(23) = CheckIndex( tower%recoNTUNC84index, tower%recoNTUNC84string) + +if (any(qc /= 0) ) then + write(string1,*)'Did not find all the required column indices.' + call error_handler(E_ERR,'decode_header',string1, source) +endif + +if (verbose) then +110 format('index for ',A20,' is ',i3) + write(*,110)tower%startstring ,tower%startindex + write(*,110)tower%endstring , tower%endindex + write(*,110)tower%neestring , tower%neeindex + write(*,110)tower%neeUNCstring , tower%neeUNCindex + write(*,110)tower%neeQCstring , tower%neeQCindex + write(*,110)tower%neeUNCstring , tower%neeUNCindex + write(*,110)tower%lestring , tower%leindex + write(*,110)tower%leQCstring , tower%leQCindex + write(*,110)tower%leUNCstring , tower%leUNCindex + write(*,110)tower%hstring , tower%hindex + write(*,110)tower%hQCstring , tower%hQCindex + write(*,110)tower%hUNCstring , tower%hUNCindex + write(*,110)tower%gppDTstring , tower%gppDTindex + write(*,110)tower%gppNTstring , tower%gppNTindex + write(*,110)tower%recoDTstring , tower%recoDTindex + write(*,110)tower%recoNTstring , tower%recoNTindex + write(*,110)tower%gppDTUNC16string , tower%gppDTUNC16index + write(*,110)tower%gppDTUNC84string , tower%gppDTUNC84index + write(*,110)tower%gppNTUNC16string , tower%gppNTUNC16index + write(*,110)tower%gppNTUNC84string , tower%gppNTUNC84index + write(*,110)tower%recoDTUNC16string, tower%recoDTUNC16index + write(*,110)tower%recoDTUNC84string, tower%recoDTUNC84index + write(*,110)tower%recoNTUNC16string, tower%recoNTUNC16index + write(*,110)tower%recoNTUNC84string, tower%recoNTUNC84index +endif + +deallocate(columns) + +end subroutine decode_header + + + +function CountChar(str1,solo) +! Count the number of instances of the single character in a character string. +! useful when parsing a comma-separated list, for example. +! Count the commas and add 1 to get the number of items in the list. + +integer :: CountChar +character(len=*), intent(in) :: str1 +character, intent(in) :: solo + +integer :: i + +CountChar = 0 +do i = 1,len_trim(str1) + if (str1(i:i) == solo) CountChar = CountChar + 1 +enddo + +end function CountChar + + + +function Match(sentence,word) +! Determine the first occurrence of the 'word' in a sentence. +! In this context, a sentence is a character array, the dimension +! of the array is the number of words in the sentence. +! This is a case-sensitive match. Trailing blanks are removed. + +integer :: Match +character(len=*), dimension(:), intent(in) :: sentence +character(len=*), intent(in) :: word + +integer :: i + +Match = 0 +WordLoop : do i = 1,size(sentence) + if (trim(sentence(i)) == trim(word)) then + Match = i + return + endif +enddo WordLoop + +end function Match + + + +function CheckIndex( myindex, context ) +! Routine to issue a warning if the index was not found. +! Returns an error code ... 0 means the index WAS found +! a negative number means the index was NOT found - an error condition. +! ALL indices checked before fatally ending. + +integer :: CheckIndex +integer, intent(in) :: myindex +character(len=*), intent(in) :: context + +if (myindex == 0) then + write(string1,*)'Did not find column header matching ',trim(context) + call error_handler(E_MSG,'decode_header:CheckIndex',string1, source) + CheckIndex = -1 ! not a good thing +else + CheckIndex = 0 ! Good to go +endif + +end function CheckIndex + + + +subroutine stringparse(str1, nwords, linenum) + +character(len=*), intent(in) :: str1 +integer , intent(in) :: nwords +integer , intent(in) :: linenum + +real(r8), allocatable, dimension(:) :: values + + +integer :: yeara, yearb, montha, monthb, daya, dayb +integer :: houra, hourb, mina, minb + +integer :: time_adjust +type(time_type) :: date_start, date_end + +if (res .eqv. .true.) then + time_adjust = 2 +else + time_adjust = 1 +endif + + + +allocate(values(nwords-time_adjust)) +values = MISSING_R8 + +! First two/one element(s) of string read in as character (time stamp) +! remainder of line read in as reals. + +if (res .eqv. .true.) then + + read(str1,*,iostat=rcio) tower%start_time,tower%end_time,values + if (rcio /= 0) then + write(string1,*)'Cannot parse line',linenum,'. Begins <',trim(str1(1:40)),'>' + call error_handler(E_ERR,'stringparse',string1, source) + endif +else + read(str1,*,iostat=rcio) tower%start_time,values + if (rcio /= 0) then + write(string1,*)'Cannot parse line',linenum,'. Begins <',trim(str1(1:40)),'>' + call error_handler(E_ERR,'stringparse',string1, source) + endif + tower%end_time=tower%start_time +endif + + +! Assign flux observations, uncertainty and QC to tower structure +! Description of FLUXNET FULLSET variables + +! start_time,end_time format YYYYMMDDHHMM +! nee units [umolCO2 m-2 s-1], VUT_REF +! neeUNC units [umolCO2 m-2 s-1], joint +! neeQC (HH,HR) no dim 0=measured, gap_filled:1=good,2=medium,3=poor,-9999=missing +! neeQC (DD-MM) no dim fraction between 0-1, % of good quality filled data +! le units [W m-2] +! leUNC units [W m-2], random +! leQC (HH,HR) no dim 0=measured, gap_filled:1=good,2=medium,3=poor,-9999=missing +! leQC (DD-MM) no dim fraction between 0-1, % of good quality filled data +! h units [W m-2] +! hUNC units [W m-2], random +! hQC no dim 0=measured, gap_filled:1=good,2=medium,3=poor,-9999=missing +! hQC (DD-MM) no dim fraction between 0-1, % of good quality filled data +! gppDT units [umolCO2 m-2 s-1] VUT_REF daytime partition +! gppNT units [umolCO2 m-2 s-1] VUT_REF nighttime partition +! gppUNCDT[xx] units [umolCO2 m-2 s-1] 16,84 percentile +! gppUNCNT[xx] units [umolCO2 m-2 s-1] 16,84 percentile +! recoDT units [umolCO2 m-2 s-1] VUT_REF daytime partition +! recoNT units [umolCO2 m-2 s-1] VUT_REF nighttime partition +! recoUNCDT[xx] units [umolCO2 m-2 s-1] 16,84 percentile +! recoUNCNT[xx] units [umolCO2 m-2 s-1] 16,84 percentile + +! Convert to 'CLM-friendly' units AFTER we determine observation error variance. +! That happens in the main routine. + +! CLM history file names and units +! (CLM) NEP,GPP,ER units [gC m-2 s-1] +! (CLM) EFLX_LH_TOT_R,FSH units [W m-2] + +tower%nee = values(tower%neeindex -time_adjust ) +tower%neeUNC = values(tower%neeUNCindex -time_adjust) +tower%neeQC = nint(values(tower%neeQCindex -time_adjust)) +tower%neeQCfrac = values(tower%neeQCindex -time_adjust) +tower%le = values(tower%leindex -time_adjust ) +tower%leUNC = values(tower%leUNCindex -time_adjust ) +tower%leQC = nint(values(tower%leQCindex -time_adjust)) +tower%h = values(tower%hindex -time_adjust ) +tower%hUNC = values(tower%hUNCindex -time_adjust ) +tower%hQC = nint(values(tower%hQCindex -time_adjust)) +tower%gppDT = values(tower%gppDTindex -time_adjust) +tower%gppNT = values(tower%gppNTindex -time_adjust) +tower%recoDT = values(tower%recoDTindex -time_adjust) +tower%recoNT = values(tower%recoNTindex -time_adjust) +tower%gppDTUNC16 = values(tower%gppDTUNC16index -time_adjust) +tower%gppDTUNC84 = values(tower%gppDTUNC84index -time_adjust) +tower%gppNTUNC16 = values(tower%gppNTUNC16index -time_adjust) +tower%gppNTUNC84 = values(tower%gppNTUNC84index -time_adjust) +tower%recoDTUNC16 = values(tower%recoDTUNC16index-time_adjust) +tower%recoDTUNC84 = values(tower%recoDTUNC84index-time_adjust) +tower%recoNTUNC16 = values(tower%recoNTUNC16index-time_adjust) +tower%recoNTUNC84 = values(tower%recoNTUNC84index-time_adjust) +deallocate(values) + + + +read(tower%start_time(1:12), fmt='(i4, 4i2)') yeara,montha,daya,houra,mina +read(tower%end_time(1:12), fmt='(i4, 4i2)') yearb,monthb,dayb,hourb,minb + + +! Certain time resolutions (MM) does not define days +if (time_resolution == 'MM') then + daya = 1 + dayb = 1 +endif + + + +write(*,*)'' +write(string1, *) 'Display tower%start_time,yeara,montha,daya,houra,mina (LTC) =', tower%start_time, yeara, montha, daya, houra, mina +write(string2, *) 'Display tower%end_time,yearb,monthb,dayb,hourb,minb (LTC) =', tower%end_time, yearb, monthb, dayb, hourb, minb +if (verbose) call error_handler(E_MSG,'fluxnetfull_to_obs',string1,text2=string2) + + + +date_start= set_date(yeara,montha,daya,houra,mina,0) +date_end= set_date(yearb,monthb,dayb,hourb,minb,0) + + + +! Assign average of flux time window to obs_seq (DART time) +tower%time_obs = (date_start+date_end) / 2 + +! Covert from Fluxnet provided LTC to UTC, UTC is standard for DART and CLM +! For example, EST = UTC-5 +if (timezoneoffset < 0.0_r8) then + tower%time_obs = tower%time_obs + offset +else + tower%time_obs = tower%time_obs - offset +endif + + +! Reject NEE data where neeQC is missing +if (tower%neeQC < 0) tower%neeQC = maxgoodqc + 100 + +! If NEE (DD,WW,MM) convert from % filled QC to integer QC good/fair/poor +if (time_resolution == 'DD' .or. time_resolution == 'MM' .or. & + time_resolution == 'WW') then + if (tower%neeQCfrac >= 0.90) tower%neeQC = 1 ! >90% fill is good + if (tower%neeQCfrac < 0.90 .and. tower%neeQC >= 0.60) tower%neeQC = 2 ! 60-90% fill is fair + if (tower%neeQCfrac < 0.60 .and. tower%neeQC >= 0.0_r8) tower%neeQC = 3 ! <60% fill is poor +endif + + +! The QC values are typically missing for le and h (-9999) +! Thus manually assign fair QC values in these cases +! When energy balance is turned off +if (energy_balance .eqv. .false.) then + if (tower%leQC < 0.0_r8) tower%leQC = 2 + if (tower%hQC < 0.0_r8) tower%hQC = 2 +else ! No QC values for energy balance corrected le and h. Assign fair QC. + tower%leQC = 2 + tower%hQC = 2 +endif + + +! No qc values for gpp/reco, thus assign fair qc unless +! the gpp/reco values are negative values, then reject. +tower%gppNTQC = 2 +tower%gppDTQC = 2 +tower%recoNTQC = 2 +tower%recoDTQC = 2 + +if (tower%gppNT < 0.0_r8) tower%gppNTQC = maxgoodqc + 100 +if (tower%gppDT < 0.0_r8) tower%gppDTQC = maxgoodqc + 100 +if (tower%recoNT < 0.0_r8) tower%recoNTQC = maxgoodqc + 100 +if (tower%recoDT < 0.0_r8) tower%recoDTQC = maxgoodqc + 100 + +! Assign very bad qc to gap_filled data if user requests it +! such that maxgoodqc threshold does not add gap_filled data to obs_seq file +! If leQC and hQC are missing values already forced to 2, thus rejected +if (gap_filled .eqv. .false.) then + if (tower%neeQC >0) tower%neeQC = maxgoodqc + 100 + if (tower%leQC >0) tower%leQC = maxgoodqc + 100 + if (tower%hQC >0) tower%hQC = maxgoodqc + 100 + ! GPP and RECO are modeled data and considered gap_filled + tower%gppNTQC = maxgoodqc + 100 + tower%gppDTQC = maxgoodqc + 100 + tower%recoNTQC = maxgoodqc + 100 + tower%recoDTQC = maxgoodqc + 100 +endif + +end subroutine stringparse + + + +end program fluxnetfull_to_obs + + diff --git a/observations/obs_converters/Ameriflux/fluxnetfull_to_obs.nml b/observations/obs_converters/Ameriflux/fluxnetfull_to_obs.nml new file mode 100644 index 0000000000..02d958176b --- /dev/null +++ b/observations/obs_converters/Ameriflux/fluxnetfull_to_obs.nml @@ -0,0 +1,15 @@ +&fluxnetfull_to_obs_nml + text_input_file = 'textdata.input', + obs_out_file = 'obs_seq.out', + timezoneoffset = -1, + latitude = -1.0, + longitude = -1.0, + elevation = -1.0, + flux_height = -1.0, + maxgoodqc = 3, + gap_filled = .true. + energy_balance = .false. + time_resolution = 'HR' + verbose = .false. + / + diff --git a/observations/obs_converters/Ameriflux/fluxnetfull_to_obs.rst b/observations/obs_converters/Ameriflux/fluxnetfull_to_obs.rst new file mode 100644 index 0000000000..825b3a27ef --- /dev/null +++ b/observations/obs_converters/Ameriflux/fluxnetfull_to_obs.rst @@ -0,0 +1,399 @@ +PROGRAM ``fluxnetfull_to_obs`` +============================== + +Overview +-------- + +FLUXNET FULLSET data to DART observation sequence file +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +| This routine is designed to convert Ameriflux FLUXNET data using the FULLSET format into a DART obs_seq file. + The data can be downloaded from the `AmeriFlux server. `__ + or `FLUXNET server. `__ The AmeriFlux network + is part of `FLUXNET `__ and the converter is suitable for either network given the similarity in + the data format. This code is not compatible with the SUBSET format of Ameriflux data, however, this code could be + modified for that purpose. The FULLSET format contains additional uncertainty information as described in the ONEflux + processing technique as described by `Pastorello et al., (2020) `__ + which is required to formally propogate error into aggregated estimates of tower flux data. + +| The AmeriFlux/FLUXNET data products use local time, therefore the code converts the time into UTC such that it is compatible + with DART and CLM. More information about the AmeriFlux data product is provided below and also within the `flux variable + description webpage `__. + +The steps required to prepare Ameriflux data for an assimilation usually include: + +#. Download the Ameriflux or FLUXNET FULLSET data for the towers and years in question (see DATA SOURCES below) +#. Record the TIME ZONE, latitude, longitude, and elevation and tower height at each site. This tower metadata can be found + `here `__ or `here `__. +#. Manually provide tower metadata information via the ``fluxnetfull_to_obs_nml`` namelist as this information is + not contained in the data file itself. +#. Build the DART executables with support for the tower observations. This is done by running ``preprocess`` with + ``obs_def_tower_mod.f90`` in the list of ``input_files`` for ``preprocess_nml``. +#. Convert each Ameriflux data file individually using ``fluxnetfull_to_obs`` +#. If necessary, combine all output files for the region and timeframe of interest into one file using + :doc:`../../../assimilation_code/programs/obs_sequence_tool/obs_sequence_tool` + + +.. Note:: + The Ameriflux data typically have all years combined into one file. However, for some models (CLM, for example), + it is required to reorganize the observation sequence files into a series of files that contains ONLY the observations + for each assimilation time step. This can be achieved with the `makedaily.sh `__ script. + +Namelist +-------- + +This namelist is read from the file ``input.nml``. Namelists start with an ampersand '&' and terminate with a slash '/'. +Character strings that contain a '/' must be enclosed in quotes to prevent them from prematurely terminating the +namelist. + +:: + + &fluxnetfull_to_obs_nml + text_input_file = 'textdata.input', + obs_out_file = 'obs_seq.out', + timezoneoffset = -1, + latitude = -1.0, + longitude = -1.0, + elevation = -1.0, + flux_height = -1.0, + maxgoodqc = 3, + gap_filled = .true. + energy_balance = .false. + time_resolution = 'HR' + verbose = .false. + / + +.. container:: + + +-----------------+--------------------+-----------------------------------------------------------------------------+ + | Contents | Type | Description | + +=================+====================+=============================================================================+ + | text_input_file | character(len=128) | Name of the Ameriflux ASCII file of comma-separated values. This may be a | + | | | relative or absolute filename. | + +-----------------+--------------------+-----------------------------------------------------------------------------+ + | obs_out_file | character(len=128) | Name of the output observation sequence file. | + +-----------------+--------------------+-----------------------------------------------------------------------------+ + | timezoneoffset | real | The local time zone difference (in hours) from UTC of the flux tower site. | + | | | The code will convert from local time to UTC. | + +-----------------+--------------------+-----------------------------------------------------------------------------+ + | latitude | real | Latitude (in degrees N) of the tower. | + +-----------------+--------------------+-----------------------------------------------------------------------------+ + | longitude | real | Longitude (in degrees E) of the tower. For internal consistency, DART uses | + | | | longitudes in the range [0,360]. An input value of -90 will be converted to | + | | | 270, for example. | + +-----------------+--------------------+-----------------------------------------------------------------------------+ + | elevation | real | Surface elevation (in meters) of the tower. | + +-----------------+--------------------+-----------------------------------------------------------------------------+ + | flux_height | real | Height (in meters) of the flux instrument on the tower. | + +-----------------+--------------------+-----------------------------------------------------------------------------+ + | maxgoodqc | real | Maximum value of the observation (obs) quality control (QC) value to pass | + | | | to the output observation sequence file. Obs that have a QC exceeding this | + | | | value will be rejected from the observation sequence file. The 'filter' step| + | | | is also capable of discriminating obs based on the QC value, therefore it | + | | | is recommended to keep all obs during the conversion unless you | + | | | are certain they are bad/missing. | + +-----------------+--------------------+-----------------------------------------------------------------------------+ + | gap_filled | logical | If set to .false. excludes all observations except for measured NEE, SH and | + | | | and LE. Must be set to .true. for time resolutions greater than hourly | + | | | (HH) or half-hourly (HR). | + +-----------------+--------------------+-----------------------------------------------------------------------------+ + | energy_balance | logical | If .true. applies energy balance closure to SH and LE fluxes, otherwise | + | | | applies standard approach. | + +-----------------+--------------------+-----------------------------------------------------------------------------+ + | time_resolution | logical | Defines the flux observation time resolution. Possible values are native | + | | | resolution of hourly (HH) or half-hourly (HR), or coarser gap-filled | + | | | resolutions including daily (DD), weekly (WW) or monthly (MM). | + +-----------------+--------------------+-----------------------------------------------------------------------------+ + | verbose | logical | If .true. outputs diagnostic data to help with troubleshooting. | + +-----------------+--------------------+-----------------------------------------------------------------------------+ + +Data sources +------------ + +| The data was acquired from https://ameriflux.lbl.gov/data/download-data/ + by choosing the data product ``Ameriflux FLUXNET`` and variable set ``FULLSET``. + The filenames are organized by site and time (HH,HR,WW,DD,MM,YY) and look like: +| ``AMF_US-Ha1_FLUXNET_FULLSET_HR_1991-2020_3-5.csv``, +| ``AMF_US-xBR_FLUXNET_FULLSET_MM_2017-2021_3-5.csv`` + +| The Ameriflux product in question are ASCII files of comma-separated values taken from all years the data is available. + The first line is a comma-separated list of column descriptors, and all subsequent lines + are comma-separated numerical values. The converter presently searches for the columns pertaining to carbon, water + and energy fluxes, corresponding quality control fields, uncertainty values, and the time of the observation. Data is available + at varying time resolutions incuding: native resolution, hourly (HR) or half-hourly (HH), and aggregated resolution, daily (DD), + weekly (WW), and monthly (MM). The source data does include yearly (YY) time resolution as well, but the coarse nature of yearly flux + observations poorly constrain fast changing ecological process, thus are not supported by this converter. The required column + headers depend upon the namelist definitions, including the ``time_resolution``, ``energy_balance`` and ``gap_filled`` settings. + These variables are defined as follows: + + + +.. container:: + + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | Ameriflux Units | Ameriflux Variable | Description | DART type | DART kind | DART units | + +=================+======================+===============================+==========================+=============================+===============+ + | YYYYMMDDHHMM | TIMESTAMP_START | start of time window | N/A | N/A | Gregorian | + | | TIMESTAMP_END | end of time window | | | | + | | | (HH,HR,WW only) | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | YYYYMMDDHHMM | TIMESTAMP | time (DD and MM only) | N/A | N/A | Gregorian | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | W/m^2 | LE_F_MDS | Latent Heat (LE) Flux | TOWER_LATENT_HEAT_FLUX | QTY_LATENT_HEAT_FLUX | W/m^2 | + | | | energy_balance = .false. | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | W/m^2 | LE_RANDUNC | Uncertainty for LE Flux | N/A | N/A | W/m^2 | + | | | energy_balance = .false. | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | [0-3] integer | LE_F_MDS_QC | QC for LE Flux | N/A | N/A | [0-3] integer | + | | | energy_balance = .false. | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | W/m^2 | LE_CORR | Latent Heat (LE) Flux | TOWER_LATENT_HEAT_FLUX | QTY_LATENT_HEAT_FLUX | W/m^2 | + | | | energy_balance = .true. | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | W/m^2 | LE_CORR_JOINTUNC | Uncertainty for LE Flux | N/A | N/A | W/m^2 | + | | | energy_balance = .true. | | | | + | | | Random and Ustar contributions| | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | W/m^2 | H_F_MDS | Sensible Heat (SH) Flux | TOWER_SENSIBLE_HEAT_FLUX | QTY_SENSIBLE_HEAT_FLUX | W/m^2 | + | | | energy_balance = .false. | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | W/m^2 | H_RANDUNC | Uncertainty for SH Flux | N/A | N/A | W/m^2 | + | | | energy_balance = .false. | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | [0-3] integer | H_F_MDS_QC | QC for SH Flux | N/A | N/A | [0-3] integer | + | | | energy_balance = .false. | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | W/m^2 | H_CORR | Sensible Heat (SH) Flux | TOWER_SENSIBLE_HEAT_FLUX | QTY_SENSIBLE_HEAT_FLUX | W/m^2 | + | | | energy_balance = .true. | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | W/m^2 | H_CORR_JOINTUNC | Uncertainty for SH Flux | N/A | N/A | W/m^2 | + | | | energy_balance = .true. | | | | + | | | Random and Ustar contributions| | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | NEE_VUT_REF | Net Ecosystem Exchange (NEE) | TOWER_NETC_ECO_EXCHANGE | QTY_NET_CARBON_PRODUCTION | gC/m^2/s | + | | | Variable Ustar, reference | | | | + | | | flux approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | NEE_VUT_REF_JOINTUNC | Uncertainty for NEE Flux | N/A | N/A | gC/m^2/s | + | | | Variable Ustar, reference | | | | + | | | flux approach. Random and | | | | + | | | and Ustar contributions | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | [0-3] integer | NEE_VUT_REF_QC | QC for NEE Flux | N/A | N/A | [0-3] integer | + | | | Variable Ustar, reference | | | | + | | | flux approach. | | | | + | | | (HH or HR only) | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | [0-1] fraction | NEE_VUT_REF_QC | QC for NEE Flux | N/A | N/A | [0-3] integer | + | | | Variable Ustar, reference | | | | + | | | flux approach. | | | | + | | | (DD, WW and MM only) | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | GPP_DT_VUT_REF | Gross Primary Production (GPP)| TOWER_GPP_FLUX | QTY_GROSS_PRIMARY_PROD_FLUX | gC/m^2/s | + | | | Day partition, Variable | | | | + | | | Ustar, reference approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | GPP_NT_VUT_REF | Gross Primary Production (GPP)| TOWER_GPP_FLUX | QTY_GROSS_PRIMARY_PROD_FLUX | gC/m^2/s | + | | | Night partition, Variable | | | | + | | | Ustar, reference approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | GPP_DT_VUT_16 | 16th percentile uncertainty | | | | + | | | estimate for GPP Flux | N/A | N/A | gC/m^2/s | + | | | Day partition, Variable | | | | + | | | Ustar, reference approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | GPP_DT_VUT_84 | 84th percentile uncertainty | | | | + | | | estimate for GPP Flux | N/A | N/A | gC/m^2/s | + | | | Day partition, Variable | | | | + | | | Ustar, reference approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | GPP_NT_VUT_16 | 16th percentile uncertainty | | | | + | | | estimate for GPP Flux | N/A | N/A | gC/m^2/s | + | | | Night partition, Variable | | | | + | | | Ustar, reference approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | GPP_NT_VUT_84 | 84th percentile uncertainty | | | | + | | | estimate for GPP Flux | N/A | N/A | gC/m^2/s | + | | | Night partition, Variable | | | | + | | | Ustar, reference approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | RECO_DT_VUT_REF | Ecosystem Respiration (ER) | TOWER_ER_FLUX | QTY_ER_FLUX | gC/m^2/s | + | | | Day partition, Variable | | | | + | | | Ustar, reference approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | RECO_NT_VUT_REF | Ecosystem Respiration (ER) | TOWER_ER_FLUX | QTY_ER_FLUX | gC/m^2/s | + | | | Night partition, Variable | | | | + | | | Ustar, reference approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | RECO_DT_VUT_16 | 16th percentile uncertainty | | | | + | | | estimate for ER Flux | N/A | N/A | gC/m^2/s | + | | | Day partition, Variable | | | | + | | | Ustar, reference approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | RECO_DT_VUT_84 | 84th percentile uncertainty | | | | + | | | estimate for ER Flux | N/A | N/A | gC/m^2/s | + | | | Day partition, Variable | | | | + | | | Ustar, reference approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | RECO_NT_VUT_16 | 16th percentile uncertainty | | | | + | | | estimate for ER Flux | N/A | N/A | gC/m^2/s | + | | | Night partition, Variable | | | | + | | | Ustar, reference approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + | umolCO2/m^2/s | RECO_NT_VUT_84 | 84th percentile uncertainty | | | | + | | | estimate for ER Flux | N/A | N/A | gC/m^2/s | + | | | Night partition, Variable | | | | + | | | Ustar, reference approach | | | | + +-----------------+----------------------+-------------------------------+--------------------------+-----------------------------+---------------+ + + +Carbon Fluxes +------------- +The flux data files come with several different approaches to estimate the carbon fluxes (NEE, GPP and ER). In the conversion code we choose the +**variable Ustar with reference flux approach (_VUT_REF)**. This choice was based on guidance from Pastorello et al., (2020) which states: + +`"The variable proposed in the SUBSET (or FULLSET) product is NEE_VUT_REF since it maintains the temporal variability (as opposed to the MEAN NEE), +it is representative of the ensemble, and the VUT method is sensitive to possible changes of the canopy (density and height) and site setup, +which can have an impact on the turbulence and consequently on the USTAR threshold. The RECO and GPP products in SUBSET (or FULLSET) are calculated +from the corresponding NEE variables filtered with the VUT method, generating RECO_NT_VUT_REF and RECO_DT_VUT_REF for ER, and +GPP_NT_VUT_REF and GPP_DT_VUT_REF for GPP. It is important to use both daytime (DT) and nighttime (NT) variables, and consider their +difference as uncertainty."` + +The reference NEE (_REF) is selected on the basis of model efficiency (MEF) thus is the single NEE data set (out of all ensemble members that sample +the Ustar uncertainty) that is most representative. The mean of the night and day partitioning methods for GPP and ER are used as the observation within +the converter. + + +Uncertainty +----------- +Multiple methods are used to estimate uncertainty within the conversion code. There are, in general, three separate sources of uncertainty +in flux data. First, **random uncertainty**, represents the random movement of eddies within the atmosphere where smaller eddies are sampled +more frequently and larger eddies less frequently. Random uncertainty (_RANDC) estimates are based on `Hollinger, D. Y. & Richardson, A. D. +Uncertainty in eddy covariance measurements and its application to physiological models. Tree Physiol. 25, 873–885 (2005)`. +Second, **Ustar uncertainty**, represents the uncertainty contribution from low turbulence conditions as calculated from the Ustar (friction velocity) +threshold. The ONEflux method uses a bootstrap sampling method to generate a 200 member ensemble from which the flux percentiles are estimated. +The third source of uncertainty, **partitioning uncertainty**, applies to GPP and ER only. The night (Reichsten et al., 2005) and day (Lasslop et al.,) partitioning methods +estimate both the contributions of photosynthesis (GPP) and ecosystem respiration (ER) as measured from the net carbon exchange (NEE). + +Within the conversion code, the flux uncertainty values (_RANDUNC) account for random uncertainty (SH and LE), where uncertainty denoted with (_JOINTUNC) indicates +combined uncertainty of random and energy balance closure uncertainty (SH and LE). The NEE uncertainty (NEE_VUT_REF_JOINTUNC) accounts for both random and Ustar contributions, +whereas the GPP and ER uncertainty combine both Ustar and partitioning method uncertainty as: + + | ``Day method Ustar uncertainty (sigma_fluxdt) = (((fluxDTUNC84-fluxDTUNC16) / 2)^2)^0.5`` + | ``Night method Ustar uncertainty (sigma_fluxnt) = (((fluxNTUNC84-fluxNTUNC16) / 2)^2)^0.5`` + + | ``Ustar and partitioning uncertainty (sigma) = (0.25 * (sigma_fluxdt)^2 + 0.25 * (sigma_fluxnt)^2)^0.5`` + +where ``flux`` stands for either GPP or ER. The 84th and 16th percentile estimates are used to generate 1 sigma estimates +for day and night Ustar uncertainty respectively. Then the contributions of the day and night Ustar uncertainty +are propogated together through standard technqiues assuming gaussian uncertainty distributions (Taylor et al, +An Introduction to Error Analysis). + +.. Note:: + + In practice the relative uncertainty reduces as the flux time resolution increases. This is likely a result of random uncertainty + decreasing with coarser time resolutions as the sample size of measurements increases. This reduced relative + uncertainty will cause a stronger impact of the observations on the prior model state (i.e. larger increments) during the ``filter`` + step. To prevent overconfident observations the ONEflux method attempts to account for as many sources of uncertainty as possible + and that is reflected in this converter code. + + If an observation has a missing uncertainty value, the code estimates an uncertainty based on an empirically-based relative uncertainty + value that reduces with increasing time resolution. The default relative uncertainty values are 20%, 10% and 5% for HH/HR, DD/WW and MM + time resolution respectively. The user can adjust these values within the source code. + + + +Quality Control +--------------- +The general QC naming convention uses an integer system (0-3) defined as the following: + + ``0 = measured``, ``1 = good quality gapfill``, ``2 = medium quality gapfill``, ``3 = poor quality gapfill``. (Refer to Pastorello et al., (2020) or + Reichstein et al. 2005 Global Change Biology for more information) + + +The QC values **do not** follow this convention for NEE, SH and LE fluxes for time resolutions coarser than the native resolution of +HH and HR. For DD, WW, and MM time resolution observations, the QC value is based on a fraction from 0-1 that indicates the fraction of the time period +that consists of measured or good quality gap-filled data. Because it is more straightforward to reject observations in DART +based on an integer value scale, the conversion code converts these fractional QC values (0-1) to integer QC values (0-3) as follows: + +#. ``QC(integer)=1 when QC(fraction) > 0.90``; +#. ``QC(integer)=2 when 0.90 >= QC(fraction) >= 0.60``; +#. ``QC(integer)=3 when 0.60 > QC(fraction) >=0``. + +.. Note:: + + The fraction QC to integer QC conversion approach was based on a qualitative assessment of flux data from Harvard Forest. + Depending on location and topography not all flux tower data will have a comparably high % of gap filled data.. The user + can change these thresholds within the source code. + +There are times when a QC value is missing or does not exist for an observation. In these cases the converter code does the following: + +#. Missing QC values (-9999) where the associated flux observation looks physically reasonable are assigned a QC = 2. +#. GPP and ER observations do not have a QC. If the observations are physically reasonable a QC = 2 is assigned. +#. There are situations where +100 is added to an existing QC value such that the observations are purposely rejected + during the conversion (assuming maxgoodQC = < 100). These situations are: + + a) When gap_filled = .false., all observations that are not measured (QC = 0) + b) When GPP or RE give non-physical negative values. + c) If NEE QC is missing. This is rare. + + + +Gap-Filling +----------- +Gap-filled data is only available for the native resolution format (HH, HR) for flux observations of +NEE, SH and LE. For all other situations choosing gap-filled data is mandatory, because measurements of fluxes at +the native resolution are frequently violated that cause the eddy covariance method to fail. These include +situations where the friction velocity (Ustar) falls below a certain value that prevents adequate land-atmosphere mixing, or when +their is instrumentation failure. In these cases, gap-filling methods (essentially models) are required to calculate daily and coarser +time resolutions. Because gap-filled data is technically modeled data (e.g. Marginal Distribution Method (MDS) which relies on +met condtions physcially and temporally similar to missing data) a user may desire only real observations. + +In general, we recommend to turn gap-filled data to ``.true.`` during the conversion process, and then use the QC value as a way to discriminate +against lower quality observations during the ``filter`` step. + + + + +Energy-Balance +-------------- + +We provide the option to use LE and SH observations that have been corrected for energy-balance closure. +In these cases the correction is based on comparing the latent and sensible heat flux against other sources +of energy loss/gain including net incoming radiation and energy radiated through the ground. + +Data Policy +----------- +It is important to recognize the flux data providers who have made their research publically available to advance +scientific research. Please see the Ameriflux data policy `here `__ +and the FLUXNET 2015 data policy provided `here `__. + + +Programs +-------- + +The ``fluxnetfull_to_obs.f90`` file is the source for the main converter program. Look at the source code where it reads the +example data file. The example code reads each text line into a character buffer and then reads from that buffer to parse up the data items. + +To compile and test, go into the work subdirectory and run the ``quickbuild.sh`` script to build the converter and a +couple of general purpose utilities. ``advance_time`` helps with calendar and time computations, and the +``obs_sequence_tool`` manipulates DART observation files once they have been created. + +To change the observation types, look in the ``DART/obs_def`` directory. If you can find an obs_def_XXX_mod.f90 file +with an appropriate set of observation types, change the 'use' lines in the converter source to include those types. +Then add that filename in the ``input.nml`` namelist file to the &preprocess_nml namelist, the 'input_files' variable. +Multiple files can be listed. Then run quickbuild.sh again. It remakes the table of supported observation types before +trying to recompile the source code. + +An example script for converting batches of files is in the ``shell_scripts`` directory. A tiny example data file is in +the ``data`` directory. These are *NOT* intended to be turnkey scripts; they will certainly need to be customized for +your use. There are comments at the top of the script saying what options they include, and should be commented enough +to indicate where changes will be likely to need to be made. + +Decisions you might need to make +-------------------------------- + +See the discussion in the :doc:`../../../guide/creating-obs-seq-real` page about what options are available +for the things you need to specify. These include setting a time, specifying an expected error, setting a location, and +an observation type. diff --git a/observations/obs_converters/Ameriflux/level4_to_obs.rst b/observations/obs_converters/Ameriflux/level4_to_obs.rst index 4b2c8b2663..b0f87d3d4a 100644 --- a/observations/obs_converters/Ameriflux/level4_to_obs.rst +++ b/observations/obs_converters/Ameriflux/level4_to_obs.rst @@ -7,12 +7,14 @@ Overview AmeriFlux level 4 data to DART observation sequence converter ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -| This routine is designed to convert the flux tower Level 4 data from the `AmeriFlux `__ - network of observations from micrometeorological tower sites. AmeriFlux is part of - `FLUXNET `__ and the converter is hoped to be a suitable starting point for the conversion of - observations from FLUXNET. As of May 2012, I have not yet tried to work with any other observations from FLUXNET. -| The AmeriFlux Level 4 products are recorded using the local time. DART observation sequence files use GMT. For more - information about AmeriFlux data products, go to http://ameriflux.lbl.gov. +| This routine was designed to convert the flux tower Level 4 data from the `AmeriFlux `__ + network of observations from micrometeorological tower sites. The download link and flux data format for this code + *has been discontinued* (e.g. ``USBar2004_L4_h.txt``). Thus if you are using flux obs converters for the first time + *PLEASE USE* the updated ``Fluxnetfull_to_obs.f90`` converter and follow the documentation there :doc:`./Fluxnetfull_to_obs` + We have kept this code available if you still use the older Ameriflux data format. Also this code uses a general approach + to calculating sensible, latent and net ecosystem exchange uncertainty, that may be helpful. +| The AmeriFlux Level 4 data products are provided in the local time of the flux tower location. DART observation sequence + files are provided in UTC, thus this routine includes a time conversion. .. warning:: diff --git a/observations/obs_converters/Ameriflux/work/input.nml b/observations/obs_converters/Ameriflux/work/input.nml index fe965b1083..2ad5583cce 100644 --- a/observations/obs_converters/Ameriflux/work/input.nml +++ b/observations/obs_converters/Ameriflux/work/input.nml @@ -27,7 +27,7 @@ text_input_file = '../data/USHa12003_L4_h.txt' obs_out_file = 'obs_seq.out' year = 2003 - timezoneoffset = -6 + timezoneoffset = -5 latitude = 42.5378 longitude = -72.1715 elevation = 353 @@ -36,6 +36,23 @@ verbose = .TRUE. / +&fluxnetfull_to_obs_nml + text_input_file = '../data/AMF_US-Ha1_FLUXNET_FULLSET_MM_1991-2020_3-5.csv' + obs_out_file = 'obs_seq.out' + timezoneoffset = -5 + latitude = 42.5378 + longitude = -72.1715 + elevation = 353 + flux_height = 29 + maxgoodqc = 3 + gap_filled = .true. + energy_balance = .false. + time_resolution = 'MM' + verbose = .true. + / + + + # This is appropriate for a days worth of flux tower observations # the obs in the file end 1 second before the time in the name. # If you are using these obs with CLM, ending 1 second before is appropriate. diff --git a/observations/obs_converters/Ameriflux/work/quickbuild.sh b/observations/obs_converters/Ameriflux/work/quickbuild.sh index 8ee2705cef..329028c0d8 100755 --- a/observations/obs_converters/Ameriflux/work/quickbuild.sh +++ b/observations/obs_converters/Ameriflux/work/quickbuild.sh @@ -14,6 +14,7 @@ LOCATION=threed_sphere programs=( +fluxnetfull_to_obs level4_to_obs obs_sequence_tool advance_time