%% reading wind speed data
clear all

CCdailywind=xlsread('dailyNQIwind summary for submission.xlsx','Sheet1','A1:C2000','basic');

%Take BB6 station for example to run the DO model;
BB6=xlsread('BB6 station information for submission.xlsx','Sheet1','A1:L165','basic');

% headers={'Date','Julian Day','Level', 'Temperature'	'Salinity',
%     'Chlorophyll a','Ammonia','N+N','ortho-P'}

DateExcel=BB6(:,1); %%Datenumber in excel
dayID=BB6(:,2);  %%% day ID from 01/01/2013;
level=BB6(:,3); %% surface or bottom 1 is surface, 0 is bottom
Temp=BB6(:,4); %%% temperature
Salinity=BB6(:,5);%%%salinity
Chlorophyll=BB6(:,6); %%%Chla
TDIN=BB6(:,7)+BB6(:,8);% total DIN, NH4+NO3+NO2
PO4=BB6(:,9);%% PO4
DOmea=BB6(:,10);%%in situ measurement for comparison;
data=[dayID,level,Temp,Salinity,Chlorophyll,TDIN,PO4,DOmea];
data_S=data(find(data(:,2)==1),:); %% surface data
data_B=data(find(data(:,2)==0),:); %% bottom data
data_ini=data(find(data(:,1)==min(dayID)),:); %%% the initial state
data_ini_S=data_ini(find(data_ini(:,2)==1),:);%%% inintial surface data
data_ini_B=data_ini(find(data_ini(:,2)==0),:);%%%in itial bottom water
DOpre(min(dayID))=data_ini_S(8);%%initial surface DO concentration
DOpre_bottom(min(dayID))=data_ini_B(8);%% initial bottom DO concentration

lastloc=1;
Sce_max=[0.75:0.25:1.25];%% the range of key coeff
Sce_ID=0;

for  i=1:3;
        for j=1:3;
            for k=1:3;
                for l=1:3
           Sce_ID=Sce_ID+1 %%model number
            for day=min(dayID):max(dayID);   
            Sce_R=Sce_max(i);
            Ver_Diff=Sce_max(j);
            Sce_Pmax=Sce_max(k);
            Sce_wind=Sce_max(l);
            loc=find(round(data_S(:,1))==day);
            if length(loc)==1  ;    % if there was measured data for this day, use it
                lastloc=loc;
                data_S_day=data_S(loc,:);
                data_B_day=data_B(loc,:);
            end
            if length(loc)==0   ;   % if there was not measured data for this day, find the next avaliable day, and use the average of these two days
                nextloc=lastloc+1;
                data_S_day=(data_S(nextloc,:)-data_S(lastloc,:))/(data_S(nextloc,1)-data_S(lastloc,1)) *(day-data_S(lastloc,1))+data_S(lastloc,:);
                data_B_day=(data_B(nextloc,:)-data_B(lastloc,:))/(data_B(nextloc,1)-data_B(lastloc,1)) *(day-data_B(lastloc,1))+data_B(lastloc,:);
            end
            
        data_S_day(1)=day ;   
        T(day)=data_S_day(3);
        Salinity(day)=data_S_day(4);

        T_bottom(day)=data_B_day(3);
        Salinity_bottom(day)=data_B_day(4);
   
        Chla(day)=data_S_day(5)/1000;
        Chla_bottom(day)=data_B_day(5)/1000;
        
        Chla_coeff=1;
        Chla(day)=Chla(day)*Chla_coeff  ;  
        Chla_bottom(day)=Chla_bottom(day)*Chla_coeff;
      
        DIN(day)=data_S_day(6);
        DIP(day)=data_S_day(7);
        
        N_coeff=1;
        DIN(day)=N_coeff*DIN(day);
        DIP(day)=N_coeff*DIP(day);
        
        DIN_bottom(day)=data_B_day(6);
        DIP_bottom(day)=data_B_day(7);
        
         N_coeff=1;
        DIN_bottom(day)=N_coeff*DIN_bottom(day);
        DIP_bottom(day)=N_coeff*DIP_bottom(day);
                
        if T(day)<=20;
           sita=1.13;
        else sita=1.047;
        end
          
           R_coeff(day)=1;
           R_coeff_bottom(day)=1;
           P_coeff(day)=1;
           P_coeff_bottom(day)=1;
          if Chla(day)>20/1000;
              P_coeff(day)=20/Chla(day)*20/Chla(day)/1000/1000;
          end
          if Chla_bottom(day)>20/1000;
              P_coeff_bottom(day)=20/Chla_bottom(day)*20/Chla_bottom(day)/1000/1000; 
          end         
          
        R(day)=Sce_R*R_coeff(day)*0.25*sita.^(T(day)-20)*Chla(day)*(50*32/12); %%% units mg/m3, the first coeff test from 0.1 to 1 to test
        
        if T_bottom(day)<=20;
           sita_bottom=1.13;
        else sita_bottom=1.047;
        end
        
        R_bottom(day)=Sce_R*R_coeff_bottom(day)*0.125*sita_bottom.^(T_bottom(day)-20)*Chla_bottom(day)*(50*32/12) ;%%% units mg/m3, the first coeff test from 0.1 to 1 to test
        S_SOD(day)=0 ; %%% sediment DO consumption is 0 for surface layer
        S_SOD_bottom(day)=exp(0.0606*(T_bottom(day)+273.15)-17.71);% unit mg/L/day, lab measuremnt; apply the temperature relation and rate
         
        %%% P was cited from Rucheng Tian 2014 and Stefan 1994
        Pmax(day)=Sce_Pmax*9.6*1.036^(T(day)-20); %units mg O2 (mg Chla)-1 h-l]
        I(day)=210*exp(-0.08-0.25*Chla(day)); %%%I0 from NCEP shortwave radiation Fennel, 2006
        MinL(day)=0.125*I(day)/(Pmax(day)^2+0.2*0.2*I(day).^2)^0.5 ;%it used to be 0.125*0.125
        N(day)=min((DIN(day)/(0.5+DIN(day))),DIP(day)/(DIP(day)+0.05));  %%%
        P(day)=P_coeff(day)*Pmax(day)*MinL(day)*Chla(day)*12*N(day);% units mg L day-1, 12 is the day time
        Pmax_bottom(day)=Sce_Pmax*9.6*1.036^(T_bottom(day)-20); %units mg O2 (mg Chla)-1 h-l]
        I_bottom(day)=25*exp(-0.08-0.25*Chla_bottom(day)); 
        MinL_bottom(day)=0.25*I_bottom(day)/(Pmax_bottom(day)^2+0.2*0.2*I_bottom(day).^2)^0.5;
        N_bottom(day)=min((DIN_bottom(day)/(0.5+DIN_bottom(day))),DIP_bottom(day)/(DIP_bottom(day)+0.05));  %%%
        P_bottom(day)=Pmax_bottom(day)*MinL_bottom(day)*Chla_bottom(day)*12*N_bottom(day);%for bottom layer
        
        CCdailywind(isnan(CCdailywind))=nanmean(CCdailywind(:,3));
        locwind=find(CCdailywind(:,2)==day); %%%wind speed match; 
        
        U10 (day)=Sce_wind*CCdailywind(locwind,3)/2 ;%unit m/s  find the wind speed in that day 
        Csat(day) = 14.6244-0.367134*T(day) + 0.0044972*T(day).^2 -0.0966*Salinity(day) + 0.00205*Salinity(day)*T(day) + 0.0002739*Salinity(day).^2;%% saturated DO concentration;
        dz=1;% assuing the surface layer depth is 1 m
        

    if U10(day)<11;
        Cd(day) = 0.0012;
    elseif U10(day)>=11 && U10(day)<=20;
            Cd(day) = (0.49+0.065.*U10(day)).*10^-3;
        elseif U10(day)>20;
            Cd(day) = 0.0018;
    end
        
