-
Notifications
You must be signed in to change notification settings - Fork 143
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Fluxnet observation converter #478
Conversation
Creating new observation converter for site level eddy covariance flux tower data (NEE,GPP,RECO,LH,SH). The existing Ameriflux_L4 converter is outdated and new database has new flux options and time aggregates
Necessary changes in syntax to allow fluxnetfull_to_obs.f90 to compile correctly. Also providing proper changes to support files input.nml and quickbuild.sh within work directory
fluxnetfull_to_obs converter compiles and runs testing on hourly flux fullset data. Header now includes timestamp which is treated differently than real values
@mjs2369 hoping you can take the lead on the review once I have pulled this PR off draft. The documentation still needs to be worked on a bunch by me, but I got the main converter (Fluxnetfull_to_obs) working today. It only functions on the hourly resolution flux data format, and I need to add the capability to convert time averaged observations, and also test it. This lays the foundation though. |
@braczka I can definitely do that. I'll keep an eye on it to review when the PR is no longer a draft. |
Adds ability to use HH,HR,DD,WW,MM,YY flux data where test file format changes
Expanded functionality such that fluxnet converter can take in HH,HR,DD,WW,MM and YY time averages
Added more fixes and warnings to make code compatiable with coarser time resolution. File formatting changes based on time resolution.
QC values change from integers to fraction depending upon time resolution. Added code to convert fraction back to integer. Provided empirical fix when uncertainty values for H,LE and NEE are missing.
The code has been tested for all relevant time resolutions: hourly, daily, weekly, and monthly. The documentation has also been completed. Taking off the draft status -- ready for review. |
Also a quicker link to get access to test data rather than the flux repository: /glade/scratch/bmraczka/Ameriflux_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. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Lines 61 and 62 mean that for the coarse time resolution (DD, WW, MM), gap_filled for NEE, latent, sensible heat, GPP and RE is true, and for the native time resolution (HH, HR), gap_filled for NEE, latent, sensible heat, GPP and RE can be true or false, correct?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, that is correct. If namelist variable gap_filled = .false.
, and time_resolution
is set to any coarse resolution like DD
, WW
or MM
, the code will error out. See lines 228-235.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I see.
rewind(iunit) | ||
call decode_header(iunit, nwords) | ||
|
||
obsloop: do iline = 2,nlines |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Could you please point me to the obs_seq.out that you generated in the test?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Sure -- I put them here: /glade/scratch/bmraczka/Ameriflux_data/
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%leQC = 2 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
tower%hQC = 2
tower%leQC = 2 | ||
tower%leQC = 2 | ||
endif | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The leQC and hQC are assigned 2 in the if-else statement from lines 1139 to 1145 if the the observation le and hQC look Okay. I am wondering why there are no statements for leQC and hQC like lines 1155-1158 for gppQC and recoQC? Is it because in le and h are never negative in the downloaded data?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Both le (latent heat) and h (sensible heat) can be positive or negative (physically reasonable), thus I don't use this as a diagnostic to adjust the QC values. On the other hand, negative gpp and reco values are physically impossible, thus I force the code to reject them by making their respective QC values big.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Makes sense.
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%recoDTQC,r8),real(tower%recoNTQC,r8)/)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
gppDTQC, gppNTQC. Or it doesn't matter since recoDTQC is always the same as gppDTQC, and recoNTQC is the same as gppNTQC?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Typo -- very good catch.
tower%gppNTQC = maxgoodqc + 100 | ||
tower%gppDTQC = maxgoodqc + 100 | ||
tower%recoNTQC = maxgoodqc + 100 | ||
tower%recoDTQC = maxgoodqc + 100 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
One last comment on lines 1168-1171. I understand the design logic you have. Say if I do would like to use the valid gpp and reco when the gap_filled is false at the HH or HR time resolution. Since the QC assigned to them in lines 1168-1171 are the same values as QC assigned when gpp, reco are negative. Do you think it is worthwhile assigning a different QC in lines 1168-1171 to distinguish the two situations so that I can pick a QC threshold to allow the valid gpp and reco when gap_filled is false at the HH or HR time resolution assimilated ?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If you want to use HH/HR GPP and RECO values, just set gap_filled=.true. Then all of the GPP and RECO values will be sent to your obs_seq file. By default they are assigned a QC=2 in the obs_seq file.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Also, to be clear, when gap_filled=.true. both gap_filled and 'true' observations are sent to the obs_seq file. 'True' observations are never excluded from the obs_seq file.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
True.
Adding better documentation on how user can adjust values that convert fractional QC to integer QC and also empirically derived variables related to uncertainty estimation for missing values
@XueliHuo Absolutely no rush on this review, but I wanted to follow up to see if there was anything else I needed to do (through documentation, code changes) to address any of your concerns? The last few commits in this PR hopefully fix inaccuracies in the documentation that you identified. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hi Brett,
Since the science review seems to be wrapping up, a few comments:
lower case for namelist and program names. (I put suggestions in for this, but ran out of steam at ~20. Just do a global find and replace). Fortran is case insensitive, various filesystems may or may not be, people may be using case-sensitive tools to generate scripts and namelists and chaos ensues. So lower case only for namelists program names.
add fluxnetfull_to_obs to .gitignore so people don't accidentally commit the executable to the repo.
You've got some duplicate code, from the obs_utilities_mod. Rather than duplicate these, just use the obs_utilities_mod versions:
use obs_utilities_mod, only: add_obs_to_seq, create_3d_obs
Cheers,
Helen
…f-rttov13 Added get_channel to obs_def_rttov13_mod.f90
use types_mod, only : r8, MISSING_R8 | ||
|
||
use utilities_mod, only : initialize_utilities, finalize_utilities, & | ||
register_module, error_handler, E_MSG, E_ERR, & |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
register_module, error_handler, E_MSG, E_ERR, & | |
error_handler, E_MSG, E_ERR, & |
|
||
use location_mod, only : VERTISHEIGHT | ||
|
||
use obs_sequence_mod, only : obs_sequence_type, obs_type, read_obs_seq, & |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
use obs_sequence_mod, only : obs_sequence_type, obs_type, read_obs_seq, & | |
use obs_sequence_mod, only : obs_sequence_type, obs_type, & |
! 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') |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@braczka this call to open_file is missing a corresponding call to close_file
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Gotcha -- took care of these in latest commits.
Hi Brett, I went through the questions I posted and the response/modifications you provided/made. I approve this pull request. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Switching my review to approve so I don't block anything.
Brett addressed my duplicate code and capital F comments.
@mjs2369 you're in charge of the release!
This is good on my end, @mjs2369 release at your convenience. Thanks all. |
bug fix: remove cloud_overlap integer from list of logicals to query
Description:
Generates a new observation converter (
Fluxnetfull_to_obs
) for eddy covariance flux tower data (carbon, water energy fluxes) which is useful for land DA and has been used in the past for CLM-DART. This converter uses the new format provided by Fluxnet/Ameriflux, and also includes new features such as GPP and RECO carbon flux observations, different time aggregations (hourly through monthly), and provides estimates of random, ustar and partitioning uncertainty. As part of this PR, documentation changes will be made to the older, deprecated ameriflux converter (level4_to_obs
) and the broken links will be fixed. New flux tower observation types will be added to accomodate the forward operator approach for time aggregated fluxes (daily through monthly).Also includes pull requests #540 (Added get_channel to obs_def_rttov13_mod.f90) and #541 (bug fix: remove cloud_overlap integer from list of logicals to query)
Fixes issue
This addresses github issue #339
Location and format of L4 Ameriflux Observations have changed
Also #537, #532
Types of changes
Documentation changes needed?
Tests
I have compiled and run a test case based on the native hourly resolution flux data from Harvard Forest (AMF_US-Ha1_FLUXNET_FULLSET_HR_1991-2020_3-5.csv). This data can be downloaded from here:
(https://ameriflux.lbl.gov/data/download-data/)
Checklist for merging
Checklist for release
Testing Datasets