/[MITgcm]/MITgcm_contrib/dyncore_ASP/mat_scripts/modif_grid.m
ViewVC logotype

Annotation of /MITgcm_contrib/dyncore_ASP/mat_scripts/modif_grid.m

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


Revision 1.1 - (hide annotations) (download)
Thu Jul 10 22:08:00 2008 UTC (15 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
add matlab scripts to generate input files & to check output

1 jmc 1.1
2     kwr=1;
3     nc=32; np=nc+1; n6x=6*nc; nPg=n6x*nc;
4     addAngle=1;
5    
6     %outpName='tile.newgrid';
7     outpName='dxC3_dXYa';
8     inpName='tile.mitgrid';
9     rDir='/home/jmc/exp/inp_CS_grid/';
10    
11     %- load angle Cos & Sin:
12     if addAngle == 1,
13     %angFil=['proj_cs',int2str(nc),'_2uEvN.bin'];
14     %angFil=['proj_cs',int2str(nc),'_2uE.bin'];
15     %anCs=rdda(angFil,[nc 6 nc],1,'real*8','b');
16     %anSn=rdda(angFil,[nc 6 nc],2,'real*8','b');
17     namR=sprintf('stdAngleCS_cs%i.bin',nc);
18     fprintf(' read files: %s ',namR);
19     anCs=rdda(namR,[nc 6 nc],1,'real*8','b');
20     namR=sprintf('stdAngleSN_cs%i.bin',nc);
21     fprintf(' , %s\n',namR);
22     anSn=rdda(namR,[nc 6 nc],1,'real*8','b');
23     %--
24     anCs=permute(anCs,[1 3 2]); anSn=permute(anSn,[1 3 2]);
25     end
26    
27     for n=1:6,
28     %-- read :
29     namF=sprintf([rDir,inpName(1:4),'%3.3i',inpName(5:end)],n);
30     fid=fopen(namF,'r','b');
31     vv1=fread(fid,'real*8');
32     fclose(fid);
33     s=size(vv1,1); k1=s/np/np;
34     fprintf(['read: ',namF,' : size: %i (%ix%ix%i)\n'],s,np,np,k1);
35     vv1=reshape(vv1,[np np k1]);
36     if addAngle == 1,
37     %-- Add angle:
38     k2=k1+2;
39     vv2=zeros(np,np,k2); vv2(:,:,[1:k1])=vv1;
40     vv2([1:nc],[1:nc],k1+1)=anCs(:,:,n);
41     vv2([1:nc],[1:nc],k1+2)=anSn(:,:,n);
42     else
43     %-- do not add angle:
44     vv2=vv1; k2=k1;
45     end
46     %- XC YC DXF DYF RA XG YG DXV DYU RAZ DXC DYC RAW RAS DXG DYG:
47     %------------------------------------------------------------------
48     %- change 1 dxC (& 1 dyC) just next to the corner:
49     % load dxC,dyC:
50     dxC=vv1(:,:,11); dyC=vv1(:,:,12);
51     %----
52     %facM=19.;
53     facM=[24.0 -2.5 2.0]; % <- changing factor in o/oo
54     mnV=[dxC(1,1) dxC(2,1) dxC(1,2)];
55     for q=1:length(facM),
56     if facM(q) ~= 0,
57     nVal=mnV(q)*(1+facM(q)/1000);
58     [I J]=find(dxC==mnV(q));
59     dxC(find(dxC==mnV(q)))=nVal;
60     fprintf('--> change %i dxC by factor= %9.6f (o/oo) (= %6.1f m)\n', ...
61     length(I),facM(q),nVal-mnV(q));
62     [I J]=find(dyC==mnV(q));
63     dyC(find(dyC==mnV(q)))=nVal;
64     fprintf('--> change %i dyC by factor= %9.6f (o/oo) (= %6.1f m)\n', ...
65     length(I),facM(q),nVal-mnV(q));
66     end
67     end
68     %------------------------------------------------------------------
69     %- replace dxC,dyC :
70     vv2(:,:,11)=dxC; vv2(:,:,12)=dyC;
71     vv1(:,:,11)=dxC; vv1(:,:,12)=dyC;
72     %---
73     dxyW=vv1(:,:,11).*vv1(:,:,16); % <- dxC*dyG
74     %- replace rAw <- dxC*dyG
75     vv2(:,:,13)=dxyW;
76     %---
77     dxyS=vv1(:,:,12).*vv1(:,:,15); % =dyC*dxG
78     %- replace rAs <- dyC*dxG
79     vv2(:,:,14)=dxyS;
80     %--
81     aWSav=zeros(np,np);
82     aWSav(1:nc,:)=vv2(1:nc,:,13)+vv2(2:np,:,13);
83     aWSav(:,1:nc)=aWSav(:,1:nc)+vv2(:,1:nc,14)+vv2(:,2:np,14);
84     aWSav=aWSav/4;
85     %- replace rAc <- rAw_i + rAs_j / 2
86     vv2(:,:,5)=aWSav;
87     %---------------------------------
88     %-- write:
89     if kwr == 1,
90     if addAngle == 1,
91     namW=sprintf([outpName,'.face%3.3i.bin'],n);
92     else
93     p=find(outpName=='.');
94     namW=sprintf([outpName(1:p-1),'%3.3i',outpName(p:end)],n);
95     end
96     fid=fopen(namW,'w','b');
97     fwrite(fid,vv2,'real*8');
98     fclose(fid);
99     fprintf([' write to file: ',namW,' %i 2D.var(%ix%i)\n'],k2,np,np);
100     end
101     end

  ViewVC Help
Powered by ViewVC 1.1.22