close all

addpath /home/stephdut/matlab/seawater_ver3_0

yearwant=2006;

 dirstr=['/nobackup1/stephdut/MONOD_SIZE_RADTRANS9/'];
 runstr=['run_size163_9spec_radtrans_interannaul360'];
 run1str=[dirstr,runstr,'/MONTHS'];
%%% read grid
 gridname=['grid_ecco.nc']; 
 X=ncread(gridname,'X');
 Y=ncread(gridname,'Y');
 Z=ncread(gridname,'Z');
 drF=ncread(gridname,'drF');
 dyG=ncread(gridname,'dyG');
 dxG=ncread(gridname,'dxG');
 rA=ncread(gridname,'rA');
 HFacC=ncread(gridname,'HFacC');
 Zp1=ncread(gridname,'Zp1');
 Depth=ncread(gridname,'Depth');
%
 drC=-diff(Z);
%
 nx=length(X); ny=length(Y); nz=length(Z);
 for k=1:nz,
   volc(k,:,:)=rA(:,:)*drF(k);
 end
 filand1=find(squeeze(HFacC(:,:,1))==0);
 land=zeros(nx,ny); land(filand1)=1;
     pmask=ones(size(HFacC));
     clear fi, fi=find(HFacC==0);
     pmask(fi)=0;

nlam=13;

%%%% READ MONTHS
 year1=1992; 
 yearnume1=43200; % for ecosystem
 yearnump1=0; % for physics
 yearnume=yearnume1+(yearwant-year1)*2880; % 83520 for 2006
 yearnump=yearnump1+(yearwant-year1)*24*360; 
 it1=0;
 mnum=1;
 for it=[1:12],
   it1=it1+1;
   time=yearnume+it*240;
   tstr=['00000',num2str(time)];
%%% read tracers
   ncid1=netcdf.open([run1str,'/ptr_tave.',tstr,'.nc']);
   [ndim,nvar,natt,unlim]=netcdf.inq(ncid1);
     for i=0:ndim-1,
      [dimname, dimlength]=netcdf.inqDim(ncid1,i);
      diml(i+1)=dimlength;
     end % for i

     for i=0:nvar-1;
      [varname, xtype,dimid,natt]=netcdf.inqVar(ncid1,i);
      if (length(dimid)==1),
        eval([varname,'=netcdf.getVar(ncid1,i,0,diml(dimid+1));']);
      end % if
      if (length(dimid)==2),
        eval([varname,'=netcdf.getVar(ncid1,i,[0 0],[diml(dimid(1)+1) diml(dimid(2)+1)] );']);
      end % if
      if (length(dimid)==3),
        eval([varname,'=netcdf.getVar(ncid1,i,[0 0 0],[diml(dimid(1)+1) diml(dimid(2)+1) diml(dimid(3)+1)] );']);
      end % if
      if (length(dimid)==4),
         eval([varname,'=netcdf.getVar(ncid1,i,[0 0 0 0],[diml(dimid(1)+1) diml(dimid(2)+1) diml(dimid(3)+1)  diml(dimid(4)+1)] );']);
       end % if
        eval([varname,'=double(',varname,');']);
     end % for i
   netcdf.close(ncid1);
%%% read dar variables etc
   ncid1=netcdf.open([run1str,'/dar_tave.',tstr,'.nc']);
   [ndim,nvar,natt,unlim]=netcdf.inq(ncid1);
     for i=0:ndim-1,
      [dimname, dimlength]=netcdf.inqDim(ncid1,i);
      diml(i+1)=dimlength;
     end % for i

     for i=0:nvar-1;
      [varname, xtype,dimid,natt]=netcdf.inqVar(ncid1,i);
      if (length(dimid)==1),
        eval([varname,'=netcdf.getVar(ncid1,i,0,diml(dimid+1));']);
      end % if
      if (length(dimid)==2),
        eval([varname,'=netcdf.getVar(ncid1,i,[0 0],[diml(dimid(1)+1) diml(dimid(2)+1)] );']);
      end % if
      if (length(dimid)==3),
        eval([varname,'=netcdf.getVar(ncid1,i,[0 0 0],[diml(dimid(1)+1) diml(dimid(2)+1) diml(dimid(3)+1)] );']);
      end % if
      if (length(dimid)==4),
         eval([varname,'=netcdf.getVar(ncid1,i,[0 0 0 0],[diml(dimid(1)+1) diml(dimid(2)+1) diml(dimid(3)+1)  diml(dimid(4)+1)] );']);
       end % if
        eval([varname,'=double(',varname,');']);
     end % for i
   netcdf.close(ncid1);
