[MITgcm-support] periodic boundary conditions in X, walls in Y
jeff polton
jpolton at ucsd.edu
Wed Sep 29 14:08:13 EDT 2004
Hi,
The temperature bounds on my model run should be between 16-18 degrees
(initially uniform T=16, no salt). Over 1 year of integration the
volume integral of heat decreases by about 0.25% (which I guess is
pretty small) because there is a persistent cooling near the lower
boundary with a lowest temp of 15.3 degrees. Also by this time approx
half the domain is less than 16 degrees by some amount. The anomalous
heating is not so pronounced and is invisible to the eye meridional
depth profiles as it is restricted to the more finely resolved boundary
layer. But the max is around 18.2 degrees.
I am using the default 2nd order centred scheme, I think that keeping
track of the diffusion is going to be a desirable thing later on.
I guess my question is whether these sorts of drifts are normal or
whether I've done something daft somewhere (which is why I've tagged on
the data and gendata.m files, just incase).
Any advice would be much appreciated,
Thanks,
Jeff
**********
** data **
***********
# ====================
# | Model parameters |
# ====================
#
# Continuous equation parameters
&PARM01
tRef=20*16.0,
sRef=20*35.0,
viscAh=5.E-4,
viscAz=5.E-4,
no_slip_sides=.TRUE.,
no_slip_bottom=.TRUE.,
viscA4=1.E9,
diffKhT=1.E-5,
diffKzT=1.E-5,
diffKhS=1.E6,
diffKzS=1.E6,
f0=1E-4,
beta=0.E-11,
tAlpha=2.0E-4,
sBeta =0.,
gravity=9.81,
rhoConst=1000.0,
# reference density of sea water, rhoNil
rhoNil=1000.0,
heatCapacity_Cp=3900.0,
rigidLid=.TRUE.,
implicitFreeSurface=.FALSE.,
saltAdvection=.FALSE.,
saltForcing=.FALSE.,
saltStepping=.FALSE.,
eosType='LINEAR',
nonHydrostatic=.FALSE.,
momAdvection=.TRUE.,
implicitViscosity=.TRUE.,
implicitDiffusion=.TRUE.,
# vert diff (m2s-1)
ivdc_kappa=8.E-3,
readBinaryPrec=64,
&
# Elliptic solver parameters
&PARM02
cg2dMaxIters=1000,
cg2dTargetResidual=1.E-9,
&
# Time stepping parameters
&PARM03
nIter0=0.,
nTimeSteps=1000000,
# decreased timestep by 10 to keep things stable with new diff
deltaT=1000,
abEps=0.1,
pChkptFreq=10000000.,
chkptFreq=0.0,
dumpFreq=10000000.,
monitorFreq=1.E-5,
# relaxation timescales (s)
tauThetaClimRelax =3600.,
cAdjFreq=0,
&
# Gridding parameters
&PARM04
usingCartesianGrid=.TRUE.,
usingSphericalPolarGrid=.FALSE.,
dXspacing=31.25E3,
dYspacing=31.25E3,
delZ=2.,2.,4.,8.,16.,32.,64.,128.,256.,512.,512.,256.,128.,64.,32.,16.,8
.,4.,2.,2.,
&
# Input datasets
&PARM05
thetaClimFile='TempSurf.sin_y',
bathyFile='topog.box',
hydrogThetaFile=,
hydrogSaltFile=,
zonalWindFile='windx.sin_y',
meridWindFile=,
&
******************
** gendata.m **
******************
% This is a matlab script that generates the input data
% Dimensions of grid
nx=32;
ny=32;
nz=20;
% Nominal depth of model (meters)
H=2048;
% Size of domain
Lx=1e6;
% Horizontal resolution (m)
dx=Lx/nx;
% Rotation
f=1e-4;
% jeff's vetical grid
delZ=[2.,2.,4.,8.,16.,32.,64.,128.,256.,512.,512.,256.,128.,64.,32.,16.,
8.,4.,2.,2.];
z=cumsum(delZ,2);
ieee='b';
accuracy='real*8';
% Flat bottom at z=-Ho
hbot=-H*ones(nx,ny);
% Walls
hbot(:,end)=0;
fid=fopen('topog.box','w',ieee); fwrite(fid,hbot,accuracy); fclose(fid);
% Wind-stress
tauMax=0.1;
x=((1:nx)-0.5)/(nx-1); % nx-1 accounts for a solid wall
y=((1:ny)-0.5)/(ny-1); % ny-1 accounts for a solid wall
[X,Y]=ndgrid(x,y);
tau=tauMax*sin(pi*Y);
fid=fopen('windx.sin_y','w',ieee); fwrite(fid,tau,accuracy);
fclose(fid);
% Heating
T0s=16.0;
xh=((1:nx)-0.5)/(nx-1); % nx-1 accounts for a solid wall
yh=((1:ny)-0.5)/(ny-1); % ny-1 accounts for a solid wall
[X,Y]=ndgrid(xh,yh);
T0=T0s+(1+cos(pi*Y));
fid=fopen('TempSurf.sin_y','w',ieee); fwrite(fid,T0,accuracy);
fclose(fid);
On Sep 29, 2004, at 10:01 AM, chris hill wrote:
> Jeff/Jane,
>
> Jeff - how big are the overshoots/undershoots and which advection
> scheme are you using? With the default centred second order scheme
> there can be numerical extrema. These are controlled by using numerical
> diffusion with that scheme. There are other advection scheme options
> where extrema do not occur, however, these have less well defined
> implied diffusion.
>
> We have a couple of more experienced users doing channel runs. I have
> mailed them to see if they would be willing to do a modest amount of
> mentoring for you guys.
>
> Chris
>
> On Tue, 2004-09-28 at 18:01, jeff polton wrote:
>> Hi Jane,
>>
>> I have been trying similar simulations; zonal periodic channel flow
>> with no slip bcs, zonal wind stress and meridional heating. I guess
>> when you say that you have "dissipating T" you mean that the
>> temperature field is evolving to values outside the initial field and
>> forcing range. If so then I have also found this problem and supposed
>> it to be a symptom of the 2nd order Adams Bashforth advection scheme
>> (http://mitgcm.org/sealion/online_documents/node80.html).
>>
>> The bad news is that I haven't resolved the problem, so I was hoping
>> that someone with more experience than me would also respond to your
>> email...
>>
>> Jeff
>>
>>
>> On Sep 24, 2004, at 11:21 AM, Jane Jimeian wrote:
>>
>>> i'm resending this message because the first message didn't get
>>> posted.
>>> i've combined the two:
>>>
>>> i'm trying to set up a channel where there are perodic boundary
>>> conditions
>>> in X and walls in Y. i doesn't seem to work. i try running the
>>> model
>>> without
>>> the using bathyFile='topog.channel' and the results look ok. when
>>> i uses topog.channel then T starts to dissipate. any suggestions
>>> would be
>>> helpful.
>>> ---------------------------------------------------------------------
>>> --
>>> ----
>>> data.pkg
>>> ---------------------------------------------------------------------
>>> --
>>> ----
>>> # Packages
>>> &PACKAGES
>>> useOBCS=.FALSE.,
>>> &
>>> ---------------------------------------------------------------------
>>> --
>>> ----
>>> data
>>> ---------------------------------------------------------------------
>>> --
>>> ----
>>> # ====================
>>> # | Model parameters |
>>> # ====================
>>> #
>>> # Continuous equation parameters
>>> &PARM01
>>> tRef= 50*0,
>>> sRef= 50*35.,
>>> viscAz=1.E-4,
>>> viscAh=2.E1,
>>> no_slip_sides=.FALSE.,
>>> no_slip_bottom=.FALSE.,
>>>
>>> viscA4=0.E12,
>>> diffKhT=2.E1,
>>> diffKzT=1.E-5,
>>> diffKhS=2.E1,
>>> diffKzS=1.E-5,
>>> f0=1.e-4,
>>> beta=0.E-11,
>>> tAlpha=2.E-4,
>>> sBeta =0.E-4,
>>> gravity=9.81,
>>> gBaro=9.81,
>>> rigidLid=.FALSE.,
>>> implicitFreeSurface=.TRUE.,
>>> eosType='LINEAR',
>>> hFacMin=0.2,
>>> nonHydrostatic=.FALSE.,
>>> readBinaryPrec=64,
>>> globalFiles=.TRUE.,
>>> &
>>>
>>> # Elliptic solver parameters
>>> &PARM02
>>> cg2dMaxIters=1000,
>>> cg2dTargetResidual=1.E-13,
>>> cg3dMaxIters=400,
>>> cg3dTargetResidual=1.E-13,
>>> &
>>>
>>> # Time stepping parameters
>>> &PARM03
>>> niter0=0,
>>> nTimeSteps=6048,
>>> deltaT=200.0,
>>> abEps=0.1,
>>> pChkptFreq=604800.0,
>>> chkptFreq=0.0,
>>> dumpFreq=86400.,
>>> monitorFreq=86400.,
>>> &
>>>
>>> # Gridding parameters
>>> &PARM04
>>> usingCartesianGrid=.TRUE.,
>>> usingSphericalPolarGrid=.FALSE.,
>>> delX=200*5.e3,
>>> delY=50*5.e3,
>>> delZ=50*30,
>>> &
>>>
>>> # Input datasets
>>> &PARM05
>>> bathyFile='topog.channel',
>>> hydrogThetaFile='theta.bin',
>>> &
>>>
>>> ---------------------------------------------------------------------
>>> --
>>> ----
>>> gendata.m
>>> ---------------------------------------------------------------------
>>> --
>>> ----
>>>
>>> % This is a matlab script that generates the input data
>>>
>>> % Dimensions of grid
>>> nx=200;
>>> ny=50;
>>> nz=50;
>>> % Nominal depth of model (meters)
>>> H=1500;
>>> % Scale of bump (m)
>>> L=25e3;
>>> % Height of bump (m)
>>> %dh=0.90*H;
>>> % Horizontal resolution (m)
>>> dx=5e3;
>>> % Rotation
>>> f=1e-4;
>>> % Stratification
>>> N=1.5 * f*L/H;
>>>
>>> % Gravity
>>> g=9.81;
>>> % E.O.S.
>>> alpha=2.e-4;
>>>
>>> z=[15:30:1500];
>>> dz=30;
>>> % dz=H/nz;
>>> sprintf('delZ = %d * %7.6g',nz,dz)
>>>
>>> x=(1:nx)*dx;x=x-mean(x);x=x+nx/3*dx;
>>> y=(1:ny)*dx;y=y-mean(y);
>>> [Y,X]=meshgrid(y,x);
>>> % z=-dz/2:-dz:-H;
>>> %z=-(cumsum(dz)-dz/2);
>>>
>>> % Temperature profile
>>> %Tz=N^2/(g*alpha)
>>> %phi=sprintf(' %8.6g',Tz*z+sum(Tz*dz)/2);
>>> %Tref=str2num(phi);
>>> %[sprintf('Tref =') sprintf(' %8.6g,',Tz*z+sum(Tz*dz)/2)]
>>>
>>> ieee='b';
>>> accuracy='real*8';
>>>
>>> % Create three dimensional temperature profile
>>> t3d=zeros(nx,ny,nz);
>>> r0=20e3;
>>> r2=X.^2 + Y.^2;
>>> z0=300;
>>> for k=1:nz
>>> t3d(:,:,k)=10-10*z(k)/H;
>>> for i=1:nx
>>> for j=1:ny
>>> if Y(i,j) >= 0
>>> if z(k) >=240
>>> zp=z(k)-240;
>>> fz=(abs(zp)*exp(-abs(zp)/z0))/z0;
>>> else
>>> fz=0;
>>> end
>>> else
>>> if z(k) < 1260
>>> zp=z(k)-1260;
>>> fz=(abs(zp)*exp(-abs(zp)/z0))/z0;
>>> else
>>> fz=0;
>>> end
>>> end
>>> % Tp=fz.*(r0.*Y(i,j)./r0^2).*exp(-r2(i,j)./r0^2).*4.*2.7*2.7;
>>> zp=z(k)-750;
>>> Tp=3.*exp(-zp*zp/300^2);
>>> if r2(i,j) > 2*r0^2
>>> rp=r2(i,j)-2*r0^2;
>>> Tp=Tp*exp(-rp^2/r0^2);
>>> end
>>> t3d(i,j,k)=t3d(i,j,k)+Tp;
>>> % t3d(i,j,k)=Tp;
>>> end
>>> end
>>> end
>>> fid=fopen('theta.bin','w',ieee); fwrite(fid,t3d,accuracy);
>>> fclose(fid);
>>>
>>> % Simple channel
>>> h=zeros(nx,ny);
>>> %h(:,1)=0;
>>> h(:,2:ny-1)=-H;
>>> %h(:,ny)=0;
>>> fid=fopen('topog.channel','w',ieee); fwrite(fid,h,accuracy);
>>> fclose(fid);
>>>
>>> --
>>> _______________________________________________
>>> MITgcm-support mailing list
>>> MITgcm-support at mitgcm.org
>>> http://dev.mitgcm.org/mailman/listinfo/mitgcm-support
>>>
>>
>> _______________________________________________
>> MITgcm-support mailing list
>> MITgcm-support at mitgcm.org
>> http://dev.mitgcm.org/mailman/listinfo/mitgcm-support
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://dev.mitgcm.org/mailman/listinfo/mitgcm-support
>
More information about the MITgcm-support
mailing list