Bathymetry.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import xgcm\n",
    "from matplotlib import pyplot as plt\n",
    "import nc_time_axis\n",
    "from pathlib import Path\n",
    "import xnemogcm as xn\n",
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path   = \"/data/dkamm/nemo_output/DINO/Stratification\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datadir = Path(path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "domcfg   = xn.open_domain_cfg(datadir=datadir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reproducing the bathymetry of DINO in python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cosh_bathy(pdistLam: float, pminlam:float, pmaxlam: float, plam: float) -> float:\n",
    "    \"\"\"        \n",
    "    pdistLam    :  Length scale of the slope   [degrees]\n",
    "    pmaxLam     :  Maximum longitude/latitude  [degrees]\n",
    "    pminLam     :  Minimum longitude/latitude  [degrees]\n",
    "    pLam        :  Value of longitude/latitude [degrees]\n",
    "    pcosh       :  Resulting tapering factor\n",
    "    \"\"\"\n",
    "    znorm = 1 + np.exp((pminlam - pmaxlam) / pdistLam)\n",
    "    pcosh = 1 - (\n",
    "        np.exp((plam - pmaxlam) / pdistLam)\n",
    "        + np.exp(-(plam - pminlam) / pdistLam)\n",
    "    ) / znorm\n",
    "    return(pcosh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss_ring(\n",
    "        pdistLam:   float,\n",
    "        plam0:      float,\n",
    "        pphi0:      float,\n",
    "        prad:       float,\n",
    "        plam:       xr.DataArray,\n",
    "        pphi:       xr.DataArray,\n",
    "        pdep_top:   float,\n",
    "        pdep_bot:   xr.DataArray\n",
    "    ) -> float:\n",
    "    \"\"\"        \n",
    "    pdistLam    : Length scale of the slope               [degrees]\n",
    "    plam0       : Longitude of the center                 [degrees]\n",
    "    pphi0       : Latitude of the center                  [degrees]\n",
    "    prad        : radius of the circle                    [degrees]\n",
    "    pLam        : Longitude                               [degrees]\n",
    "    pPhi        : Latitude                                [degrees]\n",
    "    pdep_top    : depth of the gaussian bump              [meters]\n",
    "    pdep_bot    : depth of the basin                      [meters]\n",
    "    \"\"\"\n",
    "    zx = plam - plam0\n",
    "    zy = (pphi - pphi0)\n",
    "    zexp = (- zx**2 - zy**2 + 2 * prad * np.sqrt(zx**2 + zy**2) - prad**2) / pdistLam**2\n",
    "    mask = (pdep_bot >= pdep_top) * 1\n",
    "    pring   = mask * ((pdep_top - pdep_bot) * np.exp( zexp ) + pdep_bot) + (1. - mask) * pdep_bot\n",
    "    return(pring)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def xgauss_ring(\n",
    "        pdistLam:   float,\n",
    "        plam0:      float,\n",
    "        pphi0:      float,\n",
    "        prad:       float,\n",
    "        pdep_bot:   float,\n",
    "        pdep_top:   float,\n",
    "        plam:       float,\n",
    "        pphi:       float\n",
    "    ) -> float:\n",
    "    return(xr.apply_ufunc(gauss_ring, pdistLam, plam0, pphi0, prad, pdep_bot, pdep_top, plam, pphi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cosh-shaped boundary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zminlam =   domcfg.glamt.min()\n",
    "zmaxlam =   domcfg.glamt.max()\n",
    "zminphi =   domcfg.gphit.min()\n",
    "zmaxphi =   domcfg.gphit.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zdistlam = 3.0\n",
    "zdistphi = np.cos(2 * np.pi * zmaxphi / 360.) * zdistlam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zx = cosh_bathy(zdistlam, zminlam, zmaxlam, domcfg.glamt)\n",
    "zy = cosh_bathy(zdistphi, zminphi, zmaxphi, domcfg.gphit)\n",
    "pbathy = zx * zy * ( 4000. - 2000. ) + 2000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1,3,figsize=(15,6))\n",
    "zy.plot(ax=axs[0])\n",
    "(1 - zx).plot(ax=axs[1])\n",
    "(zx * zy).plot(ax=axs[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pbathy.plot(x='glamt', y='gphit', figsize=(6,8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flattening of the bathymetry in the channel with cosh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zslp_cha = 1.5\n",
    "zcha_min = -70.\n",
    "zcha_max = -50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ztaper = cosh_bathy(zslp_cha, zcha_min, zcha_max, domcfg.gphit.where(domcfg.gphit < zcha_max)).fillna(0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pbathy_flat1 = ztaper * pbathy.max() + pbathy * (1. - ztaper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pbathy_flat1.plot(x='glamt', y='gphit', figsize=(6,8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flattening of the bathymetry in the channel with gaussian ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zds_width   = 4.\n",
    "zlam0       = (zcha_max + zcha_min) / 2\n",
    "zrad        = (zcha_max - zcha_min) / 2\n",
    "zdep_top    = 3000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pbathy_n = gauss_ring(zds_width, zminlam, zlam0, zrad, domcfg.glamt, domcfg.gphit, zdep_top, pbathy_flat1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1,3, sharey=True, figsize=(20,10))\n",
    "pbathy_n.isel(x_c=0,y_c=slice(0,70)).plot(ax=axs[0], label='0')\n",
    "pbathy_n.isel(x_c=1,y_c=slice(0,70)).plot(ax=axs[0], label='1')\n",
    "pbathy_n.isel(x_c=2,y_c=slice(0,70)).plot(ax=axs[0], label='2')\n",
    "pbathy_n.isel(x_c=3,y_c=slice(0,70)).plot(ax=axs[0], label='3')\n",
    "pbathy_n.isel(x_c=4,y_c=slice(0,70)).plot(ax=axs[0], label='4')\n",
    "axs[0].legend()\n",
    "pbathy_flat1.isel(x_c=0,y_c=slice(0,70)).plot(ax=axs[1], label='0')\n",
    "pbathy_flat1.isel(x_c=1,y_c=slice(0,70)).plot(ax=axs[1], label='1')\n",
    "pbathy_flat1.isel(x_c=2,y_c=slice(0,70)).plot(ax=axs[1], label='2')\n",
    "pbathy_flat1.isel(x_c=3,y_c=slice(0,70)).plot(ax=axs[1], label='3')\n",
    "pbathy_flat1.isel(x_c=4,y_c=slice(0,70)).plot(ax=axs[1], label='4')\n",
    "axs[1].legend()\n",
    "\n",
    "pbathy_flat1.isel(x_c=0,y_c=slice(0,70)).plot(ax=axs[2], label='0')\n",
    "pbathy_n.isel(x_c=0,y_c=slice(0,70)).plot(ax=axs[2], label='0')\n",
    "axs[2].grid('fine')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear transition from Scottia arc to western boundary (Periodic boundary conditions) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "taper = (pbathy_n.x_c > int(zdistlam))*1.0 + (pbathy_n.x_c <= int(zdistlam)) * pbathy_n.x_c / zdistlam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bathy = taper * pbathy_n + (1 - taper) * pbathy_flat1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(6,8))\n",
    "bathy.plot(x='glamt', y='gphit', ax=ax)\n",
    "ax.set_aspect('equal', adjustable='box')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tampering the bathymetry to reach the bottom in the channel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tampered_exponential(x, wall_x, wall_left, width, zdistlam, taper_dist, taper_ratio): \n",
    "    if wall_left:\n",
    "        ones    = (x >= wall_x + taper_dist) * 1.0 + (x < wall_x) * 1.0\n",
    "        sigmoid = ((x <  wall_x + taper_dist) * (x > wall_x)) * 1.0\n",
    "        sign    = 1.0\n",
    "        mask    = (x < wall_x) * 1.0\n",
    "    else:\n",
    "        ones    = (x <= wall_x - taper_dist) * 1.0 + (x > wall_x) * 1.0\n",
    "        sigmoid = ((x >  wall_x - taper_dist) * (x < wall_x)) * 1.0\n",
    "        sign    = - 1.0\n",
    "        mask    = (x > wall_x) * 1.0\n",
    "        \n",
    "    inflec  = wall_x + sign * taper_ratio * taper_dist\n",
    "    step = xr.zeros_like(x) + ones + sigmoid * 1 / (1 + np.exp( - sign * ( x - inflec) / 0.5))\n",
    "\n",
    "    znorm   = 1 + np.exp(-( width) / zdistlam)\n",
    "    exp     = 1 - (np.exp(- sign * (x - wall_x) / zdistlam)) / znorm\n",
    "\n",
    "    return (exp * (1 - step) + step) - mask\n",
    "    #return step, exp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "channel1 = tampered_exponential(domcfg.gphit, zcha_min, True, 140., zdistlam, 9., 0.8)\n",
    "channel2 = tampered_exponential(domcfg.gphit, zcha_max, False, 140., zdistlam, 9., 0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wall1 = tampered_exponential(domcfg.gphit, zminphi, True, 140., zdistlam, 9., 0.8)\n",
    "wall2 = tampered_exponential(domcfg.gphit, zmaxphi, False, 140., zdistlam, 9., 0.8)\n",
    "wall3 = tampered_exponential(domcfg.glamt, zminlam, True, 140., zdistlam, 9., 0.8)\n",
    "wall4 = tampered_exponential(domcfg.glamt, zmaxlam, False, 140., zdistlam, 9., 0.8)\n",
    "\n",
    "channel1 = tampered_exponential(domcfg.gphit, zcha_min, True, 140., zdistlam, 9., 0.8)\n",
    "channel2 = tampered_exponential(domcfg.gphit, zcha_max, False, 140., zdistlam, 9., 0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1,3, figsize=(15, 6))\n",
    "(wall1 * wall2).isel(x_c=30).plot(ax=axs[0])\n",
    "(wall3 * wall4).isel(y_c=100).plot(ax=axs[1])\n",
    "(channel2).isel(x_c=30).plot(ax=axs[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wall3_new = channel2 + wall3 * (1. - channel2)\n",
    "wall4_new = channel2 + wall4 * (1. - channel2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "taper = wall1 * wall2 * wall3_new * wall4_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bathy = taper * ( 4000. - 2000. ) + 2000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(6,8))\n",
    "bathy.plot(x='glamt', y='gphit', ax=ax)\n",
    "ax.set_aspect('equal', adjustable='box')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bathy_ring = gauss_ring(zds_width, zminlam, zlam0, zrad, domcfg.glamt, domcfg.gphit, zdep_top, bathy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "taper = (bathy_ring.x_c > int(zdistlam))*1.0 + (bathy_ring.x_c <= int(zdistlam)) * bathy_ring.x_c / zdistlam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dino = taper * bathy_ring + (1 - taper) * bathy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1, figsize=(6,8))\n",
    "dino.plot(x='glamt', y='gphit', ax=ax)\n",
    "ax.set_aspect('equal', adjustable='box')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "nemo_env",
   "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.9.15"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}