/[MITgcm]/MITgcm_contrib/high_res_cube/eddy_flux/calcUVBave.m
ViewVC logotype

Diff of /MITgcm_contrib/high_res_cube/eddy_flux/calcUVBave.m

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

revision 1.2 by edhill, Fri Aug 6 04:37:18 2004 UTC revision 1.5 by edhill, Fri Aug 13 01:45:49 2004 UTC
# Line 43  delR   = [ Line 43  delR   = [
43      341.50,364.50,387.50,410.50,433.50,456.50 ];      341.50,364.50,387.50,410.50,433.50,456.50 ];
44  R = cumsum(delR) - 0.5*delR;  R = cumsum(delR) - 0.5*delR;
45    
46    vnam = [ 'U_'; 'V_'; 'T_'; 'S_'; 'B_'; 'B1'; ...
47             'UT'; 'VT'; 'US'; 'VS'; 'UB'; 'VB' ];
48    
49  il = 1;  il = 1;
50  for il = 1:nlay  for il = 1:nlay
51        
# Line 51  for il = 1:nlay Line 54  for il = 1:nlay
54        
55    %  Zero the accumulators    %  Zero the accumulators
56    nm = 0;    nm = 0;
57    mb = zeros(ne,ne,6);    for io = 1:size(vnam,1)
58    mbm1 = zeros(ne,ne,6);      comm = sprintf('%s_acc = zeros(ne,ne,6);', vnam(io,:));
59    mbu = zeros(ne,ne,6);      eval(comm)
60    mbv = zeros(ne,ne,6);    end
61        
62    it = 1;    it = 1;
63    for it = 1:length(ilist)    for it = 1:length(ilist)
# Line 105  for il = 1:nlay Line 108  for il = 1:nlay
108        %   locPres(i,j) = -rhoConst*rC(k)*gravity*maskC(i,j,k,bi,bj)        %   locPres(i,j) = -rhoConst*rC(k)*gravity*maskC(i,j,k,bi,bj)
109        %   rhoConst = 9.997999999999999E+02  /* Ref density ( kg/m^3 ) */        %   rhoConst = 9.997999999999999E+02  /* Ref density ( kg/m^3 ) */
110        %   gravity  = 9.810000000000000E+00  /* Grav accel ( m/s^2 ) */        %   gravity  = 9.810000000000000E+00  /* Grav accel ( m/s^2 ) */
111        p = 9.998*10^2 * 9.81 * R(il) * 0.0001;        p = 0.09998 * 9.81 * R(il);
112        bcube(:,:,icf) = eh3densjmd95( squeeze(scube(:,:,icf)), ...        bcube(:,:,icf) = ...
113                                       squeeze(tcube(:,:,icf)), p );            eh3densjmd95( squeeze(scube(:,:,icf)), ...
114                            squeeze(tcube(:,:,icf)), p ) * -9.81 / 1000;
115        if il > 1        if il > 1
116          p = -9.998*10^2 * 9.81 * R(il-1) * 0.0001;          p = 0.09998 * 9.81 * R(il-1);
117          bcubem1(:,:,icf) = eh3densjmd95( scube(:,:,icf), ...          bcubem1(:,:,icf) = ...
118                                           tcube(:,:,icf), p );              eh3densjmd95( scube(:,:,icf), ...
119                              tcube(:,:,icf), p ) * -9.81 / 1000;
120        end        end
121      end      end
122      tmask   = double(tcube ~= 0.0);      tmask   = double(tcube ~= 0.0);
123        scube   = scube .* tmask;
124      bcube   = bcube .* tmask;      bcube   = bcube .* tmask;
125      bcubem1 = bcubem1 .* tmask;      bcubem1 = bcubem1 .* tmask;
126    
127      %  Fill the (ne+1)*(ne+1) grid with B,T values      %  Fill the (ne+1)*(ne+1) grid with T,S,B values
128      ni = ne;    nip1 = ni + 1;      ni = ne;    nip1 = ni + 1;
129      nj = ne;    njp1 = nj + 1;      nj = ne;    njp1 = nj + 1;
130        tcubep1 = zeros( [ nip1 njp1 6 ] );
131        scubep1 = zeros( [ nip1 njp1 6 ] );
132      bcubep1 = zeros( [ nip1 njp1 6 ] );      bcubep1 = zeros( [ nip1 njp1 6 ] );
133      for icf = 1:6      for icf = 1:6
134          tcubep1(2:nip1,2:njp1,icf) = tcube(:,:,icf);
135          scubep1(2:nip1,2:njp1,icf) = scube(:,:,icf);
136        bcubep1(2:nip1,2:njp1,icf) = bcube(:,:,icf);        bcubep1(2:nip1,2:njp1,icf) = bcube(:,:,icf);
137      end      end
138    
139        %  Do the upwind-edge T exchanges
140        tcubep1(1,2:njp1,1) = tcube(ni:-1:1,nj,5);   % -
141        tcubep1(2:nip1,1,1) = tcube(1:ni,nj,6);      % -
142        tcubep1(1,2:njp1,2) = tcube(ni,1:nj,1);      % -
143        tcubep1(2:nip1,1,2) = tcube(ni,nj:-1:1,6);   % -
144        tcubep1(1,2:njp1,3) = tcube(ni:-1:1,nj,1);   % -
145        tcubep1(2:nip1,1,3) = tcube(1:ni,nj,2);      % -
146        tcubep1(1,2:njp1,4) = tcube(ni,1:nj,3);      % -
147        tcubep1(2:nip1,1,4) = tcube(ni,nj:-1:1,2);   % -
148        tcubep1(1,2:njp1,5) = tcube(ni:-1:1,nj,3);   % -
149        tcubep1(2:nip1,1,5) = tcube(1:ni,nj,4);      % -
150        tcubep1(1,2:njp1,6) = tcube(ni,1:nj,5);      % -
151        tcubep1(2:nip1,1,6) = tcube(ni,nj:-1:1,4);   % -
152    
153        %  Do the upwind-edge S exchanges
154        scubep1(1,2:njp1,1) = scube(ni:-1:1,nj,5);   % -
155        scubep1(2:nip1,1,1) = scube(1:ni,nj,6);      % -
156        scubep1(1,2:njp1,2) = scube(ni,1:nj,1);      % -
157        scubep1(2:nip1,1,2) = scube(ni,nj:-1:1,6);   % -
158        scubep1(1,2:njp1,3) = scube(ni:-1:1,nj,1);   % -
159        scubep1(2:nip1,1,3) = scube(1:ni,nj,2);      % -
160        scubep1(1,2:njp1,4) = scube(ni,1:nj,3);      % -
161        scubep1(2:nip1,1,4) = scube(ni,nj:-1:1,2);   % -
162        scubep1(1,2:njp1,5) = scube(ni:-1:1,nj,3);   % -
163        scubep1(2:nip1,1,5) = scube(1:ni,nj,4);      % -
164        scubep1(1,2:njp1,6) = scube(ni,1:nj,5);      % -
165        scubep1(2:nip1,1,6) = scube(ni,nj:-1:1,4);   % -
166    
167      %  Do the upwind-edge B exchanges      %  Do the upwind-edge B exchanges
168      bcubep1(1,2:njp1,1) = bcube(ni:-1:1,nj,5);   % -      bcubep1(1,2:njp1,1) = bcube(ni:-1:1,nj,5);   % -
169      bcubep1(2:nip1,1,1) = bcube(1:ni,nj,6);      % -      bcubep1(2:nip1,1,1) = bcube(1:ni,nj,6);      % -
# Line 140  for il = 1:nlay Line 178  for il = 1:nlay
178      bcubep1(1,2:njp1,6) = bcube(ni,1:nj,5);      % -      bcubep1(1,2:njp1,6) = bcube(ni,1:nj,5);      % -
179      bcubep1(2:nip1,1,6) = bcube(ni,nj:-1:1,4);   % -      bcubep1(2:nip1,1,6) = bcube(ni,nj:-1:1,4);   % -
180    
181        %  Get T values on the U,V grid points
182        masku = 1.0 - abs(diff(double(tcubep1(2:nip1,:,:) == 0.0),1,2));
183        maskv = 1.0 - abs(diff(double(tcubep1(:,2:nip1,:) == 0.0),1,1));
184        diffu = 0.5*diff(tcubep1(2:nip1,:,:),1,2);
185        diffv = 0.5*diff(tcubep1(:,2:nip1,:),1,1);
186        tonu = tcube(:,:,:) + masku.*diffu;
187        tonv = tcube(:,:,:) + maskv.*diffv;
188    
189        %  Get S values on the U,V grid points
190        diffu = 0.5*diff(scubep1(2:nip1,:,:),1,2);
191        diffv = 0.5*diff(scubep1(:,2:nip1,:),1,1);
192        sonu = scube(:,:,:) + masku.*diffu;
193        sonv = scube(:,:,:) + maskv.*diffv;
194    
195      %  Get B values on the U grid points      %  Get B values on the U grid points
     masku = 1.0 - abs(diff(double(bcubep1(2:nip1,:,:) == 0.0),1,2));  
