[MITgcm-support] periodic boundary conditions in X, walls in Y
Martin Losch
mlosch at awi-bremerhaven.de
Thu Sep 30 02:31:47 EDT 2004
Jeff,
I have not followed the discussion closely, but from a quick glace at
your data and gendata.m I can see that you use temperature restoring:
&PARM03
> tauThetaClimRelax =3600.,
> &
> &PARM05
> thetaClimFile='TempSurf.sin_y',
> &
with a short restoring timescale of 1 h. For debugging, I would try
turning that off (tauThetaClimRelax=0., or comment it out) and see,
whether you lose of heat. I predict, that you won't.
Also, have you checked that you resolve the boundary layers properly?
The Munk boundary layer width for harmonic friction=pi*(Ah/beta)^(1/3)
and for biharmonic friction = pi*(A4/beta)^(1/5).
I realize that you set beta=0, but there's probably something
equivalent for that case.
Martin
On Sep 29, 2004, at 8:08 PM, jeff polton wrote:
> 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
>>
>
> _______________________________________________
> 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