PosProAX97.py

import os
import sys
import numpy as np
#np.bool = np.bool_
import xarray as xr
from scipy.io import savemat, loadmat
from scipy.interpolate import griddata
import scipy.io
from datetime import datetime
import matplotlib.pyplot as plt
import time as pytime
import subprocess
import glob
import gsw
import cmocean
import pandas as pd
#
# MOVARTEAM LIBRARIES # These files (.py) must be in the same directory
from smooth2d_loess import smooth2d_loess
from Calc_sal_Thacker_Goes_EmDr_Stom_svd_globe import make_sal


########################################################
########################################################
# BLOCK 1
########################################################
########################################################
dire = '/data/MOVAR/sumario_cruzeiros/02_proc_CH14/files_netcdf'

sys.path.append(dire)

# LIST .NC FILES
names = glob.glob(os.path.join(dire, '*.nc'))
le = len(names)
names_sorted = sorted(names, key=lambda x: os.path.basename(x).split('_')[-1])
#
########################################################
# DATA LOADING
dados_FRE = []
for ii in range(le):
    ds = xr.open_dataset(names[ii])
    dados_FRE.append({
        'long': ds['longitude'].values,
        'prof': ds['depth'].values,
        'temp': ds['temperature'].values
    })
########################################################
# GENERATE REFERENCE RADIAL - AX97, IT'S NECESSARY TO CHANGE FOR OTHER RADIALS
longitude = np.concatenate([
    np.arange(-40.9, -40.0, 0.16667),
    np.arange(-40.0, -29.6, 0.25),
    np.array([-29.65])
])[:, np.newaxis]
########################################################
p = np.polyfit([-40.9, -29.65], [-23, -20.55], 1)
latitude = np.polyval(p, longitude)
depth = np.arange(0, 810, 10)
A, B = np.meshgrid(longitude.flatten(), depth)
########################################################
# BATHYMETRY MASK FOR AX97, IT'S NECESSARY TO CHANGE FOR OTHER RADIALS
mascara = "/data2/Werner/analisesmovar/MOVAR_analises/mascara.mat"
mascara2 = scipy.io.loadmat(mascara)
mask1 = mascara2['mask1']
########################################################
results = []  # Storage for all loop results

# MAIN LOOP, READING ALL LONGITUDES, DEPTHS, AND TEMPERATURES
for i in range(len(dados_FRE)):
    long = dados_FRE[i]['long'].astype(float)
    prof = dados_FRE[i]['prof'].astype(float)
    temp = dados_FRE[i]['temp'].astype(float)

    # Set zero values to NaN
    temp[temp == 0] = np.nan
    prof[prof == 0] = np.nan

    # Check and adjust array shapes
    if temp.shape != long.shape:
        print('TEMP and LONG are not the same size, fixing it')
        long = np.tile(long, (1, temp.shape[1]))

    if temp.shape != prof.shape:
        print('TEMP and PROF are not the same size, fixing it')
        prof = np.tile(prof, (2, temp.shape[0])).T

    j = np.where(~np.isnan(temp))
    y = np.where(~np.isnan(prof))

    # Set first elements of depth to zero if they are NaN
    if np.isnan(prof[0, 0]) or np.isnan(prof[0, 3]):
        prof[:, 0] = 0
        prof[0, :] = 0

    method_interp = 'nearest'
    points = np.column_stack([long[y], prof[y]])
    values = temp[y]
    T12 = griddata(points, values, (A, B), method=method_interp)

    T12_float = T12.astype(np.float64)
    longitude_float = longitude.flatten().astype(np.float64)
    depth_float = depth.astype(np.float64)

    T, _ = smooth2d_loess(T12_float, longitude_float, depth_float,
                          np.float64(1), np.float64(50),
                          longitude_float, depth_float)

    # Apply bathymetry mask
    mask_corrigida = np.zeros_like(T)
    mask_corrigida[:, :48] = mask1
    T[mask_corrigida == 1] = np.nan

    # CALCULATE SALINITY WITH make_sal
    TT = T.copy()
    PP = np.tile(depth[:, np.newaxis], (1, TT.shape[1]))
    lat = latitude.flatten()
    lon = longitude.flatten()

    arquivo_mais_recente = names_sorted[-1]  # last file in chronological order
    nome_arquivo = os.path.basename(arquivo_mais_recente)
    
    # Extract date YYYYMM from filename
    data_str = nome_arquivo[19:-3]
    
    # Create time vector for all latitudes
    timet = np.full(len(latitude), float(data_str))    

    #time_str = os.path.basename(names[i])[19:-3]  # Ex: '201005'
    #timet = np.full(len(latitude), float(time_str))

    Pad = 0
    Pout = depth.copy()
    method_sal = 'Goes'  # or 'svd', 'Thacker', etc.

    sal = make_sal(TT, PP, lat, lon, timet, pad=Pad, Pout=Pout, method=method_sal)
    
    print(f'Type of make_sal return: {type(sal)}')
    
    if isinstance(sal, tuple):
        print(f'Size of returned tuple: {len(sal)}')
        for idx, item in enumerate(sal):
            print(f'Type of tuple item {idx}: {type(item)}')
            if hasattr(item, 'shape'):
                print(f'Shape of item {idx}: {item.shape}')
        # Assume the first element of the tuple is the salinity array:
        sal = sal[0]
    
    print(f'Calculated salinity with shape: {sal.shape}')
    
    # IN-MEMORY STORAGE
    results.append({
        'temperature': TT,
        'salinity': sal,
        'depth': depth,
        'lat': lat,
        'lon': lon,
        'time': timet
    })
  
