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

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

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

-revision 1.3 by edhill,
Thu Aug 12 15:50:38 2004 UTC
+revision 1.4 by edhill,
Thu Aug 12 18:12:02 2004 UTC
 Line 43 
 delR   = [
 .50,364.50,387.50,410.50,433.50,456.50 ];
  R = cumsum(delR) - 0.5*delR;
+ vnam = [ 'U_'; 'V_'; 'T_'; 'S_'; 'B_'; 'B1'; ...
+          'UT'; 'VT'; 'US'; 'VS'; 'UB'; 'VB' ];
  il = 1;
  for il = 1:nlay
-   a = sprintf('Slice: %g',il);
+   %  a = sprintf('Slice: %g',il);
    %  disp(a);
    %  Zero the accumulators
-   nm = 0;
+   for io = 1:size(vnam,1)
-   mb = zeros(ne,ne,6);
+     comm = sprintf('%s_acc = zeros(ne,ne,6);', vnam(io,:));
-   mbm1 = zeros(ne,ne,6);
+     eval(comm)
-   mbu = zeros(ne,ne,6);
+   end
-   mbv = zeros(ne,ne,6);
    it = 1;
    for it = 1:length(ilist)
-Line 105 
 for il = 1:nlay
+Line 107 
 for il = 1:nlay
        %   locPres(i,j) = -rhoConst*rC(k)*gravity*maskC(i,j,k,bi,bj)
        %   rhoConst = 9.997999999999999E+02  /* Ref density ( kg/m^3 ) */
        %   gravity  = 9.810000000000000E+00  /* Grav accel ( m/s^2 ) */
-       p = 9.998*10^2 * 9.81 * R(il) * 0.0001;
+       p = 0.09998 * 9.81 * R(il);
-       bcube(:,:,icf) = eh3densjmd95( squeeze(scube(:,:,icf)), ...
+       bcube(:,:,icf) = ...
-                                      squeeze(tcube(:,:,icf)), p );
+           eh3densjmd95( squeeze(scube(:,:,icf)), ...
+                         squeeze(tcube(:,:,icf)), p ) * -9.81 / 1000;
        if il > 1
-         p = 9.998*10^2 * 9.81 * R(il-1) * 0.0001;
+         p = 0.09998 * 9.81 * R(il-1);
-         bcubem1(:,:,icf) = eh3densjmd95( scube(:,:,icf), ...
+         bcubem1(:,:,icf) = ...
-                                          tcube(:,:,icf), p );
+             eh3densjmd95( scube(:,:,icf), ...
+                           tcube(:,:,icf), p ) * -9.81 / 1000;
        end
      end
      tmask   = double(tcube ~= 0.0);
+     scube   = scube .* tmask;
      bcube   = bcube .* tmask;
      bcubem1 = bcubem1 .* tmask;
-     %  Fill the (ne+1)*(ne+1) grid with B,T values
+     %  Fill the (ne+1)*(ne+1) grid with T,S,B values
      ni = ne;    nip1 = ni + 1;
      nj = ne;    njp1 = nj + 1;
+     tcubep1 = zeros( [ nip1 njp1 6 ] );
+     scubep1 = zeros( [ nip1 njp1 6 ] );
      bcubep1 = zeros( [ nip1 njp1 6 ] );
      for icf = 1:6
+       tcubep1(2:nip1,2:njp1,icf) = tcube(:,:,icf);
+       scubep1(2:nip1,2:njp1,icf) = scube(:,:,icf);
        bcubep1(2:nip1,2:njp1,icf) = bcube(:,:,icf);
      end
+     %  Do the upwind-edge T exchanges
+     tcubep1(1,2:njp1,1) = tcube(ni:-1:1,nj,5);   % -
+     tcubep1(2:nip1,1,1) = tcube(1:ni,nj,6);      % -
+     tcubep1(1,2:njp1,2) = tcube(ni,1:nj,1);      % -
+     tcubep1(2:nip1,1,2) = tcube(ni,nj:-1:1,6);   % -
+     tcubep1(1,2:njp1,3) = tcube(ni:-1:1,nj,1);   % -
+     tcubep1(2:nip1,1,3) = tcube(1:ni,nj,2);      % -
+     tcubep1(1,2:njp1,4) = tcube(ni,1:nj,3);      % -
+     tcubep1(2:nip1,1,4) = tcube(ni,nj:-1:1,2);   % -
+     tcubep1(1,2:njp1,5) = tcube(ni:-1:1,nj,3);   % -
+     tcubep1(2:nip1,1,5) = tcube(1:ni,nj,4);      % -
+     tcubep1(1,2:njp1,6) = tcube(ni,1:nj,5);      % -
+     tcubep1(2:nip1,1,6) = tcube(ni,nj:-1:1,4);   % -
+     %  Do the upwind-edge S exchanges
+     scubep1(1,2:njp1,1) = scube(ni:-1:1,nj,5);   % -
+     scubep1(2:nip1,1,1) = scube(1:ni,nj,6);      % -
+     scubep1(1,2:njp1,2) = scube(ni,1:nj,1);      % -
+     scubep1(2:nip1,1,2) = scube(ni,nj:-1:1,6);   % -
+     scubep1(1,2:njp1,3) = scube(ni:-1:1,nj,1);   % -
+     scubep1(2:nip1,1,3) = scube(1:ni,nj,2);      % -
+     scubep1(1,2:njp1,4) = scube(ni,1:nj,3);      % -
+     scubep1(2:nip1,1,4) = scube(ni,nj:-1:1,2);   % -
+     scubep1(1,2:njp1,5) = scube(ni:-1:1,nj,3);   % -
+     scubep1(2:nip1,1,5) = scube(1:ni,nj,4);      % -
+     scubep1(1,2:njp1,6) = scube(ni,1:nj,5);      % -
+     scubep1(2:nip1,1,6) = scube(ni,nj:-1:1,4);   % -
      %  Do the upwind-edge B exchanges
      bcubep1(1,2:njp1,1) = bcube(ni:-1:1,nj,5);   % -
      bcubep1(2:nip1,1,1) = bcube(1:ni,nj,6);      % -
