Logic for initializing C isotopes and maybe other variables incorrect when using init_interp

[billsacks]

Bill Sacks < sacks@ucar.edu > - 2017-07-26 06:34:15 -0600
Bugzilla Id: 2496
Bugzilla CC: cdkoven@lbl.gov, dlawren@ucar.edu, erik@ucar.edu, mvertens@ucar.edu, rfisher@ucar.edu,

Ben pointed this problem out with respect to water isotopes, but I realized that the same problem occurs in existing code for carbon isotopes and possibly other variables: For a number of variables, there is code in the Restart routine that says: If the variable isn't found on the restart file, then initialize it based on some other variable that is guaranteed to be on the restart file - e.g., initialize c13 to the bulk C value times some constant. My understanding is that this code is important when you're transitioning from a run without some option (e.g., without ciso) to a run with that option turned on.

This code works fine if use_init_interp = .false., but I'm pretty sure it won't work with use_init_interp = .true.: Imagine you're doing a run with ciso true, pointing to an finidat file that was generated with ciso false, with use_init_interp true. In this case, finidat_interp_dest.nc is written with all of the C isotope fields present, at their cold start values. In init_interp, these fields will be skipped because they're not found on the source finidat file - and thus they will remain present in the finidat_interp_dest file, with values at the cold start values. Then, back in the CLM Restart routines, the code that checks whether the variable is present on the restart file will say, "Yup, it's present", and so the code for newly-initializing the C isotope fields will never be executed.

Ben and I discussed a number of solutions yesterday, none of them completely satisfying. I won't lay those all out here. But I will lay out one more variation that I just thought of, which has the advantage of NOT requiring finding and changing all of the existing code that has patterns like the C isotope code:

When we write a restart file, add a new attribute to every field like "valid_data" (though we may want to come up with a better name). In the restartvar routine, that would be set to 1 (true) everywhere. Then, in initInterp, in the conditional that checks whether a given output variable is absent from the input file (the lines of code around the message, "variable is NOT on input file"): Set this valid_data attribute to 0 (false). (Aside: it's important that the new attribute value take up no more space than the old attribute value, according to what's allowed in netcdf when you have left 'define' mode. So, for example, it would NOT work to set it to 'true' originally and 'false' in init_interp.) Then, back in the restartvar routine, on read: it would check the value of this valid_data attribute; if 0, then it skips reading the data and sets readvar to .false.

Ben, note that, unlike the solutions we were brainstorming yesterday, this solution does NOT involve identifying the relevant variables and giving them a new interpinic_flag behavior; instead, this solution would apply to ALL variables. I believe it's safe to apply this to all variables. i.e., for these variables that aren't present on the finidat source file, the end result should be the same whether or not we read the data in the Restart routine: they should end up at their cold start values in either case. But if we go with this solution, we may want to give this a little more thought to make sure it's safe.

[bandre-ucar]

Misc notes I have on this issue from my water isotope branch:

Isotopes (water and carbon) need to be in sync with their bulk
counterparts. When reading from an initial conditions file that does
not contain isotopes, the bulk values will be updated, but the
isotpoes will remain at their cold start values, and the isotopic
ratios will be out of sync with the bulk.

Add metadata to the restart file indicating the data source as
restart, interpolated IC, or IC unavailable. The state and flux
instance variables can query the data source, and act appropriately
when reading a restart file.

Initial patch against clm4_5_16_r249 from the water isotopes branch.

[billsacks]

Here is the patch linked from the above gist, in case it gets deleted:

restart-meta-r249.diff.txt

[bandre-ucar]

It looks like the last work I did on this issue was on this svn branch from r253:

https://svn-ccsm-models.cgd.ucar.edu/clm2/branches/andre_restart_metadata/

I believe this has all the infrastructure changes needed to implement this and it was tested on the water isotopes branch. I had a set of incomplete changes on top of this branch applying the infrastructure to the carbon isotopes. I can not find that sandbox, but I believe it was stalled awaiting some changes based an email exchange with @billsacks on 2017-08-24:

