<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Atlas stabilité calottes Mars
Rappel 
Step1:  Carte en moyenne temporelle sur la totalité de l'expérience
Step2:  Carte en moyenne temporelle de $p_{sat}$ pour $H_2O$
Step3:  Carte à $L_s$ donné de $p_{sat}$ pour $H_2O$
Step4:  Carte en moyenne temporelle de $p_{sat}$ pour $CO_2$
<ASSISTANT_TASK:>
Python Code:
filename = 'resultat.nc'
import numpy as np
import matplotlib.pyplot as plt
from pylab import *
import cartopy.crs as ccrs
from netCDF4 import Dataset
%matplotlib inline
import warnings
warnings.filterwarnings('ignore')
data = Dataset(filename)
longitude=data.variables['longitude'][:]
latitude=data.variables['latitude'][:]
altitude=data.variables['altitude'][:]
Time = data.variables['Time'][:]
Ls = data.variables['Ls'][:]
dafirst = Time[0]
daint = Time[1] - dafirst
dalast = dafirst + (len(Time)-1)*daint
year = 0.
add = np.linspace(dafirst,dalast,num=len(Time)) ; add[0] = 0.
for iii in range(1,len(Ls)):
    if Ls[iii] - Ls[iii-1] < 0: year = year+1.
    add[iii] = year*360.
