Read Cell by Cell Modflow

readCVHMcbc.m

function CBC = readCVHMcbc(filename)

fid = fopen(filename,'r');
ind1 = 0;
ind2 = 1;
kper = 0;
while ~feof(fid)
    temp=[];
    KSTP = fread(fid,1,'uint');
    KPER = fread(fid,1,'uint');
    if KPER~=kper
        ind2=1;
        ind1=ind1+1;
        kper=KPER;
    end

    DESC = char(fread(fid,16,'char')');
    fprintf(['Reading ' strtrim(DESC) ' for Period ' num2str(KPER) '\n'])
    NCOL = fread(fid,1,'uint');
    NROW = fread(fid,1,'uint');
    nlay = fread(fid,1,'int');
    NLAY = abs(nlay);
    if nlay<0
        ITYPE = fread(fid,1,'uint');
        DELT = fread(fid,1,'float');
        PERTIM = fread(fid,1,'float');
        TOTIM = fread(fid,1,'float');
        if ITYPE == 0 || ITYPE == 1
            temp = nan(NROW,NCOL,NLAY);
            for k=1:NLAY
                for i=1:NROW
                    for j=1:NCOL
                        temp(i,j,k) = fread(fid,1,'float');
                    end
                end
            end

        elseif ITYPE == 2 || ITYPE == 5
            NVAL = 1;
            if ITYPE==5
                NVAL=fread(fid,1,'uint');
                if NVAL>1
                    warning('ITYPE == 5 and NVAL > 1 not implemented yet')
                    for ii = 1:NVAL-1
                        CTMP = char(fread(fid,16,'char')');
                    end
                end
            end

            NLIST=fread(fid,1,'uint32');
            if NLIST>0
                NRC=NROW*NCOL;
                temp=nan(NLIST,3+NVAL);
                for i=1:NLIST
                    ICELL = fread(fid,1,'uint');
                    tmp_icell(i,1) = ICELL;
                    [c, r, l] = ind2sub([NCOL NROW NLAY],ICELL);
                    lay = floor((ICELL-1)/NRC+1);
                    row = floor(((ICELL-(lay-1)*NRC)-1)/NCOL+1);
                    col = ICELL-(lay-1)*NRC-(row-1)*NCOL;
                    temp(i,1:3) = [lay row col];
                    for j = 1:NVAL
                        VAL = fread(fid,1,'float');
                        temp(i,3+j) = VAL;
                    end
                end
            end
        elseif ITYPE==3
            for i=1:NROW*NCOL
                temp(i,1)=fread(fid,1,'uint');
            end
            for i=1:NROW*NCOL
                temp(i,2)=fread(fid,1,'float');
            end
        elseif ITYPE == 4
            temp = nan(NCOL,NROW,1);
            for ii = 1:NROW*NCOL
                temp(ii) = fread(fid,1,'float');
            end
            temp = temp';
        end

    elseif nlay>0
        temp=nan(NROW,NCOL,NLAY);
        for k=1:NLAY
            for i=1:NROW
                for j=1:NCOL
                    temp(i,j,k)=fread(fid,1,'float');
                end
            end
        end
    end
    CBC{ind1,1}{ind2,1}=strtrim(DESC);
    CBC{ind1,1}{ind2,2}=temp;
    ind2=ind2+1;
end

fclose(fid);

end