COPSE V2.0 release

Files copied from PALEOexeter repository
branch: copseRLrelease
tag: COPSE-V2.0
SHA1 ID: 77a62885441852dd1ddcde58f91329c8a6b726b2

.gitignore

@@ -0,0 +1,5 @@
+output
+examples/copse/COPSE_test_output_full
+run_copse_tests.txt
+*~
+*asv

COPSE_setup.m

@@ -0,0 +1,63 @@
+function COPSE_setup
+% Set up COPSE paths
+
+% Add new directories to this list
+
+paleodirs={'code','code/core','code/forcings', ...
+    'code/copse', ...
+    'code/configuration','code/utils',  ...
+    'libraries/YAMLMatlab_0.4.3', ...
+    'examples/copse'};
+
+paleopath = pwd;
+
+% get current path as semicolon-separated list
+p = path;
+
+% check for any PALEO or COPSE entries in current path
+
+pentries = strsplit(p,pathsep);
+
+pcopse = {};
+for i = 1:length(pentries)
+    if ~isempty(strfind(pentries{i},'PALEO')) || ~isempty(strfind(pentries{i},'COPSE'))
+        pcopse{end+1} = pentries{i};
+    end
+end
+
+% prompt user and remove any PALEO entries
+% (eg if there are two installations ...)
+if ~isempty(pcopse)
+    fprintf('possible COPSE paths found\n');
+    for i=1:length(pcopse)
+           fprintf('         %s\n',pcopse{i})
+    end
+
+    str = input('\nRemove these paths ? Y/N [Y]\n','s');
+    if isempty(str)
+        str = 'Y';
+    end
+    if strcmpi(str,'Y')
+        for i=1:length(pcopse)
+            fprintf('removing folder %s\n',pcopse{i});
+            rmpath(pcopse{i});
+        end
+    end
+end
+
+% Add PALEO paths
+
+for i=1:length(paleodirs)
+    psepdir = strrep(paleodirs{i},'/',filesep);
+    fulldir = fullfile(paleopath,psepdir);
+    fprintf('adding folder %s\n',fulldir);
+    addpath(fulldir);
+end
+
+% Test yaml hence preload - attempt to workaround an issue on linux with crash ?
+
+yaml_file = 'libraries/YAMLMatlab_0.4.3/Tests/Data/test_import/file1.yaml';
+YamlStruct = ReadYaml(yaml_file);
+
+
+

KnownIssues.txt

