[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