femwind/coarsening_icl.m

   1 function [hzc,icl3]=coarsening_icl(dx,dy,dz,params)
   2 % [[hzc,icl3]=coarsening_icldx,dy,dz,params)
   3 % in:
   4 %   dx,dy       mesh spacings, scalars
   5 %   dz          vertical element sizes, vector
   6 %   params      structure
   7 % out:
   8 %   hzc         horizontal coarsening factors in directions 1 and 2
   9 %   icl3        coarse indices in direction 3
  10 %
  11     if ~isvector(dz),
  12           error('dz must be a vector')
  13     end
  14     dz = dz(:)';  % make sure dz is a row
  15     disp(['coarsening_icl: input dx,dy,dz=',num2str([dx,dy,dz])])
  16     dxy=min(dx,dy);  % horizontal step
  17     n3 = length(dz)+1;
  18     % decide on horizontal coarsening factor
  19     crit=(dz(1)/dxy)/params.a(3);
  20     if crit > params.minaspect
  21         hzc=[2,2]; % future proofing if they are different
  22     else
  23         hzc=[1,1];
  24     end
  25     hzcavg=sqrt(hzc(1)*hzc(2));
  26     fprintf('horizontal coarsening factor %g %g because weighted height=%g\n',...
  27         hzc, crit)
  28     icl3=zeros(1,n3); % allocate max
  29     lcl=1; % last coarse level
  30     icl3(1)=lcl;
  31     nc3=0;
  32     for i=1:n3
  33         newlcl=lcl+1; % next coarse level by 1
  34         if lcl+2 <= n3
  35             crit = ((dz(lcl)+dz(lcl+1))/(2*dxy*hzcavg/2))/params.a(3);
  36             if crit < params.maxaspect
  37                 newlcl=lcl+2; % next coarse level by 2
  38             end
  39         end
  40         lcl = newlcl;
  41         if lcl <= n3
  42             icl3(i+1)=lcl;
  43         else % at the top already
  44             nc3=i;
  45             icl3 = icl3(1:i);
  46             break
  47         end
  48     end
  49     if nc3==0
  50         error('number of coarse layers is 0')
  51     end
  52     disp(['vertical coarse layers ',num2str(icl3)])
  53     hg=[0,cumsum(dz)];
  54     hgc=hg(icl3);
  55     disp(['heights above terrain ',num2str(hg)])
  56     disp(['coarse heights above terrain ',num2str(hgc)])
  57 end