domaine_global_def; domaine;
mask=squeeze(tmask3D(:,:,1)); recl=jpi*jpj*4;

choice_plot=0;
rep_out='./';

%choiceBulkParam: 
%1->LY 2->COARE 3->LY/NCEP 4-> LY/ECMWF
%choiceDataSet:
%1->NCEP 2->ECMWF 3-> CORE (LY)

list_kcur=[1:5];
for kcur=list_kcur
switch kcur
%main set
case 1
choiceDataSet=1; choiceBulkParam=1; ycur0=1992; ycurMax=2006;
case 2
choiceDataSet=2; choiceBulkParam=1; ycur0=1992; ycurMax=2006;
case 3
choiceDataSet=3; choiceBulkParam=1; ycur0=1992; ycurMax=2006;
case 4
choiceDataSet=1; choiceBulkParam=2; ycur0=1992; ycurMax=2006;
case 5
choiceDataSet=1; choiceBulkParam=3; ycur0=1992; ycurMax=2006;
end


if choiceDataSet==1; data_cur='.ncep.'; data_heights=[10 2 2]; albedo=0.14;
elseif choiceDataSet==2; data_cur='.ecmwf.'; data_heights=[10 2 2]; albedo=0.05;
else data_cur='.core.'; data_heights=[10 10 10]; albedo=0.1;
end

if choiceBulkParam==1
suff_bulk='LY.'; fprintf([data_cur ' LY standard \n']);
elseif choiceBulkParam==2
suff_bulk='coare.'; fprintf([data_cur ' COARE standard \n']);
elseif choiceBulkParam==3
dragCur=1.4/1000; suff_bulk='LY4ncep.'; fprintf([data_cur ' LY with NCEP CST DRAG coeff \n']);
elseif choiceBulkParam==4
dragCur=1.4/1000; suff_bulk='LY4ecmwf.'; fprintf([data_cur ' LY with ECMWF CST DRAG coeff \n']);
end

cen2kel=273.15;
rhoConstFresh=999.8;
stefanBoltzmann = 5.670e-8;
ocean_emissivity=5.50e-8 / 5.670e-8;

fid_runoff=fopen('/net/ross/raid0/gforget/DATAbin/forcing/CORE/runoffannual1x1','r','b');
runoff=fread(fid_runoff,[jpi jpj],'float32').*mask;
fclose(fid_runoff);

%SST and ICE data:
fid_sst=fopen('/net/ross/raid0/gforget/DATAbin/ICECONC/HADLEY/HadISST1_SST_9206_monthly','r','b');
field_sst=fread(fid_sst,jpi*jpj*12,'float32'); fclose(fid_sst); field_sst=reshape(field_sst,jpi,jpj,12);
fid_ice=fopen('/net/ross/raid0/gforget/DATAbin/ICECONC/HADLEY/HadISST1_ICE_9206_monthly','r','b');
field_ice=fread(fid_ice,jpi*jpj*12,'float32'); fclose(fid_ice); field_ice=reshape(field_ice,jpi,jpj,12);

%for interannual sst:
rep_sst='/net/ross/raid2/king/data_1x1_92-03/obs/'; pref_sst='SST_monthly_r2_'; suff_sst='';



for ycur=ycur0:ycurMax

suff_out=[data_cur suff_bulk num2str(ycur) '.daily.bin']; list_tcur=[1:365];
doInitFiles=1; doWriteMean=0;

for tcur=list_tcur

for ttcur=1:4

hcur=(tcur-1)*4+ttcur;

%part a: load fields
if choiceDataSet==1; [dataFlds]=ncep_load_fields(ycur,hcur);
elseif choiceDataSet==2; [dataFlds]=ecmwf_load_fields(ycur,hcur);
else; [dataFlds]=core_load_fields(ycur,hcur);
end;

%part b: time interpolate sst and ice fields from monthly atlas
[mcur,mcurW]=coeff_MonthlyAtlasInterp(tcur);
sst=(1-mcurW)*field_sst(:,:,mcur(1))+mcurW*field_sst(:,:,mcur(2));
ice=(1-mcurW)*field_ice(:,:,mcur(1))+mcurW*field_ice(:,:,mcur(2));
sst=sst.*mask; ice=ice.*mask;

%for interannual sst:
sstAtlas=sst;
[sst]=sst_load_field(rep_sst,suff_sst,pref_sst,ycur,hcur);
sst(find(isnan(sst)))=sstAtlas(find(isnan(sst)));

%part c: complete data structure:
dataFlds=setfield(dataFlds,'sst',sst);
dataFlds=setfield(dataFlds,'ice',ice);
dataFlds=setfield(dataFlds,'runoff',runoff);


