%
% Lectura de archivos de Theta modelo ECCO
% Orpa Nov 2020
%

% Camino para archivos de temperatura
camino = 'E:\RESPALDO\D\orpa\Data\ECCO\THETA_nc\';
archi = dir(camino);
% Prueba leer primer archivo de la lista y le agrega el path
% archi1 = getfield(archi,{3,1},'name');
% archi1 = [camino archi1];

%Lectura_Basica_NetCDF(archi1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Abre archivo para leer lat lon prof y 
% obtener parámetros para definir región a leern
ncid = netcdf.open(archi1, 'NOWRITE');
% Lee Lat, Lon y Theta
% Variables en archivo
% 0 Time
% 1 Depth_T        Profunidades estándares
% 2 Latitude       cada 0.25
% 3 Longitude      cada 0.25
% 4 Theta   
%  FillValue -> -9.9999998e+22
%  Missing_value -> -9.9999998e+22

t = netcdf.getVar(ncid,0,'double');
z = netcdf.getVar(ncid,1,'double');
lat = netcdf.getVar(ncid,2,'double');
lon = netcdf.getVar(ncid,3,'double');

% Define región a leer
indlon = find(lon >= -120 &  lon <= -69.5);
indlat = find(lat >= -60 &  lat <= -15);
indz = find(z <= 70000 );
narchi = max(size(archi))-2;   % Número de archivos es -2 por (. y ..)
% Lee región dadda por índices [lat lon prof tiempo]  
auxlec1 = [indlon(1)-1, indlat(1)-1 , indz(1)-1 , 0];
auxlec2 = [length(indlon),length(indlat), length(indz), 1];
% Cierra el archivo 
netcdf.close(ncid)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COMIENZA ITERACIÓN PARA LEER REGIÓN EN TODOS LOS ARCHIVOS
% Prelocaliza espacio en memoria para theta
theta = zeros(length(indlon), length(indlat), length(indz), narchi);
% EXTRE REGIÓN
% Itera para todos los archivos 

for i = 1:10
    archi1 = [camino getfield(archi,{3,1},'name')];
    ncid = netcdf.open(archi1, 'NOWRITE');
    % lee archivo
    theta(:,:,:,i) = netcdf.getVar(ncid,4,auxlec1, auxlec2); % no convierte a double para ahorrar mem
    
    netcdf.close(ncid)
end

% Prueba una lectura y grafica para ver el campo
%theta = netcdf.getVar(ncid,4,auxlec1, auxlec2,'double');
%ind = (theta > 100 | theta < -10);
%theta(ind) = NaN;
% Lee Sección vertical

% % Grafico preliminar para chequeo
%pcolor(lon(indlon),lat(indlat),theta'); shading interp; colorbar; 
%colormap(jet)