I made changes to carbon isotopes related variables to change from:


if (flag=='read' .and. .not. readvar)

to 

if (flag=='read' .and. data_source%not_valid())

Where valid data is defined as:

    if (this%readvar) then
       if (this%source == data_source_restart .or. &
           this%source == data_source_interp) then
          is_valid = .true.
       end if
    end if


The remaining uses of "'read' and not readvar" is stuff like below. My best guess is that all these files and patterns are safe to change. Can you take a quick look and let me know if you see anything that should NOT be converted?

Thanks,

Ben


components/clm/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90
components/clm/src/soilbiogeochem/SoilBiogeochemNitrogenStateType.F90
components/clm/src/soilbiogeochem/SoilBiogeochemNitrogenFluxType.F90
         if (flag=='read' .and. .not. readvar) then
             call endrun(msg='ERROR:: '//trim(varname)//' is required on an initialization dataset'//&
                  errMsg(sourcefile, __LINE__))
          end if

components/clm/src/biogeophys/SoilStateType.F90
         if (flag=='read' .and. .not. readvar) then
            if (masterproc) then
               write(iulog,*) "can't find rootfr in restart (or initial) file..."
               write(iulog,*) "Initialize rootfr to default"
            end if
            call init_vegrootfr(bounds, nlevsoi, nlevgrnd, &
            this%rootfr_patch(bounds%begp:bounds%endp,1:nlevgrnd), 'water')
            call init_vegrootfr(bounds, nlevsoi, nlevgrnd, &
            this%crootfr_patch(bounds%begp:bounds%endp,1:nlevgrnd), 'carbon')
         end if


components/clm/src/biogeophys/IrrigationMod.F90
    if (flag=='read' .and. .not. readvar) then
       this%n_irrig_steps_left_patch = 0
    end if


components/clm/src/biogeophys/TemperatureType.F90
    if (flag=='read' .and. .not. readvar) then
       this%t_h2osfc_col(bounds%begc:bounds%endc) = 274.0_r8
    end if

    if (flag=='read' .and. .not. readvar) call endrun(msg=errMsg(sourcefile, __LINE__))

       if (flag=='read' .and. .not. readvar) then
          if (masterproc) write(iulog,*) "can't find t_building in initial file..."
          if (masterproc) write(iulog,*) "Initialize t_building to taf"
          this%t_building_lun(bounds%begl:bounds%endl) = this%taf_lun(bounds%begl:bounds%endl)
       end if

See some replies inline:


components/clm/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90
components/clm/src/soilbiogeochem/SoilBiogeochemNitrogenStateType.F90
components/clm/src/soilbiogeochem/SoilBiogeochemNitrogenFluxType.F90
         if (flag=='read' .and. .not. readvar) then
             call endrun(msg='ERROR:: '//trim(varname)//' is required on an initialization dataset'//&
                  errMsg(sourcefile, __LINE__))
          end if

At a glance, it looks right to convert these. It looks like at least some of these are ensuring that, if you do a run with use_nitrif_denitrif .true., then you cannot use a restart file that had use_nitrif_denitrif .false. I'm not sure if that was the intent, but changing these to the new form should maintain that check.


components/clm/src/biogeophys/SoilStateType.F90
         if (flag=='read' .and. .not. readvar) then
            if (masterproc) then
               write(iulog,*) "can't find rootfr in restart (or initial) file..."
               write(iulog,*) "Initialize rootfr to default"
            end if
            call init_vegrootfr(bounds, nlevsoi, nlevgrnd, &
            this%rootfr_patch(bounds%begp:bounds%endp,1:nlevgrnd), 'water')
            call init_vegrootfr(bounds, nlevsoi, nlevgrnd, &
            this%crootfr_patch(bounds%begp:bounds%endp,1:nlevgrnd), 'carbon')
         end if


