Vertically resolved C variables dimensioned levsoi, while N variables dimensioned levdcmp

[slevis-lmwg]

Brief summary of bug

The inconsistency described in the title should be corrected by using the levdcmp dimension for both C and N variables. Additionally we should make levdcmp = levsoi (does not include bedrock) rather than levgrnd (includes bedrock layers). For consistency we should use levdcmp rather than levgrnd when writing these variables to restart.

Extra credit: We may wish to list all the "lev" dimensions and their uses to assess whether we may consolidate some of them.

General bug information

CTSM version you are using: ctsm5.1.dev106-7-ga5b74c35b

Does this bug cause significantly incorrect results in the model's science? No

Configurations affected: bgc

Details of bug

Although this doesn't cause errors, it complicates the use of history and restart files, for example in the case of the Newton-Krylov spin-up tools. For this reason, I will open a PR to correct.

[slevis-lmwg]

@klindsay28 @wwieder Pls comment on the issue if I missed anything...

[wwieder]

my only additional questions for @ekluzek and others are:

1. How ambitious we want to clean this up throughout in the code; and
2. What it may mean for backwards compatibility?

@slevisconsulting would it be helpful to provide any examples of how/where this is done in the code or in restart and history files to facilitate the conversation, especially for anyone reading this cold?

[slevis-lmwg]

@slevisconsulting would it be helpful to provide any examples of how/where this is done in the code or in restart and history files to facilitate the conversation, especially for anyone reading this cold?

I have confirmed that dimension levdcmp corresponds to nlevdecomp_full (i.e. including bgc-inactive layers), so I'm changing it to correspond to nlevdecomp (i.e. only bgc-active layers), and we'll see if I run into trouble with that change alone.

Example from history (carbon variables using levsoi versus nitrogen using levdcmp): SoilBiogeochemCarbonStateType.F90

          if ( nlevdecomp_full > 1 ) then
             data2dptr => this%decomp_cpools_vr_col(:,1:nlevsoi,l)
             fieldname = trim(decomp_cascade_con%decomp_pool_name_history(l))//'C_vr'
             longname =  trim(decomp_cascade_con%decomp_pool_name_history(l))//' C (vertically resolved)'
             call hist_addfld2d (fname=fieldname, units='gC/m^3',  type2d='levsoi', &
                  avgflag='A', long_name=longname, &
                  ptr_col=data2dptr)
          endif

SoilBiogeochemNitrogenStateType.F90

       if ( nlevdecomp_full > 1 ) then
          data2dptr => this%decomp_npools_vr_col(:,:,l)
          fieldname = trim(decomp_cascade_con%decomp_pool_name_history(l))//'N_vr'
          longname =  trim(decomp_cascade_con%decomp_pool_name_history(l))//' N (vertically resolved)'
          call hist_addfld2d (fname=fieldname, units='gN/m^3',  type2d='levdcmp', &
               avgflag='A', long_name=longname, &
               ptr_col=data2dptr)
       endif

...and from restart (same modules, now both using dim2name='levgrnd'):

          varname=trim(decomp_cascade_con%decomp_pool_name_restart(k))//'c'
          ptr2d => this%decomp_cpools_vr_col(:,:,k)
          call restartvar(ncid=ncid, flag=flag, varname=trim(varname)//"_vr", xtype=ncd_double,  &
               dim1name='column', dim2name='levgrnd', switchdim=.true., &
               long_name='',  units='g/m3', fill_value=spval, &
               scale_by_thickness=.false., &
               interpinic_flag='interp', readvar=readvar, data=ptr2d)

       varname=trim(decomp_cascade_con%decomp_pool_name_restart(k))//'n'
       ptr2d => this%decomp_npools_vr_col(:,:,k)
       call restartvar(ncid=ncid, flag=flag, varname=trim(varname)//"_vr", xtype=ncd_double, &
            dim1name='column', dim2name='levgrnd', switchdim=.true., &
            long_name='', units='gN/m3', &
            scale_by_thickness=.false., &
            interpinic_flag='interp', readvar=readvar, data=ptr2d)

[slevis-lmwg]

Test case seems to work correctly, though I will run the full test suite to confirm.

[slevis-lmwg]

I have left the restart part of this unchanged at this time, because the model failed while writing to restart when I tried changing levgrnd to levdcmp.

So @klindsay28 I have made the C and N variables consistent in history by dimensioning them all with levdcmp (= 20), but the same variables are dimensioned with levgrnd (= 25) in restart. If this is a problem, the easier approach may be to dimension the variables in history same as in restarts, since the opposite didn't work for me. But I can also look deeper into why the opposite didn't work.

[slevis-lmwg]

I have left the restart part of this unchanged at this time, because the model failed while writing to restart when I tried changing levgrnd to levdcmp.

[...] the easier approach may be to dimension the variables in history same as in restarts, since the opposite didn't work for me.

I have gone ahead with the "easier" approach, since our immediate motivation here is to get the Newton-Krylov spin-up tool working for the ctsm. This means that the C and N variables are now consistent in history by all being dimensioned with levdcmp (=25) and the same variables are still dimensioned with levgrnd (=25) in restart.

Cheyenne test suite PASS

[slevis-lmwg]

Pushing my code changes to existing PR #1457 and the corresponding branch.

[billsacks]

I have only skimmed through the comments here, but wanted to note that, if significant changes are made to the vertical dimensions of restart variables, some thought may need to be given to the multi-level init_interp code (initInterpMultilevelContainer.F90 is a good starting point).