@@ -0,0 +1,76 @@
+COPSE known issues 2017-12-01
+-----------------------------
+
+Lenton etal (2017) COPSE Reloaded
+---------------------------------
+COPSE_reloaded_reloaded.m
+
+1) Spurious warnings reported to screen when loading some comparison output datasets from older models.
+  >>  Warning: Cannot load an object of class 'copse_model_sfbw':
+   (seen with some datasets in additional output data in COPSE_V2.0_test_output_full.zip)
+   No effect on the data loaded.
+
+Mills etal (2014) G3
+--------------------
+COPSE_mills2014g3_mills2014g3.m
+
+The below minor fixes are made to this version of the model. These changes make little difference to output and do not alter the conclusions of the original paper. The combined effect of changes relative to the original paper can be demonstrated by editing COPSE_mills2014g3_mills2014g3.m to select 
+run=copse_millsg3_millsg3_expts('g3mills2014nobugs', 'baseline');
+
+1) Land biota T limitation effectively removed (copse_landbiota_mills2014g3.m)
+
+   To reproduce, set
+    f_bug_g32014_landbiotatemp = 'Yes'
+   See https://github.com/sjdaines/PALEOexeter/issues/23 
+
+   Effect is: pCO2 is ~0.1 PAL low (200Ma-0), ~< 0.6 PAL low (250Ma-200Ma)
+              pO2 (not plotted in paper) is ~< 0.06 PAL high (200Ma-0), ~<0.1 PAL high (250Ma-200Ma)
+
+2) Sr concentration calculation was incorrect (not plotted in paper)
+    (copse_model_mills2014g3.m)
+
+    f_bug_g32014_Srconc: 'Yes' 
+
+3) Granite area calculation only used CFB area (copse_landsurfaceareas_mills2014g3.m)
+
+   f_granitearea = 'G3original'
+
+4) Degass forcing extrapolation to times earlier than 230Ma inconsistent
+   (copse_force_vandermeer.m)
+   DEGASS returned to background value (=1) at 230Ma, whereas OIB area assumed constant area for T > 230Ma
+
+   copse_force_vandermeer.extrapolate = 3
+
+5) Small (~1%) discrepancy between basalt area calculated from LIP table vs offline calculation.
+   (test_copse_load_phanlip('mills2014g3'), ie copse_load_phanlip.m  vs copse_force_revision_ba.m + either
+    bas_area_mills2014.xlsx or ggge20620-sup-0002-suppinfo2.mat)
+
+   https://github.com/sjdaines/PALEOexeter/issues/18
+
+6) Time-evolution of delta_Sr_sed
+   (copse_model_mills2014g3.m)
+
+   Updated in copse_model_reloaded.m
+
+7) Silicate weathering used 288 not 288.15 K offset, resulting in ~1% offset in pCO2
+   (copse_weathering_rates_mills2014g3.m)
+    pars.f_act_energies = 'split_bug_g32014_Toffset'
+
+    https://github.com/sjdaines/PALEOexeter/issues/21
+
+Bergman etal (2004)
+-------------------
+COPSE_bergman2004_bergman2004.m
+
+1) Organic carbon degassing rolloff at low pO2 (copse_model_bergman2004.m)
+     # COPSE 5_14 C code (and Bergman 2004) use 'O2copsecrashprevent' which rolls
+     # off organic carbon degassing at low pO2. This has a big effect at low pO2 when 
+     # oxidative weathering is oxygen-independent (ie Ordovician and earlier)
+
+     f_ocdeg        : 'O2indep'    # options 'O2indep', 'O2copsecrashprevent'
+
+2) Marine N (unlikely to affect results as N adjusts due to denitrification alone)
+   (copse_marinebiota_bergman2004.m)
+
+     f_nfix_nreplete: 'Off'       # options 'Off','Sign' (COPSE 5_14 C code has 'Sign')
+

README.txt

@@ -0,0 +1,36 @@
+
+The COPSE (Carbon, Oxygen, Phosphorus, Sulphur and Evolution) biogeochemical model predicts the coupled histories and controls on atmospheric O2, CO2 and ocean composition over Phanerozoic time.
+
+The model is described in the following publications:
+
+Bergman, N. M., Lenton, T. M., & Watson, A. J. (2004). COPSE: A new model of biogeochemical cycling over Phanerozoic time. American Journal of Science, 304(5), 397�437. http://doi.org/10.2475/ajs.304.5.397
+
+Mills, B., Daines, S. J., & Lenton, T. M. (2014). Changing tectonic controls on the long-term carbon cycle from Mesozoic to present. Geochemistry, Geophysics, Geosystems, 15(12), 4866�4884. http://doi.org/10.1002/2014GC005530
+
+Lenton, T. M., Dahl, T. W., Daines, S. J., Mills, B. J. W., Ozaki, K., Saltzman, M. R., & Porada, P. (2016). Earliest land plants created modern levels of atmospheric oxygen. Proceedings of the National Academy of Sciences, 113(35), 9704�9709. http://doi.org/10.1073/pnas.1604787113
+
+Lenton, T. M., Daines, S.J., Mills, B. J. W. (2017). COPSE reloaded: An improved model of biogeochemical cycling over Phanerozoic time. Earth Science Reviews, in revision.
+
+
+Running the model (requires Matlab version 2012 or higher):
+-----------------
+
+>> COPSE_setup                   % sets Matlab paths
+>> cd examples/copse
+>> run_copse_tests               % test against archived output (output is included for only the default set of 7 tests)
+>>
+>> COPSE_reloaded_reloaded       % runs model, plots output for Lenton etal (2017)
+>> COPSE_bergman2004_bergman2004 % Bergman etal (2004) version
+>> COPSE_millsg3_millsg3         % Mills etal (2014) version
+>> COPSE_reloaded_bergman2004    % includes results from Lenton etal (2016) - see comments in file
+
+
+Evaluation data:
+---------------
+
+Datasets are not part of the public release but are available on request from the authors.
+
+Known issues:
+------------
+
+Please see KnownIssues.txt

