{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[],"authorship_tag":"ABX9TyO7MUrXOH+ThVphA/JbMRp4"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["# Spirale d'Ekman océanique\n","$$\\left\\{\\begin{matrix} \\displaystyle {\\partial u \\over \\partial t} - f v =-  {\\partial \\over \\partial z} \\left( K {\\partial u \\over \\partial z}\\right) \\cr \\displaystyle {\\partial v \\over \\partial t} + f u= {\\partial \\over \\partial z} \\left( K {\\partial v \\over \\partial z}\\right)\\end{matrix} \\right.$$\n","$$K= \\nu + l^2 \\sqrt{\\left({\\partial \\tilde u \\over \\partial z}\\right)^2+\\left({\\partial \\tilde v \\over \\partial z}\\right)^2}\\quad \\hbox{avec} \\quad l(z) = { \\kappa z \\over 1 + \\kappa z / \\lambda}$$\n"," $\\kappa = 0,41$,  $\\displaystyle \\lambda = 2,7 \\; 10^{-4} {u_g\\over f}$ et $\\nu=0,1$ m$^2$.s$^{-1}$\n","\n"," $K(-z_0) {\\partial u\\over \\partial z}(-z_0) = u_*^2$, ${\\partial u\\over \\partial z}(-z_0)= 0$ et $u(-H)=v(-H) =0$ "],"metadata":{"id":"rr4NAG5vYdKl"}},{"cell_type":"code","execution_count":null,"metadata":{"id":"5D857EmiXmfP"},"outputs":[],"source":["import numpy as np\n","import matplotlib.pyplot as plt\n","from PIL import Image\n","import glob\n","\n","# Sous programmes \n","def inifig(xaxe=0,yaxe=0,xlab='x',ylab='y'):\n","    plt.figure(2)\n","    plt.axvline(xaxe)\n","    plt.axhline(yaxe)\n","    plt.xticks(fontsize=12)\n","    plt.yticks(fontsize=12)\n","    plt.xlabel(xlab,fontsize=16 )\n","    plt.ylabel(ylab,fontsize=16)\n","    \n","def zfi(x,le=3):\n","    miss=le-len(str(x))\n","    a='0'*miss+str(x)\n","    return a\n","    \n","def Fuv(u,v):\n","    dudz=(u[1:]-u[:-1])/dz\n","    dvdz=(v[1:]-v[:-1])/dz\n","    K=nu+l**2*np.sqrt(dudz**2+dvdz**2) \n","    Kdudz=K*dudz; Kdvdz=K*dvdz\n","    Fu=f*v[1:-1]+(Kdudz[1:]-Kdudz[:-1])/dz\n","    Fv=-f*u[1:-1]+(Kdvdz[1:]-Kdvdz[:-1])/dz\n","    return Fu,Fv,K\n","    \n","def advance(u,v):\n","    for n in range(Nt):\n","        Fu,Fv,K=Fuv(u,v)\n","        u[-1]=u[-2]+usc/K[-1]*dz\n","        v[-1]=v[-2]\n","        u[1:-1]=u[1:-1]+dt*Fu\n","        v[1:-1]=v[1:-1]+dt*Fv\n","        if (n-1)%nplot==0:\n","            t=round(n*dt/3600)\n","            plt.xlabel(r'$v$ et $u$',fontsize=16 )\n","            plt.ylabel(r'Altitude $z$',fontsize=16)    \n","            plt.xlim(xmin,xmax)\n","            plt.ylim(ymin,ymax)          \n","            plt.title(\"Vitesses $u$ et $v$ pour t=\"+zfi(t,2)+\" h\",fontsize=16)\n","            plt.plot(u,z,linewidth=3,color='green')\n","            plt.plot(v,z,linewidth=3,color='green')\n","            plt.plot(ul,z,linewidth=3,color='blue')\n","            plt.plot(vl,z,linewidth=3,color='blue')\n","            #plt.savefig(f\"{n}.png\")\n","            plt.show()\n","    return u,v\n","\n","for j in np.arange(0.05,0.16,0.01): #on fait une boucle pour faire varier us\n","    \n","    k=0.41; f=1.e-4;  nu=.1\n","    us=j; usc=us**2 \n","    H=500; z0=-10; K0=.3;\n","    gam=np.sqrt(f/(2*K0)); \n","    u0=usc/np.sqrt(f*K0)\n","    # Space\n","    Nhour=1000\n","    Nz=100; usdt=1; dt=1/usdt; Nt=Nhour*3600*usdt\n","    z=np.linspace(-H,z0,Nz+2); \n","    dz=H/Nz; \n","    zh=.5*(z[1:]+z[:-1])\n","    \n","    # Parametrization of l\n","    # Journal of the Earth Simulator, Volume 6, October 2006, 3–15\n","    la=2.7e-4*u0/f\n","    l=la*k*zh/(la-k*zh)\n","    \n","    F=False; T=True\n","    # Figure \n","    if T:\n","        plt.xlabel(r'$l$',fontsize=16 )\n","        plt.ylabel(r'$z$',fontsize=16) \n","        plt.title(r\"Longueur de mélange $l(z)$\",fontsize=16)\n","        plt.plot(l,zh,linewidth=3,color='black')\n","        #plt.savefig(\"lOcean.pdf\")\n","        #plt.savefig(\"lOcean.png\")\n","        plt.show()\n","    \n","    # Temporal\n","    nplot=5*3600*usdt;\n","    \n","    # laminar solution\n","    ul=u0*np.exp(gam*(z-z0))*np.cos(gam*(z-z0)-np.pi/4)\n","    vl=u0*np.exp(gam*(z-z0))*np.sin(gam*(z-z0)-np.pi/4)\n","    uls=ul[-1]; vls=vl[-1]\n","    xmin=-1.5*uls; xmax=1.5*uls; \n","    ymin=-H; ymax=0;\n","    \n","    # initial condition\n","    if T:\n","        u=0*z\n","        v=0*z\n","    else:\n","        u=ul.copy()\n","        v=vl.copy()\n","    \n","    # K initial\n","    if T:\n","        plt.xlabel(r'$K$',fontsize=16 )\n","        plt.ylabel(r'$z$',fontsize=16) \n","        plt.title(r\"Viscosité turbulente $K(z)$\",fontsize=16)\n","        dudz=(ul[1:]-ul[:-1])/dz\n","        dvdz=(vl[1:]-vl[:-1])/dz\n","        K=nu+l**2*np.sqrt(dudz**2+dvdz**2) \n","        plt.plot(K,zh,linewidth=3,color='magenta')\n","        plt.plot(K0+0*K,zh,linewidth=3,color='black')\n","        plt.show()\n","    \n","    # Loop\n","    if T: \n","        u,v=advance(u,v)\n","        \n","        # final \n","        t=round(Nt*dt/3600);\n","        plt.title(\"Vitesses $u$ et $v$ pour t=\"+zfi(t,2)+\" h\",fontsize=16)\n","           \n","        # avec longueur de melange\n","        plt.xlim(xmin,xmax)\n","        plt.ylim(ymin,ymax)          \n","        plt.plot(u,z,linewidth=3,color='red')\n","        plt.plot(v,z,linewidth=3,color='red')\n","        \n","        # laminar\n","        plt.xlim(xmin,xmax)\n","        plt.ylim(ymin,ymax)          \n","        plt.plot(ul,z,linewidth=3,color='blue')\n","        plt.plot(vl,z,linewidth=3,color='blue')\n","        #plt.savefig(\"uv.png\")\n","        plt.show()\n","    \n","    \n","    # spiral\n","    if T:\n","        #inifig(xlab='u',ylab='v')  \n","        plt.xlabel(r'$u$',fontsize=16 )\n","        plt.ylabel(r'$v$',fontsize=16) \n","        plt.xlim(-uls,uls*2)\n","        plt.ylim(-uls*2,uls)          \n","        plt.title(f\"Spirales d'Ekman océaniques, $u_*$ = {j}\",fontsize=16)\n","        plt.plot(ul,vl,linewidth=3,color='blue')\n","        plt.plot(u,v,linewidth=3,color='red')\n","        plt.scatter(uls,vls,marker='o',color='blue',s=100)\n","        plt.savefig(f\"spiral {j}.png\")\n","        plt.show()\n","    \n","    # K\n","    if T:\n","        plt.xlabel(r'$K$',fontsize=16 )\n","        plt.ylabel(r'$z$',fontsize=16) \n","        plt.title(r\"Viscosité turbulente $K(z)$\",fontsize=16)\n","        dudz=(u[1:]-u[:-1])/dz\n","        dvdz=(v[1:]-v[:-1])/dz\n","        K=nu+l**2*np.sqrt(dudz**2+dvdz**2) \n","        plt.plot(K,zh,linewidth=3,color='magenta')\n","        plt.plot(K0+0*K,zh,linewidth=3,color='black')\n","        #plt.savefig(\"K.png\")\n","        plt.show()\n","   \n","\n","#Pour créer le gif des vitesses   \n","#frames = []\n","#imgs = glob.glob(\"*.png\")\n","#for i in imgs:\n"," #   new_frame = Image.open(i)\n","  #  frames.append(new_frame)\n","\n","# Save the png images into a GIF file that loops forever\n","#frames[0].save('png_to_gif.gif', format='GIF',\n"," #             append_images=frames[1:],\n","  #            save_all=True,\n","   #           duration=300, loop=0)"]}]}