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

Annotation 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 - (hide 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 jmc 1.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