########################################################
########################################################
# BLOCK 2
########################################################
########################################################
# Output for Dynamic Height and Vgeo plots
output_dir = "/data2/Movar_analises/arquivos_pro_sal/plots/"

#
sa_list = []
ct_list = []
dyn_height_list = []
all_v_geo = []

# Lists to store transport results
transporte_sul_lista = []        # south transport only for negative velocities
transporte_sul_Sv_lista = []

transporte_liquido_lista = []    # net transport (south - north)
transporte_liquido_Sv_lista = []

# Climatology
ref_dep = 500
time = timet
data_str = time
time_val = int(timet[0])          
time_str = str(time_val)    
time_dt = datetime.strptime(time_str, '%Y%m') # Convert to datetime
month = time_dt.month             # extract the month correctly
print(f'Extracted month: {month}')
#  
nc_file = "/data2/tayannepf/toolbox_movar/TS/argo_CLIM_2005-2016.nc"
ds = xr.open_dataset(nc_file, decode_times=False)
ds['LONGITUDE'] = (ds['LONGITUDE'] + 360) % 360
lon_argo = ds['LONGITUDE'].values
lat_argo = ds['LATITUDE'].values
level = ds['LEVEL'].values
time_argo = ds['TIME'].values + 1
month_index = np.argmin(np.abs(time_argo - month))
level_index = np.argmin(np.abs(level - ref_dep))

########################################################
########################################################
# LOOPING THROUGH PROCESSED DATA + CLIMATOLOGY
########################################################
########################################################

print(f"Total number of results in `results`: {len(results)}")

