% bio_circadian_sol.m: Initialization of the biological model using 
% the correlation method.
clear all, close all
%-------------------------------------------------------------------------------
% Biological parameters
%-------------------------------------------------------------------------------

% Nutrient uptake

kno3 = 1.5;          % Half-saturation for NO3 uptake
knh4 = 1.5;          % Half-saturation for NH4 uptake
psi  = 3.5;          % NH4 inhibition parameter
 		  
% Conversion factors

C2chl = 80;        % mg C (mg Chl)-1
C2N   = 0.15;      % umol N (umol C)-1  
MC    = 12;        % g C (mol C)-1

% Factors to determine zoo

zoo_ini   = 1.4;
zoo_depth = 13;    % Index of oag level

% Factors to determine det

facdet  = 2;

%Fill values

fill_no3   = 0.01;
fill_nh4   = 0.01;
fill_pno3  = 0.01;
fill_pnh4  = 0.01;
fill_phyto = 0.01;
fill_zoo   = 0.01;
fill_det   = 0.01;
fill_chl   = 0.01;

%-------------------------------------------------------------------------------
% Read OA files for observable variables
%-------------------------------------------------------------------------------

OAGDIR = '/data2/AOSN2/pmoreno/ASCOT01/OAG/';

ctdfile = [OAGDIR 'INI06/oag_ini.nc']
chlfile = [OAGDIR 'BIO_INI47/oag_tchl_ascot01_june06-17.nc']
no3file = [OAGDIR 'BIO_INI53/oag_no2no3_ascot01.nc']
nh4file = [OAGDIR 'BIO_INI52/oag_nh4sio4_ascot01.nc']

[temp,chl,nh4,no3,lon,lat,z] = bio_read_oa4(ctdfile,chlfile,nh4file,no3file);

% Substitute negative (and very small) values 

ind = find(chl<fill_chl); chl(ind) = fill_chl; clear ind
ind = find(no3<fill_no3); no3(ind) = fill_no3; clear ind
ind = find(nh4<fill_nh4); nh4(ind) = fill_nh4; clear ind

%-----------------------------------------------------------------------------
% Calculation of non-observable biological variables
%-----------------------------------------------------------------------------

% Phytoplankton 

phyto = C2chl*C2N.*chl./MC;
ind = find(phyto<fill_phyto); phyto(ind) =fill_phyto;

% Zooplankton

zoo = fill_zoo.*ones(size(phyto));
zoo(:,:,1:zoo_depth) = zoo_ini;

figure, plot(squeeze(zoo(1,1,:)),z), axis ij, grid on
axis([0 1.1*zoo_ini 0 max(z)])

% Detritus

det = facdet.*phyto;

%-----------------------------------------------------------------------------
% Evaluate initial phytoplankton (NO3 and NH4 split)
%-----------------------------------------------------------------------------

Qnh4 = nh4./(knh4+nh4);
Qno3 = no3./(kno3+no3).*exp(-psi.*nh4);
Qtot = Qno3+Qnh4;

pno3 = Qno3.*phyto./Qtot;
pnh4 = Qnh4.*phyto./Qtot;

%-----------------------------------------------------------------------------
% Save the biological fields
%-----------------------------------------------------------------------------

% Substitute NaN and negative values with small value

no3  (find(isnan(no3)  |  no3<fill_no3  )) = fill_no3;
nh4  (find(isnan(nh4)  |  nh4<fill_nh4  )) = fill_nh4;
pno3 (find(isnan(pno3) | pno3<fill_pno3 )) = fill_pno3;
pnh4 (find(isnan(pnh4) | pnh4<fill_pnh4 )) = fill_pnh4;
phyto(find(isnan(phyto)|phyto<fill_phyto)) = fill_phyto;
zoo  (find(isnan(zoo)  |  zoo<fill_zoo  )) = fill_zoo;
det  (find(isnan(det)  |  det<fill_det  )) = fill_det;
chl  (find(isnan(chl)  |  chl<fill_chl  )) = fill_chl;

bio_save

disp('Done :)')
