/[MITgcm]/MITgcm/verification/MLAdjust/input.0.leith/genit.m
ViewVC logotype

Annotation of /MITgcm/verification/MLAdjust/input.0.leith/genit.m

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

Revision 1.1 - (hide annotations) (download)
Thu Mar 10 04:17:05 2005 UTC (19 years, 3 months ago) by baylor
Branch: MAIN
CVS Tags: checkpoint57t_post, checkpoint58l_post, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64t, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint57o_post, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint58e_post, checkpoint57v_post, checkpoint58u_post, checkpoint57f_post, checkpoint57s_post, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint57k_post, checkpoint57g_post, checkpoint64, checkpoint60, checkpoint61, checkpoint62, checkpoint63, checkpoint58r_post, checkpoint57i_post, checkpoint57y_post, checkpoint58g_post, checkpoint57x_post, checkpoint57m_post, checkpoint58n_post, checkpoint58x_post, checkpoint57g_pre, checkpoint58t_post, checkpoint58h_post, checkpoint58w_post, checkpoint58j_post, checkpoint57h_post, checkpoint57y_pre, checkpoint57f_pre, checkpoint58q_post, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint57r_post, checkpoint59, checkpoint58, checkpoint57h_done, checkpoint58f_post, checkpoint58d_post, checkpoint58c_post, checkpoint57w_post, checkpoint57p_post, checkpint57u_post, checkpoint58a_post, checkpoint58i_post, checkpoint57q_post, checkpoint58o_post, checkpoint57z_post, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint58y_post, checkpoint58k_post, checkpoint58v_post, checkpoint58s_post, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint58p_post, checkpoint61a, checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61q, checkpoint57j_post, checkpoint61z, checkpoint61x, checkpoint61y, checkpoint58b_post, checkpoint57h_pre, checkpoint58m_post, checkpoint57l_post 
Added verification examples of Leith, LeithD, and Smagorinsky with mnc (incl. restart from pickup.)  Model setup is Mixed-Layer Rossby Adjustment in an oceanic channel.  May be adding idealized surface diurnal forcing and restoring near walls soon.
1 baylor 1.1 %function geninit()
2     % This function generates initial condition files for the mitgcm.
3     % if init_vel=1, then the front problem is used
4     % if init_vel=0, then resting (other than small symmetry-breaking
5     % random noise) initial conditions are used.
6    
7     nx=50
8     ny=nx/2+1
9     nz=40
10     init_vel=1;
11     dxspacing=1000
12     dyspacing=dxspacing
13    
14     Lx=dxspacing*nx
15     Ly=dyspacing*ny
16    
17     %-- Params
18     g=9.81;
19     tAlpha=-2e-4;
20     f0=7.29e-5;
21     rho=1035
22     day=24*60^2;
23     prec='real*8';
24     ieee='b';
25    
26     H=200; %Max depth
27    
28     %-- Grid: x
29     dx=ones(1,nx); % uniform resolution
30     dx=dx*Lx/sum(dx);
31     xf=cumsum([0 dx]); % Face x points
32     xc=(xf(1:end-1)+xf(2:end))/2; % Centered x points
33    
34     %-- Grid: y
35     dy=ones(1,ny); % uniform resolution
36     dy=dy*Ly/sum(dy);
37     yf=cumsum([0 dy]); % Face y-points
38     yc=(yf(1:end-1)+yf(2:end))/2; %Centered y-points
39     L=yc(end)-yc(1); % this takes into account the wall of topography!!!
40    
41     %-- Grid: z
42     dh=H/nz*ones(1,nz);
43     zf=-round(cumsum([0 dh])/1)*1; % Face z points
44     dh=-diff(zf);
45     zc=(zf(1:end-1)+zf(2:end))/2; % centered z points
46     nz=length(dh);
47     H=sum(dh)
48    
49     [XT,YT,ZT]=ndgrid(xc,yc,zc); % This is the centered, temperature grid.
50     [XU,YU,ZU]=ndgrid(xf,yc,zc); % This is the U grid.
51     [XV,YV,ZV]=ndgrid(xc,yf,zc); % This is the V grid.
52     [XW,YW,ZW]=ndgrid(xc,yc,zf); % This is the W grid.
53     [XB,YB]=ndgrid(xc,yc); % This is the Bathymetry grid.
54    
55     % Stratification Params
56     Dml=200; % Mixed-layer Depth
57     D=Dml;
58     Lf=2*1e3; % Half-Width of Front
59     y0=yc(round(length(yc)/2)); % centered location of Front
60    
61     Ms=(2e-4)^2 %db/dy 2e-4 is about 1 K/50 km
62     Ns=9*Ms^2/f0^2 %db/dz
63    
64     deltatheta=Lf*Ms/g/tAlpha
65     T0=17
66     dthetadz=-Ns/g/tAlpha
67    
68     %-- Bathymetry
69     % Bathymetry is on Xc, Yc
70     hh=ones(size(XB));
71     hh(:,end)=0*hh(:,end);
72     hh=-H*hh;
73    
74     figure(1)
75     subplot(221)
76     pcolor(XB/1e3,YB/1e3,hh)
77     axis equal
78    
79     subplot(223)
80     plot(xc/1e3,dx/1e3,'.'); xlabel('X (km)');ylabel('\Delta x (km)')
81     title('\Delta x')
82    
83     subplot(222)
84     plot(dy/1e3,yc/1e3,'.'); ylabel('Y (km)');xlabel('\Delta y (km)')
85     title('\Delta y')
86    
87     subplot(224)
88     plot(dh,zc,'.'); ylabel('Z (m)');xlabel('\Delta z (m)')
89     title('\Delta z')
90    
91     fid=fopen('topo_sl.bin','w',ieee); fwrite(fid,hh,prec); fclose(fid);
92    
93     % Initial Temp Profile
94     figure(2)
95     subplot(4,2,1)
96     theta=T0+dthetadz*(ZT-ZT(1,1,end))+deltatheta*tanh((YT-y0)/Lf)/2;
97    
98     % Impose a strong initial and restoring mixed layer to Dml
99     i=1;
100     iml=1;
101    
102     while ((ZT(1,1,iml)>-Dml)&(iml<length(ZT(1,1,:))))
103     iml=iml+1;
104     end
105    
106     while (i<iml)
107     theta(:,:,i)=theta(:,:,iml);
108     i=i+1;
109     end
110    
111     [ZT(1,1,1) theta(1,1,1) theta(1,end,1) ;ZT(1,1,i) theta(1,1,i) theta(1,end,i);ZT(1,1,end) theta(1,1,end) theta(1,end,end)]
112    
113    
114     subplot(4,2,1)
115     [h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)));
116     axis([(yc(round(length(yc)/2))-1.5*Lf)/1e3 (yc(round(length(yc)/2))+1.5*Lf)/1e3 -Dml 0])
117     colorbar
118     title('Potl Temp')
119     xlabel('y (km)');ylabel('z (m)')
120    
121     subplot(4,2,2)
122     [h,c]=contourf(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(theta(1,:,:)));
123     colorbar
124     title('Potl Temp')
125     xlabel('y (km)');ylabel('z (m)')
126    
127    
128     subplot(4,2,3)
129     dens=theta*rho*tAlpha+rho;
130     [h,c]=contour(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(dens(1,:,:)));
131     title('Density')
132     xlabel('y (km)');ylabel('z (m)')
133     axis([(yc(round(length(yc)/2))-1.5*Lf)/1e3 (yc(round(length(yc)/2))+1.5*Lf)/1e3 -Dml 0])
134    
135     %Spice
136    
137     x1=xc(round(length(xc)/4)); % centered location of Front #1
138     x2=xc(round(3*length(xc)/4)); % centered location of Front #2
139     spice=T0+dthetadz*(ZT-ZT(1,1,end))+deltatheta*(tanh((XT-x1)/Lf)-tanh((XT-x2)/Lf)-1)/2;
140    
141     i=1
142     while (i<iml)
143     spice(:,:,i)=spice(:,:,iml);
144     i=i+1;
145     end
146    
147    
148     subplot(4,2,8)
149     [h,c]=contourf(squeeze(XT(:,1,:))/1e3,squeeze(ZT(:,1,:)),squeeze(spice(:,1,:)));
150     title('Spice')
151     xlabel('x (km)');ylabel('z (m)')
152    
153     subplot(4,2,4)
154     [h,c]=contour(squeeze(YT(1,:,:))/1e3,squeeze(ZT(1,:,:)),squeeze(dens(1,:,:)));
155     title('Density')
156     xlabel('y (km)');ylabel('z (m)')
157    
158     subplot(4,2,7)
159     decay=dyspacing;
160    
161     figure(2)
162     plot(YB(1,:)/1e3,hh(1,:))
163     title('topography')
164     xlabel('x (km)');ylabel('value')
165    
166     %-- Perturb Initial temperature
167     pert=rand(size(theta(:,:,1)));
168     pert=1e-5*(pert-0.5);
169     for i=1:length(theta(1,1,:))
170     theta(:,:,i)=theta(:,:,i)+pert;
171     end
172    
173     fid=fopen('thetaInitial.bin','w',ieee);
174     fwrite(fid,theta,prec); fclose(fid);
175    
176     fid=fopen('spiceInitial.bin','w',ieee);
177     fwrite(fid,spice,prec); fclose(fid);
178    
179    
  ViewVC Help
Powered by ViewVC 1.1.22