From 59ca53aaf151d6192b756b8f1905a9c06d4f2660 Mon Sep 17 00:00:00 2001
From: cfontana00 <cfontana@ogs.it>
Date: Thu, 16 May 2024 13:57:33 +0200
Subject: [PATCH] argo OK

---
 config/config_MARIA.json                      |  34 ++-
 config/config_TEMPLATE.json                   |   4 +
 datasets/ARGO/.compare_argo.py.swp            | Bin 12288 -> 0 bytes
 datasets/ARGO/compare_argo.py                 | 189 +++++++++++++-
 datasets/ARGO/get_argo.py                     |  16 +-
 datasets/ARGO/old/compare_argo.py             | 244 ++++++++++++++++++
 datasets/ARGO/old/compare_bioargo.py          | 189 ++++++++++++++
 datasets/BOUSSOLE/read_boussole.py            |  18 +-
 meteo/ECMWF/process_data.py                   |   1 +
 modules/__pycache__/fun_gen.cpython-310.pyc   | Bin 2009 -> 1796 bytes
 modules/__pycache__/fun_io.cpython-310.pyc    | Bin 4360 -> 4420 bytes
 .../fun_plot_profile.cpython-310.pyc          | Bin 0 -> 1503 bytes
 modules/fun_gen.py                            |  17 +-
 modules/fun_io.py                             |  26 +-
 modules/fun_plot_profile.py                   |  81 ++++++
 15 files changed, 758 insertions(+), 61 deletions(-)
 delete mode 100644 datasets/ARGO/.compare_argo.py.swp
 create mode 100644 datasets/ARGO/old/compare_argo.py
 create mode 100644 datasets/ARGO/old/compare_bioargo.py
 create mode 100644 modules/__pycache__/fun_plot_profile.cpython-310.pyc
 create mode 100644 modules/fun_plot_profile.py

diff --git a/config/config_MARIA.json b/config/config_MARIA.json
index 46544e6..96db83c 100644
--- a/config/config_MARIA.json
+++ b/config/config_MARIA.json
@@ -1,32 +1,32 @@
 {
 "_____":"GENERAL PARAMETERS",
 
-"date_ini":"2015-03-20",
-"date_end":"2015-03-29",
-"outdir":"/home/fontana/MODEL/OUTPUTS",
-"diagdir":"/home/fontana/MODEL/DIAGS",
+"date_ini":"2021-12-11",
+"date_end":"2021-12-12",
+"outdir":"/path/to/OUTPUT_SMP",
+"diagdir":"/path/to/diags/dir",
+"river_ini":"2015-01-03",
+"river_end":"2015-05-01",
 
 
 "_____":"RIVER PARAMETERS",
 
 "french_list":["Argens","Var"],
-"italian_list":["Magra","Arno","Ombrone","Serchio"],
-"river_ini":"2019-05-01",
-"river_end":"2019-05-31",
+"italian_list":["Arno","Magra","Ombrone","Serchio"],
 
 
 "_____":"METEO PARAMETERS",
-
-"meteo_ini":"2019-04-01",
-"meteo_end":"2019-04-30",
+"meteo_ini":"2015-03-10",
+"meteo_end":"2015-03-11",
 "day_cycle":2,
-"n_cycle":30,
+"n_cycle":6,
 "imin":"500",
 "imax":"600",
 "jmin":"350",
 "jmax":"450",
 
 
+
 "_____":"DYFAMED PARAMETERS",
 
 "some_par":"some_par",
@@ -36,7 +36,7 @@
 "some_par":"some_par",
 
 
-"_____":"PLOT PARAMETERS",
+"_____":"2D PLOT PARAMETERS",
 
 "resol":"10m",                       
 "fig_proj":"ccrs.PlateCarree()",      
@@ -53,14 +53,13 @@
 
 "_____":"SECTION PLOT PARAMETERS",
 
-
 "fig_secx":"8",                         
 "fig_secy":"4",                        
 "sec_lon1":"8.5",
 "sec_lat1":"43.",
 "sec_lon2":"10.",
 "sec_lat2":"43.5",
-"sec_dep":"200",
+"sec_dep":"600",
 "sec_resH":"0.05",  
 "sec_resV":"1",
 "cb_fraction_sec":"0.025",
@@ -68,8 +67,8 @@
 
 "_____":"PROFILE PLOT PARAMETERS",
 
-"fig_prox":"6",
-"fig_proy":"8",
+"fig_prox":"8",                         
+"fig_proy":"4",                        
 "col1":"r",
 "col2":"b",
 
@@ -77,7 +76,6 @@
 "_____":"INTERPOLATION PARAMETERS",
 
 "dump":"0.01",
-"itp_meth":"linear",
-"cmap":"viridis"
+"itp_meth":"nearest"
 
 }
diff --git a/config/config_TEMPLATE.json b/config/config_TEMPLATE.json
index 022ebe4..f1a850f 100644
--- a/config/config_TEMPLATE.json
+++ b/config/config_TEMPLATE.json
@@ -25,6 +25,10 @@
 "jmin":"350",
 "jmax":"450",
 
+"_____":"ARGO",
+
+"argo_ds":"phy",
+
 
 
 "_____":"DYFAMED PARAMETERS",
diff --git a/datasets/ARGO/.compare_argo.py.swp b/datasets/ARGO/.compare_argo.py.swp
deleted file mode 100644
index a71541abadb86555f014362886b7bb6cfb30fb0e..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001

