/[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.2 - (hide annotations) (download)
Fri Jan 23 01:35:19 2015 UTC (10 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65u, checkpoint65j, checkpoint67a, checkpoint67b, checkpoint65i, checkpoint67d, checkpoint65m, checkpoint65n, HEAD
Changes since 1.1: +75 -23 lines
update set-up: start from a 1.yr run pickup file (run with atmPhys_stepSST=F)
 using closer to present new SST file: SST_symEx3.bin ;
and then switch on atmPhys_stepSST=True and use Qflux file: Qflux_w90.bin

1 jmc 1.1
2 jmc 1.2 kwr=0;
3     %kwr=1;
4    
5     %gDir='../run_26l/';
6 jmc 1.1 %gDir='../cs_grid/';
7 jmc 1.2 gDir='../run/';
8 jmc 1.1 G=load_grid(gDir,0);
9    
10     nx=G.dims(1); ny=G.dims(2); nc=ny;
11     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;
14 jmc 1.2 y1d=[-89.5:1:90];
15 jmc 1.1 %figure(1);clf;
16     %imagesc(msk'); set(gca,'YDir','normal');
17     %colorbar
18    
19     %-- make SST field (cos shape, no noise) from grid-output file YC:
20 jmc 1.2 yy=yc*pi/90; sst1=9+19*cos(yy);
21     %sst=273.15+sst1;
22 jmc 1.1 %fname='SST_cos0.bin';
23 jmc 1.2 yy=y1d*pi/90; sst1_1d=9+19*cos(yy); %- just for doing a plot
24 jmc 1.1
25 jmc 1.2 %-- follow APE profile:
26 jmc 1.1 yy=yc*pi/180;
27 jmc 1.2 var=sin(1.5*yy); var=1.-var.*var;
28     var(find(abs(yc) > 60.))=0.; sst2=27*var;
29     %sst=273.15+sst2;
30     %fname='SST_APE_1.bin';
31     %- just for doing a plot
32     yy=y1d*pi/180;
33     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 jmc 1.1
50     figure(1);clf;
51 jmc 1.2 subplot(211)
52 jmc 1.1 var=sst;
53     grph_CS(var,xc,yc,xg,yg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
54     title('SST');
55    
56 jmc 1.2 subplot(212)
57     plot(y1d,sst1_1d,'g-'); hold on;
58     plot(y1d,sst2_1d,'b-');
59     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 jmc 1.1 %return
86    
87     %-- make initial pot-temp field by adding noise to T(iter=0) output file:
88    
89     rDir=gDir;
90     namf='T'; it=0;
91     tini=rdmds([rDir,namf],it);
92     nr=size(tini,3);
93    
94     var=rand([nx,ny]); var=var-mean(var(:));
95     yy=yc*pi/90;
96     var=var.*(2+cos(yy))/3;
97 jmc 1.2 figure(3);clf;
98 jmc 1.1 %var=sst0;
99     grph_CS(var,xc,yc,xg,yg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
100    
101     noise=1.e-3;
102     tini1=tini+noise*reshape(reshape(var,[nx*ny 1])*ones(1,nr),[nx ny nr]);
103     %size(tini1)
104    
105 jmc 1.2 if kwr > 0,
106     fname='ini_theta.bin';
107     fid=fopen(fname,'w','b'); fwrite(fid,tini1,'real*8'); fclose(fid);
108     fprintf(['write file: ',fname,'\n']);
109     end
110 jmc 1.1
111     %- spec-humid : put constant Rel-Humid in the lowest troposphere
112     relhum=0.8 ; pHum=800.e+2;
113     relhum=0.4 ; pHum=700.e+2;
114     khum=max(find(G.rC > pHum));
115    
116     %- taken from AIM -> qsat in g/kg, pIn = normalised Pressure
117     P0=1.e+5; pIn=G.rC/P0;
118     qsat=calc_Qsat(1,pIn,tini);
119     qsat=reshape(qsat,[nx ny nr]);
120 jmc 1.2 qini=qsat*1.e-3*relhum;
121 jmc 1.1 qini(:,:,khum+1:end)=0;
122    
123 jmc 1.2 figure(4);clf;
124 jmc 1.1 pax=G.rC/100; %- in mb
125     i1=1; j1=1;
126     var=squeeze(qini(i1,j1,:));
127     plot(var,pax,'k-'); hold on;
128     i1=nc/2; j1=nc/2;
129     var=squeeze(qini(i1,j1,:));
130     plot(var,pax,'r-');
131     i1=nc*2.5; j1=nc*0.5;
132     var=squeeze(qini(i1,j1,:));
133     plot(var,pax,'b-');
134     hold off
135     set(gca,'YDir','reverse');
136     grid
137     legend('mid','eq','pol');
138     title('Q-ini profile');
139    
140 jmc 1.2 if kwr > 0,
141     fname=['ini_specQ_',int2str(nr),'l.bin'];
142     fid=fopen(fname,'w','b'); fwrite(fid,qini,'real*8'); fclose(fid);
143     fprintf(['write file: ',fname,'\n']);
144     end
145 jmc 1.1
146     var=reshape(tini,[nx*ny nr]);
147     var=mean(var);
148     for n=1:ceil(nr/10)
149     is=1+(n-1)*10; ie=min(nr,n*10);
150     if n == 1, fprintf(' tRef='); else fprintf(' '); end
151     fprintf(' %5.1f,',round(10*var(is:ie))/10);
152     fprintf('\n')
153     end
154    
155     return

  ViewVC Help
Powered by ViewVC 1.1.22