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 |