Warning, /verification/atm_gray/input/gendata.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 42875574 on 2024-11-26 20:20:16 UTC
4287557481 Jean*0001
0002 kwr=0;
0003 %kwr=1;
0004
0005 %gDir='../run_26l/';
0006 %gDir='../cs_grid/';
0007 gDir='../run/';
0008 G=load_grid(gDir,0);
0009
0010 nx=G.dims(1); ny=G.dims(2); nc=ny;
0011 xc=G.xC; yc=G.yC; xg=G.xG; yg=G.yG;
0012
0013 ccB=[0 0]; shift=-1; cbV=1; AxBx=[-180 180 -90 90]; kEnv=0;
0014 y1d=[-89.5:1:90];
0015 %figure(1);clf;
0016 %imagesc(msk'); set(gca,'YDir','normal');
0017 %colorbar
0018
0019 %-- make SST field (cos shape, no noise) from grid-output file YC:
0020 yy=yc*pi/90; sst1=9+19*cos(yy);
0021 %sst=273.15+sst1;
0022 %fname='SST_cos0.bin';
0023 yy=y1d*pi/90; sst1_1d=9+19*cos(yy); %- just for doing a plot
0024
0025 %-- follow APE profile:
0026 yy=yc*pi/180;
0027 var=sin(1.5*yy); var=1.-var.*var;
0028 var(find(abs(yc) > 60.))=0.; sst2=27*var;
0029 %sst=273.15+sst2;
0030 %fname='SST_APE_1.bin';
0031 %- just for doing a plot
0032 yy=y1d*pi/180;
0033 var=sin(1.5*yy); var=1.-var.*var;
0034 var(find(abs(y1d) > 60.))=0.; sst2_1d=27*var;
0035
0036 %-- close to Symetric zonally-average current SST:
0037 yyp=abs(yc/90); phi=yyp.^3.;
0038 sst3=27.8*exp(-7*phi) - 1 ;
0039 sst=273.15+sst3;
0040 fname='SST_symEx3.bin';
0041 %- just for doing a plot
0042 yyp=abs(y1d/90); phi=yyp.^3.;
0043 sst3_1d=27.8*exp(-7*phi) - 1 ;
0044
0045 if kwr > 0,
0046 fid=fopen(fname,'w','b'); fwrite(fid,sst,'real*8'); fclose(fid);
0047 fprintf(['write file: ',fname,'\n']);
0048 end
0049
0050 figure(1);clf;
0051 subplot(211)
0052 var=sst;
0053 grph_CS(var,xc,yc,xg,yg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
0054 title('SST');
0055
0056 subplot(212)
0057 plot(y1d,sst1_1d,'g-'); hold on;
0058 plot(y1d,sst2_1d,'b-');
0059 plot(y1d,sst3_1d,'r-');
0060 hold off;
0061 axis([-90 90 -5 29]); grid
0062 legend('cos','ape','sym')
0063 title('SST ^oC');
0064 %return
0065
0066 %-- make Q-flux file (similar to the one used by Paul, hard coded
0067 % in original version of mixed_layer.f90)
0068 qflx_ampl=50. ;
0069 qflx_width=90.;
0070 yy=yc/qflx_width;
0071 yy=yy.*yy;
0072 qflx=qflx_ampl*(1-2*yy);
0073 qflx=qflx.*exp(-yy);
0074
0075 figure(2);clf;
0076 var=qflx;
0077 grph_CS(var,xc,yc,xg,yg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
0078 title('Q-flux [W/m^2]');
0079
0080 if kwr > 0,
0081 fname='Qflux_w90.bin';
0082 fid=fopen(fname,'w','b'); fwrite(fid,qflx,'real*8'); fclose(fid);
0083 fprintf(['write file: ',fname,'\n']);
0084 end
0085 %return
0086
0087 %-- make initial pot-temp field by adding noise to T(iter=0) output file:
0088
0089 rDir=gDir;
0090 namf='T'; it=0;
0091 tini=rdmds([rDir,namf],it);
0092 nr=size(tini,3);
0093
0094 var=rand([nx,ny]); var=var-mean(var(:));
0095 yy=yc*pi/90;
0096 var=var.*(2+cos(yy))/3;
0097 figure(3);clf;
0098 %var=sst0;
0099 grph_CS(var,xc,yc,xg,yg,ccB(1),ccB(2),shift,cbV,AxBx,kEnv);
0100
0101 noise=1.e-3;
0102 tini1=tini+noise*reshape(reshape(var,[nx*ny 1])*ones(1,nr),[nx ny nr]);
0103 %size(tini1)
0104
0105 if kwr > 0,
0106 fname='ini_theta.bin';
0107 fid=fopen(fname,'w','b'); fwrite(fid,tini1,'real*8'); fclose(fid);
0108 fprintf(['write file: ',fname,'\n']);
0109 end
0110
0111 %- spec-humid : put constant Rel-Humid in the lowest troposphere
0112 relhum=0.8 ; pHum=800.e+2;
0113 relhum=0.4 ; pHum=700.e+2;
0114 khum=max(find(G.rC > pHum));
0115
0116 %- taken from AIM -> qsat in g/kg, pIn = normalised Pressure
0117 P0=1.e+5; pIn=G.rC/P0;
0118 qsat=calc_Qsat(1,pIn,tini);
0119 qsat=reshape(qsat,[nx ny nr]);
0120 qini=qsat*1.e-3*relhum;
0121 qini(:,:,khum+1:end)=0;
0122
0123 figure(4);clf;
0124 pax=G.rC/100; %- in mb
0125 i1=1; j1=1;
0126 var=squeeze(qini(i1,j1,:));
0127 plot(var,pax,'k-'); hold on;
0128 i1=nc/2; j1=nc/2;
0129 var=squeeze(qini(i1,j1,:));
0130 plot(var,pax,'r-');
0131 i1=nc*2.5; j1=nc*0.5;
0132 var=squeeze(qini(i1,j1,:));
0133 plot(var,pax,'b-');
0134 hold off
0135 set(gca,'YDir','reverse');
0136 grid
0137 legend('mid','eq','pol');
0138 title('Q-ini profile');
0139
0140 if kwr > 0,
0141 fname=['ini_specQ_',int2str(nr),'l.bin'];
0142 fid=fopen(fname,'w','b'); fwrite(fid,qini,'real*8'); fclose(fid);
0143 fprintf(['write file: ',fname,'\n']);
0144 end
0145
0146 var=reshape(tini,[nx*ny nr]);
0147 var=mean(var);
0148 for n=1:ceil(nr/10)
0149 is=1+(n-1)*10; ie=min(nr,n*10);
0150 if n == 1, fprintf(' tRef='); else fprintf(' '); end
0151 fprintf(' %5.1f,',round(10*var(is:ie))/10);
0152 fprintf('\n')
0153 end
0154
0155 return
