Numerical modeling of Atmosphere and Ocean

[kkvarma]

Dear Dr. Vinu,
I could plot the height. The figure shows height of about 200 near the equator and near the northern boundaries. Could this be the expected result? I am not very sure whether the model run was correct? Is there any modification to be done in the rgmodel?
regards
kkvarma

[vinu]

Solution with the given forcing.

If the model run with the given forcing provided in the code,
(a simple zonal wind stress profile represented by a sin function
with easterlies (trade winds) in the equatorial region and
westerlies on the northern domain, the solution will be a
gyre circulation, with strong and narrow western boundary
currents (parallel to the western coast)
and broad and slow eastern boundary currents
(parallel to the eastern coast) and zonal flow everywhere else.
The corresponding thermocline (or thickness)
shape will be, larger at the center of the
domain and gradient to the domain boundaries, with strongest gradient
in the western boundary (represents the strong western boundary currents).

An example plot is given below.

(Optional)
Also, while plotting, the one grid border at the
domain boundaries can be omitted,
because that is used for applying boundary condition and is treated as
land cell. This border grid can be avoided in the plot by limiting
the output to a dimension (lon-2)*(lat-2). Following alterations
in the code may help to avoid this land cell from being printed in the output.

(with reference to the code given in the post).

line-78 & recl=(lon-2)*(lat-2)*4)
line-93 write(1,rec=lrec)((h(i,j,3),i=2,lon-1),j=2,lat-1)

Also, in the corresponding rgmodel.ctl file, the dimension values
50 and 25 has to replace with 48 and 23, respectively.





-
Thanks

[nonicad]

Do have code fortran prog fom the book :Zygmuny kowalik,

[vinu]

Sorry, I dont have the code you requested,

-Vinu

[cd264]

Dear vinu:
quite useful post for me!
long for your new lectrue.

[cd264]

hi! vinu
Is the model can be apply for the globe circulation simulation?how large can the model to simulate?what about the boundary conditions and how to prevent the unstable?

Also,happy newyear！

[vinu]

Hi Tom,

The model we were discussing in earliers posts is only a simplified
representation of thermocline without any consideration on
non-linearity and other surface ocean processes. Moreover, the
application of betaplane is limitted to only a few degrees because the
beta-value has variation with latitude (cosine function of latitude).

The above model doesnot represent the mixed layer, and thermodynamics.
Additing those two, along with non-linearity becomes a powerful tool for
tropical or midlatitude circulation. Numerous literature exists describing
the layered model with dynamics and thermodynamics and used for
several application. McCreary et al. (1994), various papers of
Jensen (1990,94,98,2003), Murtugudde et al. (1998), Qu et al. (1994,..)
etc are some examples.

The model described here, thus, is not a suitable choice for any global
circulation study. But indeed, this model could resolve some of the
fundamental structure of ocean dynamcis, such as wind driven
gyre circulation for a limitted extend of latitudilan width for which the
beta plane approximation is applicable with a given value of beta or
equatorial-mid latutude wave dynamics etc.. I have noted it in the
orignal post (second paragraph). http://www.oceanographers.net/forums...ead.php?p=1254

Literature exists using layer model with only dynamics and used for particular
application such as equatorial dynamics. Paper by Potemra et al. (2001), JGR,
Vol,106, 2407-2422 used layer model to study the propagation of
equatorlial rossby waves in the pacific and its contribution to Indian Ocean
Rossby waves.

-Thanks

[Malte]

Dear Forum,
i have tried to translate the rgmodel fortran code into a matlab file, but even though i have choosen a real small timestep compared to delta_x, the model keeps 'blowing' up, if I use a small disturbance of h as initial condition.

I tried really long finding the error in my code, but haven't been successful yet. I attach it to this post, maybe someone sees quick what goes wrong, without spending too much time on it. I am sitting on it for a few days now and still no insight...
Thank you very much in advance!
Cheers,
Malte

[Malte]

Maybe I'll post the numerical schemes here, so one doesn't have to download the file.

