/[MITgcm]/MITgcm_contrib/verification_other/atm_gray/input/gendata.m
ViewVC logotype

Contents of /MITgcm_contrib/verification_other/atm_gray/input/gendata.m

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


Revision 1.1 - (show annotations) (download)
Sun Jul 7 23:58:23 2013 UTC (12 years ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64k, checkpoint64x, checkpoint64z, checkpoint64o, checkpoint64p, checkpoint64r, checkpoint64w, checkpoint64v, checkpoint64m, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65e, checkpoint65h, checkpoint65
a 26 levels version for CS-32 grid

1
2 gDir='../run_26l/';
3 %gDir='../cs_grid/';
4 G=load_grid(gDir,0);
5
6 nx=G.dims(1); ny=G.dims(2); nc=ny;
7 xc=G.xC; yc=G.yC; xg=G.xG; yg=G.yG;
8
9 ccB=[0 0]; shift=-1; cbV=1; AxBx=[-180 180 -90 90]; kEnv=0;
10 %figure(1);clf;
11 %imagesc(msk'); set(gca,'YDir','normal');
12 %colorbar
13
14 %-- make SST field (cos shape, no noise) from grid-output file YC:
15 %yy=yc*pi/90;
16 %sst=273.15+9+19*cos(yy);
17 %fname='SST_cos0.bin';
18
19 yy=yc*pi/180;
20 sst=sin(1.5*yy); sst=1.-sst.*sst;
21 sst(find(abs(yc) > 60.))=0.;
22 sst=273.15+27*sst;
23 fname='SST_APE_1.bin';
24
25 %fid=fopen(fname,'w','b'); fwrite(fid,sst,'real*8'); fclose(fid);
26 %fprintf(['write file: ',fname,'\n']);
27
28 figure(1);clf;
29 var=sst;
30 grph_CS(var,xc,yc,xg,yg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
31 title('SST');
32
33 %figure(2);clf;
34 %yy=[-90:1:90]*pi/180;
35 %var=sin(1.5*yy); var=1.-var.*var; var(find(abs(yy) > pi/3))=0.; var=27*var;
36 %plot(yy*180/pi,var,'r-'); axis([-90 90 -1 28]); grid
37 %return
38
39 %-- make initial pot-temp field by adding noise to T(iter=0) output file:
40
41 rDir=gDir;
42 namf='T'; it=0;
43 tini=rdmds([rDir,namf],it);
44 nr=size(tini,3);
45
46 var=rand([nx,ny]); var=var-mean(var(:));
47 yy=yc*pi/90;
48 var=var.*(2+cos(yy))/3;
49 figure(2);clf;
50 %var=sst0;
51 grph_CS(var,xc,yc,xg,yg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
52
53 noise=1.e-3;
54 tini1=tini+noise*reshape(reshape(var,[nx*ny 1])*ones(1,nr),[nx ny nr]);
55 %size(tini1)
56
57 %fname='ini_theta.bin';
58 %fid=fopen(fname,'w','b'); fwrite(fid,tini1,'real*8'); fclose(fid);
59 %fprintf(['write file: ',fname,'\n']);
60
61 %- spec-humid : put constant Rel-Humid in the lowest troposphere
62 relhum=0.8 ; pHum=800.e+2;
63 relhum=0.4 ; pHum=700.e+2;
64 khum=max(find(G.rC > pHum));
65
66 %- taken from AIM -> qsat in g/kg, pIn = normalised Pressure
67 P0=1.e+5; pIn=G.rC/P0;
68 qsat=calc_Qsat(1,pIn,tini);
69 qsat=reshape(qsat,[nx ny nr]);
70 qini=qsat*1.e-3*relhum;
71 qini(:,:,khum+1:end)=0;
72
73 figure(3);clf;
74 pax=G.rC/100; %- in mb
75 i1=1; j1=1;
76 var=squeeze(qini(i1,j1,:));
77 plot(var,pax,'k-'); hold on;
78 i1=nc/2; j1=nc/2;
79 var=squeeze(qini(i1,j1,:));
80 plot(var,pax,'r-');
81 i1=nc*2.5; j1=nc*0.5;
82 var=squeeze(qini(i1,j1,:));
83 plot(var,pax,'b-');
84 hold off
85 set(gca,'YDir','reverse');
86 grid
87 legend('mid','eq','pol');
88 title('Q-ini profile');
89
90 %fname=['ini_specQ_',int2str(nr),'l.bin'];
91 %fid=fopen(fname,'w','b'); fwrite(fid,qini,'real*8'); fclose(fid);
92 %fprintf(['write file: ',fname,'\n']);
93
94 var=reshape(tini,[nx*ny nr]);
95 var=mean(var);
96 for n=1:ceil(nr/10)
97 is=1+(n-1)*10; ie=min(nr,n*10);
98 if n == 1, fprintf(' tRef='); else fprintf(' '); end
99 fprintf(' %5.1f,',round(10*var(is:ie))/10);
100 fprintf('\n')
101 end
102
103 return

  ViewVC Help
Powered by ViewVC 1.1.22