for i, result in enumerate(results):
      
    try:
        
        print(f"Processing file {i+1} of {len(results)}")
    
        time_value = result['time'][0]  # get the first value of the time array
        data_str = str(time_value)      # convert to string
        
        s2 = np.array(result['salinity'], dtype=np.float64)
        tt = np.array(result['temperature'], dtype=np.float64)
        pp = np.array(result['depth'], dtype=np.float64)
        lon = np.array(result['lon'], dtype=np.float64).flatten()
        lat = np.array(result['lat'], dtype=np.float64).flatten()
        timet = result['time']
    
        #p = pp.T
        p = pp.reshape(-1, 1)    
        lon2d = lon.reshape(1, -1)
        lat2d = lat.reshape(1, -1)
    
        # Save original NaN masks
        nan_mask_s2 = np.isnan(s2)
        nan_mask_tt = np.isnan(tt)
    
        sa = gsw.SA_from_SP(s2, p, lon2d, lat2d)
        ct = gsw.CT_from_t(sa, tt, p)
        sigma_theta = gsw.sigma0(sa, ct)
    
        dyn_height_100 = gsw.geo_strf_dyn_height(sa, ct, p, 100)
        dyn_height_180 = gsw.geo_strf_dyn_height(sa, ct, p, 180)
        dyn_height_500 = gsw.geo_strf_dyn_height(sa, ct, p, 500)
    
        dyn_height_all = np.zeros((81, 49))  # adjust if dimensions are different
        dyn_height_all[:, 0] = dyn_height_100[:, 0]
        dyn_height_all[:, 1] = dyn_height_180[:, 1]
        dyn_height_all[:, 2:] = dyn_height_500[:, 2:]
    
        # CALCULATE GEOPOTENTIAL ANOMALY
        F = dyn_height_all.shape[1]
        if 'lon' in locals():
            dxx = np.abs(np.gradient(lon)) * 111320  # in meters
        else:
            # If there's no longitude, assume a fictitious uniform spacing (e.g., 10 km)
            dxx = np.ones(F) * 10000  # 10 km between stations
        
        gpan = dyn_height_all.copy()
        fi, fj = np.where(gpan == 0.0)
        f = np.where(gpan == 0.0)
        
        # Mark missing values as NaN
        gpan[f] = np.nan
        
        kp = 50
        # Zero the reference level in all columns
        gpan[kp, :] = 0.0
        
        lf = len(fi)
        as_ = np.abs(fj - F / 2)
        saa = np.argsort(as_)
        fi = fi[saa]
        fj = fj[saa]
        
        for l in range(lf):
            lj = fj[l]
            lj1 = min(F - 1, lj + 1)
            lj2 = min(F - 1, lj + 2)
            ljm1 = max(0, lj - 1)
            ljm2 = max(0, lj - 2)
        
            dgpan = gpan[fi[l]:kp + 1, lj2] - gpan[fi[l]:kp + 1, lj1]
            dgpan = 1 / (2 * dxx[lj1]) * dgpan
        
            signo = -1
            gpan[fi[l]:kp + 1, lj] = gpan[fi[l]:kp + 1, lj1] + signo * dgpan * dxx[lj]
        
            if np.isnan(gpan[fi[l], lj]):
                signo = 1
                dgpan = gpan[fi[l]:kp + 1, ljm1] - gpan[fi[l]:kp + 1, ljm2]
                dgpan = 1 / (2 * dxx[ljm1]) * dgpan
                gpan[fi[l]:kp + 1, lj] = gpan[fi[l]:kp + 1, ljm1] + signo * dgpan * dxx[lj]
        
            gpan[:fi[l], lj] = gpan[:fi[l], lj] + gpan[fi[l], lj]
    
        # Argo Climatology
        lon = (lon + 360) % 360
        lonmin_idx = np.argmin(np.abs(lon_argo - lon.min()))
        lonmax_idx = np.argmin(np.abs(lon_argo - lon.max()))
        latmin_idx = np.argmin(np.abs(lat_argo - lat.min()))
        latmax_idx = np.argmin(np.abs(lat_argo - lat.max()))
    
        if lonmin_idx > lonmax_idx:
            lonmin_idx, lonmax_idx = lonmax_idx, lonmin_idx
        if latmin_idx > latmax_idx:
            latmin_idx, latmax_idx = latmax_idx, latmin_idx
    
        addep_sub = ds['ADDEP'].isel(
            LONGITUDE=slice(lonmin_idx, lonmax_idx + 1),
            LATITUDE=slice(latmin_idx, latmax_idx + 1),
            LEVEL=level_index,
            TIME=month_index
        ).values
    
        addep_sub = np.nan_to_num(addep_sub, nan=np.nan)
        lon_sub = lon_argo[lonmin_idx:lonmax_idx + 1]
        lat_sub = lat_argo[latmin_idx:latmax_idx + 1]
        lon_grid, lat_grid = np.meshgrid(lon_sub, lat_sub)
        lon_flat = lon_grid.ravel()
        lat_flat = lat_grid.ravel()
        addep_flat = addep_sub.T.ravel()
        valid_mask = ~np.isnan(addep_flat) & (addep_flat != 0)
    
        addep_interp = griddata(
            points=np.array([lon_flat[valid_mask], lat_flat[valid_mask]]).T,
            values=addep_flat[valid_mask],
            xi=np.array([lon, lat]).T,
            method='linear'
        )
    
        if np.any(np.isnan(addep_interp)):
            addep_interp_nearest = griddata(
                points=np.array([lon_flat[valid_mask], lat_flat[valid_mask]]).T,
                values=addep_flat[valid_mask],
                xi=np.array([lon, lat]).T,
                method='nearest'
            )
            addep_interp[np.isnan(addep_interp)] = addep_interp_nearest[np.isnan(addep_interp)]
    
        addep_interp = addep_interp #print(ds['ADDEP'].attrs) / LEAVING 'units': 'dynamic meters' =  = dm²/s²  TO --> (m²/s²)
        dyn_height_total = gpan + addep_interp[np.newaxis, :]
        dyn_height_total[nan_mask_s2 | nan_mask_tt] = np.nan
        #####
        dyn_height_total_m = dyn_height_total/9.81
        ###################################################################################### 
        dyn_height_all_m = (dyn_height_all / 9.81)      
        #dyn_height_total_cm = (dyn_height_total_m) * 100
        #dyn_height_m = (dyn_height_all / 9.81)
        #dyn_height_cm = (dyn_height_all / 9.81) * 100
        ######################################################################################   
        v_geo = gsw.geostrophic_velocity(dyn_height_total, lon, lat)[0] # It cuts one dimension 81x48      
        v_geo_cm = v_geo * 100
        ######################################################################################   
        lon_plot = (lon + 180) % 360 - 180    
        ######################################################################################
        # Append variables to the list
        sa_list.append(sa)
        ct_list.append(ct)
        dyn_height_list.append(dyn_height_total)
        all_v_geo.append(v_geo)
        ######################################################################################
        ######################################################################################
        #################### CALCULATION OF GEOSTROPHIC TRANSPORT ############################
        ######################################################################################
        ######################################################################################
        dx = gsw.distance(lon, lat)                # vector with distances between each point
        dx = np.array(dx)[np.newaxis, :]         # reshape to (1, n_profiles)
        
        # 2. Calculate vertical intervals (dz) between depth levels [m]
        dz = np.diff(pp)                           # difference between depths (pp)
        dz = np.append(dz, dz[-1])                 # keeps the last interval to close
        dz = dz[:, np.newaxis]                     # reshape to (n_depth, 1)
        dz = dz[:v_geo.shape[0]]                 # adjust to the correct dimension
        
        # --- Southward transport (negative velocities only) ---
        # 3a. Apply mask: keep only negative velocities (southward)
        vg_sul = np.where(v_geo < 0, v_geo, 0) # zero out positive velocities
        
        # 4a. Expand dx and dz for broadcasting to dimensions (n_depth, n_profiles)
        dx_exp = np.repeat(dx, vg_sul.shape[0], axis=0)
        dz_exp = np.repeat(dz, vg_sul.shape[1], axis=1)
        
        # 5a. Calculate local flux (m³/s) for southward transport
        fluxo_local_sul = vg_sul * dx_exp * dz_exp
        
        # 6a. Sum all local fluxes up to 500 m (50 levels) and 10 horizontal profiles
        transporte_sul_m3s = np.nansum(fluxo_local_sul[0:50, 0:10])
        
        # 7a. Convert to Sverdrups (1 Sv = 1e6 m³/s)
        transporte_sul_Sv = transporte_sul_m3s / 1e6
        
        # 8a. Store results
        transporte_sul_lista.append(transporte_sul_m3s)
        transporte_sul_Sv_lista.append(transporte_sul_Sv)
        
        print(f"Southward transport: {transporte_sul_m3s:.2e} m³/s ({transporte_sul_Sv:.2f} Sv)")
        
        
        # Net transport (considering positive and negative velocities) ---
        # 3b. Use all velocities, no filtering
        vg_liquido = v_geo
        
        # 4b. Expand dx and dz for broadcasting
        dx_exp = np.repeat(dx, vg_liquido.shape[0], axis=0)
        dz_exp = np.repeat(dz, vg_liquido.shape[1], axis=1)
        
        # 5b. Calculate net local flux (m³/s)
        fluxo_local_liquido = vg_liquido * dx_exp * dz_exp
        
        # 6b. Sum all local fluxes up to 500 m and 10 profiles
        transporte_liquido_m3s = np.nansum(fluxo_local_liquido[0:50, 0:10])
        
        # 7b. Convert to Sverdrups
        transporte_liquido_Sv = transporte_liquido_m3s / 1e6
        
        # 8b. Store results
        transporte_liquido_lista.append(transporte_liquido_m3s)
        transporte_liquido_Sv_lista.append(transporte_liquido_Sv)
        
        print(f"Net transport (south-north balance): {transporte_liquido_m3s:.2e} m³/s ({transporte_liquido_Sv:.2f} Sv)")
        
          ######################################################################################
        # PLOTTING!
        ## Plot 1 - Dynamic Height
        plt.figure(figsize=(10, 3))
        plt.plot(lon_plot, dyn_height_all_m[0, :], label='Relative', color='blue')
        plt.plot(lon_plot, dyn_height_total_m[0, :], label='Total (with climatology)', color='red')
        plt.xlabel("Longitude")
        plt.ylabel("Dynamic Height (cm)")
        plt.title(f"Dynamic Height: {data_str}")
        plt.grid(True)
        plt.legend()
        plt.tight_layout()
        plt.savefig(os.path.join(output_dir,f"altura_dinamica_{data_str}_{i+1}.png"))
        #plt.show()
        plt.close()
        
        ## Plot 4 - Contour + mask
        mask_gray = np.where(np.isnan(v_geo_cm), 1, np.nan)
        plt.figure(figsize=(8, 3))
        
        plt.imshow(mask_gray,extent=[lon_plot[:-1].min(), lon_plot[:-1].max(), p[:, 0].max(), p[:, 0].min()],cmap='gray',
            alpha=0.5,origin='upper',aspect='auto')
        
        cs = plt.contourf(lon_plot[:-1], p[:, 0], v_geo_cm,vmin=-75,vmax=75,levels=np.arange(-75,80,5), cmap='cmo.balance')
        
        #cs = plt.contourf(lon_plot[:-1], p[:, 0], v_geo_cm,levels=20, cmap='cmo.balance')
        contours = plt.contour(lon_plot[:-1], p[:, 0], v_geo_cm, levels=10, colors='k', linewidths=0.5)
        plt.colorbar(cs, label='v_geo (cm/s)')
        plt.xlabel('Station')
        plt.ylabel('Depth (m)')
        plt.title(f'VGEO: {data_str}')
        plt.tight_layout()
        plt.savefig(os.path.join(output_dir, f"vgeo_contour_masked_{data_str}_{i+1}.png"))
        #plt.show()
        plt.close()
        
    ######################################################################################
    ######################################################################################
    except Exception as e:
        print(f"Error in file {i+1}: {e}")