%%%
%%% load physics
   timestep=24*30; timestep1=24;

   atheta=zeros(nx,ny,nz); asalt=zeros(nx,ny,nz); 
   auvel=zeros(nx,ny,nz); avvel=zeros(nx,ny,nz); 
   for itt=1:30,
    time1=yearnump+(it-1)*timestep+(itt-1)*timestep1;
    if (time1<1000000); timestr=['0000',num2str(time1)]; end
    if (time1<100000); timestr=['00000',num2str(time1)]; end
    if (time1<10000); timestr=['000000',num2str(time1)]; end
    if (time1<1000); timestr=['0000000',num2str(time1)]; end
    if (time1<100); timestr=['00000000',num2str(time1)]; end
%
    fid=fopen(['/nobackup1/benw/ECCO_Iter3.73daily/DDtheta/DDtheta.',timestr,'.data'],'r','b');
    field=fread(fid,'float32'); fclose(fid);
    field=reshape(field,nx,ny,nz);
    atheta=atheta+field/30;
%
   end %itt
%%% finish reading physics
%%%
   sst(it,:,:)=atheta(:,:,1);
%%%
%%%%%%%%%%%%%%%
     for ilam=1:nlam,
       if (ilam<10),
         eval(['upwell(:,:,ilam)=squeeze(dar_Eutop0',num2str(ilam),'(:,:,1));']);
         eval(['etot(:,:,',num2str(ilam),')=squeeze(dar_Ed0',num2str(ilam),'(:,:,1))+squeeze(dar_Es0',num2str(ilam),'(:,:,1));']);
       else
         eval(['upwell(:,:,ilam)=squeeze(dar_Eutop',num2str(ilam),'(:,:,1));']);
         eval(['etot(:,:,',num2str(ilam),')=squeeze(dar_Ed',num2str(ilam),'(:,:,1))+squeeze(dar_Es',num2str(ilam),'(:,:,1));']);
       end % if
     end % for ilam
     etot=double(etot);

     reflect(it,:,:,:)=upwell./etot;
%%
end % for it
%%%
%%%%%%%%%%%%
%%
 darwin_waves=[400.0, 425.0, 450.0, 475.0, 500.0, 525.0, 550.0, 575.0, 600.0, 625.0, 650.0, 675.0, 700.0];
 darwin_wavebands=[400.0, 412.5, 437.5, 462.5, 487.5, 512.5, 537.5, 562.5, 587.5, 612.5, 637.5, 662.5, 687.5, 700.0];
 month=[1:12]; nt=length(month);
%%
 ncid=netcdf.create(['ecco_radtrans9_size163_interannaul360_2006_formatt.nc'],'NC_WRITE');
 dimnx=netcdf.defDim(ncid,'longitude',nx);
 dimny=netcdf.defDim(ncid,'latitude',ny);
 dimnlam=netcdf.defDim(ncid,'waveband',nlam);
 dimnlamo=netcdf.defDim(ncid,'wavebandlimits',nlam+1);
 dimnt=netcdf.defDim(ncid,'month',nt);

 long_ID=netcdf.defVar(ncid,'longitude','double',[dimnx]);
 lat_ID=netcdf.defVar(ncid,'latitude','double',[dimny]);
 wvband_ID=netcdf.defVar(ncid,'wavebandcenter','double',[dimnlam]);
 wvbando_ID=netcdf.defVar(ncid,'wavebandlimits','double',[dimnlamo]);
 month_ID=netcdf.defVar(ncid,'month','double',[dimnt]);

 sst_ID=netcdf.defVar(ncid,'SST','double',[dimnt,dimnx,dimny]);
 reflect_ID=netcdf.defVar(ncid,'Reflectance','double',[dimnt,dimnx,dimny,dimnlam]);

 netcdf.endDef(ncid);

 netcdf.putVar(ncid,long_ID,X);
 netcdf.putVar(ncid,lat_ID,Y);
 netcdf.putVar(ncid,wvband_ID, darwin_waves);
 netcdf.putVar(ncid,wvbando_ID, darwin_wavebands);
 netcdf.putVar(ncid,month_ID,month);

 netcdf.putVar(ncid,sst_ID,sst);
 netcdf.putVar(ncid,reflect_ID,reflect);

 netcdf.close(ncid);
