Adding history variables for use with Newton-Krylov #1457
Conversation
| @@ -708,6 +727,7 @@ subroutine SetValues ( this, num_column, filter_column, value_column) | |||
| i = filter_column(fi) | |||
| this%decomp_cpools_transport_tendency_col(i,j,k) = value_column | |||
| this%decomp_cpools_sourcesink_col(i,j,k) = value_column | |||
| this%decomp_k_col(i,j,k) = value_column | |||
There was a problem hiding this comment.
I moved decomp_k from the "transitions" loop to the "pools" loop as part of this variable's corrected 3rd dimension.
| integer, private :: i_s2s3 | ||
| integer, private :: i_s3s1 | ||
| integer, private :: i_cwdl2 | ||
| integer, private :: i_cwdl3 |
There was a problem hiding this comment.
I made these variables available to this entire module to accommodate my changes to pathfrac_decomp_cascade. Similarly in deprecated SoilBiogeochemDecompCascadeCNMod.F90.
There was a problem hiding this comment.
are the parallel changes needed in CNMod?
There was a problem hiding this comment.
I made the changes in CNMod to avoid errors when building the code. Otherwise no.
| if (.not. use_fates) then | ||
| pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = cwd_fcel | ||
| pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = cwd_flig | ||
| end if |
There was a problem hiding this comment.
I now calc. pathfrac_decomp_cascade right after decomp_k in subr. decomp_rate_constants_bgc (and in deprecated subr. decomp_rate_constants_cn) because otherwise the order of initializations did not permit adding pathfrac_decomp_cascade to history.
In particular, leaving pathfrac_decomp_cascade in subr. init_decompcascade_bgc (i.e. unchanged):
- as a
soilbiogeochem_state_instvariable (i.e. unchanged) and adding to history resulted in messed up bounds forcascade_donor_poolin SoilBiogeochemStateType.F90 - as a
soilbiogeochem_carbonflux_instvariable and adding to history resulted in messed up bounds forpathfrac_decomp_cascadein subr. init_decompcascade_bgc
I'm happy to explain in greater detail if requested.
There was a problem hiding this comment.
I'll take your word on this :)
|
Cheyenne test-suite OK |
| if (.not. use_fates) then | ||
| pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = cwd_fcel | ||
| pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = cwd_flig | ||
| end if |
There was a problem hiding this comment.
I'll take your word on this :)
| @@ -39,6 +39,18 @@ module SoilBiogeochemDecompCascadeCNMod | |||
| integer, private :: i_soil3 = -9 ! SOM third pool | |||
| integer, private :: i_soil4 = -9 ! SOM fourth pool | |||
|
|
|||
| real(r8), private :: cwd_fcel | |||
There was a problem hiding this comment.
we don't need this in CNMod, since the CN model is being deprecated.
There was a problem hiding this comment.
True, but I think that the case would have stopped in the build phase, so I made the changes anyway.
| @@ -110,7 +110,7 @@ subroutine SoilBiogeochemPotential (bounds, num_soilc, filter_soilc, & | |||
| initial_cn_ratio => decomp_cascade_con%initial_cn_ratio , & ! Input: [real(r8) (:) ] c:n ratio for initialization of pools | |||
|
|
|||
| rf_decomp_cascade => soilbiogeochem_state_inst%rf_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] respired fraction in decomposition step (frac) | |||
| pathfrac_decomp_cascade => soilbiogeochem_state_inst%pathfrac_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] what fraction of C leaving a given pool passes through a given transition (frac) | |||
| pathfrac_decomp_cascade => soilbiogeochem_carbonflux_inst%pathfrac_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] how much C from the donor pool passes through a given transition (frac) | |||
There was a problem hiding this comment.
detail should still be what fraction of C, as this is a fraction, not a flux (which seems implied by 'how much C').
There was a problem hiding this comment.
Changing to
what fraction of C passes from donor to receiver pool through a given transition
Does this seem ok?
| @@ -34,7 +34,6 @@ module SoilBiogeochemStateType | |||
| real(r8) , pointer :: fpi_col (:) ! (no units) fraction of potential immobilization | |||
| real(r8), pointer :: fpg_col (:) ! (no units) fraction of potential gpp | |||
| real(r8) , pointer :: rf_decomp_cascade_col (:,:,:) ! (frac) respired fraction in decomposition step | |||
| real(r8) , pointer :: pathfrac_decomp_cascade_col (:,:,:) ! (frac) what fraction of C leaving a given pool passes through a given transition | |||
There was a problem hiding this comment.
why is this being removed from this module?
Seeing this code it also makes me wonder if we need to:
- write out
fpi_col(N limitation of decomp), or if it's brought intodecomp_k? - write out
rf_decomp_cascade_col, or is this brought intopathfracsomewhere?
There was a problem hiding this comment.
pathfrac got moved from here to SoilBiogeochemCarbonFluxType for reasons outlined in this comment
FPI is on the master list as an inactive field.
rf_decomp_cascade is not. Let me know if I should add it.
There was a problem hiding this comment.
Let's go ahead and write out fpi and rf_decomp_cascade to history. I think we'll want these for the precondiotnal
doc/ChangeLog
Outdated
| --------------------------------------- | ||
| Substantial timing or memory changes: | ||
| [For timing changes, you should at least check the PFS test in the test | ||
|
|
There was a problem hiding this comment.
are these history fields going to be written out for each column or grid cell average? I'd assume the former, but wondered if N-K solver can handle column level history files @klindsay28?
@ekluzek I also wondered how we easily turn on this output for a N-K spinup? Should this feature come into this PR, or just handled more manually at this stage?
There was a problem hiding this comment.
Regarding the second question:
One can activate so-called "inactive" history fields using the namelist. In this testing phase I would rely on that. If after the testing we decide to go with Newton-Krylov, then we can create a corresponding compset. Would you respond differently, @ekluzek ?
There was a problem hiding this comment.
are these history fields going to be written out for each column or grid cell average? I'd assume the former, but wondered if N-K solver can handle column level history files @klindsay28?
When added to the primary history file, history fields appear as grid cell average by default. With appropriate namelist selections, they can also be sent to history as column variables.
There was a problem hiding this comment.
@klindsay28, grid cell average variables are dimensioned by lat, lon, while column variables are written as a column vector. The latter complicates a bit the task of finding where you are on the planet, which I think is why @wwieder asked the question.
There was a problem hiding this comment.
The solver will be reading and modifying restart files. Are these column based? (I don't have the CLM/CTSM diagram of grid cell/PFT/column handy.) If so, then the solver will be dealing with column level variables, and column level history output then seems natural. If there is history output needed for the preconditioner that is cell averages, then I think there will be an extra step of mapping the column index to the proper grid cell.
| integer, private :: i_s2s3 | ||
| integer, private :: i_s3s1 | ||
| integer, private :: i_cwdl2 | ||
| integer, private :: i_cwdl3 |
There was a problem hiding this comment.
are the parallel changes needed in CNMod?
|
Izumi test-suite OK |
There are more of these in store; however, I need to merge my work from ESCOMP#1457 while ESCOMP#1457 is pending, in order to make progress here
|
Running tests that compare against dev053: I will update this post as tests complete. |
|
Thanks for these changes @slevisconsulting . I just wanted to confirm that in this PR are the additional variables still set to inactive by default with the plan of manually turning them on as we're testing the N-K output? |
@wwieder that is correct. The last test listed a few posts up will confirm that namelist variable |
|
Notes from today's mtg:
|
|
@slevisconsulting just add a subdirectory to cime_config/usermods_dirs/ with a name something like newton_krylov_spinup. You can add the user_nl_clm there, as well as a shell_commands file if you need to do any xml settings (you should set it to a cold start, so you'll want that). |
| @@ -0,0 +1,15 @@ | |||
| hist_dov2xy = .true.,.false. | |||
| hist_nhtfrq = 0,-175200 | |||
There was a problem hiding this comment.
@wwieder the above two lines set the auxiliary history file to 20-yr average vector output.
nhtfrq is in hours: 24 hrs/day * 365 days/yr * 20 yrs = 175200
Two lines down is the list of variables, and I think I got all of them.
|
Thanks @slevisconsulting dumb question here, but what do I need to do when creating a case to get these user_mods to be active? Or will they just be included by default? |
|
@wwieder this is done similar to the user-mods used for NEON. It's an option to create_newcase. Use ./create_newcase --help to remind yourself of the options. But, it would be using the "--user-mods-dirs newton_krylov_spinup" option. |
|
OK, questions here are moving outside of this issue / PR, but what's the most effective way to create a case for a random point and populate the datm.streams file correctly? you can see my initial attempts here /glade/work/wwieder/ctsm/ctsm5.1_N-K_test/cime/scripts/RandomBoreal_0 |
I didn't look at your case directory, but the quickest/simplest approach that I think should work is to point to the existing global atm files rather than generate new ones for this boreal point. The model will run slower due to interpolation calculations, but you won't have to generate a whole new atmosphere dataset. But @ekluzek may have another suggestion. |
|
sorry to pester here, @ekluzek, but can you help make suggestions for how to set up a generic single point (58N, 55E) with less headache than I did here /glade/work/wwieder/ctsm/ctsm5.1_N-K_test/cime/scripts/RandomBoreal_0 |
|
Hey @wwieder OK so I looked at your case, and also created a case to model from. Two issues in your case is that you need to set CLM_USRDAT_NAME to something that refers to the site you are working on. Using your nomenclature I used "RANDOMBOREAL1". The other thing is the compset you are using which is configured to get tower site data (that's what the 1Pt refers to in the name). So you want a regular compset name so that you get the global forcing data. You also set the site location with PTS_LAT and PTS_LON xml variables. So my example case is here: /glade/work/erik/ctsm_worktrees/branch1/cime/scripts/cases/RandomBoreal_0 With createnew case command:
I did these xml commands (added them to shell_commands) NOTE: CASESTR is not required, but makes sure a few things are set with something besides UNSET Added this to user_nl_clm fsurdat = '/glade/scratch/wwieder/single_point/surfdata_hist_16pfts_Irrig_CMIP6_simyr2000_RandomBoreal_c210823.nc' and this to user_nl_mosart frivinp_rtm = '/dev/null' (This is something that is required with latest updates in externals, when MOSART_MODE=null (i.e. when compset has MOSART rather than SROF). This is a little glitch that we should fix in CTSM so you don't have to do this (There's a MOSART issue about this) |
|
Thanks Erik. After submitting, now I'm getting this: Do you have any suggestions here? looks like this is true in your case too /glade/work/erik/ctsm_worktrees/branch1/cime/scripts/cases/RandomBoreal_0/CaseDocs/datm_in |
|
Hi @wwieder. Hmmm. Actually there's a simple step that I left out, that you need to do. My case did run even with that being UNSET. I think that's actually normal. Although it probably means we should clean this up and give it some value like "UNUSED" to make it more clear that it's not needed. From comparing my case to yours, the differences I see are: PTS_LON,PTS_LAT, and CLM_USRDAT_NAME aren't set in your case, which are the critical things to make sure are set correctly. Oh, and I see that you added these settings to your ./shell_commands -- but they weren't invoked by the case. So that just means you need to run ./shell_commands manually. That's an important step that I neglected to say. The case does run shell_commands, itself, but it must be either before case.setup or with case.setup.... ./shell_commands |
|
Thanks @ekluzek that's what I needed to do! We'll have to work on some new documentation for this somewhere. I'll discuss with @danicalombardozzi where this should go. |
|
Last question, how should can I share the build from the first case with another nearly identical case? |
I haven't used this in a long while, but try: |
Happy to discuss -- just let me know when! |
|
@wwieder the important argument to create_clone is "--keepexe" which will then shared the built case exe with the new case. There's some notes on it and help with ./create_clone --help |
|
the shell_commands that were created with this case didn't have all the needed information. I copied this from a different case and things seems to be running now. |
|
@klindsay28 you should be able to clone this case for the next step in your testing. @slevisconsulting and @ekluzek thanks for helping work this out. |
|
@ekluzek I should have been more clear. This is a single point case, but it's just reading in the global GSWP3 datm files. This is something we can work to better document in the updated nuopc workflow for running generic single point cases with global input data. |
|
Linking @ekluzek 's comment |
Correct perched water table calculation Modifies the calculation of the frost table and perched water table layers. This brings in the answer-changing aspects of the Hillslope Hydrology branch (ESCOMP#1715). Specific changes are: * PerchedWaterTable 1) Frost table depth a) original frost table determination looped from the top of the soil column downward to the index of the first layer whose temperature was <= freezing, and whose neighbor above had a temperature above freezing. The frost table depth is then given by the node of that soil layer, i.e. z(k_frz). b) in the new method, the same index is found, but the depth of the frost table is given by the depth of the top of the frozen layer, i.e. zi(k_frz-1). Note zi(k_frz) would be the bottom of layer k_frz. 2) Perched water table depth a) in the original formulation, a loop from k_frz to layer 1 was used to identify the deepest layer in 1:k_frz whose volumetric soil moisture was greater than a threshold given by sat_lev (e.g. sat_lev = 0.9). b) in the new method, the search is only done if k_frz is greater than 1. The rationale is that if k_frz = 1, then zwt_perched has already been initialized to the frost_table depth (which is equal to the top of the uppermost soil layer), and therefore no search is required. 3) Determining perched water table depth within layer identified by index k_perch in 2) a) in the original formulation, the perched water table depth was calculated by linearly interpolating between layers k_perch and k_perch+1, with no consideration of their relative values. In the case where the deeper layer was drier than the layer above it, this could result in values far outside the soil layer. b) in the new formulation, if the deeper layer is drier than the layer above it (s1 > s2), then the perched water table depth is simply given by the depth of the upper surface of layer k_perch, i.e. zi(c,k_perch-1). * PerchedLateralFlow 1) Removal of icefrac calculation a) in the original calculation, the frozen layer was included, so an ice impedance factor was calculated using icefrac. If only unfrozen layers are used, no ice impedance factor is needed. b) in the new formulation, the icefrac variable and loop are removed. 2) Move loop calculating frost and perched water table depths a) in the original formulation, the frost and perched water table depths were calculated in the same loop as the calculation of the lateral flow from the perched saturated zone. b) in preparation for the hillslope hydrology branch, this calculation is moved into its own loop. 3) q_perch calculation a)in the original formulation, q_perch was calculated by summing over layers k_perch to k_frost. However, because k_frost is now identified as the frozen layer, and its depth the top of the frozen layer, it should not be included in the calculation; only the unfrozen layers above it should be included. b) in the new formulation, the loop is bounded by k_frost-1 instead of k_frost. 4) Removal of water from perched saturated zone a) the in the original method, the frozen layer (k_frost) was included in the loop. Also, the drainage was defined to be negative, which was confusing. b) in the new method, the frozen layer is not included in the loop; water is only removed from the unfrozen layers above k_frost. Calculate drainage as positive values, which are then subtracted from the soil moisture in each layer. Resolved conflicts: doc/ChangeLog doc/ChangeSum
|
I'm running a copy of this case For the same purpose I am also running a global 10x15 case. |
Reuse some files generated in initialization when rerunning A new directory is now in your run directory init_generated_files. These two run time generated files - finidat_interp_dest.nc and the land fraction file ctsm_landfrac.nc (which is new in this PR) - are put in that directory when they are created. Then, if they already exist, they are reused rather than being regenerated. (Technically, it looks for a status file that flags successful creation of this file.) This applies when you rerun a startup (non-restart) run a second or subsequent time, as is common when doing development and testing. For tests, theinit_generated_files directory is removed before running, so that rerunning a test will do the same thing the second time. There are occasional, rare cases when a user would need to manually remove the init_generated_files directory – e.g., when changing things about N dominant landunits or PFTs. Also, output was cleaned up in a few cases so that it did not span multiple lines. Resolved conflicts: doc/ChangeLog doc/ChangeSum
|
@slevis-lmwg is this something that was close enough that it come in as it is? Or is there more work to do on it? We kind of abandoned it a while back, so if it takes too long to get spun up on it again, don't worry about it. But, if it's close enough that it could come in as is, it would be worth bringing it in. |
|
@wwieder also has a say here. |
|
I think this suggestion makes sense, @slevis-lmwg hopefully it's not hard to bring this draft up to main if the N-K is effectively implimented. |
|
Sounds good. We'll leave it as it is... |
samsrabin
added
good first issue
|
Closing as won't fix. Anderson Acceleration is the currently planned method for mimics spinups. |
Description of changes
decomp_k already had the necessary infrastructure to appear in ctsm history. It's an "inactive" variable and appears in the documentation's master list as K_*. It was dimensioned wrong in the code, so I corrected that in my first commit to this PR, thankfully without changing answers.
Next I will add
pathfrac_decomp_cascadeto history.Specific notes
Contributors other than yourself, if any:
@wwieder @ekluzek
CTSM Issues Fixed (include github issue #):
#1455 #1825
Are answers expected to change (and if so in what way)?
No. Variable
pathfrac_decomp_cascadewill appear in history optionally.Any User Interface Changes (namelist or namelist defaults changes)?
No. Later we may consider interface changes to add these two inactive variables to history when spinning up with Newton-Krylov.
Testing performed, if any:
The following test PASSES with the mods of the first commit to this PR:
ERP_P36x2_D_Ld3.f10_f10_mg37.I1850Clm50BgcCrop.cheyenne_gnu.clm-extra_outputs -c /glade/p/cgd/tss/ctsm_baselines/ctsm5.1.dev051