This original code looks buggy for a number of reasons:

(1) These initializations were already done in SoilStateInitTimeConst, and so I don't think they should be needed here.

(2) The initializations in SoilStateInitTimeConst are within a use_fates conditional:

    ! Initialize root fraction 
    ! Note that fates has its own root fraction root fraction routine and should not
    ! use the following since it depends on patch%itype - which fates should not use

    if (.not. use_fates) then
        call init_vegrootfr(bounds, nlevsoi, nlevgrnd, &
             soilstate_inst%rootfr_patch(begp:endp,1:nlevgrnd),'water')
        call init_vegrootfr(bounds, nlevsoi, nlevgrnd, &
             soilstate_inst%crootfr_patch(begp:endp,1:nlevgrnd),'carbon')
     end if

whereas the ones here are not. Presumably there is some problem if those things actually end up getting called when use_fates is .true.

(3) Looking at some larger code context, I see:

    call restartvar(ncid=ncid, flag=flag, varname='HK', xtype=ncd_double,  &
         dim1name='column', dim2name='levgrnd', switchdim=.true., &
         long_name='hydraulic conductivity', units='mm/s', &
         interpinic_flag='interp', readvar=readvar, data=this%hk_l_col)

     if( use_dynroot ) then
         call restartvar(ncid=ncid, flag=flag, varname='root_depth', xtype=ncd_double,  &
              dim1name='pft', &
              long_name='root depth', units='m', &
              interpinic_flag='interp', readvar=readvar, data=this%root_depth_patch)
 
         call restartvar(ncid=ncid, flag=flag, varname='rootfr', xtype=ncd_double,  &
              dim1name='pft', dim2name='levgrnd', switchdim=.true., &
              long_name='root fraction', units='', &
              interpinic_flag='interp', readvar=readvar, data=this%rootfr_patch)
     end if
         if (flag=='read' .and. .not. readvar) then
            if (masterproc) then
               write(iulog,*) "can't find rootfr in restart (or initial) file..."
               write(iulog,*) "Initialize rootfr to default"
            end if
            call init_vegrootfr(bounds, nlevsoi, nlevgrnd, &
            this%rootfr_patch(bounds%begp:bounds%endp,1:nlevgrnd), 'water')
            call init_vegrootfr(bounds, nlevsoi, nlevgrnd, &
            this%crootfr_patch(bounds%begp:bounds%endp,1:nlevgrnd), 'carbon')
         end if

It looks like the final block was probably within the use_dynroot conditional at one point, but then was pulled outside of it for some reason. But now it doesn't make any sense: if use_dynroot is false, then the 'readvar' is going to refer to the HK read, not the rootfr read.

So my sense is that this should probably just be removed, but I haven't given this a close inspection, so someone else should confirm that.



components/clm/src/biogeophys/IrrigationMod.F90
    if (flag=='read' .and. .not. readvar) then
       this%n_irrig_steps_left_patch = 0
    end if

The two conditionals here (one for n_irrig_steps_left_patch and one for irrig_rate_patch) can simply be removed: setting n_irrig_steps_left_patch = 0 doesn't accomplish anything because it is already set to 0 in initialization; and setting irrig_rate_patch isn't important in this case, because it's unused if n_irrig_steps_left_patch = 0.

