{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import netCDF4 as nc\n",
    "import matplotlib.pylab as plt\n",
    "import imp\n",
    "import csv\n",
    "import pandas as pd\n",
    "import random as rnd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "651\n"
     ]
    }
   ],
   "source": [
    "# Read atmospheric forcing\n",
    "\n",
    "NumTensemble = 600\n",
    "Tlen = 651\n",
    "\n",
    "fname = \"../MAGICC/CMIP5_RCP2500/Temp_CMIP5_RCP85_HistConstraint.dat\" # File to read\n",
    "df = pd.read_csv(fname,sep='\\s+',index_col=0,header=None)\n",
    "df.columns.name = \"ensemble member\"\n",
    "df.index.name = \"Time\"\n",
    "SAT = np.array(df.values)\n",
    "\n",
    "print(len(SAT[:,1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Normalize and crop temperature series\n",
    "Temp = []\n",
    "Tavebeg = 0\n",
    "Taveend = 80\n",
    "tbeg = 51\n",
    "tend = 251\n",
    "for i in range(len(SAT[1,:])):\n",
    "    SATave = np.mean(SAT[Tavebeg:Taveend,i])\n",
    "    SAT[:,i] = SAT[:,i]-SATave\n",
    "SAT = SAT[tbeg:tend,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read ocean scaling\n",
    "\n",
    "NumOmodel = 19\n",
    "\n",
    "fname = \"../ScalingCoefficients/OceanScaling/OS_R1.dat\" # File to read\n",
    "with open(fname) as f:\n",
    "    OS_NoDelay_R1 = np.array([float(row.split(\"\\t\")[0]) for row in f])\n",
    "with open(fname) as f:\n",
    "    OS_WiDelay_R1 = np.array([float(row.split(\"\\t\")[3]) for row in f])\n",
    "with open(fname) as f:\n",
    "    OS_Delay_R1 = np.array([float(row.split(\"\\t\")[2]) for row in f])\n",
    "\n",
    "fname = \"../ScalingCoefficients/OceanScaling/OS_R2.dat\" # File to read\n",
    "with open(fname) as f:\n",
    "    OS_NoDelay_R2 = np.array([float(row.split(\"\\t\")[0]) for row in f])\n",
    "with open(fname) as f:\n",
    "    OS_WiDelay_R2 = np.array([float(row.split(\"\\t\")[3]) for row in f])\n",
    "with open(fname) as f:\n",
    "    OS_Delay_R2 = np.array([float(row.split(\"\\t\")[2]) for row in f])\n",
    "\n",
    "fname = \"../ScalingCoefficients/OceanScaling/OS_R3.dat\" # File to read\n",
    "with open(fname) as f:\n",
    "    OS_NoDelay_R3 = np.array([float(row.split(\"\\t\")[0]) for row in f])\n",
    "with open(fname) as f:\n",
    "    OS_WiDelay_R3 = np.array([float(row.split(\"\\t\")[3]) for row in f])\n",
    "with open(fname) as f:\n",
    "    OS_Delay_R3 = np.array([float(row.split(\"\\t\")[2]) for row in f])\n",
    "\n",
    "fname = \"../ScalingCoefficients/OceanScaling/OS_R4.dat\" # File to read\n",
    "with open(fname) as f:\n",
    "    OS_NoDelay_R4 = np.array([float(row.split(\"\\t\")[0]) for row in f])\n",
    "with open(fname) as f:\n",
    "    OS_WiDelay_R4 = np.array([float(row.split(\"\\t\")[3]) for row in f])\n",
    "with open(fname) as f:\n",
    "    OS_Delay_R4 = np.array([float(row.split(\"\\t\")[2]) for row in f])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read melting sensitivity\n",
    "fname = \"../ScalingCoefficients/MeltSensitivity/MeltSensitivity.dat\" # File to read\n",
    "with open(fname) as f:\n",
    "    MeltSensitivity = np.array([float(row) for row in f])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ISSM_JPL = ISSM_JPL\n",
    "# Read ice scaling\n",
    "fname = \"../ScalingCoefficients/IceScaling/ES_ISSM_JPL_scaledto08.dat\" # File to read\n",
    "with open(fname) as f:\n",
    "    FS_ISSM_JPL_BM08 = np.array([float(row) for row in f])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ISSM_JPL = ISSM_JPL\n",
    "# Read response functions\n",
    "\n",
    "fname = \"../RFunctions/RF_ISSM_JPL_BM08_R1.dat\" # File to read\n",
    "with open(fname) as f:\n",
    "    RF_ISSM_JPL_BM08_R1 = np.array([float(row) for row in f])\n",
    "\n",
    "fname = \"../RFunctions/RF_ISSM_JPL_BM08_R2.dat\" # File to read\n",
    "with open(fname) as f:\n",
    "    RF_ISSM_JPL_BM08_R2 = np.array([float(row) for row in f])\n",
    "\n",
    "fname = \"../RFunctions/RF_ISSM_JPL_BM08_R3.dat\" # File to read\n",
    "with open(fname) as f:\n",
    "    RF_ISSM_JPL_BM08_R3 = np.array([float(row) for row in f])\n",
    "\n",
    "fname = \"../RFunctions/RF_ISSM_JPL_BM08_R4.dat\" # File to read\n",
    "with open(fname) as f:\n",
    "    RF_ISSM_JPL_BM08_R4 = np.array([float(row) for row in f])\n",
    "\n",
    "fname = \"../RFunctions/RF_ISSM_JPL_BM08_R5.dat\" # File to read\n",
    "with open(fname) as f:\n",
    "    RF_ISSM_JPL_BM08_R5 = np.array([float(row) for row in f])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.53455 0.5555555555555556\n",
      "0.89375 0.6666666666666666\n",
      "0.5288 0.5555555555555556\n",
      "0.4707 1.0\n"
     ]
    }
   ],
   "source": [
    "EnsembleSize = 20000\n",
    "scaled_forcing = False\n",
    "\n",
    "countR1 = 0\n",
    "countR2 = 0\n",
    "countR3 = 0\n",
    "countR4 = 0\n",
    "\n",
    "SL_wTd_nos_base_ISSM_JPL_R1_RCP85 = [0] * (tend-tbeg)\n",
    "SL_wTd_nos_base_ISSM_JPL_R2_RCP85 = [0] * (tend-tbeg)\n",
    "SL_wTd_nos_base_ISSM_JPL_R3_RCP85 = [0] * (tend-tbeg)\n",
    "SL_wTd_nos_base_ISSM_JPL_R4_RCP85 = [0] * (tend-tbeg)\n",
    "SL_wTd_nos_base_ISSM_JPL_R5_RCP85 = [0] * (tend-tbeg)\n",
    "\n",
    "for i in range(EnsembleSize):\n",
    "\n",
    "    # Select forcing randomly\n",
    "\n",
    "    # select global warming path\n",
    "    iTens = rnd.randint(0,NumTensemble-1)\n",
    "    Temp = np.array(SAT[:,iTens])\n",
    "\n",
    "    # select ocean model\n",
    "    iOmod = rnd.randint(0,NumOmodel-1)\n",
    "    OS_R1 = OS_WiDelay_R1[iOmod]\n",
    "    OS_R2 = OS_WiDelay_R2[iOmod]\n",
    "    OS_R3 = OS_WiDelay_R3[iOmod]\n",
    "    OS_R4 = OS_WiDelay_R4[iOmod]\n",
    "    OS_R5 = OS_WiDelay_R4[iOmod]\n",
    "\n",
    "    tau_R1 = int(OS_Delay_R1[iOmod])\n",
    "    tau_R2 = int(OS_Delay_R2[iOmod])\n",
    "    tau_R3 = int(OS_Delay_R3[iOmod])\n",
    "    tau_R4 = int(OS_Delay_R4[iOmod])\n",
    "    tau_R5 = int(OS_Delay_R4[iOmod])\n",
    "\n",
    "    if tau_R1>0:\n",
    "        countR1 = countR1+1\n",
    "    if tau_R2>0:\n",
    "        countR2 = countR2+1\n",
    "    if tau_R3>0:\n",
    "        countR3 = countR3+1\n",
    "    if tau_R4>0:\n",
    "        countR4 = countR4+1\n",
    "    \n",
    "    Temp_R1 = np.append(np.zeros(tau_R1),Temp[:tend-tbeg-tau_R1])\n",
    "    Temp_R2 = np.append(np.zeros(tau_R2),Temp[:tend-tbeg-tau_R2])\n",
    "    Temp_R3 = np.append(np.zeros(tau_R3),Temp[:tend-tbeg-tau_R3])\n",
    "    Temp_R4 = np.append(np.zeros(tau_R4),Temp[:tend-tbeg-tau_R4])\n",
    "    Temp_R5 = np.append(np.zeros(tau_R5),Temp[:tend-tbeg-tau_R5])\n",
    "    \n",
    "    # select melting sensitivity\n",
    "    MS_R1 = rnd.uniform(MeltSensitivity[0],MeltSensitivity[1])\n",
    "    MS_R2 = rnd.uniform(MeltSensitivity[0],MeltSensitivity[1])\n",
    "    MS_R3 = rnd.uniform(MeltSensitivity[0],MeltSensitivity[1])\n",
    "    MS_R4 = rnd.uniform(MeltSensitivity[0],MeltSensitivity[1])\n",
    "    MS_R5 = rnd.uniform(MeltSensitivity[0],MeltSensitivity[1])\n",
    "\n",
    "    # Compose forcing time series\n",
    "    M_R1 = MS_R1*OS_R1*Temp_R1\n",
    "    M_R2 = MS_R2*OS_R2*Temp_R2\n",
    "    M_R3 = MS_R3*OS_R3*Temp_R3\n",
    "    M_R4 = MS_R4*OS_R4*Temp_R4\n",
    "    M_R5 = MS_R5*OS_R5*Temp_R5\n",
    "\n",
    "    M_R1[M_R1 < 0.0] = 0.0\n",
    "    M_R2[M_R2 < 0.0] = 0.0\n",
    "    M_R3[M_R3 < 0.0] = 0.0\n",
    "    M_R4[M_R4 < 0.0] = 0.0\n",
    "    M_R5[M_R5 < 0.0] = 0.0\n",
    "    \n",
    "    # Scaling of forcing\n",
    "    if (scaled_forcing == True):\n",
    "        for i in range(len(M_R1)):\n",
    "            if M_R1[i] > 0.0:\n",
    "                dump = np.log(M_R1[i]/8)/np.log(2.0)\n",
    "                M_R1[i] = M_R1[i] * FS_ISSM_JPL_BM08[1]**dump\n",
    "            if M_R2[i] > 0.0:\n",
    "                dump = np.log(M_R2[i]/8)/np.log(2.0)\n",
    "                M_R2[i] = M_R2[i] * FS_ISSM_JPL_BM08[2]**dump\n",
    "            if M_R3[i] > 0.0:\n",
    "                dump = np.log(M_R3[i]/8)/np.log(2.0)\n",
    "                M_R3[i] = M_R3[i] * FS_ISSM_JPL_BM08[3]**dump\n",
    "            if M_R4[i] > 0.0:\n",
    "                dump = np.log(M_R4[i]/8)/np.log(2.0)\n",
    "                M_R4[i] = M_R4[i] * FS_ISSM_JPL_BM08[4]**dump\n",
    "            if M_R5[i] > 0.0:\n",
    "                dump = np.log(M_R5[i]/8)/np.log(2.0)\n",
    "                M_R5[i] = M_R5[i] * FS_ISSM_JPL_BM08[5]**dump\n",
    "\n",
    "    # Linear response\n",
    "    SL = []\n",
    "    SL.append(0)\n",
    "    for t in range(1,tend-tbeg):\n",
    "        dSL = 0\n",
    "        for tp in range(0,t):\n",
    "            dSL = dSL + M_R1[tp]*RF_ISSM_JPL_BM08_R1[t-tp]\n",
    "        SL.append(dSL)\n",
    "    SL_wTd_nos_base_ISSM_JPL_R1_RCP85=np.vstack([SL_wTd_nos_base_ISSM_JPL_R1_RCP85, SL])\n",
    "\n",
    "    SL = []\n",
    "    SL.append(0)\n",
    "    for t in range(1,tend-tbeg):\n",
    "        dSL = 0\n",
    "        for tp in range(0,t):\n",
    "            dSL = dSL + M_R2[tp]*RF_ISSM_JPL_BM08_R2[t-tp]\n",
    "        SL.append(dSL)\n",
    "    SL_wTd_nos_base_ISSM_JPL_R2_RCP85=np.vstack([SL_wTd_nos_base_ISSM_JPL_R2_RCP85, SL])\n",
    "\n",
    "    SL = []\n",
    "    SL.append(0)\n",
    "    for t in range(1,tend-tbeg):\n",
    "        dSL = 0\n",
    "        for tp in range(0,t):\n",
    "            dSL = dSL + M_R3[tp]*RF_ISSM_JPL_BM08_R3[t-tp]\n",
    "        SL.append(dSL)\n",
    "    SL_wTd_nos_base_ISSM_JPL_R3_RCP85=np.vstack([SL_wTd_nos_base_ISSM_JPL_R3_RCP85, SL])\n",
    "\n",
    "    SL = []\n",
    "    SL.append(0)\n",
    "    for t in range(1,tend-tbeg):\n",
    "        dSL = 0\n",
    "        for tp in range(0,t):\n",
    "            dSL = dSL + M_R4[tp]*RF_ISSM_JPL_BM08_R4[t-tp]\n",
    "        SL.append(dSL)\n",
    "    SL_wTd_nos_base_ISSM_JPL_R4_RCP85=np.vstack([SL_wTd_nos_base_ISSM_JPL_R4_RCP85, SL])\n",
    "\n",
    "    SL = []\n",
    "    SL.append(0)\n",
    "    for t in range(1,tend-tbeg):\n",
    "        dSL = 0\n",
    "        for tp in range(0,t):\n",
    "            dSL = dSL + M_R5[tp]*RF_ISSM_JPL_BM08_R5[t-tp]\n",
    "        SL.append(dSL)\n",
    "    SL_wTd_nos_base_ISSM_JPL_R5_RCP85=np.vstack([SL_wTd_nos_base_ISSM_JPL_R5_RCP85, SL])\n",
    "\n",
    "SL_wTd_nos_base_ISSM_JPL_SU_RCP85 = SL_wTd_nos_base_ISSM_JPL_R1_RCP85+SL_wTd_nos_base_ISSM_JPL_R2_RCP85+SL_wTd_nos_base_ISSM_JPL_R3_RCP85+SL_wTd_nos_base_ISSM_JPL_R4_RCP85+SL_wTd_nos_base_ISSM_JPL_R5_RCP85\n",
    "\n",
    "print(countR1/EnsembleSize,10/19)\n",
    "print(countR2/EnsembleSize,17/19)\n",
    "print(countR3/EnsembleSize,10/19)\n",
    "print(countR4/EnsembleSize,9/19)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Time = range(1900,2100)\n",
    "ncfile = nc.Dataset('EnsembleSingleModelProjections/SL_wTd_nos_base_ISSM_JPL_RCP85.nc','w', format='NETCDF4')\n",
    "ncfile.createDimension('Time', None)\n",
    "ncfile.createDimension('Emember', None)\n",
    "\n",
    "SL_wTd_nos_base_R1 = ncfile.createVariable('EAIS', 'f4', ('Emember', 'Time'))\n",
    "SL_wTd_nos_base_R2 = ncfile.createVariable('Ross', 'f4', ('Emember', 'Time'))\n",
    "SL_wTd_nos_base_R3 = ncfile.createVariable('Amundsen', 'f4', ('Emember', 'Time'))\n",
    "SL_wTd_nos_base_R4 = ncfile.createVariable('Weddell', 'f4', ('Emember', 'Time'))\n",
    "SL_wTd_nos_base_R5 = ncfile.createVariable('Peninsula', 'f4', ('Emember', 'Time'))\n",
    "SL_wTd_nos_base_SU = ncfile.createVariable('Antarctica', 'f4', ('Emember', 'Time'))\n",
    "t = ncfile.createVariable('Time', 'i4', 'Time')\n",
    "\n",
    "SL_wTd_nos_base_R1[:,:] = SL_wTd_nos_base_ISSM_JPL_R1_RCP85\n",
    "SL_wTd_nos_base_R2[:,:] = SL_wTd_nos_base_ISSM_JPL_R2_RCP85\n",
    "SL_wTd_nos_base_R3[:,:] = SL_wTd_nos_base_ISSM_JPL_R3_RCP85\n",
    "SL_wTd_nos_base_R4[:,:] = SL_wTd_nos_base_ISSM_JPL_R4_RCP85\n",
    "SL_wTd_nos_base_R5[:,:] = SL_wTd_nos_base_ISSM_JPL_R5_RCP85\n",
    "SL_wTd_nos_base_SU[:,:] = SL_wTd_nos_base_ISSM_JPL_SU_RCP85\n",
    "t[:] = Time\n",
    "\n",
    "SL_wTd_nos_base_R1.units = 'meter'\n",
    "SL_wTd_nos_base_R2.units = 'meter'\n",
    "SL_wTd_nos_base_R3.units = 'meter'\n",
    "SL_wTd_nos_base_R4.units = 'meter'\n",
    "SL_wTd_nos_base_R5.units = 'meter'\n",
    "SL_wTd_nos_base_SU.units = 'meter'\n",
    "\n",
    "t.units = 'years'\n",
    "\n",
    "ncfile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "SL_wTd_nos_base_ISSM_JPL_SU_RCP85_50pc = np.percentile(SL_wTd_nos_base_ISSM_JPL_SU_RCP85, 50, axis=0, out=None, overwrite_input=False, interpolation='linear', keepdims=False)\n",
    "SL_wTd_nos_base_ISSM_JPL_SU_RCP85_83pc = np.percentile(SL_wTd_nos_base_ISSM_JPL_SU_RCP85, 83.33, axis=0, out=None, overwrite_input=False, interpolation='linear', keepdims=False)\n",
    "SL_wTd_nos_base_ISSM_JPL_SU_RCP85_17pc = np.percentile(SL_wTd_nos_base_ISSM_JPL_SU_RCP85, 16.66, axis=0, out=None, overwrite_input=False, interpolation='linear', keepdims=False)\n",
    "SL_wTd_nos_base_ISSM_JPL_SU_RCP85_95pc = np.percentile(SL_wTd_nos_base_ISSM_JPL_SU_RCP85, 5, axis=0, out=None, overwrite_input=False, interpolation='linear', keepdims=False)\n",
    "SL_wTd_nos_base_ISSM_JPL_SU_RCP85_05pc = np.percentile(SL_wTd_nos_base_ISSM_JPL_SU_RCP85, 95, axis=0, out=None, overwrite_input=False, interpolation='linear', keepdims=False)\n",
    "SL_wTd_nos_base_ISSM_JPL_SU_RCP85_99pc = np.percentile(SL_wTd_nos_base_ISSM_JPL_SU_RCP85, 99, axis=0, out=None, overwrite_input=False, interpolation='linear', keepdims=False)\n",
    "SL_wTd_nos_base_ISSM_JPL_SU_RCP85_01pc = np.percentile(SL_wTd_nos_base_ISSM_JPL_SU_RCP85, 1, axis=0, out=None, overwrite_input=False, interpolation='linear', keepdims=False)\n",
    "\n",
    "np.savetxt(\"PercentilesSingleModelProjections/SL_wTd_nos_base_ISSM_JPL_SU_RCP85_50pc.dat\", SL_wTd_nos_base_ISSM_JPL_SU_RCP85_50pc, delimiter=\",\")\n",
    "np.savetxt(\"PercentilesSingleModelProjections/SL_wTd_nos_base_ISSM_JPL_SU_RCP85_83pc.dat\", SL_wTd_nos_base_ISSM_JPL_SU_RCP85_83pc, delimiter=\",\")\n",
    "np.savetxt(\"PercentilesSingleModelProjections/SL_wTd_nos_base_ISSM_JPL_SU_RCP85_17pc.dat\", SL_wTd_nos_base_ISSM_JPL_SU_RCP85_17pc, delimiter=\",\")\n",
    "np.savetxt(\"PercentilesSingleModelProjections/SL_wTd_nos_base_ISSM_JPL_SU_RCP85_05pc.dat\", SL_wTd_nos_base_ISSM_JPL_SU_RCP85_05pc, delimiter=\",\")\n",
    "np.savetxt(\"PercentilesSingleModelProjections/SL_wTd_nos_base_ISSM_JPL_SU_RCP85_95pc.dat\", SL_wTd_nos_base_ISSM_JPL_SU_RCP85_95pc, delimiter=\",\")\n",
    "np.savetxt(\"PercentilesSingleModelProjections/SL_wTd_nos_base_ISSM_JPL_SU_RCP85_01pc.dat\", SL_wTd_nos_base_ISSM_JPL_SU_RCP85_01pc, delimiter=\",\")\n",
    "np.savetxt(\"PercentilesSingleModelProjections/SL_wTd_nos_base_ISSM_JPL_SU_RCP85_99pc.dat\", SL_wTd_nos_base_ISSM_JPL_SU_RCP85_99pc, delimiter=\",\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp3=plt.figure()\n",
    "for i in range(EnsembleSize):\n",
    "    plt.plot(Time,SL_wTd_nos_base_ISSM_JPL_SU_RCP85[i,:])\n",
    "\n",
    "plt.show()\n",
    "fp3.savefig(\"Figures/SL_wTd_nos_base_ISSM_JPL_SU_RCP85_alllines.pdf\", bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}