-Line 140 
 for il = 1:nlay
+Line 177 
 for il = 1:nlay
      bcubep1(1,2:njp1,6) = bcube(ni,1:nj,5);      % -
      bcubep1(2:nip1,1,6) = bcube(ni,nj:-1:1,4);   % -
+     %  Get T values on the U,V grid points
+     masku = 1.0 - abs(diff(double(tcubep1(2:nip1,:,:) == 0.0),1,2));
+     maskv = 1.0 - abs(diff(double(tcubep1(:,2:nip1,:) == 0.0),1,1));
+     diffu = 0.5*diff(tcubep1(2:nip1,:,:),1,2);
+     diffv = 0.5*diff(tcubep1(:,2:nip1,:),1,1);
+     tonu = tcube(:,:,:) + masku.*diffu;
+     tonv = tcube(:,:,:) + maskv.*diffv;
+     %  Get S values on the U,V grid points
+     diffu = 0.5*diff(scubep1(2:nip1,:,:),1,2);
+     diffv = 0.5*diff(scubep1(:,2:nip1,:),1,1);
+     sonu = scube(:,:,:) + masku.*diffu;
+     sonv = scube(:,:,:) + maskv.*diffv;
      %  Get B values on the U grid points
-     masku = 1.0 - abs(diff(double(bcubep1(2:nip1,:,:) == 0.0),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));
      diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1);
+     bonu = bcube(:,:,:) + masku.*diffu;
      bonv = bcube(:,:,:) + maskv.*diffv;
      if ((newyear ~= oldyear) && (nm > 0)) || (it == length(ilist))
        %  Write if its a new year or the last of all entries
-       bout = sprintf('b-ave-%d', oldyear);
+       disp(['  ==>  Writing files for ==> ' sprintf('%d',oldyear)]);
-       mout = sprintf('bm1-ave-%d', oldyear);
+       for io = 1:size(vnam,1)
-       uout = sprintf('ub-ave-%d', oldyear);
+         comm = sprintf( ...
-       vout = sprintf('vb-ave-%d', oldyear);
+             '%s_ido = fopen(''%s_tave_%d'', ''a'', ''ieee-be'');', ...
-       disp(['  ==>  Writing files "' bout '" "' mout '" "' uout '" "' vout '"']);
+             vnam(io,:), vnam(io,:), oldyear );
-       bido = fopen(bout, 'a', 'ieee-be');
+         eval(comm)
-       mido = fopen(mout, 'a', 'ieee-be');
+         comm = sprintf( '%s_acc = %s_acc ./ nm;',vnam(io,:), ...
-       uido = fopen(uout, 'a', 'ieee-be');
+                         vnam(io,:));
-       vido = fopen(vout, 'a', 'ieee-be');
+         eval(comm)
-       mb = mb / nm;
+         comm = sprintf('fwrite(%s_ido, %s_acc, ''real*4'');', ...
-       mbm1 = mbm1 / nm;
+                        vnam(io,:), vnam(io,:) );
-       mbu = mbu / nm;
+         eval(comm)
-       mbv = mbv / nm;
+         comm = sprintf( 'fclose(%s_ido);', vnam(io,:) );
-       fwrite(bido, mb,   'real*4');
+         eval(comm)
-       fwrite(mido, mbm1, 'real*4');
-       fwrite(uido, mbu,  'real*4');
+         %  Re-zero the accumulators
-       fwrite(vido, mbv,  'real*4');
+         comm = sprintf('%s_acc = zeros(ne,ne,6);', vnam(io,:));
-       fclose(bido);
+         eval(comm)
-       fclose(mido);
+       end
-       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);
        nm = 0;
      end
      %  Accumulate partial sums
-     mb = mb + bcube;
+     U__acc = U__acc + ucube;
-     mbm1 = mbm1 + bcubem1;
+     V__acc = V__acc + vcube;
-     mbu = mbu + bonu.*ucube;
+     T__acc = T__acc + tcube;
-     mbv = mbv + bonv.*vcube;
+     S__acc = S__acc + scube;
+     B__acc = B__acc + bcube;
+     B1_acc = B1_acc + bcubem1;
+     UT_acc = UT_acc + tonu.*ucube;
+     VT_acc = VT_acc + tonv.*vcube;
+     US_acc = US_acc + sonu.*ucube;
+     VS_acc = VS_acc + sonv.*vcube;
+     UB_acc = UB_acc + bonu.*ucube;
+     VB_acc = VB_acc + bonv.*vcube;
      nm = nm + 1;
      oldyear = newyear;
    end
  end

 Legend:



Removed from v.1.3
 


changed lines


 
Added in v.1.4
 Legend:



Removed from v.1.3
 


changed lines


 
Added in v.1.4
-Removed from v.1.3
+Added in v.1.4

	ViewVC Help
Powered by ViewVC 1.1.22