gael/bulkMatlab/loop_1x1_bulk.m


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


1	gforget	1.1
2			domaine_global_def; domaine;
3			mask=squeeze(tmask3D(:,:,1)); recl=jpijpj4;
4
5			choice_plot=0;
6			rep_out='./';
7
8			%choiceBulkParam:
9			%1->LY 2->COARE 3->LY/NCEP 4-> LY/ECMWF
10			%choiceDataSet:
11			%1->NCEP 2->ECMWF 3-> CORE (LY)
12
13			list_kcur=[1:5];
14			for kcur=list_kcur
15			switch kcur
16			%main set
17			case 1
18			choiceDataSet=1; choiceBulkParam=1; ycur0=1992; ycurMax=2006;
19			case 2
20			choiceDataSet=2; choiceBulkParam=1; ycur0=1992; ycurMax=2006;
21			case 3
22			choiceDataSet=3; choiceBulkParam=1; ycur0=1992; ycurMax=2006;
23			case 4
24			choiceDataSet=1; choiceBulkParam=2; ycur0=1992; ycurMax=2006;
25			case 5
26			choiceDataSet=1; choiceBulkParam=3; ycur0=1992; ycurMax=2006;
27			end
28
29
30			if choiceDataSet==1; data_cur='.ncep.'; data_heights=[10 2 2]; albedo=0.14;
31			elseif choiceDataSet==2; data_cur='.ecmwf.'; data_heights=[10 2 2]; albedo=0.05;
32			else data_cur='.core.'; data_heights=[10 10 10]; albedo=0.1;
33			end
34
35			if choiceBulkParam==1
36			suff_bulk='LY.'; fprintf([data_cur ' LY standard \n']);
37			elseif choiceBulkParam==2
38			suff_bulk='coare.'; fprintf([data_cur ' COARE standard \n']);
39			elseif choiceBulkParam==3
40			dragCur=1.4/1000; suff_bulk='LY4ncep.'; fprintf([data_cur ' LY with NCEP CST DRAG coeff \n']);
41			elseif choiceBulkParam==4
42			dragCur=1.4/1000; suff_bulk='LY4ecmwf.'; fprintf([data_cur ' LY with ECMWF CST DRAG coeff \n']);
43			end
44
45			cen2kel=273.15;
46			rhoConstFresh=999.8;
47			stefanBoltzmann = 5.670e-8;
48			ocean_emissivity=5.50e-8 / 5.670e-8;
49
50			fid_runoff=fopen('/net/ross/raid0/gforget/DATAbin/forcing/CORE/runoffannual1x1','r','b');
51			runoff=fread(fid_runoff,[jpi jpj],'float32').*mask;
52			fclose(fid_runoff);
53
54			%SST and ICE data:
55			fid_sst=fopen('/net/ross/raid0/gforget/DATAbin/ICECONC/HADLEY/HadISST1_SST_9206_monthly','r','b');
56			field_sst=fread(fid_sst,jpijpj12,'float32'); fclose(fid_sst); field_sst=reshape(field_sst,jpi,jpj,12);
57			fid_ice=fopen('/net/ross/raid0/gforget/DATAbin/ICECONC/HADLEY/HadISST1_ICE_9206_monthly','r','b');
58			field_ice=fread(fid_ice,jpijpj12,'float32'); fclose(fid_ice); field_ice=reshape(field_ice,jpi,jpj,12);
59
60			%for interannual sst:
61			rep_sst='/net/ross/raid2/king/data_1x1_92-03/obs/'; pref_sst='SST_monthly_r2_'; suff_sst='';
62
63
64
65			for ycur=ycur0:ycurMax
66
67			suff_out=[data_cur suff_bulk num2str(ycur) '.daily.bin']; list_tcur=[1:365];
68			doInitFiles=1; doWriteMean=0;
69
70			for tcur=list_tcur
71
72			for ttcur=1:4
73
74			hcur=(tcur-1)*4+ttcur;
75
76			%part a: load fields
77			if choiceDataSet==1; [dataFlds]=ncep_load_fields(ycur,hcur);
78			elseif choiceDataSet==2; [dataFlds]=ecmwf_load_fields(ycur,hcur);
79			else; [dataFlds]=core_load_fields(ycur,hcur);
80			end;
81
82			%part b: time interpolate sst and ice fields from monthly atlas
83			[mcur,mcurW]=coeff_MonthlyAtlasInterp(tcur);
84			sst=(1-mcurW)field_sst(:,:,mcur(1))+mcurWfield_sst(:,:,mcur(2));
85			ice=(1-mcurW)field_ice(:,:,mcur(1))+mcurWfield_ice(:,:,mcur(2));
86			sst=sst.mask; ice=ice.mask;
87
88			%for interannual sst:
89			sstAtlas=sst;
90			[sst]=sst_load_field(rep_sst,suff_sst,pref_sst,ycur,hcur);
91			sst(find(isnan(sst)))=sstAtlas(find(isnan(sst)));
92
93			%part c: complete data structure:
94			dataFlds=setfield(dataFlds,'sst',sst);
95			dataFlds=setfield(dataFlds,'ice',ice);
96			dataFlds=setfield(dataFlds,'runoff',runoff);
97
98
99			%part 1: bulk formulae
100
101			if choiceBulkParam==1; [bulkFlds]=exf_bulk_largeyeager04(dataFlds,data_heights);
102			elseif choiceBulkParam==2; [bulkFlds]=gmaze_bulk_coare(dataFlds,data_heights,albedo);
103			else; [bulkFlds]=exf_bulk_largeyeager04(dataFlds,data_heights,dragCur);
104			end;
105
106			%[bulkFlds]=gfdl_largeyeager04(dataFlds);
107
108
109			%part 2: net fluxes
110
111			%-- Compute net from downward and downward from net longwave and
112			% shortwave radiation, if needed.
113			% lwflux = Stefan-Boltzmann constant * emissivity * SST - lwdown
114			% swflux = - ( 1 - albedo ) * swdown
115			lwflux=stefanBoltzmannocean_emissivity((dataFlds.sst+cen2kel).^4)-dataFlds.lwdn;
116			swflux=-(1-albedo)*dataFlds.swdn;
117			% Net surface heat flux.
118			% qnet = ( - hs - hl + lwflux + swflux );
119			qnet = (1-dataFlds.ice).*( - bulkFlds.hs - bulkFlds.hl + lwflux + swflux );
120			% Salt flux from Precipitation and Evaporation.
121			% empmr = ( evap - precip - runoff );
122			evap=bulkFlds.evap*rhoConstFresh; %m/s->kg/m^2/s
123			empmr = (1-dataFlds.ice).*( evap - dataFlds.precip ) - dataFlds.runoff;
124			% put in structure as well:
125			netFlds=struct('lwflux',lwflux,'swflux',swflux,'qnet',qnet,'empmr',empmr);
126
127			%part 3: time average
128			swflux = (1-dataFlds.ice).( swflux ); turbflux = (1-dataFlds.ice).( - bulkFlds.hs - bulkFlds.hl + lwflux );
129			writeFlds=struct('ustress',bulkFlds.ustress,'vstress',bulkFlds.vstress,'qnet',netFlds.qnet,'empmr',netFlds.empmr,'turbflux',turbflux,'swflux',swflux,'precip',dataFlds.precip);
130
131			averagesFields(doInitFiles,doWriteMean,rep_out,suff_out,writeFlds);
132			doInitFiles=0;
133
134			%part 4: plot results
135			if choice_plot==1; fig0=(choiceDataSet-1)*3;
136			%sign conventions are not the same in exf and largyeager ... I use largeyeager
137			%convention for plot, positive=into the ocean
138			figure(fig0+1);
139			subplot(2,2,1); pcolor(lon2D_t,lat2D_t,bulkFlds.ustress); shading flat; caxis([-1 1]*0.2); colorbar;
140			subplot(2,2,2); pcolor(lon2D_t,lat2D_t,bulkFlds.vstress); shading flat; caxis([-1 1]*0.1); colorbar;
141			subplot(2,2,3); pcolor(lon2D_t,lat2D_t,-qnet); shading flat; caxis([-1 1]*200); colorbar;
142			subplot(2,2,4); pcolor(lon2D_t,lat2D_t,-empmr); shading flat; caxis([-1 1]*1e-4); colorbar;
143			figure(fig0+2);
144			subplot(2,2,1); pcolor(lon2D_t,lat2D_t,-swflux); shading flat; caxis([-1 1]*300); colorbar;
145			subplot(2,2,2); pcolor(lon2D_t,lat2D_t,-lwflux); shading flat; caxis([-1 1]*100); colorbar;
146			subplot(2,2,3); pcolor(lon2D_t,lat2D_t,bulkFlds.hl); shading flat; caxis([-1 0.5]*200); colorbar;
147			subplot(2,2,4); pcolor(lon2D_t,lat2D_t,bulkFlds.hs); shading flat; caxis([-1 1]*50); colorbar;
148			figure(fig0+3);
149			subplot(3,1,1); pcolor(lon2D_t,lat2D_t,dataFlds.precip); shading flat; caxis([-1 1]*1e-4); colorbar;
150			subplot(3,1,2); pcolor(lon2D_t,lat2D_t,-evap); shading flat; caxis([-1 1]*1e-4); colorbar;
151			subplot(3,1,3); pcolor(lon2D_t,lat2D_t,dataFlds.runoff); shading flat; caxis([-1 1]*1e-4); colorbar;
152			return;
153			end%if choice_plot==1
154
155			end%for ttcur=1:4
156			doWriteMean=1;
157			averagesFields(doInitFiles,doWriteMean,rep_out,suff_out,writeFlds);
158			doWriteMean=0;
159			end%for tcur=1:365
160
161			doInitFiles=-1; doWriteMean=0;
162			averagesFields(doInitFiles,doWriteMean,rep_out,suff_out,writeFlds);
163
164			end%for ycur=ycur0:ycurMax
165			end%for kcur=1:3
166
167