(As an aside: Probably this isn't worth doing now, but it actually looks like the backwards compatibility stuff in this module (avoiding a read if the variables have the wrong size) could mostly be removed at this point: The backwards compatibility for n_irrig_steps_left and irrig_rate is 3 years old, and I think this wouldn't work right with init_interp (triggering a problem similar to http://bugs.cgd.ucar.edu/show_bug.cgi?id=2240). And the backwards compatibility for irrig_rate_demand could now be reworked to use your new infrastructure, but I don't think it's worth the effort, since we'll probably just remove this backwards compatibility stuff after the clm5 release.)





components/clm/src/biogeophys/TemperatureType.F90
    if (flag=='read' .and. .not. readvar) then
       this%t_h2osfc_col(bounds%begc:bounds%endc) = 274.0_r8
    end if

This block is unnecessary because it was already set to this value in initCold.


    if (flag=='read' .and. .not. readvar) call endrun(msg=errMsg(sourcefile, __LINE__))

       if (flag=='read' .and. .not. readvar) then
          if (masterproc) write(iulog,*) "can't find t_building in initial file..."
          if (masterproc) write(iulog,*) "Initialize t_building to taf"
          this%t_building_lun(bounds%begl:bounds%endl) = this%taf_lun(bounds%begl:bounds%endl)
       end if


Yes, looks right to replace this with your new infrastructure.

[billsacks]

(Not having looked back at the notes in this issue for a while.) I think that resolving this issue will also set water isotope values correctly when not using init_interp, when they aren't on the finidat_file. If not, we need a separate issue for that.

[billsacks]

I'm adding an error-check to the code that prevents running init_interp from a case without C isotopes to a case with C isotopes. Once we have fixed this issue, search the code for references to "BUG... #67" for things to remove related to this error-check.

[czarakas]

I'm planning to run some coupled CESM (CAM+CLM) experiments with 13C on, and I'd ideally like to use a finidat file generated from a run where isotopes were off. I'm only going to be looking at 13C GPP, so I don't think it should matter whether 13C is spun up. I am considering bypassing the error check with for_testing_allow_interp_non_ciso_to_ciso = .true. so that I can use the initial conditions file without isotopes and turn on the isotopes when I start running the new simulations. Do you think this workaround will work OK (for 13C GPP only) given this issue?

[dlawrenncar]

I can't remember how this works. We will discuss it at our software meeting tomorrow and get back to you.

…

On Mon, Jan 4, 2021 at 2:22 PM Claire Zarakas ***@***.***> wrote: I'm planning to run some coupled CESM (CAM+CLM) experiments with 13C on, and I'd ideally like to use a finidat file generated from a run where isotopes were off. I'm only going to be looking at 13C GPP, so I don't think it should matter whether 13C is spun up. I am considering bypassing the error check with for_testing_allow_interp_non_ciso_to_ciso = .true. so that I can use the initial conditions file without isotopes and turn on the isotopes when I start running the new simulations. Do you think this workaround will work OK (for 13C GPP only) given this issue? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#67 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AFABYVGGZKFBKVT7MFAWDQDSYIWQBANCNFSM4EIRZDMA> .

[billsacks]

@czarakas - From reading back through this issue, and from what I remember about it: This issue only applies if you are doing a run with C isotopes where the finidat file is from a run without C isotopes AND you are using init_interp.

So I think a workaround is:

• Do a short (e.g., 1-day) run with the desired configuration, but with isotopes turned off. If you are using init_interp in this run, then you will get a finidat_interp_dest.nc file in your run directory.
• Do the actual run (with C isotopes on), using that finidat_interp_dest.nc file as your finidat file (in user_nl_clm), with use_init_interp = .false.

If that doesn't work for some reason, then we think that it should be okay to just bypass this error check if all you care about is GPP - since it shouldn't matter that your C isotope values are nowhere near spunup for the sake of GPP C13. However, we aren't 100% sure about that, which is why I'd suggest first trying the above workaround.

[czarakas]

Thanks @billsacks! That makes sense to me. I'll plan on using the workaround you suggested.

[billsacks]

There may be more significant issues with going from an initial conditions file without C isotopes to a run with C isotopes, leading to garbage values rather than just starting with cold start C isotope values -- though maybe the former can arise as a result of the latter? See recent discussion in #1358

I guess the important thing is that, once we address the known issue documented here, we should do some careful sanity checks to ensure that it's really working correctly to interpolate from a non-ciso finidat to a ciso case.