max daylength is hard-coded for present-day orbital parameters

[billsacks]

Bill Sacks < [email protected] > - 2013-10-16 09:49:44 -0600
Bugzilla Id: 1843
Bugzilla CC: [email protected], [email protected], [email protected], [email protected], [email protected],

In doing my daylength refactor (https://svn-ccsm-models.cgd.ucar.edu/clm2/branches/daylength_refactor/), I discovered that max daylength is hard-coded for present-day orbital parameters. In talking with Nan, this could be a problem for present-day runs.

It is easiest to see the problem on the above-referenced branch (which will come to the trunk soon).

I believe that the daylength calculation (in models/lnd/clm/src/clm4_5/biogeophys/DaylengthMod.F90 - see the daylength function) is agnostic to orbital parameters, so that should be fine - as long as it is getting the true, current solar declination angle, which I believe is the case.

The problem is in src/clm4_5/main/iniTimeConst.F90 - see the setting of max_decl there. There are actually two related problems here:

(1)

[billsacks]

Bill Sacks < [email protected] > - 2013-10-16 10:06:16 -0600

Oops, didn't finish that; to finish:

The two related problems are:

(1) max_dayl is computed using an orbital parameter that is hard-coded for present-day

(2) max_dayl is set in initialization, then doesn't change... this isn't itself a problem currently, but would be a problem if we set it based on the actual, current orbital parameters

So I propose the following solution:

(1) Create a new function in DaylengthMod that computes the max daylength for a given grid cell. This would look like:

elemental real(r8) function max_daylength(lat, obliquity)
   max_decl = obliquity
   if (lat < 0) max_decl = -max_decl
   max_daylength = daylength(lat, max_decl)
end function

(2) Do NOT set max_dayl in iniTimeConst. Instead, call the above function wherever it is needed - currently, this is just in one place, in CanopyFluxesMod.F90. So here we could have a local variable, max_dayl, and calculate it as:

max_dayl(begg:endg) = max_daylength(lat, obliquity)

(3) We then have to get the obliquity from the coupler, which I believe is already done in this call in lnd_comp_mct:

   call shr_orb_decl( calday     , eccen, mvelpp, lambm0, obliqr, declin  , eccf )

and then pass obliquity as an argument to the driver, which will then pass it along to CanopyFluxes.

I would like confirmation from Nan and Erik on these points:

(a) Nan (or Erik): Is obliquity the correct parameter to use here? What we want is the maximum solar declination angle; in iniTime const this is set as:

  ! +/- 23.4667 degrees = +/- 0.409571 radians, use negative value for S. Hem

So, is that the obliquity?

(b) Mostly for Erik: Does this seem like a reasonable solution? And do you agree with the overall design I am proposing? In particular, does it make sense to get obliquity in lnd_comp_mct and pass it down the call chain? The alternative would be to call shr_orb_decl from within the max_daylength function, but from our brief shuttle conversation, it sounded like there is some design reason to get these sort of parameters at the top level and then pass them down?

[billsacks]

Bill Sacks < [email protected] > - 2013-10-16 10:11:01 -0600

This change will almost certainly change answers slightly, even for present-day runs (because I imagine that the hard-coded max declination angle in iniTimeConst doesn't agree completely with CESM's obliquity).

So for testing purposes, this should be done in two steps:

(1) Simply change the hard-coded max declination angle in iniTimeConst to agree with CESM's obliquity. This will change answers. (If this differs for different time periods, e.g. 1850 vs 2000, I may want to introduce a kludge so that the hard-coded value keys off of the time period of the run (or maybe introduce a call in iniTimeConst to shr_orb_decl) - so that this hard-coded value always agrees with CESM's obliquity at initialization.

(2) Do my proposed refactoring. For runs in which the orbital parameters are fixed throughout the run, this should be bfb with #1.

[billsacks]

Bill Sacks < [email protected] > - 2013-10-16 10:13:16 -0600

Oops, when I said "or maybe introduce a call in iniTimeConst to shr_orb_decl" I meant shr_orb_params

[billsacks]

Bill Sacks < [email protected] > - 2013-12-17 11:18:05 -0700

Un-assigning myself from this... I was only assigned to this because it came up in the course of my work a few months ago, but it's not currently relevant for me.

[billsacks]

Bill Sacks < [email protected] > - 2015-05-08 09:44:02 -0600

Looking back at this bug report, I noticed a typo in the initial bug report:

In talking with Nan, this could be a problem for present-day runs.

Should read:

In talking with Nan, this could be a problem for PALEO runs.

[billsacks]

Erik Kluzek < [email protected] > - 2017-11-02 15:01:17 -0600

The simplest change to accomplish this is to just do...

Index: main/clm_initializeMod.F90
===================================================================
--- main/clm_initializeMod.F90	(revision 87306)
+++ main/clm_initializeMod.F90	(working copy)
@@ -385,11 +385,10 @@
     call InitDaylength(bounds_proc, declin=declin, declinm1=declinm1)
              
     ! Initialize maximum daylength, based on latitude and maximum declination
-    ! maximum declination hardwired for present-day orbital parameters, 
-    ! +/- 23.4667 degrees = +/- 0.409571 radians, use negative value for S. Hem
+    ! given by the obliquity use negative value for S. Hem
 
     do g = bounds_proc%begg,bounds_proc%endg
-       max_decl = 0.409571
+       max_decl = obliqr
        if (grc%lat(g) < 0._r8) max_decl = -max_decl
        grc%max_dayl(g) = daylength(grc%lat(g), max_decl)
     end do

This looks correct to me, since in shr_orb_calc

declin = asin(sin(obliqr)*sin(lamb))

lamb represents the earth going around the sun, so we just want the point where sin(lamb)==1. So you have

declin = asin(sin(obliqr)*1) = obliqr

So the above change is the quickest solution to the right answer. Other refactoring could be done later.

[billsacks]

Bill Sacks < [email protected] > - 2017-11-02 16:41:17 -0600

Erik: At a glance, your solution looks like a simple and reasonable solution for the first problem:

(1) max_dayl is computed using an orbital parameter that is hard-coded for present-day

but doesn't address

(2) max_dayl is set in initialization, then doesn't change... this isn't itself a problem currently, but would be a problem if we set it based on the actual, current orbital parameters

It probably still makes sense to bring your solution to the trunk, since that solves the first-order problem. But then we should open a separate issue to eventually solve (2), since my understanding is that that's needed if you have time-evolving orbital parameters.

[billsacks]

Erik Kluzek < [email protected] > - 2017-11-02 17:32:31 -0600

A few years ago there was only the option to set orbital parameters at initialization. It looks like there is now the hardware in place to actually change it each year. But, this isn't normally exercised, and I couldn't find any tests that actually turn it on. Compsets don't seem to trigger it either. Since it's untested I wonder if it actually works? But, in any case, the simple change is an improvement, since it'll be aligned for the first year, and just slowly be off which each passing year. The simple change will be right for everything we test, and just off for a feature that isn't tested.

[billsacks]

Bill Sacks < [email protected] > - 2017-11-02 17:53:41 -0600

Okay, thanks for those details.

[billsacks]

Erik Kluzek < [email protected] > - 2017-11-08 16:06:15 -0700

clm4_5_17_r263 does a partial fix to this by correcting it at initialization. variable_years will be wrong though.

[billsacks]

@ekluzek I'm closing this issue, and opening #260 to document the remaining issue.