# Final stacking for NetCDF
SA_final = np.stack(sa_list, axis=-1)
CT_final = np.stack(ct_list, axis=-1)
dyn_height_final = np.stack(dyn_height_list, axis=-1)
v_geo_final = np.stack(all_v_geo, axis=-1)

# Check how many files were read:
print(f"Number of files: {len(results)}")
print("Processing and final saving completed for all files.")

########################################################
########################################################
# BLOCK 4 - STATISTICS
########################################################
########################################################
# Calculate the max and min values for SA_list
SA_max = np.nanmax(sa_list)
SA_min = np.nanmin(sa_list)

# Calculate the max and min values for CT_list
CT_max = np.nanmax(ct_list)
CT_min = np.nanmin(ct_list)

# Calculate the max and min values for dyn_height_list
dyn_height_max = np.nanmax(dyn_height_list)
dyn_height_min = np.nanmin(dyn_height_list)

# Calculate the max and min values for all_v_geo
v_geo_max = np.nanmax(all_v_geo)
v_geo_min = np.nanmin(all_v_geo)

# Display the results
print(f"SA - Max value: {SA_max}, Min value: {SA_min}")
print(f"CT - Max value: {CT_max}, Min value: {CT_min}")
print(f"Dyn Height - Max value: {dyn_height_max}, Min value: {dyn_height_min}")
print(f"VGeo - Max value: {v_geo_max}, Min value: {v_geo_min}")
########################################################
# ALL VARIABLES FOR NETCDF / Stacking along the last dimension (time)
SA_final = np.stack(sa_list, axis=-1)  
CT_final = np.stack(ct_list, axis=-1)
dyn_height_final = np.stack(dyn_height_list, axis=-1)
v_geo_final = np.stack(all_v_geo, axis=-1)
########################################################
v_geo_stack = np.array(all_v_geo)  
v_geo_mean = np.nanmean(v_geo_stack, axis=0)  
v_geo_mean_max = np.nanmax(v_geo_mean)
v_geo_mean_min = np.nanmin(v_geo_mean)
print(f"VGeo - Max value: {v_geo_mean_max}, Min value: {v_geo_mean_min}")

