%modifies the "reference" core normal year windfield
%Lovenduski et al. (2013) apply a 14% increase in windspeed in the southern ocean
%                         between 1958 to 2007, i.e. in 50 years
%                         this drives, in their model, an outgassing of
%                         0.007+-0.002 PgC/yr/m^2 (!units unclear), equivalent to a 0.368 Pg/C/yr
%                         in the souther ocean (the definition here is unclear - their pictures
%                         are consistent with assuming that the SO starts south of 50S.

%define input files
 infil_u = '/Users/hdietze/mom4p1_pubrel_18dec2009/exp/mom_sis_bling_template/INPUT/u_10_mod.clim.nc'
 infil_v = '/Users/hdietze/mom4p1_pubrel_18dec2009/exp/mom_sis_bling_template/INPUT/v_10_mod.clim.nc'

%some ini
 addpath('/usr/local/mbib/netcdf/matlab_netCDF_OPeNDAP/')
 !rm -f ../dat/*

%get grid
 x = ncread(infil_u,'LON1');
 y = ncread(infil_v,'LAT1');

%define region of increase
 southof = -44; %south of this Latitude wind stresss will be increased
 Iy      = find(y<southof);

%define amplitude of yearly increase
 incpye = 14./50 %[percent/year]

 
for i=1:50
 %create output file
  outfil_u = ['../dat/u_10_mod.clim.' num2str(i+1000) '.nc'];
  outfil_v = ['../dat/v_10_mod.clim.' num2str(i+1000) '.nc'];
  eval(['!ncks -v U_10_MOD ' infil_u ' ' outfil_u]);
  eval(['!ncks -v V_10_MOD ' infil_v ' ' outfil_v]);
 %get original wind field (read in everytime to check copy)
  u = ncread(outfil_u,'U_10_MOD');
  v = ncread(outfil_v,'V_10_MOD');
 %increase the wind speed
  u(:,Iy,:) = u(:,Iy,:).* (i*incpye+100)/100;
  v(:,Iy,:) = v(:,Iy,:).* (i*incpye+100)/100;
 %save increased wind speed
  ncwrite(outfil_u,'U_10_MOD',u);
  ncwrite(outfil_v,'V_10_MOD',v);
end

