%% Script to solve/plot 1-layer equations for Hadley Circulation, based on: % * Sobel and Schenider, 2009, "Single Layer Axisymmetric Model for a Hadley % Circulation with Parameterized Eddy Momentum Forcing", Journal of % Advances in Modeling Earth Systems, Volume 1, Article #10, 11 pp. % * Source hereafter referred to as "SS09" % % * Coded by Tim Cronin (twc) 6/21/2012-7/6/2012 % % * notes: % * 1) EMFD parameterization is difficult to stabilize % * 2) To reproduce SS09 results, the conversion from theta to T in the % meridional momentum equation must be turned off; otherwise their % thermal forcing gives zonal winds that are too weak % * 3) At present this conversion factor th_T is set to unity (effectively % commented out) % * 4) To aid with convergence, the numerics are smoothed both with an % Asselin Filter (modification of the leapfrog time-scheme) and an % empirical smoothing method that periodically replaces the % prognostic variables with their mean values over the past iplot % timesteps -- both of these should not affect the equilibrium % solution -- they should only suppress transients % % * Notes by Malte Jansen 6/24/2012: % * 1) Added a 2nd order Shapiro filter (essentially the same as hyperdiff) % This makes EMFD somewhat more stable and causes it to equilibrate better % It still tends to develop shock-like features - which conveniently % generally seem to travel outside of the plotted domain before equilibrating % * >> twc - 6/25/2012: added option for 4th order Shapiro filter % * 2) Probably would be good to have some kind of implicit scheme % Even without the EMFD issues, this should allow MUCH longer timestep, % which seems to be limited by gravity waves (with speeds on the order of 100m/s % - compare to meridional advection which typically has at most 1 m/s) % * >> twc - 7/6/2012 - Malte also recognized that initialization from an % isothermal state produces far fewer transients than initialization % from the radiative equilibrium temperature %% Parameter initialization %clear all; % Parameters for meridional grid (as in SS09 text, top of p. 4) % name value units description ngu = 800; % number of grid points for u and theta ngv = ngu+1; % number of grid points for v dy = 39.3e3; % m grid spacing in y dyi = 1/dy; % 1/m inverse grid spacing (multiplication clearer than division in some cases) y1 = 9.439e6; % m meridional location beyond which differential forcing vanishes y0 = 1e6; % m meridional location of peak RCE temperature y0arr = (0:0.2e6:2e6);% m array of y0 values to use for contour plots % if only one simulation is desired, set y0arr to y0 %y0arr = y0; % m % % % construct grid yu = dy*((1:ngu)-(1+ngu)/2)'; % y-grid for u, theta yv = dy*((1:ngv)-(1+ngv)/2)'; % y-grid for v % output arrays uout = zeros(length(y0arr),ngu); % output array for u vout = zeros(length(y0arr),ngv); % output array for v % Parameters from Table 1 of SS09 % name value units description tau = 37*86400; % days*s/day=s -- relaxation timescale for heating/cooling towards th_equil H = 16e3; % m -- Tropopause height delta = 4e3; % m -- depth of layers in which merdional flow occurs T0 = 300; % K -- reference temperature delthz = 60; % K -- vertical potential temperature stratification over depth H delthy = 50; % K -- RCE equator-to-pole temperature gradient th00 = 330; % K -- Background tropospheric-mean potential temperature epsu = 1e-8; % 1/s -- Background Rayleigh Drag coefficient vd = 2.5; % m/s -- Parameterized Eddy Momentum Flux coefficient kv = 7786; % m^2/s -- Diffusivity on v (meridional velocity) beta = 2e-11; % m^-1 s^-1 -- Meridional Gradient of Coriolis Parameter % Parameters in SS09 equations left undefined % name value units description g = 9.81; % m/s^2 -- gravitational acceleration (OK) ptrop = 200; % hPa -- tropopause pressure (guess) psurf = 1000; % hPa -- surface pressure (guess) th_T = 0.63;%(ptrop/psurf)^(2/7); % conversion factor from theta to T, % from non-numbered equation following 2.3 on p.2, this factor should be used % as a multiplier on d(theta)/dy in eqn. 2.2 to get dT/dy -- but it appears % to be unity (or very nearly so) based on their model results % Parameters for time-stepping scheme (--twc guesses/empirical choices for stable performance) % name value units description dt = 200; % s -- model timestep tmax = 5e7; % s -- end time of integration oshap = 2; % -- Order of Shapiro filter to be used (valid options: none, 2, 4) nshap = 20; % -- number of timesteps for Shapiro Filter smoothing nt = round(tmax/dt); % number of model timesteps iplot = 500; % -- plotting interval (in timesteps) afilt = 0.05; % -- Asselin filter coefficient fsmooth = 0.05; % -- every iplot steps, set prognostic variables % equal to the weighted mean of their current values and their time-mean values for the % last iplot steps, using u(new) = (1-fsmooth)*u(old) + fsmooth*ubar % set to zero to remove time-smoothing %% Loop over y0 values for iy0=1:length(y0arr) %% Variable initialization % 3 prognostic variables: theta, u, v y0=y0arr(iy0); % Define RCE theta profile (th_equil) based on SS09 eqn 3.2 th_equil = th00 - delthy*(sin(pi*yu/(2*y1)).^2-2*sin(pi*yu/(2*y1))*sin(pi*y0/(2*y1))); % set to a flat profile for |y|>y1 using array logic iythmin = find(yu<-y1); % find indices corresponding to values of yu less than -y1 iythmin = iythmin(end)+1; % minimum yu-index for which RCE theta profile is valid iythmax = find(yu>y1); % find indices corresponding to values of yu greater than +y1 iythmax = iythmax(1)-1; % maximum yu-index for which RCE theta profile is valid th_equil=th_equil.*(abs(yu)<=y1)+... % if |yu|<=y1 then use the RCE profile th_equil(iythmin)*(yu<-y1)+... % else if yu<-y1: use the RCE profile at its southernmost valid point th_equil(iythmax)*(yu>y1); % else if yu>y1: use the RCE profile at its northernmost valid point % theta (th*): MJ -- initialize to isothermal initial state th1=0*th_equil+th00-delthy; th2=th1; % th1, th2, th3 are values for use in leapfrogging/Asselin time scheme th3=th1; thext=[th1(1); th1; th1(end)]; % thext has size length(th1)+2; for use in difference eqns % zonal wind (u): initialize to zero everywhere %u1=-g*H/T0*th_T*0.5*dyi*(thext(3:end)-thext(1:end-2))./(beta*yu); % twc commented out initialization to thermal wind value 7/6 u1=zeros(ngu,1); u2=u1; u3=u1; uext=[u1(1); u1; u1(end)]; % meridional wind (v): initialize to zero v1=zeros(ngv,1); v2=v1; v3=v1; vext=[v1(1); v1; v1(end)]; % total tendency arrays for three variables fth=zeros(ngu,1); % K/s fu=zeros(ngu,1); % m/s^2 fv=zeros(ngv,1); % m/s^2 % theta-tendencies (K/s) fthadv=zeros(ngu,1); % z-advection: delta/H*delthz*dyi*diff(v+,v-) fthrad=zeros(ngu,1); % radiative relaxation: (th_equil-th)/tau % u-tendencies (m/s^2) fucorl=zeros(ngu,1); % coriolis: beta*yu*mean(v+,v-) fuadvy=zeros(ngu,1); % y-advection: if v>0, dyi*(v-)*diff(u,u-); if v<0, dyi*(v+)*diff(u+,u) fuadvz=zeros(ngu,1); % z-advection: dyi*heaviside{diff(v+,v-)}*u*diff(v+,v-) frayd=zeros(ngu,1); % rayleigh drag: epsu*u femfd=zeros(ngu,1); % eddy momentum flux divergence array % v-tendencies (m/s^2) fvcorl=zeros(ngv,1); % coriolis: 0.5*beta*yv*mean(u+,u-) fvadv=zeros(ngv,1); % y-advection: if v>0, dyi*(v)*diff(v,v-); if v<0, dyi*(v)*diff(v+,v) fvpgrad=zeros(ngv,1); % pressure gradient: 0.5*g*H/T0*dyi*diff(th+,th-) fvdiff=zeros(ngv,1); % diffusive: kv*dyi*dyi*((v+)-2v+(v-))/2 % time-averages for plotting and/or using to accelerate convergence to % steady state solutions ubar=0*u1; vbar=0*v1; thbar=0*th1; figure(1); %huline=plot(yu,u2,'b'); %hold on; %hvline=plot(yv,v2,'g'); %hthline=plot(yu,(th2-th_equil),'r'); plot(yu,u2,yv,v2,yu,(th2-th_equil)); xlim([-1.5e7 1.5e7]); ylim([-30 100]); legend('u','v','\theta-\theta_E') drawnow; pause(1.0); %% Integration in time for it=1:nt % extended arrays for finite differencing uext=[0; u2; 0]; % assume u->0 at both ends %uext=[u2(1); u2; u2(end)]; % assume du/dy->0 at both ends - causes crazy oszillations in v... vext=[0; v2; 0]; % assume v->0 at both ends thext=[th2(1); th2; th2(end)]; % assume dtheta/dy is zero at both ends % CALCULATE THETA-TENDENCIES fthadv=delta/H*delthz*dyi*(v2(2:end)-v2(1:end-1)); fthrad=(th_equil-th2)/tau; % total theta-tendency is sum of components fth=-fthadv+fthrad; % CALCULATE U-TENDENCIES fucorl=beta*yu.*(v2(1:end-1)+(v2(2:end)))/2; fuadvy=dyi*(v2(1:end-1).*((v2(1:end-1)+v2(2:end))>0).*(uext(2:end-1)-uext(1:end-2))+... (v2(2:end).*((v2(1:end-1)+v2(2:end))<=0).*(uext(3:end)-uext(2:end-1)))); fuadvz=dyi*((v2(2:end)-v2(1:end-1))>0).*(u2.*(v2(2:end)-v2(1:end-1))); frayd=epsu*u2; %Leapfrog or not seem about similarly (un)stable for this term femfd=vd.*(u2>0).*sign(yu).*dyi.*[0;u2(3:end)-u2(1:end-2);0]./2; %femfd=vd.*(u1>0).*sign(yu).*dyi.*[0;u1(3:end)-u1(1:end-2);0]./2; %femfd=vd.*(uext(3:end)>0).*(u2>0).*(uext(1:end-2)>0).*sign(yu).*dyi.*... % ((uext(3:end)-u2).*(yu>0)+(u2-uext(1:end-2)).*(yu<0)); %femfd=femfd.*((uext(3:end)-uext(2:end-1)).*(uext(2:end-1)-uext(1:end-2))>0); %femfd=femfd.*(abs(femfd)=emfdmax); % tried different defs, setting to 0 due to stability issues %femfd=0; % total u-tendency is sum of components fu=fucorl-fuadvy-fuadvz-frayd-femfd; % CALCULATE V-TENDENCIES fvcorl=0.5*beta*yv.*(uext(1:end-1)+uext(2:end))/2; fvadv=dyi*(v2.*(v2>0).*(vext(2:end-1)-vext(1:end-2))+... (v2.*(v2<=0).*(vext(3:end)-vext(2:end-1)))); fvpgrad=0.5*g*H/T0*th_T*dyi*(thext(2:end)-thext(1:end-1)); %use leapfrog for diff? - doesn't seem to make an obvious difference fvdiff=kv*dyi*dyi*(vext(3:end)-2*v2+vext(1:end-2))/2; %fvdiff=kv*dyi*dyi*(v1ext(3:end)-2*v1+v1ext(1:end-2))/2; % total v-tendency is sum of components fv=-fvcorl-fvadv-fvpgrad+fvdiff; % ADVANCE QUANTITIES ONE STEP u3=u1+2*dt*fu; v3=v1+2*dt*fv; th3=th1+2*dt*fth; % APPLY SHAPIRO FILTER to eliminate gridscale noise every nshap % timesteps % % Use filter of order 'oshap' % % See: Shapiro, 1975, "Linear Filtering," Mathematics of Computation, % Volume 29, #132, pp. 1094-1097 % % Table 2: Stencils for different-order filters: % filter order x(i) x(i+/-1) x(i+/-2) x(i+/-3) x(i+/-4) % 1 2 1 % 2 10 4 -1 % 3 44 15 -6 1 % 4 186 56 -28 8 -1 % if mod(it,nshap)==0 if oshap==2 uext4=[0;0;u3;0;0]; vext4=[0;0;v3;0;0]; thex4=[th3(1);th3(1);th3;th3(end);th3(end)]; u3=1/16*(-uext4(1:end-4)+4*uext4(2:end-3)+10*uext4(3:end-2)+4*uext4(4:end-1)-uext4(5:end)); v3=1/16*(-vext4(1:end-4)+4*vext4(2:end-3)+10*vext4(3:end-2)+4*vext4(4:end-1)-vext4(5:end)); th3=1/16*(-thex4(1:end-4)+4*thex4(2:end-3)+10*thex4(3:end-2)+4*thex4(4:end-1)-thex4(5:end)); elseif oshap==4 uext8=[0;0;0;0;u3;0;0;0;0]; vext8=[0;0;0;0;v3;0;0;0;0]; thex8=[th3(1);th3(1);th3(1);th3(1);th3;th3(end);th3(end);th3(end);th3(end)]; u3=1/256*(186*u3+56*(uext8(4:end-5)+uext8(6:end-3))-28*(uext8(3:end-6)+uext8(7:end-2))+8*(uext8(2:end-7)+uext8(8:end-1))-1*(uext8(1:end-8)+uext8(9:end))); v3=1/256*(186*v3+56*(vext8(4:end-5)+vext8(6:end-3))-28*(vext8(3:end-6)+vext8(7:end-2))+8*(vext8(2:end-7)+vext8(8:end-1))-1*(vext8(1:end-8)+vext8(9:end))); th3=1/256*(186*th3+56*(thex8(4:end-5)+thex8(6:end-3))-28*(thex8(3:end-6)+thex8(7:end-2))+8*(thex8(2:end-7)+thex8(8:end-1))-1*(thex8(1:end-8)+thex8(9:end))); end end % APPLY ASSELIN TIME-FILTER u1=u2+afilt*(u1-2*u2+u3); v1=v2+afilt*(v1-2*v2+v3); th1=th2+afilt*(th1-2*th2+th3); u2=u3; v2=v3; th2=th3; % accumulate average values over iplot timesteps ubar=ubar+1/iplot*u2; vbar=vbar+1/iplot*v2; thbar=thbar+1/iplot*th2; if (mod(it,iplot)==0) % plot the three prognostic variables every iplot timesteps figure(1) plot(yu,0.1*u2,yv,10*v2,yu,(th2-th_equil)); xlim([-1e7 1e7]); title({['day = ' num2str(it*dt/86400+0.0001,4) ]; '{\color{blue} 0.1*u (m/s)},{\color{green} 10*v (m/s)},{\color{red} \theta-\theta_E (K)}'}); xlabel('meridional distance y, m'); ylabel({'{\color{blue} 0.1*u (m/s)},{\color{green} 10*v (m/s)},{\color{red} \theta-\theta_E (K)}'}); %plot(yu,ubar,yv,vbar,yu,(thbar-th_equil)); drawnow; % every iplot timesteps, smooth prognostic variables towards their % mean values over the past averaging period th1=(1-fsmooth)*th1+fsmooth*thbar; th2=th1; th3=th1; u1=(1-fsmooth)*u1+fsmooth*ubar; u2=u1; u3=u1; v1=(1-fsmooth)*v1+fsmooth*vbar; v2=v1; v3=v1; % reset average values to zero thbar=0*thbar; ubar=0*ubar; vbar=0*vbar; % if the model goes unstable, decrease the timestep and try again if (sum(isnan(u2))>0) dt=dt/1.5; fprintf(1, 'WARNING: NaNs found in zonal wind; model probably unstable\n'); fprintf(1, '\tDecreasing timestep by factor of 1.5, to %4.2f s \n', dt); fprintf(1, '\tWill attempt to re-initialize with th=th_equil, u=v=0 \n'); % reset prognostic variables so condition doesn't immediately activate % again th1=0*th_equil+th00-delthy; th2=th1; th3=th1; thbar=0*th1; u1=zeros(ngu,1); u2=u1; u3=u2; ubar=0*u1; v1=zeros(ngv,1); v2=v1; v3=v2; vbar=0*v1; end end end % time integration loop % set u and v output arrays based on results at end of simulation uout(iy0,:) = u3; vout(iy0,:) = v3; end % loop over y0 values %% Plot output as in SS09 Figures 6 or 7 if length(y0arr) is greater than 1 if length(y0arr)>1 % make plot like 7a (u) in SS09 figure('Position',[1 1 335 276],'Color','w'); contourf(yu,y0arr,uout,(-50:10:80)); shading flat; caxis([-40 80]); colorbar; xlim([-1e7 1e7]); set(gca,'XTick',[-1e7;0;1e7]); xlabel('y (m)'); ylim([0.2e6 2e6]); set(gca,'YTick',[0.5e6; 1e6; 1.5e6; 2e6]); ylabel('y_0 (m)'); title('u (m/s)'); % twc - 7/6/2012 -- print commands are specialized for my file structure; % they can be modified % %figname='C:\DOCUMENT\12.980\HadCirc\SobelSchneider7a-u-reproduced.emf'; %print('-dmeta','-painters',figname); %figname='C:\DOCUMENT\12.980\HadCirc\SobelSchneider7a-u-reproduced.pdf'; %print('-dpdf',figname); % make plot like 7b (v) in SS09 figure('Position',[1 1 335 276],'Color','w'); contourf(yv,y0arr,vout,(-3:0.5:3)); shading flat; caxis([-3 3]); colorbar; xlim([-1e7 1e7]); set(gca,'XTick',[-1e7;0;1e7]); xlabel('y (m)'); ylim([0.2e6 2e6]); set(gca,'YTick',[0.5e6; 1e6; 1.5e6; 2e6]); ylabel('y_0 (m)'); title('v (m/s)'); % twc - 7/6/2012 -- print commands are specialized for my file structure; % they can be modified % %figname='C:\DOCUMENT\12.980\HadCirc\SobelSchneider7b-v-reproduced.emf'; %print('-dmeta','-painters',figname); %figname='C:\DOCUMENT\12.980\HadCirc\SobelSchneider7b-v-reproduced.pdf'; %print('-dpdf',figname); end