v_geo_mean_cm = v_geo_mean*100
v_geo_mean_cm_max = np.nanmax(v_geo_mean_cm)
v_geo_mean_cm_min = np.nanmin(v_geo_mean_cm)
print(f"VGeo - Max value: {v_geo_mean_cm_max}, Min value: {v_geo_mean_cm_min}")

########################################################
########################################################
# TRANSPORT STATISTICS
########################################################
########################################################

# Southward transport statistics
print(f"\nNumber of files (southward transport): {len(transporte_sul_Sv_lista)}")
transporte_sul_array = np.array(transporte_sul_Sv_lista)
media_sul = np.mean(transporte_sul_array)
desvio_sul = np.std(transporte_sul_array)
print(f"Mean Southward Transport (Sv): {media_sul:.3f}")
print(f"Standard Deviation Southward Transport (Sv): {desvio_sul:.3f}")

# Net transport statistics
print(f"\nNumber of files (net transport): {len(transporte_liquido_Sv_lista)}")
transporte_liquido_array = np.array(transporte_liquido_Sv_lista)
media_liquido = np.mean(transporte_liquido_array)
desvio_liquido = np.std(transporte_liquido_array)
print(f"Mean Net Transport (Sv): {media_liquido:.3f}")
print(f"Standard Deviation Net Transport (Sv): {desvio_liquido:.3f}")


