/[MITgcm]/MITgcm_contrib/gael/profilesMatlabProcessing/profiles_misc/MITprof_resample.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/profilesMatlabProcessing/profiles_misc/MITprof_resample.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.4 - (hide annotations) (download)
Fri Oct 9 23:43:34 2015 UTC (9 years, 9 months ago) by gforget
Branch: MAIN
Changes since 1.3: +3 -3 lines
- bug fix (avoid variance loss).

1 gforget 1.1 function [profOut]=MITprof_resample(profIn,fldIn,filOut,method,varargin);
2     %object: resample a set of fields in file filFldIn with specified time
3     % line timeIn to the positions of profIn and add to file filOut
4     %inputs: profIn is a gcmfaces field (nan-masked; up to N3,N4 dimensions)
5     % fldIn is a description of the fields being resampled including
6     % the corresponding file name and additional specs :
7     % fldIn.name, fldIn.long_name, fldIn.units, fldIn.fil (file
8     % name) and fldIn.tim (time axis specification). Supported
9     % fldIn.tim spec: 'const' (for time invariant climatology),
10     % 'monclim' (for monthly climatology), 'monser' (for monthly
11     % time series)
12     % filOut is the output MITprof file name (if un-specified
13     % the resul may only be returned as a function argument)
14     % method may be 'polygons' (or 'TriScatteredInterp' ... via
15     % gcmfaces_interp_2d in a loop ... to be implemented later)
16     %outputs: profOut is the MITprof structure where the interpolated values
17     % were appended to profIn (if un-specified the result
18     % may only be returned to output file)
19     %
20     %filFldIn is assumed to be 3D and binary at this point
21    
22     gcmfaces_global;
23    
24     doOut=~isempty(who('filOut')); doOutInit=false;
25     if doOut; doOut=~isempty(filOut); end;
26     if doOut; doOutInit=isempty(dir(filOut)); end;
27    
28     if isempty(who('method')); method='polygons'; end;
29    
30     %1) deal with time line
31     if strcmp(fldIn.tim,'monclim');
32     tim_fld=[-0.5:12.5]; rec_fld=[12 1:12 1];
33     tmp1=datevec(profIn.prof_date);
34     tmp2=datenum([tmp1(:,1) ones(profIn.np,2) zeros(profIn.np,3)]);
35     tim_prof=(profIn.prof_date-tmp2);
36     tim_prof(tim_prof>365)=365;
37     tim_prof=tim_prof/365*12;%neglecting differences in months length
38 gforget 1.3 elseif strcmp(fldIn.tim,'const')|strcmp(fldIn.tim,'std');
39     tim_fld=[1 2]; rec_fld=[1 1];
40 gforget 1.1 tim_prof=1.5*ones(profIn.np,1);
41     else;
42     error('this case remains to be implemented');
43     end;
44    
45     lon=profIn.prof_lon; lat=profIn.prof_lat;
46     depIn=-mygrid.RC; depOut=profIn.prof_depth;
47     profOut=NaN*ones(profIn.np,profIn.nr);
48    
49     %2) loop over record pairs
50     for tt=1:length(rec_fld)-1;
51 gforget 1.2 fld0=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt));
52     fld1=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt+1));
53 gforget 1.1 ndim=length(size(fld0{1}));
54     fld=cat(ndim+1,fld0,fld1);
55     %
56     ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));
57     if ~isempty(ii);
58 gforget 1.3 arr=gcmfaces_interp(fld,lon(ii),lat(ii),method);
59 gforget 1.1 arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
60     end;
61     %
62     k0=floor(tim_prof(ii)); k1=k0+1;
63     a0=tim_prof(ii)-k0; a0=a0*ones(1,profIn.nr);
64     profOut(ii,:)=(1-a0).*arr2(:,:,1)+a0.*arr2(:,:,2);
65 gforget 1.4 if strcmp(fldIn.tim,'std');
66     profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:)));
67     end;
68 gforget 1.1 end;
69    
70     %3) deal with file output
71     if doOutInit;
72     %create a header only file to later append resampled fields
73     prof=profIn;
74     tmp1=fieldnames(prof);
75     nt=length(prof.prof_date);
76     nr=length(prof.prof_depth);
77     for ii=1:length(tmp1);
78     tmp2=prod(size(getfield(prof,tmp1{ii})));
79     if tmp2==nt*nr; prof=rmfield(prof,tmp1{ii}); end;
80     end;
81     MITprof_write(filOut,prof);
82     end;
83    
84     if doOut;
85     %add the array itelf
86     MITprof_addVar(filOut,fldIn.name,'double',{'iDEPTH','iPROF'},profOut);
87    
88     %add its attributes
89     nc=ncopen(filOut,'write');
90     ncaddAtt(nc,fldIn.name,'long_name',fldIn.long_name);
91     ncaddAtt(nc,fldIn.name,'units',fldIn.units);
92     ncaddAtt(nc,fldIn.name,'missing_value',fldIn.missing_value);
93     ncaddAtt(nc,fldIn.name,'_FillValue',fldIn.FillValue);
94     ncclose(nc);
95     end;
96    
97     %4) deal with argument output
98     if nargout>0; profOut=setfield(profIn,fldIn.name,profOut); end;
99    
100    

  ViewVC Help
Powered by ViewVC 1.1.22