It depends on your data structure, ( NEMO data, projected coordinates, data gird). Anyway, here is a simple example where ur lon and lat are just arrays not matrices, your main variable temp is 4D ( lon,lat,depth, time). You can play with the code to fit it to your data structure.
%% read the older file
file = 'global-reanalysis-phy_2012_Jan_Apr_fill_domain_SST_currents_dep.nc';
lon = ncread(file,'longitude');
lat = ncread(file,'latitude');
time = ncread(file,'time');
depth = ncread(file,'depth');
temp = ncread(file,'thetao');
%% create the new file
ncnc=netcdf.create(['new','.nc'] ,'NC_WRITE')
lonID=netcdf.defDim(ncnc,'lon', length(lon));
vlongID=netcdf.defVar(ncnc,'lon','double',lonID);
netcdf.putAtt(ncnc,vlongID,'units','degrees');
netcdf.putAtt(ncnc,vlongID,'axis','long');
depID=netcdf.defDim(ncnc,'depth',length(depth));
vdepID=netcdf.defVar(ncnc,'depth','double',depID);
netcdf.putAtt(ncnc,vdepID,'units','m');
latID=netcdf.defDim(ncnc,'lat', length(lat));
vlatID=netcdf.defVar(ncnc,'lat','double',latID);
netcdf.putAtt(ncnc,vlatID,'units','degrees');
netcdf.putAtt(ncnc,vlatID,'axis','lat');
tID=netcdf.defDim(ncnc,'time',length(time));
vtID=netcdf.defVar(ncnc,'time','double',tID);
netcdf.putAtt(ncnc,vtID,'units','days since 1950-1-1');
veddyID=netcdf.defVar(ncnc,'temp','double',[lonID,latID,depID, tID]);
netcdf.putAtt(ncnc,veddyID,'units','C');
netcdf.endDef(ncnc)
% start put data
netcdf.putVar(ncnc,vdepID,depth);
netcdf.putVar(ncnc,vtID,time);
netcdf.putVar(ncnc,vlatID,lat);
netcdf.putVar(ncnc,vlongID,lon);
netcdf.putVar(ncnc,veddyID,temp);
netcdf.close(ncnc)