########################################################
########################################################
# BLOCK 5 - NETCDF CREATION
########################################################
########################################################
# Output directory setup
output_dir = "/data2/Movar_analises/arquivos_pro_sal/NETCDF/"
os.makedirs(output_dir, exist_ok=True)
####
# Depth
depth_netcdf = np.arange(0, 810, 10)  

print("SA_final shape:", SA_final.shape)
print("CT_final shape:", CT_final.shape)
print("dyn_final shape:", dyn_height_final.shape)
print("vgeo_final shape:", v_geo_final.shape)
print("Depth:", len(depth_netcdf), "Station:", len(lon_plot), "Time:", len(time))
print("Time length:", len(time))
########################################################
# Block for NetCDF creation
output_nc_file = os.path.join(output_dir, "AX97_MOVAR.nc")

# Creating coordinates for v_geo (between stations)
station_vgeo = np.arange(len(lon_plot)-1)  # 48 stations for v_geo
lon_vgeo = (lon_plot[:-1] + lon_plot[1:]) / 2  # Mean longitude between stations
lat_vgeo = (lat[:-1] + lat[1:]) / 2  # Mean latitude between stations

# Building the time dimension
time_list = [int(r['time'][0]) for r in results]
time = pd.to_datetime([f"{t}01" for t in time_list], format="%Y%m%d")
ordem = np.argsort(time)
time = time[ordem]
time_str = np.array([str(t) for t in time_list], dtype='U6')[ordem]


SA_final = SA_final[:, :, ordem]  
CT_final = CT_final[:, :, ordem]
dyn_height_final = dyn_height_final[:, :, ordem]
v_geo_final = v_geo_final[:, :, ordem]  


