drw_ncl.rar_fvcom_ncl

%%%计算纳潮量 clc clear casename='tfw'; n_element=5825; %三角形单元数 n_id =3126; %节点数 n_obc =14; %开边界点数 grd_file ='D:\陵水\土福湾\57_yuan.grd'; %grd 文件名 单位 米(转成平面坐标系)，未删除网格前的原始文件 t_beg = 118; t_end = 144; dir = 'D:\陵水\土福湾\sms'; depth=load('D:\陵水\土福湾\dep.dat'); %地形 avail_element=load('D:\陵水\土福湾\element.dat'); %%%三角网格面积计算 fid=fopen(grd_file,'r'); temp=fgetl(fid); temp=fgetl(fid); for i=1:n_id str2=fgetl(fid); A(i,:)=strread(str2,'%f','delimiter',' '); %节点坐标 end A(:,4)=[]; for i=1:n_element str2=fgetl(fid); B(i,:)=strread(str2,'%f','delimiter',' ');%三角形单元 end B(:,2)=[]; fclose(fid); temp_a=((A(B(:,2),2)-A(B(:,3),2)).^2+(A(B(:,2),3)-A(B(:,3),3)).^2).^(0.5); %a边 temp_b=((A(B(:,2),2)-A(B(:,4),2)).^2+(A(B(:,2),3)-A(B(:,4),3)).^2).^(0.5); %b边 temp_c=((A(B(:,3),2)-A(B(:,4),2)).^2+(A(B(:,3),3)-A(B(:,4),3)).^2).^(0.5); %c边 temp_p=(temp_a+temp_b+temp_c)./2; temp_s=((temp_p-temp_a).*(temp_p-temp_b).*(temp_p-temp_c).*temp_p).^0.5; %面积 clear temp_a temp_b temp_c temp_p; sss=sum(temp_s(avail_element)); for i=t_beg:t_end file=[dir,'\',casename,'el_',sprintf('%04d',i),'.xy']; ele=load(file); %水位 ele(:,4:end)=[]; ele_temp=ele+depth; temp_d=(ele_temp(B(:,2),3)+ele_temp(B(:,3),3)+ele_temp(B(:,4),3))/3; %ele+dep 水位加深度=实际水深 temp_d(temp_d<=0.05)=0; temp_v=temp_s.*temp_d; %纳水体积 ncl(i)=sum(temp_v(avail_element)); disp(i); end fid=fopen('D:\陵水\土福湾\ncl_h.dat','w+'); for i=t_beg:t_end; fprintf(fid,'%20.4f',ncl(i)); fprintf(fid,'\n'); end fclose(fid);