Merge pull request #233 from sophmaca/main

Adding water isotope exercise to paleo challenge

notebooks/challenge/paleo/exercise_3.ipynb

@@ -0,0 +1,373 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "f406f992-92bd-4b17-9bd3-b99c5c8abaf3",
+   "metadata": {},
+   "source": [
+    "# 3: Water isotope tracers in CESM\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b90d4773-7ca0-4131-ab07-517608a3e976",
+   "metadata": {},
+   "source": [
+    "\n",
+    "<div class=\"alert alert-info\">\n",
+    "<strong>Exercise: Run a preindustrial simulation with water isotope tracers</strong><br><br>\n",
+    "\n",
+    "Download isotope-enabled CESM1.3 (iCESM1.3) code (the version of CESM used in this tutorial does not include water isotope capabilities). \n",
+    "\n",
+    "Create, configure, build and run a fully coupled preindustrial case called ``b.e13.B1850.f19_g17.piControl.001`` following [CESM naming conventions](https://www.cesm.ucar.edu/models/cesm2/naming-conventions) including water isotope tracers. \n",
+    "\n",
+    "Run for 1 year. \n",
+    "\n",
+    "</div>\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "65b2cbda-2d54-48ee-898b-4c391f16ca79",
+   "metadata": {},
+   "source": [
+    "\n",
+    "<div class=\"alert alert-warning\">  \n",
+    "<details>\n",
+    "\n",
+    "<summary> <font face=\"Times New Roman\" color='blue'>Click here for hints</font> </summary>\n",
+    "<br>\n",
+    "\n",
+    "**Where is the code for iCESM1.3 located?**\n",
+    "\n",
+    "- https://github.com/NCAR/iCESM1.3_iHESP_hires \n",
+    "\n",
+    "**What is the resolution for B1850?**\n",
+    "\n",
+    "- Use resolution ``f19_g17`` for fast throughput \n",
+    "\n",
+    "**Which XML variable should you change to tell the model to run for one year?**\n",
+    "\n",
+    "- Use ``STOP_OPTION`` and ``STOP_N`` \n",
+    "\n",
+    "**How to check if each XML variable is modified correctly?**\n",
+    "\n",
+    "- Use ``xmlquery -p`` \n",
+    "\n",
+    "</details>\n",
+    "</div>\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7dd602b7-372d-4f36-b6d1-df8e22ba1646",
+   "metadata": {},
+   "source": [
+    "\n",
+    "<div class=\"alert alert-success\">   \n",
+    "<details>\n",
+    "<summary><font face=\"Times New Roman\" color='blue'>Click here for the solution</font></summary><br>\n",
+    " \n",
+    "**# Download iCESM1.3 code** \n",
+    "\n",
+    "Set environment variables with the commands:\n",
+    "\n",
+    "```\n",
+    "cd /glade/work/$USER/code \n",
+    "git clone https://github.com/NCAR/iCESM1.3_iHESP_hires iCESM1.3_iHESP_hires \n",
+    "cd iCESM1.3_iHESP_hires \n",
+    "./manage_externals/checkout_externals \n",
+    "```\n",
+    "\n",
+    " \n",
+    "**# Set environment variables** \n",
+    "\n",
+    "Set environment variables with the commands:\n",
+    "    \n",
+    "**For tcsh users** \n",
+    "    \n",
+    "```\n",
+    "set CASENAME=b.e13.B1850C5.f19_g16.piControl.001\n",
+    "set CASEDIR=/glade/u/home/$USER/cases/$CASENAME\n",
+    "set RUNDIR=/glade/derecho/scratch/$USER/$CASENAME/run\n",
+    "set COMPSET=B1850C5\n",
+    "set RESOLUTION=f19_g16\n",
+    "set PROJECT=UESM0013\n",
+    "```\n",
+    "\n",
+    "Note: You should use the project number given for this tutorial.\n",
+    "\n",
+    "**For bash users** \n",
+    "    \n",
+    "```\n",
+    "export CASENAME=b.e13.B1850C5.f19_g16.piControl.001\n",
+    "export CASEDIR=/glade/u/home/$USER/cases/$CASENAME\n",
+    "export RUNDIR=/glade/derecho/scratch/$USER/$CASENAME/run\n",
+    "export COMPSET=B1850C5\n",
+    "export RESOLUTION=f19_g16\n",
+    "export PROJECT=UESM0013\n",
+    "```\n",
+    "\n",
+    "Note: You should use the project number given for this tutorial.\n",
+    "\n",
+    "**# Make a case directory**\n",
+    "\n",
+    "If needed create a directory `cases` into your home directory:\n",
+    "    \n",
+    "```\n",
+    "mkdir /glade/u/home/$USER/cases/\n",
+    "```\n",
+    "    \n",
+    "\n",
+    "**# Create a new case**\n",
+    "\n",
+    "Create a new case with the command ``create_newcase``:\n",
+    "```\n",
+    "cd /glade/work/$USER/code/iCESM1.3_iHESP_hires/cime/scripts/\n",
+    "./create_newcase --case $CASEDIR --res $RESOLUTION --compset $COMPSET --project $PROJECT --run-unsupported \n",
+    "```\n",
+    "\n",
+    "**# Change the job queue**\n",
+    "\n",
+    "If needed, change ``job queue``.<br>\n",
+    "For instance, to run in the queue ``main``.\n",
+    "```  \n",
+    "cd $CASEDIR\n",
+    "./xmlchange JOB_QUEUE=main\n",
+    "```\n",
+    "This step can be redone at anytime in the process. \n",
+    "\n",
+    "**# Setup**\n",
+    "\n",
+    "Invoke ``case.setup`` with the command:\n",
+    "```    \n",
+    "cd $CASEDIR\n",
+    "./case.setup    \n",
+    "```  \n",
+    "\n",
+    "You build the namelists with the command:\n",
+    "```\n",
+    "./preview_namelists\n",
+    "```\n",
+    "This step is optional as the script ``preview_namelists`` is automatically called by ``case.build`` and ``case.submit``. But it is nice to check that your changes made their way into:\n",
+    "```\n",
+    "$CASEDIR/CaseDocs/atm_in\n",
+    "```\n",
+    "\n",
+    "\n",
+    "**# Set run length**\n",
+    "\n", 
+    "```\n",
+    "./xmlchange STOP_N=1,STOP_OPTION=nyears\n",
+    "```\n",
+    "\n",
+    "**# Build the run**\n",
+    "\n",
+    "```\n",
+    "qcmd -A $PROJECT -- ./case.build\n",
+    "```\n",
+    "\n",
+    "**# Which namelist variables enable water isotope tracers?**\n",
+    "\n",
+    "- Notice that the steps to set up this isotope-enabled preindustrial simulation are very similar to a preindustrial simulation without isotopes (e.g., Paleo Exercise 1) \n",
+    "- In iCESM1.3, it is assumed you will run with water isotope tracers so each compset include isotope settings by default \n",
+    "- Use ``./xmlquery`` to explore how ``FLDS_WISO``, ``CAM_CONFIG_OPTS``, and ``OCN_TRACER_MODULES`` differ between this iCESM1.3 case and that of Exercise 1 \n",
+    "- Also, take a look at namelist settings for each CESM component with variables contain ``wiso`` in ``$CASEDIR/Buildconf`` \n",
+    "\n",
+    "\n",
+    "**# Submit the run**\n",
+    "\n",
+    "```\n",
+    "./case.submit\n",
+    "```\n",
+    "------------\n",
+    "\n",
+    "**# Check on your run**\n",
+    "\n",
+    "\n",
+    "After submitting the job, use ``qstat -u $USER`` to check the status of your job. \n",
+    "It may take ~16 minutes to finish the one-year simulation. \n",
+    "\n",
+    "**# Check your solution**\n",
+    "\n",
+    "When the run is completed, look at the history files into the archive directory.  \n",
+    "    \n",
+    "(1) Check that your archive directory on derecho (The path will be different on other machines): \n",
+    "```\n",
+    "cd /glade/derecho/scratch/$USER/archive/$CASENAME/atm/hist\n",
+    "ls \n",
+    "```\n",
+    "\n",
+    "As your run is one-year, there should be 12 monthly files (``h0``) for each model component. \n",
+    "\n",
+    "\n",
+    "Success! Let's plot the results. \n",
+    "\n",
+    "</details>\n",
+    "</div>\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "65b2cbda-2d54-48ee-898b-4c391f16ca79",
+   "metadata": {},
+   "source": [
+    "\n",
+    "<div class=\"alert alert-warning\">  \n",
+    "<details>\n",
+    "\n",
+    "<summary> <font face=\"Times New Roman\" color='blue'>Click here to visualize results</font> </summary>\n",
+    "<br>\n",
+    "\n",
+    "##Option 1: \n",
+    "\n",
+    "**# Use NCO to calculate the oxygen isotopic composition of precipitation**\n",
+    "\n",
+    "The ratio of heavy (<sup>18</sup>O) to light (<sup>16</sup>O) isotopes are most commonly expressed relative to a standard in delta (δ) notation: \n",
+    "\n",
+    "####δ<sup>18</sup>O = (R<sub>sample</sub> - R<sub>std</sub>)/(R<sub>std</sub>) * 1000‰ \n",
+    "\n",
+    "where \n",
+    "- R<sub>sample</sub> = ratio of <sup>18</sup>O to <sup>16</sup>O in sample \n",
+    "- R<sub>std</sub> = ratio of <sup>18</sup>O to <sup>16</sup>O in a standard \n",
+    "\n",
+    "Thus, the δ<sup>18</sup>O of a sample which is identical to the standard would be 0‰, positive values indicate a greater proportion of <sup>18</sup>O than the standard, and negative values indicate a lower proportion of <sup>18</sup>O. \n",
+    "\n",
+    "In isotope-enabled CESM, the relative abundances of <sup>16</sup>O and <sup>18</sup>O are already adjusted to their naturally occurring global abundances (99.757% and 0.205%, respectively), so we do not include R<sub>std</sub> in the calculation of δ<sup>18</sup>O. Rather, isotope variables in CESM are expressed in delta (δ) notation as: \n",
+    "\n",
+    "\n",
+    "####δ<sup>18</sup>O = ((PRECRC_H218Or + PRECSC_H218Os + PRECRL_H218OR + PRECSL_H218OS)/ \n",
+    "####(PRECRC_H216Or + PRECSC_H216Os + PRECRL_H216OR + PRECSL_H216OS) - 1) * 1000‰ \n",
+    "\n",
+    "\n",
+    "- Use ``ncdump /glade/derecho/scratch/$USER/$CASENAME/atm/hist/$CASENAME.cam.h0.0001-01.nc | less`` to define each isotope variable above \n",
+    "- For example, search for ``PRECRC_H218Or`` using ``/PRECRC_H218Or + <enter>``) \n",
+    "\n",
+    "To calculate the δ<sup>18</sup>O of precipitation from the simulation using NCO, \n",
+    "```\n",
+    "cd /glade/derecho/scratch/$USER/archive/$CASENAME/atm/hist\n",
+    "ncap2 -s 'd18Op=((PRECRC_H218Or+PRECSC_H218Os+PRECRL_H218OR+PRECSL_H218OS)/(PRECRC_H216Or+PRECSC_H216Os+PRECRL_H216OR+PRECSL_H216OS) - 1)*1000.' -v $CASENAME.cam.h0.0001-12.nc d18Op.$CASENAME.cam.h0.0001-12.nc \n",
+    "```\n",
+    "\n",
+    "**# Use Ncview to visualize precipitation δ<sup>18</sup>O**\n",
+    "\n",
+    "Earth's orbital configuration influences incoming solar insolation.\n",
+    "Take a look at the ``d18Op`` variable we calculated for 1 month in the pre-industrial run.\n",
+    "```  \n",
+    "module load ncview\n",
+    "ncview d18Op.$CASENAME.cam.h0.0001-12.nc \n",
+    "```\n",
+    "\n",
+    "##Option 2: \n",
+    "\n",
+    "**# Use Python to calculate and plot the oxygen isotopic composition of precipitation**\n",
+    "\n",
+    "The following Python code will produce a plot of precipitation δ<sup>18</sup>O for 1 month. \n",
+    "\n",
+    "```  \n",
+    "import xarray as xr \n",
+    "import numpy as np \n",
+    "import matplotlib.pyplot as plt \n",
+    "import cartopy.crs as ccrs \n",
+    "from cartopy.util import add_cyclic_point \n",
+    "\n",
+    "def calculate_d18Op(ds): \n",
+    "    # Compute precipitation δ18O with iCESM output \n",
+    "    \n",
+    "    # Parameters \n",
+    "    # ds: xarray.Dataset contains necessary variables \n",
+    "    \n",
+    "    # Returns \n",
+    "    # ds: xarray.Dataset with δ18O added \n",
+    "    \n",
+    "    # convective & large-scale rain and snow, respectively \n",
+    "    p16O = ds.PRECRC_H216Or + ds.PRECSC_H216Os + ds.PRECRL_H216OR + ds.PRECSL_H216OS \n",
+    "    p18O = ds.PRECRC_H218Or + ds.PRECSC_H218Os + ds.PRECRL_H218OR + ds.PRECSL_H218OS \n",
+    "    \n",
+    "    # avoid dividing by small number here \n",
+    "    p18O = p18O.where(p16O > 1.E-18, 1.E-18) \n",
+    "    p16O = p16O.where(p16O > 1.E-18, 1.E-18) \n",
+    "    d18O = (p18O / p16O - 1.0) * 1000.0 \n",
+    "    \n",
+    "    ds['p16O'] = p16O \n",
+    "    ds['p18O'] = p18O \n",
+    "    ds['d18O'] = d18O \n",
+    "    return ds \n",
+    "\n",
+    "# Read in monthly file \n",
+    "case = 'b.e13.B1850C5.f19_g16.piControl.001' \n",
+    "file = '/glade/derecho/scratch/macarew/archive/'+case+'/atm/hist/'+case+'.cam.h0.0001-12.nc' \n",
+    "ds = xr.open_mfdataset(file, parallel=True, \n",
+    "                       data_vars='minimal', \n",
+    "                       coords='minimal', \n",
+    "                       compat='override') \n",
+    "\n",
+    "# Call function to add preciptation d18O to dataset \n",
+    "ds = calculate_d18Op(ds) \n",
+    "\n",
+    "fig, ax = plt.subplots( \n",
+    "    nrows=1, ncols=1, \n",
+    "    figsize=(6, 2), \n",
+    "    subplot_kw={'projection': ccrs.Robinson(central_longitude=210)}, \n",
+    "    constrained_layout=True) \n",
+    "\n",
+    "d18O_new, lon_new = add_cyclic_point(ds.d18O[0,:,:], ds.lon) \n",
+    "\n",
+    "# Plot model results using contourf \n",
+    "p0 = ax.contourf(lon_new, ds.lat, d18O_new, \n",
+    "                 levels=np.linspace(-30, 0, 16), extend='both', \n",
+    "                 transform=ccrs.PlateCarree()) \n",
+    "\n",
+    "plt.colorbar(p0, ax=ax) \n",
+    "ax.set_title('Dec δ18Op of PI') \n",
+    "ax.coastlines(linewidth=0.5) \n",
+    "```\n",
+    "\n",
+    "\n",
+    "**# Questions for reflection:**\n",
+    "- Do you notice any spatial patterns in precipitation δ<sup>18</sup>O? \n",
+    "\n",
+    "</details>\n",
+    "</div>\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "472131c7-88f9-4863-a2bc-d7364333542d",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "815be0bc-515a-474b-a3dd-b7ba02831b9a",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "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.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}