%part 1: bulk formulae 

if choiceBulkParam==1; [bulkFlds]=exf_bulk_largeyeager04(dataFlds,data_heights);
elseif choiceBulkParam==2; [bulkFlds]=gmaze_bulk_coare(dataFlds,data_heights,albedo);
else; [bulkFlds]=exf_bulk_largeyeager04(dataFlds,data_heights,dragCur);
end;

%[bulkFlds]=gfdl_largeyeager04(dataFlds);


%part 2: net fluxes 

%--   Compute net from downward and downward from net longwave and
%     shortwave radiation, if needed.
%     lwflux = Stefan-Boltzmann constant * emissivity * SST - lwdown
%     swflux = - ( 1 - albedo ) * swdown
lwflux=stefanBoltzmann*ocean_emissivity*((dataFlds.sst+cen2kel).^4)-dataFlds.lwdn;
swflux=-(1-albedo)*dataFlds.swdn;
%     Net surface heat flux.
%     qnet = ( - hs - hl + lwflux + swflux );
qnet = (1-dataFlds.ice).*( - bulkFlds.hs - bulkFlds.hl + lwflux + swflux );
%     Salt flux from Precipitation and Evaporation.
%     empmr = ( evap - precip - runoff );
evap=bulkFlds.evap*rhoConstFresh; %m/s->kg/m^2/s
empmr = (1-dataFlds.ice).*( evap - dataFlds.precip ) - dataFlds.runoff;
%	put in structure as well:
netFlds=struct('lwflux',lwflux,'swflux',swflux,'qnet',qnet,'empmr',empmr);

%part 3: time average
swflux = (1-dataFlds.ice).*( swflux ); turbflux = (1-dataFlds.ice).*( - bulkFlds.hs - bulkFlds.hl + lwflux );
writeFlds=struct('ustress',bulkFlds.ustress,'vstress',bulkFlds.vstress,'qnet',netFlds.qnet,'empmr',netFlds.empmr,'turbflux',turbflux,'swflux',swflux,'precip',dataFlds.precip);

averagesFields(doInitFiles,doWriteMean,rep_out,suff_out,writeFlds);
doInitFiles=0; 

%part 4: plot results
if choice_plot==1; fig0=(choiceDataSet-1)*3;
%sign conventions are not the same in exf and largyeager ... I use largeyeager 
%convention for plot, positive=into the ocean
figure(fig0+1);
subplot(2,2,1); pcolor(lon2D_t,lat2D_t,bulkFlds.ustress); shading flat; caxis([-1 1]*0.2); colorbar;
subplot(2,2,2); pcolor(lon2D_t,lat2D_t,bulkFlds.vstress); shading flat; caxis([-1 1]*0.1); colorbar;
subplot(2,2,3); pcolor(lon2D_t,lat2D_t,-qnet); shading flat; caxis([-1 1]*200); colorbar;
subplot(2,2,4); pcolor(lon2D_t,lat2D_t,-empmr); shading flat; caxis([-1 1]*1e-4); colorbar;
figure(fig0+2);
subplot(2,2,1); pcolor(lon2D_t,lat2D_t,-swflux); shading flat; caxis([-1 1]*300); colorbar;
subplot(2,2,2); pcolor(lon2D_t,lat2D_t,-lwflux); shading flat; caxis([-1 1]*100); colorbar;
subplot(2,2,3); pcolor(lon2D_t,lat2D_t,bulkFlds.hl); shading flat; caxis([-1 0.5]*200); colorbar;
subplot(2,2,4); pcolor(lon2D_t,lat2D_t,bulkFlds.hs); shading flat; caxis([-1 1]*50); colorbar;
figure(fig0+3);
subplot(3,1,1); pcolor(lon2D_t,lat2D_t,dataFlds.precip); shading flat; caxis([-1 1]*1e-4); colorbar;
subplot(3,1,2); pcolor(lon2D_t,lat2D_t,-evap); shading flat; caxis([-1 1]*1e-4); colorbar;
subplot(3,1,3); pcolor(lon2D_t,lat2D_t,dataFlds.runoff); shading flat; caxis([-1 1]*1e-4); colorbar;
return;
end%if choice_plot==1

end%for ttcur=1:4
doWriteMean=1;
averagesFields(doInitFiles,doWriteMean,rep_out,suff_out,writeFlds);
doWriteMean=0;
end%for tcur=1:365

doInitFiles=-1; doWriteMean=0; 
averagesFields(doInitFiles,doWriteMean,rep_out,suff_out,writeFlds);

end%for ycur=ycur0:ycurMax
end%for kcur=1:3