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

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

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

revision 1.1 by jmc, Sun Jul 7 23:58:23 2013 UTC revision 1.2 by jmc, Fri Jan 23 01:35:19 2015 UTC
# Line 1  Line 1 
1    
2  gDir='../run_26l/';  kwr=0;
3    %kwr=1;
4    
5    %gDir='../run_26l/';
6  %gDir='../cs_grid/';  %gDir='../cs_grid/';
7    gDir='../run/';
8  G=load_grid(gDir,0);  G=load_grid(gDir,0);
9    
10  nx=G.dims(1); ny=G.dims(2); nc=ny;  nx=G.dims(1); ny=G.dims(2); nc=ny;
11  xc=G.xC; yc=G.yC; xg=G.xG; yg=G.yG;  xc=G.xC; yc=G.yC; xg=G.xG; yg=G.yG;
12    
13  ccB=[0 0]; shift=-1; cbV=1; AxBx=[-180 180 -90 90]; kEnv=0;  ccB=[0 0]; shift=-1; cbV=1; AxBx=[-180 180 -90 90]; kEnv=0;
14    y1d=[-89.5:1:90];
15  %figure(1);clf;  %figure(1);clf;
16  %imagesc(msk'); set(gca,'YDir','normal');  %imagesc(msk'); set(gca,'YDir','normal');
17  %colorbar  %colorbar
18    
19  %-- make SST field (cos shape, no noise) from grid-output file YC:  %-- make SST field (cos shape, no noise) from grid-output file YC:
20  %yy=yc*pi/90;  yy=yc*pi/90; sst1=9+19*cos(yy);
21  %sst=273.15+9+19*cos(yy);  %sst=273.15+sst1;
22  %fname='SST_cos0.bin';  %fname='SST_cos0.bin';
23    yy=y1d*pi/90; sst1_1d=9+19*cos(yy); %- just for doing a plot
24    
25    %-- follow APE profile:
26  yy=yc*pi/180;  yy=yc*pi/180;
27  sst=sin(1.5*yy); sst=1.-sst.*sst;  var=sin(1.5*yy); var=1.-var.*var;
28  sst(find(abs(yc) > 60.))=0.;  var(find(abs(yc) > 60.))=0.; sst2=27*var;
29  sst=273.15+27*sst;  %sst=273.15+sst2;
30  fname='SST_APE_1.bin';  %fname='SST_APE_1.bin';
31    %- just for doing a plot
32  %fid=fopen(fname,'w','b'); fwrite(fid,sst,'real*8'); fclose(fid);  yy=y1d*pi/180;
33  %fprintf(['write file: ',fname,'\n']);  var=sin(1.5*yy); var=1.-var.*var;
34    var(find(abs(y1d) > 60.))=0.; sst2_1d=27*var;
35    
36    %-- close to Symetric zonally-average current SST:
37    yyp=abs(yc/90); phi=yyp.^3.;
38    sst3=27.8*exp(-7*phi) - 1 ;
39    sst=273.15+sst3;
40    fname='SST_symEx3.bin';
41    %- just for doing a plot
42    yyp=abs(y1d/90); phi=yyp.^3.;
43    sst3_1d=27.8*exp(-7*phi) - 1 ;
44    
45    if kwr > 0,
46     fid=fopen(fname,'w','b'); fwrite(fid,sst,'real*8'); fclose(fid);
47     fprintf(['write file: ',fname,'\n']);
48    end
49    
50  figure(1);clf;  figure(1);clf;
51    subplot(211)
52  var=sst;  var=sst;
53  grph_CS(var,xc,yc,xg,yg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);  grph_CS(var,xc,yc,xg,yg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
54  title('SST');  title('SST');
55    
56  %figure(2);clf;  subplot(212)
57  %yy=[-90:1:90]*pi/180;  plot(y1d,sst1_1d,'g-'); hold on;
58  %var=sin(1.5*yy); var=1.-var.*var; var(find(abs(yy) > pi/3))=0.; var=27*var;  plot(y1d,sst2_1d,'b-');
59  %plot(yy*180/pi,var,'r-'); axis([-90 90 -1 28]); grid  plot(y1d,sst3_1d,'r-');
60    hold off;
61    axis([-90 90 -5 29]); grid
62    legend('cos','ape','sym')
63    title('SST ^oC');
64    %return
65    
66    %-- make Q-flux file (similar to the one used by Paul, hard coded
67    %                     in original version of mixed_layer.f90)
68    qflx_ampl=50. ;
69    qflx_width=90.;
70    yy=yc/qflx_width;
71    yy=yy.*yy;
72    qflx=qflx_ampl*(1-2*yy);
73    qflx=qflx.*exp(-yy);
74    
75    figure(2);clf;
76    var=qflx;
77    grph_CS(var,xc,yc,xg,yg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
78    title('Q-flux [W/m^2]');
79    
80    if kwr > 0,
81     fname='Qflux_w90.bin';
82     fid=fopen(fname,'w','b'); fwrite(fid,qflx,'real*8'); fclose(fid);
83     fprintf(['write file: ',fname,'\n']);
84    end
85  %return  %return
86    
87  %-- make initial pot-temp field by adding noise to T(iter=0) output file:  %-- make initial pot-temp field by adding noise to T(iter=0) output file:
# Line 46  nr=size(tini,3); Line 94  nr=size(tini,3);
94  var=rand([nx,ny]); var=var-mean(var(:));  var=rand([nx,ny]); var=var-mean(var(:));
95  yy=yc*pi/90;  yy=yc*pi/90;
96  var=var.*(2+cos(yy))/3;  var=var.*(2+cos(yy))/3;
97  figure(2);clf;  figure(3);clf;
98  %var=sst0;  %var=sst0;
99  grph_CS(var,xc,yc,xg,yg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);  grph_CS(var,xc,yc,xg,yg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
100    
# Line 54  noise=1.e-3; Line 102  noise=1.e-3;
102  tini1=tini+noise*reshape(reshape(var,[nx*ny 1])*ones(1,nr),[nx ny nr]);  tini1=tini+noise*reshape(reshape(var,[nx*ny 1])*ones(1,nr),[nx ny nr]);
103  %size(tini1)  %size(tini1)
104    
105  %fname='ini_theta.bin';  if kwr > 0,
106  %fid=fopen(fname,'w','b'); fwrite(fid,tini1,'real*8'); fclose(fid);   fname='ini_theta.bin';
107  %fprintf(['write file: ',fname,'\n']);   fid=fopen(fname,'w','b'); fwrite(fid,tini1,'real*8'); fclose(fid);
108     fprintf(['write file: ',fname,'\n']);
109    end
110    
111  %- spec-humid : put constant Rel-Humid in the lowest troposphere  %- spec-humid : put constant Rel-Humid in the lowest troposphere
112  relhum=0.8 ; pHum=800.e+2;  relhum=0.8 ; pHum=800.e+2;
# Line 67  khum=max(find(G.rC > pHum)); Line 117  khum=max(find(G.rC > pHum));
117  P0=1.e+5; pIn=G.rC/P0;  P0=1.e+5; pIn=G.rC/P0;
118  qsat=calc_Qsat(1,pIn,tini);  qsat=calc_Qsat(1,pIn,tini);
119  qsat=reshape(qsat,[nx ny nr]);  qsat=reshape(qsat,[nx ny nr]);
120  qini=qsat*1.e-3*relhum;  qini=qsat*1.e-3*relhum;
121  qini(:,:,khum+1:end)=0;  qini(:,:,khum+1:end)=0;
122    
123  figure(3);clf;  figure(4);clf;
124  pax=G.rC/100; %- in mb  pax=G.rC/100; %- in mb
125  i1=1; j1=1;  i1=1; j1=1;
126  var=squeeze(qini(i1,j1,:));  var=squeeze(qini(i1,j1,:));
# Line 87  grid Line 137  grid
137  legend('mid','eq','pol');  legend('mid','eq','pol');
138  title('Q-ini profile');  title('Q-ini profile');
139    
140  %fname=['ini_specQ_',int2str(nr),'l.bin'];  if kwr > 0,
141  %fid=fopen(fname,'w','b'); fwrite(fid,qini,'real*8'); fclose(fid);   fname=['ini_specQ_',int2str(nr),'l.bin'];
142  %fprintf(['write file: ',fname,'\n']);   fid=fopen(fname,'w','b'); fwrite(fid,qini,'real*8'); fclose(fid);
143     fprintf(['write file: ',fname,'\n']);
144    end
145    
146  var=reshape(tini,[nx*ny nr]);  var=reshape(tini,[nx*ny nr]);
147  var=mean(var);  var=mean(var);
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22