ga_full_GArun.m

function ga_full_GArun(Nsupply_struct,name_curr,varargin)


%% GA_FULL_GA_RUN: for a given forcing (mostly nitrate supply and forcing), runs the plankton model alongside Lagrangian trajectories
% then concatene the result into monthly maps.
% 
% Use:
% ga_full_GArun(Nsupply_struct,name_curr,varargin)
%
% Required inputs:
% 	Nsupply_struct 		structure containing .lat, .time, .Nsupply (in mmolC/m3/d), .name (used to save outputs)
% 	name_curr 			prefix of the current dataset to be used (netcdf files as generated by utils/ga_write_ariane_currents into dir_ariane_global/currents_data)
%
% Optional inputs: 
% 	'options_plankton_model' 	options for the plankton model, default {} (e.g. for krill parameterization: {'gmax_big',0.6*0.6,'eZ',0.1*0.6,'mZ',0.05*16/106*0.6})
%
% For instance, to use a krill parameterization:
%	ga_full_GArun(Nsupply_struct,name_curr,'gmax_big',0.6*0.6,'eZ',0.1*0.6,'mZ',0.05*16/106*0.6)
%
% Outputs: (saved in outputs/)
% 	daily trajectories computed by ga_growthadvection (saved yearly)
%	monthly maps (one file)
%
% Monique Messié, 2021 for public version
% Reference: Messié, M., D. A. Sancho-Gallegos, J. Fiechter, J. A. Santora, and F. P. Chavez (2022). 
%			Satellite-based Lagrangian model reveals how upwelling and oceanic circulation shape krill hotspots in the California Current System.
%			Frontiers in Marine Science, in press, https://doi.org/10.3389/fmars.2022.835813


		
global dir_ariane_global dir_output_global
dir_curr=[dir_ariane_global,'currents_data/'];
if nargin<2, error('Need to provide both Nsupply_struct and name_curr'), end
arg=ga_read_varargin(varargin,{'options_plankton_model',{}});
GArun_name=[Nsupply_struct.name,'_',name_curr];

		
% Get available time for currents		
curr_list=ga_dir2filenames(dir_curr,[name_curr,'_0*']);
time_curr=nan(length(curr_list),1);
for itime=1:length(curr_list)
	time_curr(itime)=ncread([dir_curr,curr_list{itime}],'time'); 
end
	
% Set Nsupply and currents on the same time grid. Here, regrids monthly Nsupply onto daily time steps.
Nsupply_struct_daily=struct();
Nsupply_struct_daily.lat=Nsupply_struct.lat;
Nsupply_struct_daily.time=time_curr;
Nsupply_struct_daily.Nsupply=nan(length(Nsupply_struct_daily.lat),length(Nsupply_struct_daily.time));
for ilat=1:length(Nsupply_struct_daily.lat)
	Nsupply_struct_daily.Nsupply(ilat,:)=interp1(Nsupply_struct.time,Nsupply_struct.Nsupply(ilat,:),Nsupply_struct_daily.time);
end
		
% In case Nsupply was monthly (otherwise comment out): fill the beginning of the first month and the end of the last month with constant Nsupply
[y,m,~]=datevec(Nsupply_struct_daily.time); 
hasdata=max(~isnan(Nsupply_struct_daily.Nsupply),[],1);
ifirst=find(hasdata,1,'first'); ilast=find(hasdata,1,'last');
ifirstmonth=find(y==y(ifirst) & m==m(ifirst),1,'first');	% first day of the first month with data
ilastmonth=find(y==y(ilast) & m==m(ilast),1,'last');		% last day of the last month with data
Nsupply_struct_daily.Nsupply(:,ifirstmonth:ifirst)=repmat(Nsupply_struct_daily.Nsupply(:,ifirst),1,ifirst-ifirstmonth+1);
Nsupply_struct_daily.Nsupply(:,ilast:ilastmonth)=repmat(Nsupply_struct_daily.Nsupply(:,ilast),1,ilastmonth-ilast+1);

% Retains only days for which we have Nsupply data (currents are already taken care of by regridding Nsupply on current time steps)
hasdata=max(~isnan(Nsupply_struct_daily.Nsupply),[],1);		 % recompute now that beginning/end months are filled
Nsupply_struct_daily.time=Nsupply_struct_daily.time(hasdata);
Nsupply_struct_daily.Nsupply=Nsupply_struct_daily.Nsupply(:,hasdata);
[y,~,~]=datevec(Nsupply_struct_daily.time); liste_years=unique(y)';
		
% Load coastline
load('inputs/coastline_California.mat','coast_x','coast_y')

%% -------- loop on years
for year=liste_years

	% sets inputs and outputs for that year
	time_all=Nsupply_struct_daily.time(Nsupply_struct_daily.time>=datenum(year,1,1) & Nsupply_struct_daily.time<=datenum(year,12,31));
	zoo_all=cell(length(time_all),1); 
		

	%% -------- loop on days
	for itime=1:length(time_all)
		time=time_all(itime); disp(datestr(time))

		% inputs Nsupply
		itime_Nsupply = Nsupply_struct_daily.time==time;
		inputs=struct();
		inputs.Nsupply=Nsupply_struct_daily.Nsupply(:,itime_Nsupply); 
		inputs.time=time;
		inputs.lat=Nsupply_struct_daily.lat; inputs.lon=inputs.lat*NaN; 
		igood_lat=~isnan(inputs.Nsupply); 
		% original positions
		for ilat=find(igood_lat)'
			icoast=ga_find_index(coast_y,inputs.lat(ilat)); 
			inputs.lon(ilat)=min(interp1(1:length(coast_x),coast_x,icoast));
		end

		% calcul growthadvection
		zoo_all{itime}=ga_growthadvection(inputs,name_curr,time,'options_plankton_model',arg.options_plankton_model);

	end

	% saving runs per year
	save([dir_output_global,'zoo_Lagrangian_',GArun_name,'_',num2str(year),'.mat'],'zoo_all','time_all','-v7.3')

end


%% -------- concatenation into monthly maps

output_grid=struct('lon',-135:0.125:-110,'lat',28:0.125:48);		% California grid at 1/8° resolution
ga_concatenation(GArun_name,output_grid,{'Z_big'});


return