index script ships with the rest

functions/pism/pism.functions

@@ -417,7 +417,7 @@ pism_prepare_forcing() {
            echo "Index values at $RUN_NUMBER_pism"
            index_values=$(sed "${RUN_NUMBER_pism}q;d" $index_file)
            echo $index_values | read kyr index
-           index_script=${index_script:-${HOME}/first.py}
+           index_script=${index_script:-${FUNCTION_PATH}/pism_index_script.py}
 
            echo "Using PI_file=$PI_file"
            echo "Using LGM_file=$LGM_file"

functions/pism/pism_index_script.py

+#!/usr/bin/env python
+
+# Import from Standard Library
+import argparse
+import logging
+import sys
+
+# Import Third-Party Packages
+
+from climlab.solar.insolation import daily_insolation
+from climlab.solar.orbital import OrbitalTable
+import climlab
+import numpy as np
+import xarray as xr
+
+
+def parse_arguments():
+    parser = argparse.ArgumentParser(description='Generates a pseudo ECHAM6 file given an index and two climate endmembers')
+    parser.add_argument("warm_file")
+    parser.add_argument("cold_file")
+    parser.add_argument("index", type=float)
+    parser.add_argument("kyr", type=int)
+    parser.add_argument("outfile_name")
+    parser.add_argument("-v", action="store_true", dest="debug")
+    return parser.parse_args()
+
+
+def index_forcing(kyr, index, infile1, infile2, outfile):
+    """
+    # synthesizes atm. forcing files for PISM, based on one index, two climatologies 
+    # (a "warm" and a "cold" climate in infile1 and infile2 respectively) and a specific time in kyrs
+    # (LGM is -21, does not seem to work for future!!!)
+    #
+    # input file is expected to have climatologic monthly temp2, precip and geosp on a lonxlatx12 grid
+    #
+    # output forcing goes to netcdf outfile and includes 
+    # temperature T = index*warm.temp2+(1-index)*cold.temp2
+    # precipitation PP analogue to temp2
+    # surface elevation z analogue to temp2
+    # and solar radiation is calculated from climlab for each month on the same grid  
+    """
+
+    logging.debug("kyr=%s" % str(kyr))
+    logging.debug("type(kyr)=%s" % str(type(kyr)))
+    logging.debug("index=%s" % str(index))
+    logging.debug("infile1=%s" % str(infile1))
+    logging.debug("infile2=%s" % str(infile2))
+    logging.debug("outfile=%s" % str(outfile))
+
+    warm = xr.open_dataset(infile1)
+    lat = warm.lat
+    lon = warm.lon
+    nx = lon.shape[0]
+    ny = lat.shape[0]
+    dtype = 1
+    if dtype == 1:
+        SOLm = np.zeros((365, ny))
+        daya = np.arange(1.78, 366.8)  # (2.3,367.3)
+    else:
+        SOLm = np.zeros((360, ny))
+        daya = np.arange(0, 360)
+
+    srad0dI = np.zeros((12, ny, nx))
+    sradtest = np.zeros((12, ny, nx))
+    w_temp = warm.temp2.values
+    w_aprc = warm.aprc.values
+    w_aprl = warm.aprl.values
+    w_z = warm.geosp.values
+    w_srad0d = warm.srad0d.values
+    w_sradsu = warm.sradsu.values
+    w_srads = warm.srads.values
+    w_alb = warm.albedo.values  # albedo
+    w_sradsd = w_srads - w_sradsu  # downward solar radiation at the surface
+    w_tao = w_sradsd / w_srad0d  # else surface/TOA downward solar radiation
+
+    cold = xr.open_dataset(infile2)
+    c_temp = cold.temp2.values
+    c_aprc = cold.aprc.values
+    c_aprl = cold.aprl.values
+    c_z = cold.geosp.values
+    c_srad0d = cold.srad0d.values
+    c_sradsu = cold.sradsu.values
+    c_srads = cold.srads.values
+    c_alb = cold.albedo.values  # albedo
+    c_sradsd = c_srads - c_sradsu  # downward solar radiation at the surface
+    c_tao = c_sradsd / c_srad0d  # else surface/TOA downward solar radiation
+
+    # subset of orbital parameters for a specified year in kyr (LGM is -21, careful: not yet defined for kyear>0!)
+    orb = OrbitalTable.interp(kyear=kyr)
+    orb2 = OrbitalTable.interp(kyear=0)
+    orb_PI = {"ecc": 0.016724, "long_peri": 282.157, "obliquity": 23.4468}
+    orb_LGM = {"ecc": 0.018994, "long_peri": 294.42, "obliquity": 22.949}
+    # insolation values for the specific year for all days of a year and all latitudes of the atmospheric grid
+
+    SOL = daily_insolation(lat=lat, day=daya, orb=orb, S0=1365.0, day_type=dtype)
+    # compute monthly mean insolation from daily values
+    if dtype == 1:
+        days_per_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
+        firstday = 0
+    else:
+        days_per_month = 30
+        firstday = 360 - 80
+        # firstday=0
+
+    for i in range(0, 12):
+        if dtype == 1:
+            days = np.arange(firstday, firstday + days_per_month[i])  # dtype =1
+            firstday = firstday + days_per_month[i]
+        else:
+            if (firstday + days_per_month) < 360:
+                days = np.arange(firstday, firstday + days_per_month)
+                firstday = firstday + days_per_month
+            else:
+                first_day_of_next_month = firstday + days_per_month - 360
+                days = np.concatenate((range(firstday, 360), range(0, first_day_of_next_month)))
+                firstday = first_day_of_next_month
+        SOLm = SOL[days, :].mean(axis=0)
+        for j in range(0, ny):
+            srad0dI[i, j, :] = SOLm[j]
+
+    # combine "warm" and "cold" fields according to index
+    T2I = (index * w_temp) + ((1 - index) * c_temp)
+    PPCI = index * w_aprc + (1 - index) * c_aprc
+    PPLI = index * w_aprl + (1 - index) * c_aprl
+    ZI = index * w_z + (1 - index) * c_z
+    taoI = index * w_tao + (1 - index) * c_tao  # scaled transmissivity taoI
+    albI = index * w_alb + (1 - index) * c_alb  # scaled albedo
+    sradsdI = (
+        taoI * srad0dI
+    )  # downward surf. radiation from taoI and TOA insolation from climlab
+    sradsuI = -albI * sradsdI  # but we actually only need UPWARD surf. radiation
+    sradsI = (
+        1 - albI
+    ) * sradsdI  # and net surface radiation... I hope we can use albedo here...
+
+    out = xr.open_dataset(infile1)
+    out.temp2.values = T2I
+    out.aprc.values = PPCI
+    out.aprl.values = PPLI
+    out.geosp.values = ZI
+    out.srads.values = sradsI
+    out.sradsu.values = sradsuI
+    out.srad0d.values = srad0dI
+    out.albedo_vis.values = albI
+    out.to_netcdf(path=outfile, mode="w", compute=True)
+
+
+if __name__ == "__main__":
+    args = parse_arguments()
+    if args.debug:
+        logging.basicConfig(level=logging.DEBUG)
+    index_forcing(args.kyr, args.index, args.warm_file, args.cold_file, args.outfile_name)