Ls_true = add + Ls
Explanation: Atlas stabilité calottes Mars
Rappel : pour enregistrer une figure, placer cette ligne après la figure en question :
<code> fig.savefig('figure.png') </code>
Générer ensuite la figure avec la cellule correspondante, puis exécuter cette ligne. La figure sera sauvegardée dans le dossier de la simulation correspondante. Vous pouvez bien sûr changer le nom ou l'extension (JPG, PDF).
End of explanation
# Paramètres utilisateurs -----------------------------------------
earthtopo = False # ajouter les traits de côte actuels
varname = 'tsurf'
vmin = 120
vmax = 280
# Code ------------------------------------------------------------
dataplt = data.variables[varname][:,:,:]
fig = plt.figure(figsize=(12,8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
if (earthtopo): ax.coastlines(resolution="110m",linewidth=1)
gl = ax.gridlines(linestyle='--',color='black',
             draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False
clevs = np.linspace(vmin,vmax,29)
plt.contourf(longitude, latitude, np.mean(dataplt[:,:,:],axis=0),
             clevs, transform=ccrs.PlateCarree(),cmap="jet")
plt.title(r"Température de surface moyenne", size=14)
cb = plt.colorbar(ax=ax, orientation="vertical", pad=0.02, aspect=16, shrink=0.8)
cb.set_label(r'K',size=12,rotation=0,labelpad=15)
cb.ax.tick_params(labelsize=10)
plt.show()
Explanation: Carte en moyenne temporelle sur la totalité de l'expérience
End of explanation
def psatw(temp):
  # METHOD GOFF GRATCH (HygroLP) - OVER WATER
  # -----------------------------------------
  log10ew = -7.90298*(373.16/temp-1) \
    + 5.02808 * np.log10(373.16/temp) \
    - 1.3816e-7 * (10**(11.344 * (1-temp/373.16))-1) \
    + 8.1328e-3 * (10**(-3.49149 *(373.16/temp-1))-1) \
    + np.log10(1013.246)
  return 100 * (10**(log10ew))
def psati(temp):
  # METHOD GOFF GRATCH (HygroLP) - OVER ICE
  # ---------------------------------------
  log10ei = -9.09718*(273.16/temp-1) \
    - 3.56654*np.log10(273.16/temp) \
    + 0.876793*(1-temp/273.16) \
    + np.log10(6.1071)
  return  100 * (10**(log10ei))
tzero = 273.15
temp = np.linspace(-80+tzero,tzero,81)
plt.yscale('log')
plt.plot(temp,psatw(temp))
plt.plot(temp,psati(temp))
plt.show()
# Paramètres utilisateurs -----------------------------------------
earthtopo = False # ajouter les traits de côte actuels
ph2oatmo = 0.05e-2*610. # assumed mean water vapor partial pressure
vmin = 0.
vmax = 1.
# Code ------------------------------------------------------------
tsurfnc = data.variables['tsurf'][:,:,:]
dataplt = ph2oatmo/psati(tsurfnc)
fig = plt.figure(figsize=(12,8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
if (earthtopo): ax.coastlines(resolution="110m",linewidth=1)
gl = ax.gridlines(linestyle='--',color='black',
             draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False
clevs = np.linspace(vmin,vmax,21)
plt.contourf(longitude, latitude, np.mean(dataplt[:,:,:],axis=0),
             clevs, transform=ccrs.PlateCarree(),cmap="jet")
plt.title(r"Saturation ratio", size=14)
cb = plt.colorbar(ax=ax, orientation="vertical", pad=0.02, aspect=16, shrink=0.8)
cb.set_label(r'Pa',size=12,rotation=0,labelpad=15)
cb.ax.tick_params(labelsize=10)
plt.show()
Explanation: Carte en moyenne temporelle de $p_{sat}$ pour $H_2O$
End of explanation
# Paramètres utilisateurs -----------------------------------------
earthtopo = False # ajouter les traits de côte actuels
year_user = 1 # année de simulation à regarder
Ls_user = 90. # longitude solaire choisie
varname = 'tsurf'
vmin = 80.
vmax = 280.
# Code ------------------------------------------------------------
Ls_true_user = year_user*360. + Ls_user
Ls_ind = np.where(abs(Ls_true-Ls_true_user)==
                  abs(Ls_true-Ls_true_user).min())[0]
print("La valeur la plus proche trouvée est Ls = "
      + str(Ls_true[Ls_ind]-year_user*360.)
      + " pour l'année " + str(year_user))
# Code ------------------------------------------------------------
var = data.variables[varname][:,:,:]
dataplt = var
fig = plt.figure(figsize=(12,8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
if (earthtopo): ax.coastlines(resolution="110m",linewidth=1)
gl = ax.gridlines(linestyle='--',color='black',
             draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False
clevs = np.linspace(vmin,vmax,29)
plt.contourf(longitude, latitude, np.squeeze(dataplt[Ls_ind,:,:]),
             clevs, transform=ccrs.PlateCarree(),cmap="jet")
plt.title(r"Température de surface", size=14)
cb = plt.colorbar(ax=ax, orientation="vertical", pad=0.02, aspect=16, shrink=0.8)
cb.set_label(r'K',size=12,rotation=0,labelpad=15)
cb.ax.tick_params(labelsize=10)
plt.show()
# Paramètres utilisateurs -----------------------------------------
earthtopo = False # ajouter les traits de côte actuels
year_user = 1 # année de simulation à regarder
Ls_user = 270. # longitude solaire choisie
ph2oatmo = 0.05e-2*610. # assumed mean water vapor partial pressure
vmin = 0.
vmax = 1.
# Code ------------------------------------------------------------
Ls_true_user = year_user*360. + Ls_user
Ls_ind = np.where(abs(Ls_true-Ls_true_user)==
                  abs(Ls_true-Ls_true_user).min())[0]
print("La valeur la plus proche trouvée est Ls = "
      + str(Ls_true[Ls_ind]-year_user*360.)
      + " pour l'année " + str(year_user))
# Code ------------------------------------------------------------
tsurfnc = data.variables['tsurf'][:,:,:]
dataplt = ph2oatmo/psati(tsurfnc)
fig = plt.figure(figsize=(12,8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
if (earthtopo): ax.coastlines(resolution="110m",linewidth=1)
gl = ax.gridlines(linestyle='--',color='black',
             draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False
clevs = np.linspace(vmin,vmax,29)
plt.contourf(longitude, latitude, np.squeeze(dataplt[Ls_ind,:,:]),
             clevs, transform=ccrs.PlateCarree(),cmap="jet")
plt.title(r"Saturation ratio", size=14)
cb = plt.colorbar(ax=ax, orientation="vertical", pad=0.02, aspect=16, shrink=0.8)
cb.set_label(r'NU',size=12,rotation=0,labelpad=15)
cb.ax.tick_params(labelsize=10)
plt.show()
Explanation: Carte à $L_s$ donné de $p_{sat}$ pour $H_2O$
End of explanation
def psatco2(temp):
  return 1.382 * 1e12 * np.exp(-3182.48/temp)
temp = np.linspace(100,200,81)
plt.yscale('log')
plt.plot(temp,psatco2(temp))
plt.show()
# Paramètres utilisateurs -----------------------------------------
earthtopo = False # ajouter les traits de côte actuels
pco2atmo = 610. # CO2 pressure
vmin = 0.
vmax = 1.
# Code ------------------------------------------------------------
tsurfnc = data.variables['tsurf'][:,:,:]
dataplt = pco2atmo/psatco2(tsurfnc)
fig = plt.figure(figsize=(12,8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_global()
if (earthtopo): ax.coastlines(resolution="110m",linewidth=1)
gl = ax.gridlines(linestyle='--',color='black',
             draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False
clevs = np.linspace(vmin,vmax,21)
plt.contourf(longitude, latitude, np.mean(dataplt[:,:,:],axis=0),
             clevs, transform=ccrs.PlateCarree(),cmap="jet")
plt.title(r"Saturation ratio", size=14)
cb = plt.colorbar(ax=ax, orientation="vertical", pad=0.02, aspect=16, shrink=0.8)
cb.set_label(r'Pa',size=12,rotation=0,labelpad=15)
cb.ax.tick_params(labelsize=10)
plt.show()
Explanation: Carte en moyenne temporelle de $p_{sat}$ pour $CO_2$
End of explanation
<END_TASK>