code/configuration/copse_modify_struct.m

@@ -0,0 +1,11 @@
+function strct = copse_modify_struct( strct, fld, newval )
+%Modify a struct, error if field doesn't already exist
+
+if ~isfield(strct,fld)
+    error ('struct "%s" no field "%s"',inputname(1),fld); 
+end
+
+strct.(fld)=newval;
+
+end
+

code/configuration/paleo_modelbuilder.m

@@ -0,0 +1,45 @@
+classdef paleo_modelbuilder
+    % Create and configure a model from configuration file.
+    %
+    %
+    methods(Static)
+        function tm  = createModel(configfile, configname)
+
+            LN = 'paleo_modelbuilder.createModel';  L = paleo_log.getLogger();
+
+            % Read configuration from yaml file
+            cfgset = paleo_parameterset(configfile, configname);
+
+            % Create model:     
+            % locate the relevant section in the cfg structure
+            modelcfg = paleo_parameterset.findElem(cfgset.configdata, 'model');                   
+            % create the model
+            L.log(L.DEBUG, LN, sprintf('creating tm ctorstr=''%s''\n', modelcfg.class));                
+            tm                                          = eval(modelcfg.class);
+
+            %%%%%%%% Set forcing functions for this run
+            tm.force = {};            
+            forcecfg = paleo_parameterset.findElem(cfgset.configdata, 'model.force');           
+            for i = 1:length(forcecfg)
+                ctorstr = forcecfg{i};
+                newforce = eval(ctorstr);
+                tm.force{end+1} = newforce;
+            end
+
+            tm.perturb = {};
+            perturbcfg = paleo_parameterset.findElem(cfgset.configdata, 'model.perturb');
+            for i = 1:length(perturbcfg)
+                ctorstr = perturbcfg{i};
+                newperturb = eval(ctorstr);
+                tm.perturb{end+1} = newperturb;
+            end
+
+            %%%%%%%%%% Set top-level parameters 
+            tm.pars = cfgset.configureobject(tm.pars, 'model.pars', 'pars');
+
+        end
+
+    end
+end
+
+

code/configuration/paleo_parameterset.m