literal 12288
zcmeI2O>7%Q6vrp6fbvmj(+gbC5Gz`nZtS!@G>T*iW1A3UH<9BYRjDx(@2s=YdUn|z
zyY{8rK&ZWNDc=p0(kmcwfFlTrUXkF&i31-81mea8_|JYMPFf{gfm!M2wP)U&_uhZs
zj1*<pXP;fH)4BX<g6kL|lRp&gZ;#j$R~{wA_Jm}DshNGRe`<FrF)S=Mmdfq=Ql(le
zR7%ZKyI!pnEGAjVWmqUR7M2PfujjLXw^`8j^8RSw@On}!pcObsfxfh7=8lpxr%%l)
z)6CRk^u*-qL0)=ZE1(t73TOqi0$KsBfL1^&pcVKJDIom`@*?Uvl+^TGa-Z09Pk!o~
zRzNGD70?Q31+)TM0j+>mKr5gX&<bb;v;zM^1(-$1;Uk2+b}y30|Nqm!|G$5TkT1ba
za1{(d2dsk{I16&%1o-PgLVgB6fgizl;A?Obd;o5MSHUH)37TL9JOO6FJ>YI|7x?2S
zA-{lIU<+IWuYeFRPz5G94vvA}9w6jb@D2D3d<wRJ1R;1H)WIA$3o_s^I0U}9pODYN
zb#NJMf+jcv9sv)7-|r*j6L1Aw218H>li)b`6T003UxAOoN8m$n1H2EegZIF@U<+IW
z?|`?#Ti{Lb26!E)XK)p~3|<2I(h6t=v;tZIt$<eG|Dyn<yKAMm_|TEr9LCiTKX8Q1
z+Jg1CX-U?_xEkIWpheo{vTZvqS4LT7F8sXhdR~y-v~oGqj@*=2hAwbPQ`2~*O$#en
zGLM%sxpaCwAIiX2U$Pdyna1?_nd#olv}NS-5=&Z+U@o##%P9jm<1?!q)*NftTAOX<
z@}cjd1fyx>@G~;a(6T26Jdjkyi6g5eiq!R()gG|GVe2l>c06G_U9;y!V{w>4R6JE>
zn6%SlzDb?X_3&c=JN4-Z^~r`Nb=f+1P1+ZZ3@;>+Rn47#5O5)B8i?YRQ>Z0*5N;=v
z1Vv7U+~AUVMzjET9b5U)A|{&E`f|Hed1iH`S*;lJkpapWa}Dq0w!BiR?eEsvaN}@e
z!OBn9O5<LLW_v-DizDcov3uS$c1JPyB-@%-mEswpKs3n34Ic0;<QJ+iU)`15af>k5
zrH-)pkV767Gnz+Gn6)A9Jju|RJr{;*p6636<OiN=GWf>rn|S5J4d!#SglVJ}D5x*2
z3&O~e#3;@VF<$8NEXL6Ax5r9DW2(Y!NAT^C)b~P1I-bDa1?b6mH5LSH6!j?f(8%XS
zv@F*310KS_3fGL~Ms;OJ>K5sYI>}7=e7RO_pIfazS#20CBKp15rB-S_wbZEZ%q0QU
zdd<ap71qc<v92whTUczaR-#cR)6r5h9VC{@Wf<h64UxPeoekN5TTXD&C`7MEiw=ZJ
z#?%rH>&COj$$00$P`y-bjB6G>7-p%Rr(W<@I;(x-G8y0bK;4T`+>YAfA@=v}ny1S|
zCO8FT)Ej}<qjq1kyIfGG=X==Uk0S*3bevGKqxWA>?fK|fLVB?qiVH_!9z_q<PPJ59
zmZOuQA29WBsX7wIa>Q4)%VCb5BS{|ViAGt711=pMM$MqchCw=ZXG<htdLgmX5*=@+
z!tKFbb*Jm0N&JxtA$)R?1;P>CFdrf<7tQM8(8*3YUEu|MDwn#XW@AOAO}>4n8O%1k
WzElk_J$3HW65g7Ro-8yv!1r%?f<f&7

diff --git a/datasets/ARGO/compare_argo.py b/datasets/ARGO/compare_argo.py
index 613d104..d901523 100644
--- a/datasets/ARGO/compare_argo.py
+++ b/datasets/ARGO/compare_argo.py
@@ -7,10 +7,18 @@
 
 from fun_gen import *
 from fun_io import *
+from fun_plot_2D import *
+from fun_plot_profile import *
 import sys,os
-import datetime as dt
 import xarray as xr
 import datetime as dt
+from datetime import datetime
+
+from mpl_toolkits.axes_grid1 import host_subplot
+import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.use("Agg")
+
 
 
 # Get args
@@ -23,6 +31,17 @@
 load_config(config)
 from fun_gen import *
 
+tag = ''
+if argo_ds == 'bgc':
+  tag = 'bio'
+
+# Create arborescence
+# -------------------
+os.system('mkdir -p '+diagdir)
+os.system('mkdir -p '+diagdir+'/'+config)
+savedir = diagdir+'/'+config+'/ARGO'
+os.system('mkdir -p '+savedir)
+
 # Load coordinates
 # ----------------
 lon_mod,lat_mod,lev_mod=load_coords()
@@ -31,54 +50,196 @@
 # Load data
 # ---------
 savedir = diagdir+'/'+config+'/ARGO'
-ds = xr.load_dataset(savedir+'/argo.nc')
+ds = xr.open_dataset(savedir+'/'+tag+'argo.nc')
+
+lat = np.array( ds['LATITUDE'])
+lon = np.array(ds['LONGITUDE'])
+time = np.array(ds['TIME'])
+num = np.array(ds['PLATFORM_NUMBER'])
+cyc = np.array(ds['CYCLE_NUMBER'])
+pres = np.array(ds['PRES_ADJUSTED'])
 
-lat = ds['LATITUDE']
-lon = ds['LONGITUDE']
-time = ds['TIME']
-num = ds['PLATFORM_NUMBER']
-cyc = ds['CYCLE_NUMBER']
-pres = ds['PRES_ADJUSTED']
 
 time = np.array(time,dtype=str)
 
+
 # Define profiles position
 # ------------------------
 
 lon_uni = np.unique(lon)
 print(lon_uni.shape[0],'profiles found\n')
 
+
 # Loop on profiles
 # ----------------
+print('Processing\n')
 for plon in lon_uni:
 
   # Get all indexes
   idx = np.where( lon == plon )
 
-  for var in ['CHLA_ADJUSTED','PSAL_ADJUSTED','TEMP_ADJUSTED']:
+  # Get data lon/lat
+  data_lon = lon[idx]
+  data_lat = lat[idx]
+  data_pres = pres[idx]
+
+  pnum = str(num[idx][0])
+  pcyc = str(cyc[idx][0])
+  ptime = str(time[idx][0])
+
+
+  print('Float',pnum,'cycle',pcyc)
+  print('-----')
+
+
+  full_data = []
+  full_val = []
 
-    if var == 'CHLA_ADJUSTED':
+  varlist = ['PSAL_ADJUSTED','TEMP_ADJUSTED']
+
+  if argo_ds == 'bgc':
+     varlist = np.concatenate((['CHLA'],varlist))
+
+  for var in varlist:
+
+    if var == 'CHLA':
        mod_var = 'chl'
     elif var == 'PSAL_ADJUSTED':
        mod_var = 'so'
     elif var == 'TEMP_ADJUSTED':
        mod_var = 'thetao'
 
+    # Get data
+    data = ds[var][idx]
+
     # Get current variable parameters
     vname, ftag, cmap, islog, vmod, vmin, vmax, label, units\
         = load_variable(config,mod_var)
 
     # Convert date
-    dstr = str(time[idx][0]).split('T')[0]
-    jd = dt.datetime.strptime(dstr,'%Y-%m-%d').toordinal()
+    dstr = str(time[idx][0]).replace('.000000000','')
+    date = dt.datetime.strptime(dstr,'%Y-%m-%dT%H:%M:%S')
+    jd = date.toordinal()
+ 
+    hours,mnt = date.strftime('%H'),date.strftime('%M')
+    hours = int(np.round(float(hours)+float(mnt)/60.))
+    if hours == 24:
+      hours = 0
+
 
     # Get filename
     fname,dtag = get_filename(int(np.floor(jd)),ftag)
 
-    print(fname,dtag)
+
+    # Interpolate model on data
+    try:
+      val  = get_model_val_3d(fname,hours,vname,\
+               lon_mod,lat_mod,lev_mod,\
+                 data_lon,data_lat,data_pres)
+    except Exception as e:
+      print('Interpolation failed')
+      val = data_pres
+      val[:] = np.nan
+
+    full_val.append(val)
+    full_data.append(data)
+
+
+  full_val = np.array(full_val)
+  full_data = np.array(full_data)
+
+  # Plot profile
+  # ------------
+
+  os.system('mkdir -p '+savedir+'/'+tag+str(pnum))
+
+  fout = savedir+'/'+tag+str(pnum)+'/profile_'+pcyc+'.'+fig_fmt
+  ptime = ptime.replace('.000000000','')
+
+  if argo_ds == 'bgc':
+    mchl,mpsal,mtemp = full_val[0],full_val[1],full_val[2]
+    dchl,dpsal,dtemp = full_data[0],full_data[1],full_data[2]
+
+    plot_profiles(fout,ptime,-data_pres,\
+               mchl,dchl,'Chlorophyll (mg.m$^{-3}$)',\
+               mpsal,dpsal,'Salinity',\
+               mtemp,dtemp,'Temperature (°C)')
+  else:
+    mpsal,mtemp = full_val[0],full_val[1]
+    dpsal,dtemp = full_data[0],full_data[1]
+
+    plot_profiles(fout,ptime,-data_pres,\
+               mpsal,dpsal,'Salinity',\
+               mtemp,dtemp,'Temperature (°C)')
+
+
+
+
+# Plot trajectories
+# -----------------
+
+print('Plot trajectories ...')
+
+# Load proj
+exec('proj = ' + fig_proj)
+
+pos = []
+prof = []
+
+# Loop on Argo float
+for n in np.unique(num):
+
+   print(n)
+
+
+   # Initialize map
+   extent = [lon_mod.min(),lon_mod.max(),lat_mod.min(),lat_mod.max()]
+   fig, ax = plt.subplots(1,1,figsize=(float(fig_sx), float(fig_sy)), subplot_kw={'projection': proj})
+   init_fig(ax,extent,proj)
+
+   # Get all indexes
+   idx = np.where( num == n )
+
+   plon = np.unique( lon[idx])
+   plat = np.unique( lat[idx])
+   txt = np.unique(cyc[idx])
+
+   pos.append([plon,plat,n])
+
+
+   plt.plot(plon,plat,marker='+',linestyle='-')
+
+   for i in [0,plon.shape[0]-1]:
+     plt.text(plon[i],plat[i],str(txt[i]),fontsize=8)
+
+   # Save figure
+   fout = savedir+'/'+tag+str(n)+'/map'+'.'+fig_fmt
+   savefig(fout)
+
+   plt.close()
+
+
+# Figure for all trajectories
+fig, ax = plt.subplots(1,1,figsize=(float(fig_sx), float(fig_sy)), subplot_kw={'projection': proj})
+init_fig(ax,extent,proj)
+
+
+for p in pos:
+   
+   lon = np.array(p[0])
+   lat = np.array(p[1])
+   
+   plt.plot(lon,lat,marker='+',linestyle='none')
+   plt.text(lon[0],lat[0],str(p[2]))
+
+# Save figure
+fout = savedir+'/map'+tag+'.'+fig_fmt
+savefig(fout)
+
+plt.close()
+
 
 
-    exit()
 
 
 
diff --git a/datasets/ARGO/get_argo.py b/datasets/ARGO/get_argo.py
index 4a8012e..82663dd 100644
--- a/datasets/ARGO/get_argo.py
+++ b/datasets/ARGO/get_argo.py
@@ -1,6 +1,8 @@
 # ------------------------- #
 # Routine to get Argo data  #
 # ------------------------- #
+import warnings
+warnings.simplefilter("ignore")
 
 import numpy as np
 
@@ -52,11 +54,14 @@
 # -------------
 box = [lon_min, lon_max, lat_min, lat_max, 0, 300, date_ini, date_end]
 
-f = DataFetcher(ds='bgc', mode='expert', params='all',
+f = DataFetcher(ds=argo_ds, mode='expert', params='all',
                 parallel=True, progress=True, cache=False,
                 chunks_maxsize={'time': 30},
                )
 
+
+
+
 f = f.region(box).load()
 
 ds = f.data
@@ -67,6 +72,13 @@
 
 # Save data
 # ---------
-ds.to_netcdf(path=savedir+'/argo.nc', mode='w', format=None, group=None, engine=None, encoding=None, unlimited_dims=None, compute=True, invalid_netcdf=False)
+tag=''
+if argo_ds == 'bgc':
+  tag = 'bio'
+
+fname = savedir+'/'+tag+'argo.nc'
+ds.to_netcdf(path=fname, mode='w', format=None, group=None, engine=None, encoding=None, unlimited_dims=None, compute=True, invalid_netcdf=False)
+
+print('[SAVED]',fname)
 
 
diff --git a/datasets/ARGO/old/compare_argo.py b/datasets/ARGO/old/compare_argo.py
new file mode 100644
index 0000000..35dc884
--- /dev/null
+++ b/datasets/ARGO/old/compare_argo.py
@@ -0,0 +1,244 @@
+# ------------------------- #
+# Routine to get Argo data  #
+# ------------------------- #
+import warnings
+warnings.simplefilter("ignore")
+import numpy as np
+
+from fun_gen import *
+from fun_io import *
+from fun_plot_2D import *
+import sys,os
+import datetime as dt
+import xarray as xr
+import datetime as dt
+
+from mpl_toolkits.axes_grid1 import host_subplot
+import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.use("Agg")
+
+
+
+# Get args
+# --------
+config = sys.argv[1]  # Configuration name
+
+
+# Load config file
+# ----------------
+load_config(config)
+from fun_gen import *
+
+# Create arborescence
+# -------------------
+os.system('mkdir -p '+diagdir)
+os.system('mkdir -p '+diagdir+'/'+config)
+savedir = diagdir+'/'+config+'/ARGO'
+os.system('mkdir -p '+savedir)
+
+# Load coordinates
+# ----------------
+lon_mod,lat_mod,lev_mod=load_coords()
+
+
+# Load data
+# ---------
+savedir = diagdir+'/'+config+'/ARGO'
+ds = xr.open_dataset(savedir+'/argo.nc')
+
+lat = np.array( ds['LATITUDE'])
+lon = np.array(ds['LONGITUDE'])
+time = np.array(ds['TIME'])
+num = np.array(ds['PLATFORM_NUMBER'])
+cyc = np.array(ds['CYCLE_NUMBER'])
+pres = np.array(ds['PRES_ADJUSTED'])
+
+
+
+time = np.array(time,dtype=str)
+
+# Define profiles position
+# ------------------------
+
+lon_uni = np.unique(lon)
+print(lon_uni.shape[0],'profiles found\n')
+
+# Loop on profiles
+# ----------------
+print('Processing\n')
+for plon in lon_uni:
+
+
+  # Get all indexes
+  idx = np.where( lon == plon )
+
+  # Get data lon/lat
+  data_lon = lon[idx]
+  data_lat = lat[idx]
+  data_pres = pres[idx]
+
+  pnum = str(num[idx][0])
+  pcyc = str(cyc[idx][0])
+  ptime = str(time[idx][0])
+
+
+
+  full_data = []
+  full_val = []
+
+  for var in ['PSAL_ADJUSTED','TEMP_ADJUSTED']:
+
+    if var == 'PSAL_ADJUSTED':
+       mod_var = 'so'
+    elif var == 'TEMP_ADJUSTED':
+       mod_var = 'thetao'
+
+    # Get data
+    data = ds[var][idx]
+
+    # Get current variable parameters
+    vname, ftag, cmap, islog, vmod, vmin, vmax, label, units\
+        = load_variable(config,mod_var)
+
+    # Convert date
+    dstr = str(time[idx][0])
+
+    print('Float',pnum,'cycle',pcyc)
+
+    date = dt.datetime.strptime(dstr,'%Y-%m-%dT%H:%M:%S.000000000')
+    jd = date.toordinal()
+
+    hours,mnt = date.strftime('%H'),date.strftime('%M')
+    hours = int(np.round(float(hours)+float(mnt)/60.))
+
+
+    # Get filename
+    fname,dtag = get_filename(int(np.floor(jd)),ftag)
+
+
+    # Interpolate model on data
+    try :
+      val  = get_model_val_3d(fname,hours,vname,\
+               lon_mod,lat_mod,lev_mod,\
+                 data_lon,data_lat,data_pres)
+    except Exception as e :
+      print('Interp failed')
+      print(e)
+      val = data_lon
+      val[:] = np.nan
+
+   
+
+    full_val.append(val)
+    full_data.append(data)
+
+
+  full_val = np.array(full_val)
+  full_data = np.array(full_data)
+
+  # Plot profile
+  # ------------
+
+  fig = plt.figure(figsize=(float(fig_prox), float(fig_proy)))
+  ax1 = fig.add_subplot(111)
+  ax2 = ax1.twiny()
+
+  psal,temp = full_val[0],full_val[1]
+
+  # Plot model
+  ax1.plot(psal,-data_pres,c=col2,label='Salinity')
+  ax2.plot(temp,-data_pres,c=col3,label='Temperature')
+
+  psal,temp = full_data[0],full_data[1]
+
+  # Plot data
+  ax1.plot(psal,-data_pres,c=col2,marker='.',linestyle='None',alpha=0.5)
+  ax2.plot(temp,-data_pres,c=col3,marker='.',linestyle='None',alpha=0.5)
+
+
+  ax1.set_xlabel('Salinity')
+  ax1.set_ylabel('Depth (m)')
+  ax1.xaxis.label.set_color(col2)
+  ax1.spines["bottom"].set_color(col2)
+  ax1.tick_params(axis='x',colors=col2, which='both')
+  ax1.patch.set_alpha(0.0)
+
+  ax2.set_xlabel('Temperature (°C)')
+  ax2.xaxis.label.set_color(col2)
+  ax2.spines["top"].set_color(col3)
+  ax2.xaxis.label.set_color(col3)
+  ax2.tick_params(colors=col3, which='both')
+
+
+  plt.title(ptime)
+
+
+
+  # Save figure
+  # -----------
+  os.system('mkdir -p '+savedir+'/'+str(pnum))
+
+  fout = savedir+'/'+pnum+'/profile_'+pcyc+'.'+fig_fmt
+  savefig(fout)
+  print('')
+
+  plt.close()
+
+  exit()
+
+
+# Plot trajectories
+# -----------------
+
+print('Plot trajectories ...')
+
+
+# Loop on Argo float
+for n in np.unique(num):
+
+   print(n)
+
+   # Load proj
+   exec('proj = ' + fig_proj)
+
+   # Initialize map
+   extent = [lon_mod.min(),lon_mod.max(),lat_mod.min(),lat_mod.max()]
+   fig, ax = plt.subplots(1,1,figsize=(float(fig_sx), float(fig_sy)), subplot_kw={'projection': proj})
+   init_fig(ax,extent,proj)
+
+   # Get all indexes
+   idx = np.where( num == n )
+
+   plon = np.unique( lon[idx])
+   plat = np.unique( lat[idx])
+   txt = np.unique(cyc[idx])
+
+   plt.plot(plon,plat,marker='+',linestyle='-')
+
+   for i in [0,plon.shape[0]-1]:
+     plt.text(plon[i],plat[i],str(txt[i]),fontsize=8)
+
+   # Save figure
+   fout = savedir+'/'+str(n)+'/map'+'.'+fig_fmt
+   savefig(fout)
+
+   plt.close()
+
+
+  
+  
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/datasets/ARGO/old/compare_bioargo.py b/datasets/ARGO/old/compare_bioargo.py
new file mode 100644
index 0000000..6c6250e
--- /dev/null
+++ b/datasets/ARGO/old/compare_bioargo.py
@@ -0,0 +1,189 @@
+# ------------------------- #
+# Routine to get Argo data  #
+# ------------------------- #
+import warnings
+warnings.simplefilter("ignore")
+import numpy as np
+
+from fun_gen import *
+from fun_io import *
+from fun_plot_2D import *
+from fun_plot_profile import *
+import sys,os
+import xarray as xr
+import datetime as dt
+from datetime import datetime
+
+from mpl_toolkits.axes_grid1 import host_subplot
+import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.use("Agg")
+
+
+
+# Get args
+# --------
+config = sys.argv[1]  # Configuration name
+
+
+# Load config file
+# ----------------
+load_config(config)
+from fun_gen import *
+
+# Create arborescence
+# -------------------
+os.system('mkdir -p '+diagdir)
+os.system('mkdir -p '+diagdir+'/'+config)
+savedir = diagdir+'/'+config+'/ARGO'
+os.system('mkdir -p '+savedir)
+
+# Load coordinates
+# ----------------
+lon_mod,lat_mod,lev_mod=load_coords()
+
+
+# Load data
+# ---------
+savedir = diagdir+'/'+config+'/ARGO'
+ds = xr.open_dataset(savedir+'/bioargo.nc')
+
+lat = np.array( ds['LATITUDE'])
+lon = np.array(ds['LONGITUDE'])
+time = np.array(ds['TIME'])
+num = np.array(ds['PLATFORM_NUMBER'])
+cyc = np.array(ds['CYCLE_NUMBER'])
+pres = np.array(ds['PRES_ADJUSTED'])
+
+
+time = np.array(time,dtype=str)
+
+
+# Define profiles position
+# ------------------------
+
+lon_uni = np.unique(lon)
+print(lon_uni.shape[0],'profiles found\n')
+
+
+# Loop on profiles
+# ----------------
+print('Processing\n')
+for plon in lon_uni:
+
+  # Get all indexes
+  idx = np.where( lon == plon )
+
+  # Get data lon/lat
+  data_lon = lon[idx]
+  data_lat = lat[idx]
+  data_pres = pres[idx]
+
+  pnum = str(num[idx][0])
+  pcyc = str(cyc[idx][0])
+  ptime = str(time[idx][0])
+
+
+  print('Float',pnum,'cycle',pcyc)
+
+
+  full_data = []
+  full_val = []
+
+  for var in ['CHLA','PSAL_ADJUSTED','TEMP_ADJUSTED']:
+
+    if var == 'CHLA':
+       mod_var = 'chl'
+    elif var == 'PSAL_ADJUSTED':
+       mod_var = 'so'
+    elif var == 'TEMP_ADJUSTED':
+       mod_var = 'thetao'
+
+    # Get data
+    data = ds[var][idx]
+
+    # Get current variable parameters
+    vname, ftag, cmap, islog, vmod, vmin, vmax, label, units\
+        = load_variable(config,mod_var)
+
+    # Convert date
+    dstr = str(time[idx][0]).replace(':00.000000000','')
+    date = dt.datetime.strptime(dstr,'%Y-%m-%dT%H:%M')
+    jd = date.toordinal()
+ 
+    hours,mnt = date.strftime('%H'),date.strftime('%M')
+    hours = int(np.round(float(hours)+float(mnt)/60.))
+
+
+    # Get filename
+    fname,dtag = get_filename(int(np.floor(jd)),ftag)
+
+
+    # Interpolate model on data
+    val  = get_model_val_3d(fname,hours,vname,\
+               lon_mod,lat_mod,lev_mod,\
+                 data_lon,data_lat,data_pres)
+
+    full_val.append(val)
+    full_data.append(data)
+
+
+  full_val = np.array(full_val)
+  full_data = np.array(full_data)
+
+  # Plot profile
+  # ------------
+  mchl,mpsal,mtemp = full_val[0],full_val[1],full_val[2]
+  dchl,dpsal,dtemp = full_data[0],full_data[1],full_data[2]
+
+  os.system('mkdir -p '+savedir+'/bio_'+str(pnum))
+  fout = savedir+'/bio_'+str(pnum)+'/profile_'+pcyc+'.'+fig_fmt
+
+  ptime = ptime.replace(':00.000000000','')
+  plot_profiles(fout,ptime,-data_pres,\
+               mchl,dchl,'Chlorophyll (mg.m$^{-3}$)',\
+               mpsal,dpsal,'Salinity',\
+               mtemp,dtemp,'Temperature (°C)')
+
+
+
+
+# Plot trajectories
+# -----------------
+
+print('Plot trajectories ...')
+
+# Loop on Argo float
+for n in np.unique(num):
+
+   print(n)
+
+   # Load proj
+   exec('proj = ' + fig_proj)
+
+   # Initialize map
+   extent = [lon_mod.min(),lon_mod.max(),lat_mod.min(),lat_mod.max()]
+   fig, ax = plt.subplots(1,1,figsize=(float(fig_sx), float(fig_sy)), subplot_kw={'projection': proj})
+   init_fig(ax,extent,proj)
+
+   # Get all indexes
+   idx = np.where( num == n )
+
+   plon = np.unique( lon[idx])
+   plat = np.unique( lat[idx])
+   txt = np.unique(cyc[idx])
+
+   plt.plot(plon,plat,marker='+',linestyle='-')
+
+   for i in [0,plon.shape[0]-1]:
+     plt.text(plon[i],plat[i],str(txt[i]),fontsize=8)
+
+   # Save figure
+   fout = savedir+'/bio_'+str(n)+'/map'+'.'+fig_fmt
+   savefig(fout)
+
+   plt.close()
+
+
+
+
diff --git a/datasets/BOUSSOLE/read_boussole.py b/datasets/BOUSSOLE/read_boussole.py
index 13f3262..898f006 100644
--- a/datasets/BOUSSOLE/read_boussole.py
+++ b/datasets/BOUSSOLE/read_boussole.py
@@ -2,7 +2,8 @@
 # Routine reading BOUSSOLE data     #
 # and interpolating model variables #
 #####################################
-
+import warnings
+warnings.simplefilter("ignore")
 import numpy as np
 
 from fun_gen import *
@@ -49,9 +50,6 @@
      year = line[4]
      hour = str(int(float(line[5])))
 
-     if year == '2022':
-       year = '2021' # !!!!!!!!!!!!
-
      try:
        jd = dt.datetime.strptime(year+mon+day,'%Y%m%d').toordinal()
        jd = jd + float(line[5])/24.
@@ -75,9 +73,6 @@
 print('Processing day:\n')
 for jd in range(jdini,jdend+1):
 
-   # Get filename
-   fname,dtag = get_filename(int(np.floor(jd)),ftag)
-   print(dtag)
     
    profile=[]
                                                                
@@ -90,6 +85,12 @@
    # ----------------
    if idx.any() : 
 
+     d = dt.datetime.fromordinal(int(np.floor(jd)))
+     print('Data found on',d.strftime('%Y-%m-%d'))
+
+     # Get filename
+     fname,dtag = get_filename(int(np.floor(jd)),ftag)
+
      # Select data for this layer
      slc = data[idx].squeeze()
      slc = np.flipud(slc)
@@ -105,7 +106,6 @@
 
 
      # Store profile
-     d = dt.datetime.fromordinal(int(np.floor(jd)))
 
      # Write data 
      np.savetxt(savedir+'/boussole_profile_'+d.strftime('%Y-%m-%d')+'.dat',profile)
@@ -122,7 +122,7 @@
      plt.xlabel(label+' ('+units+')')
      plt.ylabel("Depth (m)")
 
-     #savefig(savedir+'/boussole_profile_'+d.strftime('%Y-%m-%d')+'.'+fig_fmt)
+     savefig(savedir+'/boussole_profile_'+d.strftime('%Y-%m-%d')+'.'+fig_fmt)
      plt.show()
 
      print('')
diff --git a/meteo/ECMWF/process_data.py b/meteo/ECMWF/process_data.py
index 299147e..86efe8c 100644
--- a/meteo/ECMWF/process_data.py
+++ b/meteo/ECMWF/process_data.py
@@ -28,6 +28,7 @@
 load_config(config)
 from fun_gen import *
 
+
 # Try to read date as parameters
 # ------------------------------
 try : 
diff --git a/modules/__pycache__/fun_gen.cpython-310.pyc b/modules/__pycache__/fun_gen.cpython-310.pyc
index 3ff56f4eedee33c5ea64baa012fc7ce81f7175aa..fe59a5a502df50fa194555328e09210aad8fc13d 100644
GIT binary patch
delta 206
zcmcb~-@?b6&&$ij00fR-Y}14{^1fzfWSabwS(CAAvNDUKUkytNt0V)EWCN4zk_;uR
z3)p~s4oQXu><bwdG6LD0Kv51L%_Ygu%m`+2gIGWj9<VN6Nrp7eU<OUT$wI7#S}z%a
zw!LHm5zHWh1w^od2sRMG4k9>!geK?Yv#hd{-?Ijb+W@)6azKKMS&ET^nTL^uk?B7h
T@8nFjuZ+o)<JrrZ^?1|)B|a+T

delta 421
zcmZ9H%}T>S6ouzzGLz{vleYd6K?t}gS;XQtg5WbKuH2-?N}(~LK_wKrs4KVQ1GJBz
zZ{s6G*DhS?olvsq44iwuduPrtFWzTI<w{9{>*P|6<h%N&0Q#4P@TCvz(^<h%J_jIJ
zohi(%+KH*1wtTW8iPUb$MN$oUNS%hdNWF$OOfGhqZ*;R(OfC4rLB_;|SPsY%5HPT4
zVKq@uI}>9jTlKYVpNv{(B;3&-_R*~DMr*5-7s>rN8BVT)xa3z^au>vvg`Y0jB)>sV
zjrX$g?I25@gDg++U}c*of5a{VTLF&aiW}&yJrc<ZuwS(uqXH$ATCsh57p)H2dH)#A
h!y$^KdmyL_HdudxO7NO3fK&a(eqg9y`4f%Y{y&cXR(AjZ

diff --git a/modules/__pycache__/fun_io.cpython-310.pyc b/modules/__pycache__/fun_io.cpython-310.pyc
index 78b7051a0234edaabaaa582943e2f84a534da0b7..4d59a8ff597f3dd51e56a7d18879871e1642f338 100644
GIT binary patch
delta 2310
zcmaJ?OK%)S5T5Ry*Urw)KKxkQv0}$@OjeEwc_=^<K_DiEpoD|NBjMpJ+q3rC`!GG@
zgtdCIkqm-M5Mlnnni~fWNIAig0Dl4I!a0Hi5=c2kP}OTYmIcgezpCo0sjlgI^rsW|
zC%jfB<7n{Q{%*FK%ir}z$-TSN*E<^5>$yh$TApa!;O0l#wF0-e4YZG^xC6AGr@0Gs
zfM>V|bdYCx4rr0*c>(AU@8kVIhxq^>1UkaWyjChcz|Ag~NTM&+Tkiw47e$Tpa##l5
zngts8VP*4ALd^It(gtkONiq`0^a2@(6MCRC3JxmhW&5#!!EZPZ|DX%N38!_fK_ZQ4
z>xZSrSy78fNV_c5CD|aH3~QStqM^|>CFQ!3*p_4{N#Bx8B^i=-txfRurq;YPs%xR$
zzDTsl1oycxB~3{}r<>l?mfhvd21B!}8VCI?vuX6`O>$F{OzP68GB^h7vLFzsFyimo
z5!(|s3TFI}&5%rdR4;m>04pU+n5Zbm7xj^g!^n>8z&_xTbiGyP3zb$&@L(Di0@S;d
zQkxWrr+zN+=s2-&)4v@iN}#T=n6RPvzWzyP97J&dAYrv;m{`FFYrenki%DG1Ab1Ex
zge*c1A&;QqGsO__3A<4i39I`z#33b?g%JN1t+VKd3t+o&eD)j|2AzIDN*yvu9HJZE
zqy0pu5a}cF6eubTVEah}gE@W`{=v@xl3^AjJ=D5H>YNU1_=l_-TPtj5#oStfJa5@p
z>qoH~>JneDr7liMTf*y-NCUEIRlzh$$y8V01llZYYdK`~W4asBvS%gjuF)$InVQ$2
zM$GmK)QAH$k`FEE#NU|3OU^9^wNMklVNH~b%ZbTbjdHE2D&hf_#RCdZa>WD?aS&k~
z0TU&TAoTKsg)-1fM7oy<W2v<!0x=UmZ;f>@7vd<w|K{Sv>}y~X;9fQ@&;k?=q=f3M
zACfZ3fZLsz(_mIf0c<}5z@V&}?@5Rr4qL-h>!>ay4H`KO7Nv;>y_en~QCf8tlP2a*
zriZi@9U5axDrzg*IZTSHtaNGZvf7b$tTM{VkT%C`>$nc@da!DR9a~~Psq|z<x{Hh(
zcXUu@L79R+%)sV`^Z=Q3p>WtcEXqass33D)XH$~}g?U8ADC!DPUS?4SA8848GGApf
zi}fOH;2dt=(WAc5Sa$IYroc&RNW-alherKjhQpv(_V9bSBVl@O7(|nO@guvGI|c^`
z<IimP;Zlo>6Y;6ki-*-rI0G7S79cUI(CeY9iNv_M<O{!78mHrK>R4t#6@WO2<}>jx
zsj=xPWS>A7Q?CbMS%hN8pLiBk&mcS<A9N<?_Dq%>_G+WW{W^4Pec=Ref@y%W*QPmY
zvnet~2FWNLXQRZXqim24Qk{zV_%r7VawR^H9_hRV#_R6e^Dn&d((L$!SLWUwhk>TP
z>I!cAIzYl!>o|AHH~eaCQQ!y@a{wil=+)L*D0VFOOg?c9tUiZ-fKveiqJyNsa#XGN
zCP(~j`YZc5Iy1yvyyEVkRDGZ(mR?<JR`$BYS@&i9wOb@!{Db?`^>L7U)6s+@So4Ek
z->LYlL~r{r=Vs6bU8>P1UO-Nrw~klZo*fcbjceYnmK#357yq6aCM&V&{Xr`6L+|Y4
zaV(b54|aKCh;nn$7cb+&<9V>CRKaS!7A9tSb=7Zjfn6x@G!i=q#cEh<_(>|ntD)8`
z*I}I0(BQuU7E>BjM6Y0$(TTat>$M<ENG;JhINbX!j{25Hs}=Y^#=m6`O}F8#(!Rv}
z5RRqXRwFAhnrn^Kb}w)Ytrz5b@tNEtxf);39a>V?3LaZxHU04X?1hs;M?<VMbzP~u
zK;YaF^9Y#u-nF!X9PYX+G*!jwYP8Qn^{&<Z=U^sj!QuGVT(RFJwq+0619sLPvnT9=
IU2qHk0N@PXE&u=k

delta 2327
zcmZ`)OK%)S5T5RyN6+l+dc1bl4_1le*p3$`KoTGUg?Ao=M^IuZ;gQG;?%B;c`!GFY
zJ8Si3kxV3n$RW%D4mqqP4#0s6H%=TPg&)9a4jhpXr`(C4YS!LZR+!Oz-Bn%PUDNed
z^}We2CUebf*3sbk_P2xOyWX8#k$it=?&FFsjL)?XJz)w9+CjmE4Q*aHA_MJ^a77l{
zVUZI9(C!hQ7=*SU@?r?u5iu<GKszePyjCjwg&UkNkwmZ5nx8<+E3zttANfDq1Rf*K
z&=%k=Iz`6fAL$t~6hEN*J5Rv0MV!0+n9xCQ-+(W?3m}M~HLXtC8qwBgBHCq<t~9|$
zv<)K2n6^pUBr>|DB0^6kx+NK@#I__em2^eA)+VfUS!)z^Ejq_rZxXF-z|vl1D<g8E
zOxNAima<E^4TfePXadG>8V#e@u8}oOQN@&AX7J<Ic}c)YsmFh?LwrD5C>U{GKb&!~
zn?)$Zuj(Uy$&?022^C?&YRv|KA0=ijSPg2S97V&iK8&Mfrq=Yu!eX;2MOXqwh(<1@
zv_Nv{>yjLuApA4>ufybQd{6&wb{{S=0g$k2BTB6B)0H4t599&V?{cXFE{&w7tG)zx
zf>k+{ioTSOb;zj%(E`{m96j|C7>2_DlsaUJI7Byc5BCzCf;SK4(=f0|0k$6uDSmIm
z7yb^Qs34s&t<55>OB54yOv4YNsc&gmUt?@(Aii5R=K5h|f<k*sG5M^r6mUTpN>_$5
z%V62nH?=lz<hMrz1o9T$iRJmenX<Y@j~0s0Yd{RMRRvzyz>C4iQnu0}-nD@j&QfLv
zFYKE(+5<D7hs-FM@jvG1L5a~SxpF@=auT5*7CDDvFB<U?tJo<yX$<5c&?m-qUsiw~
zq7L+^igkn*!i_XC@+qABII8p`bL`YRU=vONaGIxik|z)$s*`at#h~38r4wM(S>%A*
z4+a$FH2V$#v*F}5JjXN(Nu9RsI%_*bgB)}=NZU!Xiz)*%Q8}a9bsdspOB&IxYp-I2
zGN~0)<_@cEspIs`+|*20mcY>KsE5f}&{`2!nj$KzT$QOXVccdgl7kT&GSP*tb(ICA
zDgzW^S6F+X?X?HhK-byS)L@DpB4ZTC2HT#>q1?{B6=f6{DzD5`4!NT&=-I-&t+(@$
zvE=qnK=DycSc2cC?V%_upj0g7@Na~phD2rq3W&<bQ+COF22KykpIHy0YfT}K##imv
zXVa?i0*uHP0TQDOi5;cvNQ|{>fed;iIT7EtkMQ9wKE#HzZ*CeV%^>VeXNHk4BiUa{
z9><~Q5uS_RbEeMjsvsUFY`rOh8l-J);h4AvreP6)Q;%{sLuN>UjME8LB%BsmfflGv
z<%Rec=WB8?em67Txdg`R?nm>lT{?Yg;>@|TS0_qVugqP<ZRY_Jwp_zf<F5wgYDGSY
z`qKa<mgwc?N+kQ1yDFS~AFRHFFT~ma4nbwhGmobFMFu)8{xkCfKLVYSAy3EOxO-2g
z^iRu{US4S|?p_irF5Vkq4#1%<lqMHS6lX4#+4zBbzthjf#Dqgx3Bq1FrXE{~-U^`n
zVK(;SfE7``f{;c)$2)D;o{5_lYSAkD^+5atwxI{$!M44Z6PGzeRzbN)<Kwx%i60*r
zIQhsES;h^pDkO&V8<jwwLB*r`5u`)ma;*wihQGWVG=#)Nl{OCWFp|qrwH_pPgttYt
z;n$%2RM6o6{zOh|K%*YFkQIr!Bx=<#N=P-)1uVGxC!O>s>&<2u{2GsV({pPut2CCF
zH{iVdR$Af`qp?z7ZuJ~T*LqG~i7Va|xg3AxO?NOXz4OGoU(O=nbs@1D$qNX`ir&qM
yCnj;99nPc_Pwz*55~#gW3tq-6IK(XBj1zp_^zHSz$qRg#5AY(N<awRv-TZ&*XW3%_

diff --git a/modules/__pycache__/fun_plot_profile.cpython-310.pyc b/modules/__pycache__/fun_plot_profile.cpython-310.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..0de2f9a58d8e47fe8c8e6db66daba617b386e3a1
GIT binary patch
literal 1503
zcmZ`(&2QsG6rUN79p}4E)34pK0*M3RuvzXvSvJ*jU|R%xsszy_V<&EHk9Ry_(=wJ2
z<=X!tmAG)>Z{RQB#A&!J;zDkm5eMEl*`->EvHaefH}Adqy$@S7n>7UM*WZr9e_VwA
zu$$Y9z~&=(^9MjU;y6Z6R~(m|#N^ZgJ&9@Jo_ZK@htu!TX@!$7P|sb07wTby-#s2%
z@#X4*P=)s~y!oFPLSqPl5?tVf6eMv9C!qy|EL@DV%bjx~w>2%ilxU~$l&2~hU13i5
z5T|!RDqzp?J}P|RF7PUF4|q-aXH`|x_zD9hS2!fRa*d1n1QiWkU(~sOgNvqaE*iXg
zgW<Nu0(MqXE$-^d6&au}t3_+2U$k$KJ$rTnr+9t(ojuhpxUJ2pZgJX0yQtXENN?;T
zOS!v?){<u>$01I!S#*jm-9FpUo4Tew>#w(X?V17xtL?nh@;YzuCU0Fsd3AfSskcg;
z?X%9=)}p)eTx+&@d$s2qeDm4?dq4@b`xzqc@XkJZ&3<>W$5+v;9hJY7cHP(Ppbl{M
z;?>zc-{w0IdDrF{lCmCODSvTCmlafZwSQ&ncfj{b?l0cBK`<H!Mq>Y1Lc3u!j2}jm
zVJ5}=m2_b$V-XLt=~@_#C#J&WL=1UuoY7>wbij3K{g-wyVX!I?B#$nGrGWiqbr57~
z`sO=HuSX0X8#iX7AU62Ep~oT(CSU{W#bJ-i76`^mm^=@psm4(n<mw_0jLYI|$_$<x
zIudHya|~tkC^sG;6EZjM^Jz4m8V?SYNY)V{l%qsuO@n-``|EJ>P4}nAzhC}vXi)DS
zAf2VdNsvmLgQS+g2Zy9*Nd<^;qEs0t4pQS}u`(W{Hj{uTj0IDs&Nv_Dvr!fcWomh#
zh8L^sJeFla%8e^oIstR_JW4MN{m*W`vUBG)5!tp^Rjq%Sk*TX_d_K&WWC>)MXSNc?
z&6t8LOS4r$8rCBdd88td!X=S{VpVeXbpSBIURv<(I2L)(>lz2T_P>eBJqy82q4%Mo
zp_nN{Rg?g9Zrf)Zh$s)rim>^621kalx$#680tVwIf<sR!XNr}{{D|boVEp$7{i#TT
zekfAKQr7?M^OHxzCntlWPy2(1Pah7S91Z&5I*Wt6Z(CvOW5+Ly;^2OEVVbuUxgC%E
z6Cek~cRfr=m-^T~dxSbQY@Zg{z%9qe4e~bLA$M?%bZ{5@q>CFR^~nJy-`0p<+Q<O}
z>-q8?Fla3Rhy>JNmcZ2{&CF)PltmCnBS2F&&s%Z%7Pyq|WuCXlBOA+qf0fmT_Or<C
M-`p>uHgLoH7g|@2!~g&Q

literal 0
HcmV?d00001

diff --git a/modules/fun_gen.py b/modules/fun_gen.py
index 723b5a3..7c62ed6 100644
--- a/modules/fun_gen.py
+++ b/modules/fun_gen.py
@@ -39,6 +39,7 @@ def load_config(config):
    print('----------------------------------')
    print('=> Config file '+config+' loaded')
    print('----------------------------------\n')
+
  except Exception as e :
   print('Error in argument or config')
   print(e)
@@ -67,14 +68,14 @@ def load_variable(config,var):
      if par[0] == var:
        break
    
-   print('----------------')
-   print('Variable :',var)
-   print('cmap :',par[2])
-   print('log :',par[3])
-   print('min/max mode :',par[4])
-   print('vmin :',par[5])
-   print('vmax :',par[6])
-   print('----------------\n')
+   #print('----------------')
+   #print('Variable :',var)
+   #print('cmap :',par[2])
+   #print('log :',par[3])
+   #print('min/max mode :',par[4])
+   #print('vmin :',par[5])
+   #print('vmax :',par[6])
+   #print('----------------\n')
 
 
    return par[0],\
diff --git a/modules/fun_io.py b/modules/fun_io.py
index b97c1c3..32015ec 100644
--- a/modules/fun_io.py
+++ b/modules/fun_io.py
@@ -11,6 +11,8 @@
 import keyword
 import inspect
 import time
+from netCDF4 import Dataset
+
 
 
 # ------------ #
@@ -43,10 +45,11 @@ def load_coords():
   from fun_gen import outdir
 
   try :
-    ds=xr.open_dataset(glob(outdir+'/*TEMP*.nc')[0])
-    lon=np.array(ds['longitude'])
-    lat=np.array(ds['latitude'])
-    levels=np.array(ds['depth'])
+    search = outdir+'/*TEMP*.nc'
+    ds = xr.open_dataset(glob(search)[0])
+    lon = np.array(ds['longitude'])
+    lat = np.array(ds['latitude'])
+    levels = np.array(ds['depth'])
     ds.close()
 
   except Exception as e :
@@ -73,7 +76,7 @@ def get_var_2D(fname,var,lev):
 # ----------------- #
 # Load 3D variables #
 # ----------------- #
-def get_var_3D(fname,var,**kargs):
+def get_var_3D(fname,hours,var,**kargs):
 
   dmn = kargs['domain']
 
@@ -84,11 +87,14 @@ def get_var_3D(fname,var,**kargs):
   try:
 
     ds = xr.open_dataset(fname)
+
+    # Reduce domain
     if dmn:
-      arr = ds[var][23,idz,idy,idx].squeeze()
+      arr = ds[var][hours,idz,idy,idx]
+      arr = arr.squeeze()
 
     else:
-      arr = ds[var][23,:,:,:].squeeze()
+      arr = ds[var][hours,:,:,:].squeeze()
 
     ds.close()
 
@@ -101,7 +107,7 @@ def get_var_3D(fname,var,**kargs):
 # ------------------------ #
 # Get value from 3D fields #
 # ------------------------ #
-def get_model_val_3d(fname,var,lon_mod,lat_mod,lev_mod,lon,lat,depth):
+def get_model_val_3d(fname,hours,var,lon_mod,lat_mod,lev_mod,lon,lat,depth):
 
   start = time.time()
 
@@ -113,8 +119,8 @@ def get_model_val_3d(fname,var,lon_mod,lat_mod,lev_mod,lon,lat,depth):
   idy = np.where( (lat_mod > lat[0]-dump) & (lat_mod < lat[0]+dump) )
   idz = np.where( (lev_mod < depth[-1]+800))
 
-  var3d = np.array(get_var_3D(fname,var,domain=[idz,idy,idx]))
-  
+  var3d = np.array(get_var_3D(fname,hours,var,domain=[idz,idy,idx]))
+
   # Prepare interpolation
   LAT,LEV,LON = np.meshgrid(lat_mod[idy],lev_mod[idz],lon_mod[idx]) 
   LAT,LON,LEV = LAT.flatten(),LON.flatten(),LEV.flatten()
diff --git a/modules/fun_plot_profile.py b/modules/fun_plot_profile.py
new file mode 100644
index 0000000..ad3dfc6
--- /dev/null
+++ b/modules/fun_plot_profile.py
@@ -0,0 +1,81 @@
+##################################
+# Routine plotting profile plots #
+##################################
+
+from fun_gen import *
+from fun_io import *
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+# -------------
+# Plot profiles
+# -------------
+def plot_profiles(*arg):
+
+  from fun_gen import fig_prox,\
+                      fig_proy,\
+                      col_phy,col_bcg,\
+                      argo_ds
+  
+  if argo_ds == 'bgc':
+    colors = col_bcg
+  else: 
+    colors = col_phy
+
+
+  # Get number of profiles
+  nprof = int((len(arg)-3)/3)
+  fout = arg[0]
+  time = arg[1]
+  depth = arg[2]
+
+  fig = plt.figure(figsize=(float(fig_prox), float(fig_proy)))
+
+  # Create axes
+  axes = []
+  ax1 = fig.add_subplot(111)
+  ax1.set_ylabel('Depth (m)')
+  axes.append(ax1)
+
+  for i in range(1,nprof):
+    axes.append(ax1.twiny())
+
+  # Plot
+  for i in range(0,nprof):
+    ax = axes[i]
+
+    offset = 3+i*3
+    model,data,label = arg[offset],arg[offset+1],arg[offset+2]
+
+    ax.plot(model,depth,c=colors[i],label=label)
+    ax.plot(data,depth,c=colors[i],marker='.',linestyle='None',alpha=0.5)
+
+    ax.set_xlabel(label)
+    ax.xaxis.label.set_color(colors[i])
+    ax.tick_params(axis='x',colors=colors[i], which='both')
+
+    if i == 0:
+      ax.spines["bottom"].set_color(colors[i])
+      ax.patch.set_alpha(0.0)
+
+    if i > 0:
+      ax.spines["top"].set_color(colors[i])
+
+    if i > 1:
+      ax.spines["top"].set_position(("axes", 1.0+i*0.04))
+     
+  plt.title(time)  
+
+  # Save figure
+  savefig(fout)
+  print('')
+
+  plt.close()
+
+
+  
+
+
+
+