########################################################
########################################################

# Creating the xarray dataset
ds_out = xr.Dataset(
    {
        "transporte_sul": (["time"], transporte_sul_Sv_lista),
        "transporte_liquido": (["time"], transporte_liquido_Sv_lista),  
        "SA": (["depth", "station", "time"], SA_final),
        "CT": (["depth", "station", "time"], CT_final),
        "dyn_height": (["depth", "station", "time"], dyn_height_final),
        "v_geo": (["depth", "station_vgeo", "time"], v_geo_final),
    },
    coords={
        "time": time,
        "time_label": ("time", time_str),
        "depth": depth_netcdf,
        "station": np.arange(len(lon_plot)),  
        "station_vgeo": station_vgeo,  
        "lon": (["station"], lon_plot),
        "lat": (["station"], lat),
        "lon_vgeo": (["station_vgeo"], lon_vgeo),
        "lat_vgeo": (["station_vgeo"], lat_vgeo),
    }
)

# Adding attributes (keeping the same as the previous code)
ds_out.transporte_sul.attrs = {
    "long_name": "Total geostrophic south transport",
    "units": "Sv",
    "description": "Integrated volume transport calculated from geostrophic velocities"
}

ds_out.transporte_liquido.attrs = {
    "long_name": "Total geostrophic transport liquid",
    "units": "Sv",
    "description": "Integrated volume transport calculated from geostrophic velocities"
}


ds_out.SA.attrs = {
    "long_name": "Absolute Salinity",
    "units": "g/kg",
    "standard_name": "sea_water_absolute_salinity"
}

ds_out.CT.attrs = {
    "long_name": "Conservative Temperature",
    "units": "°C",
    "standard_name": "sea_water_conservative_temperature"
}

ds_out.dyn_height.attrs = {
    "long_name": "Dynamic height",
    "units": "cm",
    "description": "Dynamic height relative to 500 dbar"
}

ds_out.v_geo.attrs = {
    "long_name": "Geostrophic velocity",
    "units": "m/s",
    "description": "Geostrophic velocity perpendicular to the section, calculated between stations"
}

# Adding attributes for the new coordinates
ds_out.lon_vgeo.attrs = {
    "long_name": "Longitude at v_geo points",
    "units": "degrees_east",
    "description": "Longitude at midpoint between stations where v_geo is calculated"
}

ds_out.lat_vgeo.attrs = {
    "long_name": "Latitude at v_geo points",
    "units": "degrees_north",
    "description": "Latitude at midpoint between stations where v_geo is calculated"
}

# Replacing possible NaNs with filled values (fill_value)
fill_value = -9999.0
ds_out = ds_out.fillna(fill_value)

# Saving to NetCDF
encoding = {
    var: {"_FillValue": fill_value} for var in ds_out.data_vars
}

# Add time variable encoding
encoding["time"] = {
    "units": "days since 2004-08-01 00:00:00", # START DATE FOR AX97, SET THE START DATE OF YOUR CRUISE!
    "calendar": "gregorian",
    "dtype": "float64",
    "_FillValue": None
}

# Adding global attributes (signature)
ds_out.attrs["title"] = "Geostrophic Transport Dataset"
ds_out.attrs["summary"] = "Geostrophic transport and oceanographic variables from section analysis"
ds_out.attrs["institution"] = "- LOF/IGEO/UFRJ"
ds_out.attrs["creator_name"] = "Werner de Barros (PhD Student), Tayanne Pires (PhD Student), Mauro Cirano (UFRJ Teacher)"
ds_out.attrs["creator_email"] = "wbarros@oceanica.ufrj.br"
ds_out.attrs["date_created"] = np.datetime_as_string(np.datetime64("now"), unit="s")
ds_out.attrs["Conventions"] = "CF-1.8"
ds_out.attrs["history"] = f"Created on {ds_out.attrs['date_created']} with xarray"
ds_out.to_netcdf(output_nc_file, encoding=encoding)

print(f"NetCDF file saved to: {output_nc_file}")

# END