Sc(day)=(1953.4-128.*T(day)+3.9918.*T(day).^2-0.050091.*T(day).^3).*(1+Salinity(day).*3.14.*10.^-3);% calculate Sc
%Cs=GGo2_units(t,S,'umol'); %calculate Cs, saturated [O2]
u_star_a(day)=U10(day).*sqrt(Cd(day));
u_star_w(day)=0.034.*u_star_a(day);
ks(day)=0.00013.*u_star_a(day).*((Sc(day)/660).^-0.5); %%A19
kc(day)=5.56.*(u_star_w(day).^3.86); %%A20
kp(day)=5.5.*(u_star_w(day).^2.76).*((Sc(day)/660).^-0.67);%%A21
d_P(day)=1.5244.*(u_star_w(day).^1.06);
kc_factor=1;
kp_factor=1;
ks_factor=1;
Fc(day)=-kc_factor.*kc(day).*0.20946*86400*32; %unit  mg L-1 day-1  
Fp(day)=-kp_factor.*kp(day).*((1+d_P(day)).*Csat(day)-DOpre(day)).*86400;%mg L-1 day-1  
Fs(day)=-ks(day).*ks_factor.*(DOpre(day)-Csat(day)).*86400;%convert umol/kg to mol/m3, final unit mol/m2/day
 
