istart=1;
close all

addpath /home/stephdut/matlab/seawater_ver3_0

if (istart==1)
 clear all

 load colorgam.mat

 dirstr=['/nobackup1/stephdut/MONOD_SIZE_RADTRANS9/'];
 runstr=['run_size163_9spec_radtrans_newwab'];
 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;
end % for istart

 kn2=8;

%%%% READ MONTHS
 it1=0;
 mnum=1;
 for it=[1:12],
   it1=it1+1;
   time=25920+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
   istart=78936;
   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=istart+(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/stephdut/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;
%
    fid=fopen(['/nobackup1/stephdut/ECCO_Iter3.73daily/DDsalt/DDsalt.',timestr,'.data'],'r','b');
    field=fread(fid,'float32'); fclose(fid);
    field=reshape(field,nx,ny,nz);
    asalt=asalt+field/30;
%
    fid=fopen(['/nobackup1/stephdut/ECCO_Iter3.73daily/DDuvel/DDuvel.',timestr,'.data'],'r','b');
    field=fread(fid,'float32'); fclose(fid);
    field=reshape(field,nx,ny,nz);
    auvel=auvel+field/30;
%
    fid=fopen(['/nobackup1/stephdut/ECCO_Iter3.73daily/DDvvel/DDvvel.',timestr,'.data'],'r','b');
    field=fread(fid,'float32'); fclose(fid);
    field=reshape(field,nx,ny,nz);
    avvel=avvel+field/30;
%%
   end %itt
%%% finish reading physics
%%%
   theta=atheta;
   salt=asalt;
%%% calculate MLD
    depthgrid=zeros(size(theta));
    for k=1:nz
      depthgrid(:,:,k)=Z(k);
    end % for k
    latgrid=zeros(size(theta));
    for j=1:ny,
      latgrid(:,j,:)=Y(j);
    end % for j
    press=sw_pres(depthgrid,latgrid);
    pressref=zeros(size(theta));
    press0=zeros(size(theta));
    insitut=sw_temp(salt,theta,press,pressref);
    dens=sw_dens(salt,insitut,press);
    sigma=dens-1000;
    sigmat=sw_dens(salt,insitut,press0)-1000;
    sigmatheta=sw_dens(salt,theta,press0)-1000;
    dens=dens.*pmask;
    sigma=sigma.*pmask;
    sigmat=sigmat.*pmask;
    sigmatheta=sigmatheta.*pmask;
%% linearly interpolate
    sigint=zeros(nx,ny,5000);
    ks=ceil(drF(1)/2);
    for k=1:nz-1,
      clear fi, fi=find(squeeze(pmask(:,:,k+1))==0);
      cdz=ceil(drF(k)/2)+floor(drF(k+1)/2);
      rcdz=(drF(k)/2)+(drF(k+1)/2);
      cdz=floor(rcdz);
      for kk=1:cdz,
         tmp=(rcdz-kk+1)*squeeze(sigmatheta(:,:,k))+(kk-1)*squeeze(sigmatheta(:,:,k+1));
         tmp=tmp/rcdz;
         if (kk>(drF(k)/2)), tmp(fi)=0; end
         sigint(:,:,ks)=tmp;
         ks=ks+1;
      end % for kk 
    end % for k
    for i=1:nx, for j=1:ny,
       clear kk
       kk=find(sigint(i,j,:)>sigint(i,j,5)+0.125);
       if (isempty(kk)==0),
          mld(i,j)=kk(1);
       else
          mld(i,j)=0;
       end
     end, end % for i j
     amld=mld;
%
%%%
%% integrate diaz over 1:kn1
      totphy(it1,:,:,1:kn2)=squeeze(Phy01(:,:,1:kn2)+Phy02(:,:,1:kn2)+Phy03(:,:,1:kn2)+ ...
			    Phy04(:,:,1:kn2)+Phy05(:,:,1:kn2)+Phy06(:,:,1:kn2)+ ...
			    Phy07(:,:,1:kn2)+Phy08(:,:,1:kn2)+Phy09(:,:,1:kn2))*120*12/1e3;
      cocco(it1,:,:,1:kn2)=squeeze(Phy09(:,:,1:kn2))*120*12/1e3;
      diatom(it1,:,:,1:kn2)=squeeze(Phy01(:,:,1:kn2))*120*12/1e3;
      otherlarge(it1,:,:,1:kn2)=squeeze(Phy02(:,:,1:kn2))*120*12/1e3;
      syn(it1,:,:,1:kn2)=squeeze(Phy03(:,:,1:kn2))*120*12/1e3;
      pro(it1,:,:,1:kn2)=squeeze((Phy04(:,:,1:kn2)+Phy05(:,:,1:kn2)))*120*12/1e3;
      unidiaz(it1,:,:,1:kn2)=squeeze((Phy06(:,:,1:kn2)))*120*12/1e3;
      tricho(it1,:,:,1:kn2)=squeeze((Phy07(:,:,1:kn2)))*120*12/1e3;
      picoeuk(it1,:,:,1:kn2)=squeeze(Phy08(:,:,1:kn2))*120*12/1e3;
%%
      totchl(it1,:,:,1:kn2)=squeeze(dar_Chl_ave(:,:,1:kn2));
      coccochl(it1,:,:,1:kn2)=squeeze(Chl09(:,:,1:kn2));
      diatomchl(it1,:,:,1:kn2)=squeeze(Chl01(:,:,1:kn2));
      otherlargechl(it1,:,:,1:kn2)=squeeze(Chl02(:,:,1:kn2));
      synchl(it1,:,:,1:kn2)=squeeze(Chl03(:,:,1:kn2));
      prochl(it1,:,:,1:kn2)=squeeze((Chl04(:,:,1:kn2)+Chl05(:,:,1:kn2)));
      unidiazchl(it1,:,:,1:kn2)=squeeze((Chl06(:,:,1:kn2)));
      trichochl(it1,:,:,1:kn2)=squeeze((Chl06(:,:,1:kn2)));
      picoeukchl(it1,:,:,1:kn2)=squeeze(Chl08(:,:,1:kn2));

      din(it,:,:,1:kn2)=squeeze(NO3(:,:,1:kn2)+NO2(:,:,1:kn2)+NH4(:,:,1:kn2));
      po4(it,:,:,1:kn2)=squeeze(PO4(:,:,1:kn2));
      sil(it,:,:,1:kn2)=squeeze(SiO2(:,:,1:kn2));
      fet(it,:,:,1:kn2)=squeeze(FeT(:,:,1:kn2));
%
      spar(it,:,:,1:kn2)=squeeze(dar_PAR_ave(:,:,1:kn2));
%
      stheta(it,:,:,1:kn2)=squeeze(atheta(:,:,1:kn2));
      ssalt(it,:,:,1:kn2)=squeeze(asalt(:,:,1:kn2));
      suvel(it,:,:,1:kn2)=squeeze(auvel(:,:,1:kn2));
      svvel(it,:,:,1:kn2)=squeeze(avvel(:,:,1:kn2));
      smld(it,:,:)=amld;
%%
end % for it
%%%
%%%%%%%%%%%%
%%
 month=[1:12]; nt=length(month);
 ncid=netcdf.create(['ecco_radtrans9_size163_forbcodmo_phy.nc'],'NC_WRITE');
 dimnx=netcdf.defDim(ncid,'longitude',nx);
 dimny=netcdf.defDim(ncid,'latitude',ny);
 dimnz=netcdf.defDim(ncid,'depth',kn2);
 dimnt=netcdf.defDim(ncid,'month',nt);

 long_ID=netcdf.defVar(ncid,'longitude','double',[dimnx]);
 lat_ID=netcdf.defVar(ncid,'latitude','double',[dimny]);
 month_ID=netcdf.defVar(ncid,'month','double',[dimnt]);
 depth_ID=netcdf.defVar(ncid,'depth','double',[dimnz]);
 layer_ID=netcdf.defVar(ncid,'layerthickness','double',[dimnz]);

 cocco_ID=netcdf.defVar(ncid,'Coccolithophore','double',[dimnt,dimnx,dimny,dimnz]);
 totphy_ID=netcdf.defVar(ncid,'TotalPhytoplanktonBiomass','double',[dimnt,dimnx,dimny,dimnz]);
 diatom_ID=netcdf.defVar(ncid,'Diatom','double',[dimnt,dimnx,dimny,dimnz]);
 otherlarge_ID=netcdf.defVar(ncid,'OtherLarge','double',[dimnt,dimnx,dimny,dimnz]);
 picoeuk_ID=netcdf.defVar(ncid,'PicoEukaryote','double',[dimnt,dimnx,dimny,dimnz]);
 pro_ID=netcdf.defVar(ncid,'Prochlorococcus','double',[dimnt,dimnx,dimny,dimnz]);
 syn_ID=netcdf.defVar(ncid,'Synecococcus','double',[dimnt,dimnx,dimny,dimnz]);
 unidiaz_ID=netcdf.defVar(ncid,'UnicellDiaz','double',[dimnt,dimnx,dimny,dimnz]);
 tricho_ID=netcdf.defVar(ncid,'Tricho','double',[dimnt,dimnx,dimny,dimnz]);

 coccochl_ID=netcdf.defVar(ncid,'CoccolithophoreChl','double',[dimnt,dimnx,dimny,dimnz]);
 totchl_ID=netcdf.defVar(ncid,'TotalChl','double',[dimnt,dimnx,dimny,dimnz]);
 diatomchl_ID=netcdf.defVar(ncid,'DiatomChl','double',[dimnt,dimnx,dimny,dimnz]);
 otherlargechl_ID=netcdf.defVar(ncid,'OtherLargeChl','double',[dimnt,dimnx,dimny,dimnz]);
 picoeukchl_ID=netcdf.defVar(ncid,'PicoEukaryoteChl','double',[dimnt,dimnx,dimny,dimnz]);
 prochl_ID=netcdf.defVar(ncid,'ProchlorococcusChl','double',[dimnt,dimnx,dimny,dimnz]);
 synchl_ID=netcdf.defVar(ncid,'SynecococcusChl','double',[dimnt,dimnx,dimny,dimnz]);
 unidiazchl_ID=netcdf.defVar(ncid,'UnicellDiazChl','double',[dimnt,dimnx,dimny,dimnz]);
 trichochl_ID=netcdf.defVar(ncid,'TrichoCell','double',[dimnt,dimnx,dimny,dimnz]);
%%%
 netcdf.endDef(ncid);

 netcdf.putVar(ncid,long_ID,X);
 netcdf.putVar(ncid,lat_ID,Y);
 netcdf.putVar(ncid,month_ID,month);
 netcdf.putVar(ncid,depth_ID,squeeze(Z(1:kn2)));
 netcdf.putVar(ncid,layer_ID,squeeze(drF(1:kn2)));

 netcdf.putVar(ncid,totphy_ID,totphy);
 netcdf.putVar(ncid,cocco_ID,cocco);
 netcdf.putVar(ncid,diatom_ID,diatom);
 netcdf.putVar(ncid,otherlarge_ID,otherlarge);
 netcdf.putVar(ncid,pro_ID,pro);
 netcdf.putVar(ncid,syn_ID,syn);
 netcdf.putVar(ncid,picoeuk_ID,picoeuk);
 netcdf.putVar(ncid,unidiaz_ID,unidiaz);
 netcdf.putVar(ncid,tricho_ID,tricho);


 netcdf.putVar(ncid,totchl_ID,totchl);
 netcdf.putVar(ncid,coccochl_ID,coccochl);
 netcdf.putVar(ncid,diatomchl_ID,diatomchl);
 netcdf.putVar(ncid,otherlargechl_ID,otherlargechl);
 netcdf.putVar(ncid,prochl_ID,prochl);
 netcdf.putVar(ncid,synchl_ID,synchl);
 netcdf.putVar(ncid,picoeukchl_ID,picoeukchl);
 netcdf.putVar(ncid,unidiazchl_ID,unidiazchl);
 netcdf.putVar(ncid,trichochl_ID,trichochl);

 netcdf.close(ncid);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 ncid=netcdf.create(['ecco_radtrans9_size163_forbcodmo_nutrients.nc'],'NC_WRITE');
 dimnx=netcdf.defDim(ncid,'longitude',nx);
 dimny=netcdf.defDim(ncid,'latitude',ny);
 dimnz=netcdf.defDim(ncid,'depth',kn2);
 dimnt=netcdf.defDim(ncid,'month',nt);

 long_ID=netcdf.defVar(ncid,'longitude','double',[dimnx]);
 lat_ID=netcdf.defVar(ncid,'latitude','double',[dimny]);
 month_ID=netcdf.defVar(ncid,'month','double',[dimnt]);
 depth_ID=netcdf.defVar(ncid,'depth','double',[dimnz]);
 layer_ID=netcdf.defVar(ncid,'layerthickness','double',[dimnz]);

 din_ID=netcdf.defVar(ncid,'DIN','double',[dimnt,dimnx,dimny,dimnz]);
 po4_ID=netcdf.defVar(ncid,'PO4','double',[dimnt,dimnx,dimny,dimnz]);
 sil_ID=netcdf.defVar(ncid,'SiO2','double',[dimnt,dimnx,dimny,dimnz]);
 fet_ID=netcdf.defVar(ncid,'FeT','double',[dimnt,dimnx,dimny,dimnz]);
%%%
 netcdf.endDef(ncid);

 netcdf.putVar(ncid,long_ID,X);
 netcdf.putVar(ncid,lat_ID,Y);
 netcdf.putVar(ncid,month_ID,month);
 netcdf.putVar(ncid,depth_ID,squeeze(Z(1:kn2)));
 netcdf.putVar(ncid,layer_ID,squeeze(drF(1:kn2)));

  netcdf.putVar(ncid,din_ID,din);
  netcdf.putVar(ncid,po4_ID,po4);
  netcdf.putVar(ncid,sil_ID,sil);
  netcdf.putVar(ncid,fet_ID,fet);


 netcdf.close(ncid);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 ncid=netcdf.create(['ecco_radtrans9_size163_forbcodmo_physics.nc'],'NC_WRITE');
 dimnx=netcdf.defDim(ncid,'longitude',nx);
 dimny=netcdf.defDim(ncid,'latitude',ny);
 dimnz=netcdf.defDim(ncid,'depth',kn2);
 dimnt=netcdf.defDim(ncid,'month',nt);

 long_ID=netcdf.defVar(ncid,'longitude','double',[dimnx]);
 lat_ID=netcdf.defVar(ncid,'latitude','double',[dimny]);
 month_ID=netcdf.defVar(ncid,'month','double',[dimnt]);
 depth_ID=netcdf.defVar(ncid,'depth','double',[dimnz]);
 layer_ID=netcdf.defVar(ncid,'layerthickness','double',[dimnz]);


 par_ID=netcdf.defVar(ncid,'PAR','double',[dimnt,dimnx,dimny,dimnz]);
 theta_ID=netcdf.defVar(ncid,'Temp','double',[dimnt,dimnx,dimny,dimnz]);
 salt_ID=netcdf.defVar(ncid,'Salt','double',[dimnt,dimnx,dimny,dimnz]);
 uvel_ID=netcdf.defVar(ncid,'Uvel','double',[dimnt,dimnx,dimny,dimnz]);
 vvel_ID=netcdf.defVar(ncid,'Vvel','double',[dimnt,dimnx,dimny,dimnz]);
 mld_ID=netcdf.defVar(ncid,'MLD','double',[dimnt,dimnx,dimny]);

 netcdf.endDef(ncid);

 netcdf.putVar(ncid,long_ID,X);
 netcdf.putVar(ncid,lat_ID,Y);
 netcdf.putVar(ncid,month_ID,month);
 netcdf.putVar(ncid,depth_ID,squeeze(Z(1:kn2)));
 netcdf.putVar(ncid,layer_ID,squeeze(drF(1:kn2)));


  netcdf.putVar(ncid,par_ID,spar);

  netcdf.putVar(ncid,theta_ID,stheta);
  netcdf.putVar(ncid,salt_ID,ssalt);
  netcdf.putVar(ncid,uvel_ID,suvel);
  netcdf.putVar(ncid,vvel_ID,svvel);
  netcdf.putVar(ncid,mld_ID,smld);

 netcdf.close(ncid);