The integration scheme (centered in time and space):

    for j=2:jmax-1
        for i=2:imax-1
            
        % calculate the coriolis force at the same position, one step in the    past (km): 
        coriolis = f(j)*v(i,j,km);
        
        % calculate the horizontal pressure gradient 
        % (finite diff. centered in y space of the first derivative)::
        horiz_pressure = -g_reduced*(h(i+1,j,k) - h(i-1,j,k))/d2x; 
        
        % Diffusion term: A_H * Laplace*U
        % Laplace*U = (d^2 U / d x^2) + (d^2 U / d y^2)
        % Centered finite difference of the second derivative
        diffx = (u(i-1,j,k) - 2*u(i,j,k) + u(i+1,j,k))/dx2;
        diffy = (u(i-1,j,k) - 2*u(i,j,k) + u(i+1,j,k))/dy2;
        diffusion = A_H*(diffx + diffy);
        
        % windstress
        windstress = 0;
        
        % calculate U:
        u(i,j,kp)=u(i,j,km)+delta_t*(coriolis-horiz_pressure+diffusion+windstress);
        end
    end 
    
 
    
    % Calculate meridional transport (y-direction) 
    % Finite Difference Formulation of dV/dt:
    for j=2:jmax-1
        for i=2:imax-1

        % calculate the coriolis force at the same position, one step in the past (km): 
        coriolis = -f(j)*u(i,j,km);
        
         % calculate the horizontal pressure gradient 
         % (finite diff. centered in y space of the first derivative):
        horiz_pressure = -g_reduced*(h(i,j+1,k) - h(i,j-1,k))/d2y; 
        
        
        % Diffusion term: A_H * Laplace*V
        % Laplace*V = (d^2 V / d x^2) + (d^2 V / d y^2)
        % Centered finite difference of the second derivative
        diffx = (v(i-1,j,k) - 2*v(i,j,k) + v(i+1,j,k))/dx2;
        diffy = (v(i-1,j,k) - 2*v(i,j,k) + v(i+1,j,k))/dx2;
        diffusion = A_H*(diffx + diffy);
        
        % windstress
        windstress = 0;
        
        % calculate V:
        v(i,j,kp)=v(i,j,km)+delta_t*(coriolis+horiz_pressure+diffusion+windstress);
        end
    end
        
    
    
    
% Continuity equation.
% Using forward finite difference instead of centered finite difference!
    for j=2:jmax-1
        for i=2:imax-1
          
          % hdiv = dU/dx + dV/dy  
          % use centered finite difference:
         % hdiv = (-H*(u(i+1,j,k) - u(i-1,j,k))/d2x) + ((v(i,j+1,k) - v(i,j-1,k))/d2y);
          hdiv = (-H*((u(i+1,j,k)-u(i-1,j,k))/d2x+(v(i,j+1,k)-v(i,j-1,k))/d2y)+ A_H*((h(i+1,j,k)+h(i-1,j,k)-2.0*h(i,j,k))/dx2 + (h(i,j+1,k)+h(i,j-1,k)-2.0*h(i,j,k))/dy2));
          % new elevation of the thermocline h:  
          h(i,j,kp) = h(i,j,km) + delta_t*hdiv;  
          
          
        end
    end

with the boundary cond. applied every timestep:

 for j=1:jmax                  
     
     % West:
     u(1,j,kp) = -u(3,j,kp);     % no normal flow across the rigid boundary
     u(2,j,kp) =  0.0;           % no-slip at the boundary
     
     v(1,j,kp) = -v(3,j,kp);     % same as above but for v-velocity
     v(2,j,kp) =  0.0;           %
     
     h(1,j,kp) = h(3,j,kp);      % thickness restored to neighbouring value
     
     % East
     u(imax,j,kp) = -u(imax-2,j,kp); 
     u(imax-1,j,kp) = 0.0;
     
     v(imax,j,kp) = -v(imax-2,j,kp);
     v(imax-1,j,kp) = 0.0;
     
     h(imax,j,kp) = h(imax-2,j,kp);
    end

    for i=1:imax
     
     % North:
     u(i,1,kp) = -u(i,3,kp);     
     u(i,2,kp) =  0.0;           
     
     v(i,1,kp) = -v(3,j,kp);    
     v(i,2,kp) =  0.0;           
     
     h(i,j,kp) = h(i,3,kp);    
    

     % South:
     u(i,jmax,kp) = -u(i,jmax-2,kp);     
     u(i,jmax-1,kp) =  0.0;           
     
     v(i,jmax,kp) = -v(i,jmax-2,kp);    
     v(i,jmax-1,kp) =  0.0;           
     
     h(i,jmax,kp) = h(i,jmax-2,kp);  
    
   end

[vinu]

HI Malte,

I have seen one of your (other post?) asking about C.F.L. limits.
I think it is this post only.

I havent checked your code in detail, but a quick look found that

diffx = (u(i-1,j,k) - 2*u(i,j,k) + u(i+1,j,k))/dx2;
diffy = (u(i-1,j,k) - 2*u(i,j,k) + u(i+1,j,k))/dy2;

your diffx and diffy are doing same operation. The first should be derivative
in x and second should be in Y. that means. u(i,j-1,k) - 2*u + u(i,j+1,k)

If I could find any other fixes I will report,

Also please see the "Numerical ocean model" post for C.F.L. details.




-Vinu