Fb(day)=(Fc(day)+Fp(day)); %%%total air-sea flux
     
         Kz(day)=Ver_Diff*0.00000132*5; %%% diffusion coeff between surface and bottom water
          
         %%% average value from Jianhong Xue's paper and Justin 1996
         %Kz is the vertical diffusion, the bigger Kz, the larger diffusion
         %coefficients, but weaker stratification, the smaller density
         %gredients,which means the lower surface temperature warming
         %rate.the smaller Kz, the smaller diffusion coefficents,and
         %stronger stratification, the bigger density gredients, which
         %means the higher surface warming rate.

        Diffusion(day)=Kz(day)*3600*24*(DOpre(day)-DOpre_bottom(day));% units mg/L/day  Kz is from Justic 1996        
        S(day)=(P(day)-R(day)+Fs(day)+Fb(day)-S_SOD(day)- Diffusion(day));%-surface layer      
        S_bottom(day)=(P_bottom(day)-R_bottom(day)-S_SOD_bottom(day)+ Diffusion(day));%- bottom layer
        DOpre(day+1)=DOpre(day)+S(day); %%the DO in next day is euqaion the present DO plus delta DO      
        DOpre_bottom(day+1)=DOpre_bottom(day)+S_bottom(day);
        
    if DOpre(day+1)<0;
        DOpre(day+1)=0;
    end
    if DOpre_bottom(day+1)<0;
        DOpre_bottom(day+1)=0;
    end
        end
       
        DOcell_surface{Sce_ID}=[DOpre];%% put the surface DO predication from difference models in cell format
        DOcell_bottom{Sce_ID}=[DOpre_bottom];%% put the surface DO predication from difference models in cell format
        %process_surface{Sce_ID}=[P;R;Fs;Fb;S_SOD;Diffusion];
        %process_bottom{Sce_ID}=[P_bottom;R_bottom;S_SOD_bottom;Diffusion];
        end
    end
        end
end




%%%%output data organization
 cell_DO_pre=[];
for m=1:i*j*k*l;
 %%%% summary all the predicted daily value, and do the moving average of predicted DO in both surface and bottom
 ID=min(dayID):max(dayID);
 cell_DO_pre{m}(:,1)=ID';
 cell_DO_pre{m}(:,2)=DOcell_surface{m}(min(dayID):max(dayID))';
 cell_DO_pre{m}(:,3)=DOcell_bottom{m}(min(dayID):max(dayID))';
 cell_DO_pre{m}(:,4)=tsmovavg(cell_DO_pre{m}(:,2),'s',30,1); %%moving average of surface
 cell_DO_pre{m}(:,5)=tsmovavg(cell_DO_pre{m}(:,3),'s',30,1); %%moving average of bottom
end
 DO_pre_summary=cell2mat(cell_DO_pre); %%change cell into matrix
 Average_pre_surface=mean(DO_pre_summary(:,4:5:i*j*k*l*5),2); %%%average surface prediction from 81 models;
 Std_pre_surface=std(DO_pre_summary(:,4:5:i*j*k*l*5),0,2);%%%DO std surface prediction from 81 models;
 Average_pre_bottom=mean(DO_pre_summary(:,5:5:i*j*k*l*5),2);%%%average bottom prediction from 81 models;
 Std_pre_bottom=std(DO_pre_summary(:,5:5:i*j*k*l*5),0,2);%%%std bottom prediction from 81 models;
 %%%%%%%%%%%%%%%%%%%%%%%%
 
 
%%%%%%%%%%output figure 
figure
A=shadedErrorBar(x2mdate(DO_pre_summary(:,1)+41275),Average_pre_surface,Std_pre_surface,'lineProps','r-','transparent',1,'patchSaturation',0.1);hold on;
B=shadedErrorBar(x2mdate(DO_pre_summary(:,1)+41275),Average_pre_bottom,Std_pre_bottom,'lineProps','b-','transparent',1,'patchSaturation',0.1)

C=scatter(x2mdate(data_S(:,1)+41275), data_S(:,8),20,'filled','MarkerfaceColor','r');hold on;
D=scatter(x2mdate(data_B(:,1)+41275), data_B(:,8),20,'filled','MarkerfaceColor','b');hold on;

box on;
legend([[A.patch,B.patch],C,D],{'Surf. Pre.','Bott. Pre.','Surf. Mea','Bott. Mea.'})
ylim([0,12]);
datetick('x','mmmyyyy');
set(gcf, 'PaperUnits', 'centimeters');
 x_width=50 ;y_width=25;
 set(gcf, 'PaperPosition', [0 0 x_width y_width]); %
 title('BB6 DO prediction','FontWeight','bold','FontSize',25);
xlabel ('Date','FontWeight','bold','LineWidth',1,'FontSize',25);
ylabel ('DO conc (mg/L)','FontWeight','bold','LineWidth',1,'FontSize',25);
set(gca,'FontSize',25,'FontWeight','bold','LineWidth',1);
saveas(gcf,'Example BB6 DO distribution with coeff change from 0_75 to 1_25.jpg');

%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%output excel

model_S_B=[x2mdate(DO_pre_summary(:,1)+41275),Average_pre_surface,Average_pre_bottom];
model_S_B(isnan(model_S_B(:,2)),:)=[];
mea_S_B=[x2mdate(data_S(:,1)+41275),data_S(:,5),data_B(:,5)];
header={'Date','Predicted_Surface_DO','Predicted_Bottom_DO'};
csvwrite_with_headers('Example Predicted DO for station BB6.csv',model_S_B,header); %%%for output file