@@ -0,0 +1,151 @@
+classdef paleo_parameterset
+    % A parameter set contains 'configdata' read from 'configname' in yaml 'configfile'.
+    %
+    % The yaml configfile should contain configuration data at the top level, eg
+    %
+    % copse_reloaded:
+    %       <configuration> 
+    %
+    % copse_bergman2004:
+    %      <another configuration>
+    %
+
+    properties(SetAccess=private)
+        configfile   = '';   % configfile read
+        configname   = '';   % named configuration from that file
+
+        configdata;                   % struct with data from configfile.(configname)
+    end
+
+    properties(Constant)
+        reservedFields = {'class'};    % fields in config file to ignore;
+    end
+
+    methods
+        function obj = paleo_parameterset(configfile, configname)
+            % Read configuration data from yaml parameter file (which may have multiple configs within)
+
+            % Read entire parameter file (which may contain multiple configurations)
+            obj.configfile = configfile;
+            filects = ReadYaml(configfile); 
+
+            % Get the data for the specified configname
+            obj.configname = configname;
+            if  isfield(filects, obj.configname)
+                obj.configdata = filects.(obj.configname);
+            else
+                error('no configname %s in configfile %s', obj.configname, obj.configfile);
+            end
+        end
+
+
+        function target = configureobject(obj, target, configpath, targetname)
+            % Set properties/fields on a single object 'target' from location 'configpath' in 'configdata'
+            %
+            % Input:
+            % target    - object to be configured 
+            % configpath - path within obj.configdata containing configuration for 'target'
+            % targetname - [optional] name of target object. This makes the target object fields
+            %              available when evaluating a parameter, to allow eg par1 = (2*pars.par2)        
+            %
+            % Returns:
+            % target    - configured object
+
+            LN = 'paleo_parameterset.configureobject'; L = paleo_log.getLogger(LN); 
+
+
+            % Find the configuration data at location configpath
+            % (configdata has been read from file, configpath is specified as an argument)
+            targetconfig = paleo_parameterset.findElem(obj.configdata, configpath);
+
+            if isempty(targetconfig)
+                error('no config for configfile %s configname %s configpath %s', obj.configfile, obj.configname, configarg);
+            else
+                % configure target ...
+
+                rpars = fields(targetconfig);
+
+                % apply config to target object properties, ignoring reserved fields
+                for rp = 1:length(rpars)
+                    if ~any(strcmp(rp, obj.reservedFields)) % ignore reserved fields
+                        subval = paleo_parameterset.subValue(targetconfig.(rpars{rp}), targetname, target);
+                        L.log(L.TRACE, LN, sprintf('<target>.%s=%s\n', rpars{rp}, subval)); 
+                        target.(rpars{rp}) = subval;
+                    end
+                end
+            end
+        end
+    end
+
+
+    methods(Static)
+        function subval = subValue(val, ctxtname, ctxtdata)
+            % parameter substitution.
+            %
+            % Input:
+            % val      - parameter value either as a string (which will be evaluated if an expression in round brackets),
+            %            or a string without brackets or non-string (which will just be returned directly)
+            % ctxtname - [optional] name of a struct or object (eg 'pars') to provide as context for evaluation.
+            % ctxtdata - [optional] value of this struct or object
+            %
+            % Returns: 
+            % subval   - parameter value, either as-supplied if non-string, or evaluated from a supplied string
+
+            if nargin < 2
+                ctxtname = '';
+            end
+
+            if ischar(val)
+                % evaluate string-valued parameter if of from '(expr)'
+
+                % guard against inadvertent use of 'false' or 'true' which are parsed as strings, not logicals
+                % use true and false (no quotes)
+                if any(strcmp(val,{'false','true'}))
+                    error('''false'' and ''true'' strings illegal in parameter value - remove quotes');
+                end
+                % TODO guard against strings that look like they might be supposed to be numeric expressions
+                % or attempts to set numeric fields to strings?
+
+                % If string is of form '(expr)', evaluate 
+                if ~isempty(val) && strcmp(val(1),'(') && strcmp(val(end),')')
+                    % evaluate contents of brackets so (2+2) = 4 etc
+                    if ~isempty(ctxtname)
+                        % provide optional context for evaluation
+                        eval([ctxtname '= ctxtdata;']);
+                    end
+                    subval = eval(val);
+                else
+                    % no brackets - just return the raw string
+                    subval = val;
+                end
+            else
+                % return unmodified non-string value
+                subval = val;
+            end
+        end
+
+
+        function elem = findElem(data, path)
+            % find element within structure-tree 'data' at location 'path'
+            % eg
+            % data.x.y = elem
+            % path = 'x.y' 
+
+            % split path on .              
+            psplit = strsplit(path,'.');
+
+            % navigate into data to find elem
+            elem = data;
+            if ~isempty(psplit{1}) % NB: strsplit('','.') returns {''} 
+                for i=1:length(psplit)
+                    if isfield(elem, psplit{i})
+                        elem = elem.(psplit{i});
+                    else
+                        error('path %s not present',  path);                       
+                    end
+                end
+            end           
+        end
+    end
+end
+