196      diffu = 0.5*diff(bcubep1(2:nip1,:,:),1,2);      diffu = 0.5*diff(bcubep1(2:nip1,:,:),1,2);
     bonu = bcube(:,:,:) + masku.*diffu;  
       
     %  Get B values on the V grid points  
     maskv = 1.0 - abs(diff(double(bcubep1(:,2:nip1,:) == 0.0),1,1));  
197      diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1);      diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1);
198        bonu = bcube(:,:,:) + masku.*diffu;
199      bonv = bcube(:,:,:) + maskv.*diffv;      bonv = bcube(:,:,:) + maskv.*diffv;
200            
201      if ((newyear ~= oldyear) && (nm > 0)) || (it == length(ilist))      if ((newyear ~= oldyear) && (nm > 0)) || (it == length(ilist))
202    
203        %  Write if its a new year or the last of all entries        %  Write if its a new year or the last of all entries
204        bout = sprintf('b-ave-%d', oldyear);        disp(['  ==>  Writing files for ==> ' sprintf('%d',oldyear)]);
205        mout = sprintf('bm1-ave-%d', oldyear);        for io = 1:size(vnam,1)
206        uout = sprintf('ub-ave-%d', oldyear);          comm = sprintf( ...
207        vout = sprintf('vb-ave-%d', oldyear);              '%s_ido = fopen(''%s_tave_%d'', ''a'', ''ieee-be'');', ...
208        disp(['  ==>  Writing files "' bout '" "' mout '" "' uout '" "' vout '"']);              vnam(io,:), vnam(io,:), oldyear );
209        bido = fopen(bout, 'a', 'ieee-be');          eval(comm)
210        mido = fopen(mout, 'a', 'ieee-be');          comm = sprintf( '%s_acc = %s_acc ./ nm;',vnam(io,:), ...
211        uido = fopen(uout, 'a', 'ieee-be');                          vnam(io,:));
212        vido = fopen(vout, 'a', 'ieee-be');          eval(comm)
213        mb = mb / nm;          comm = sprintf('fwrite(%s_ido, %s_acc, ''real*4'');', ...
214        mbm1 = mbm1 / nm;                         vnam(io,:), vnam(io,:) );
215        mbu = mbu / nm;          eval(comm)
216        mbv = mbv / nm;          comm = sprintf( 'fclose(%s_ido);', vnam(io,:) );
217        fwrite(bido, mb,   'real*4');          eval(comm)
218        fwrite(mido, mbm1, 'real*4');          
219        fwrite(uido, mbu,  'real*4');          %  Re-zero the accumulators
220        fwrite(vido, mbv,  'real*4');          comm = sprintf('%s_acc = zeros(ne,ne,6);', vnam(io,:));
221        fclose(bido);          eval(comm)
222        fclose(mido);        end
223        fclose(uido);  
       fclose(vido);  
         
       %  Re-zero the accumulators  
       mb = zeros(ne,ne,6);  
       mbm1 = zeros(ne,ne,6);  
       mbu = zeros(ne,ne,6);  
       mbv = zeros(ne,ne,6);  
224        nm = 0;        nm = 0;
225                
226      end      end
227            
228      %  Accumulate partial sums      %  Accumulate partial sums
229      mb = mb + bcube;      U__acc = U__acc + ucube;
230      mbm1 = mbm1 + bcubem1;      V__acc = V__acc + vcube;
231      mbu = mbu + bonu.*ucube;      T__acc = T__acc + tcube;
232      mbv = mbv + bonv.*vcube;      S__acc = S__acc + scube;
233        B__acc = B__acc + bcube;
234        B1_acc = B1_acc + bcubem1;
235        
236        UT_acc = UT_acc + tonu.*ucube;
237        VT_acc = VT_acc + tonv.*vcube;
238        US_acc = US_acc + sonu.*ucube;
239        VS_acc = VS_acc + sonv.*vcube;
240        UB_acc = UB_acc + bonu.*ucube;
241        VB_acc = VB_acc + bonv.*vcube;
242    
243      nm = nm + 1;      nm = nm + 1;
244            
245      oldyear = newyear;      oldyear = newyear;
246            
247    end    end
248  end  end
249    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22