[MITgcm-support] Temperature flux through solid boundaries?

Christopher L. Wolfe clwolfe at ucsd.edu
Mon Jan 22 16:46:57 EST 2007


I am running MITgcm in a semi-enclosed basin with a zonal channel 
spanning the lowest 1/8 of the basin and solid vertical walls 
elsewhere. I have surface temperature forcing, but no salinity. The 
temperature is advected via the flux limiting scheme 33 and dissipated 
through Laplacian and biharmonic diffusion. My problem is that the 
normal derivative of the temperature (estimated via centered 
differences or just "by eye) does not appear to go to zero at the 
vertical walls, implying a heat flux into (or through) the boundaries. 
I can figure out if there's something wrong with my configuration or if 
I'm just misinterpreting how the boundary condition is applied.

Any help understanding this would be greatly appreciated. The relevant 
configuration files follow,
Christopher

-----------------------------------------------------------
Dr. Christopher L. Wolfe             		 858-534-4560
Physical Oceanography Research Division    OAR 357
Scripps Institution of Oceanography, UCSD  clwolfe at ucsd.edu
-----------------------------------------------------------

data:
# ====================
# | Model parameters |
# ====================
#
# Continuous equation parameters
  &PARM01
  tRef=-9.9, -9.9, -9.9, -9.9, -9.9, -9.9, -9.9
       -9.9, -9.9, -9.9, -9.9, -9.9, -9.9, -9.9
       -9.9, -9.9, -9.9, -9.9, -9.9, -9.9
  sRef=20*35.,
# Prandtl number P = 25
  viscAh=10E-0,
  viscAz=2E-3,
  no_slip_sides=.TRUE.,
  no_slip_bottom=.FALSE.,
  viscA4=50.E7,
  diffK4T=0.E9,
  diffKhT=8E-5,
  diffKzT=8E-5,
  diffKhS=1.E0,
  diffKzS=1.E0,
  f0=-1.3683e-04,
  beta=3.4208e-11,
  tAlpha=2.E-4,
  sBeta =0.,
  gravity=10.,
  rhoConst=1000.,
  rhoNil=1000.,
  heatCapacity_Cp=3900.,
  rigidLid=.TRUE.,
  implicitFreeSurface=.FALSE.,
  saltAdvection=.FALSE.,
  saltForcing=.FALSE.,
  saltStepping=.FALSE.,
  eosType='LINEAR',
  nonHydrostatic=.FALSE.,
  momAdvection=.TRUE.,
  implicitViscosity=.TRUE.,
  implicitDiffusion=.FALSE.,
  readBinaryPrec=64,
  tempAdvScheme=33,
  staggerTimeStep=.TRUE.,
  bottomDragLinear=4.4E-4,
  debugLevel=-1,
  useJamartWetPoints=.TRUE.,
  &

# Elliptic solver parameters
  &PARM02
  cg2dMaxIters=40,
  cg2dTargetResidual=4.E-7,
  cg3dMaxIters=40,
  cg3dTargetResidual=4.E-7,
  &

# Time stepping parameters
  &PARM03
  nIter0=0007650000,
  nTimeSteps=225000,
  deltaT=500.,
  abEps=0.1,
  pChkptFreq=112.5E+06,
  chkptFreq=11.25E6,
  dumpFreq=112.5E6,
  monitorFreq=11.25E6,
  tauThetaClimRelax=764400.,
  cAdjFreq=0,
&

# Gridding parameters
  &PARM04
  usingCartesianGrid=.TRUE.,
  usingSphericalPolarGrid=.FALSE.,
  dXspacing=4.4643E3,
  dYspacing=4.4643E3,
  delZ=  12.0505,  14.9159,  18.4328,  22.7337,  27.9692,
         34.3067,  41.9245,  51.0025,  61.7057,  74.1598,
         88.4181, 104.4189, 121.9376, 140.5412, 159.5568,
        178.0710, 194.9769, 209.0770, 219.2372, 224.5642
  &

# Input datasets
  &PARM05
  thetaClimFile='TempSurf.two.sin_y',
  bathyFile='topo_two.box',
  hydrogThetaFile=,
  hydrogSaltFile=,
  zonalWindFile='windx.two.sin_y_clean',
  meridWindFile=,
  &



gendata.m:
% This is a matlab script that generates the input data

% Dimensions of grid
nx=448;
ny=1792;
nz=20;
% Nominal depth of model (meters)
H=2000;
% Size of domain
Lx=2e6;
Ly=8e6;
% Horizontal resolution (m)
dx=Lx/nx;
dy=Ly/ny;

ieee='b';
accuracy='real*8';

% Topog
% - semi-enclosed
hbot=-H*ones(nx,ny);
% Walls
hbot(:,end)=0;
hbot(end,:)=0;

% open channel for the most southern 1000 km.
Nc = round(1e6/dy);
hbot(end,1:Nc)=-H;
fid=fopen('topo_two.box','w',ieee); fwrite(fid,hbot,accuracy); 
fclose(fid);


% Wind-stress
tauMax=0.1;
sig = 0.18;
a = 0.8;

x=((1:nx)-0.5)/nx;
y=((1:ny)-0.5)/ny;
y = 4*y - 2;
[X,Y]=ndgrid(x,y);
%tau=tauMax*sin(3*pi*Y) + 0.00*rand(nx,ny);
tau=tauMax*(-cos(3*pi*y/4)+a*exp(-y.^2/(2*sig^2)));
fid=fopen('windx_two.weakeq','w',ieee); fwrite(fid,tau,accuracy); 
fclose(fid);


% Heating
Ts = -10; % coldest SH temp
Tn = -8.5;  % coldest NH temp
sigt=.5;
at=2.5;

%T0=.5*(Ts-Tn)*cos(pi*Y) + .5*(Ts+Tn)*cos(2*pi*Y);
T0=-.5*(Ts-Tn)*sin(pi*Y/4) - .5*(Ts+Tn)*cos(pi*Y/2);
T1=T0(1,:)+at*exp(-(y-2).^2/(2*sigt)^2);
fid=fopen('TempSurf.two_valveoff','w',ieee);
fwrite(fid,T0,accuracy); fclose(fid);



SIZE.h:
C $Header: /u/gcmpack/MITgcm/model/inc/SIZE.h,v 1.26 2001/09/21 
15:13:31 cnh Exp $
C $Name:  $

C
CBOP
C    !ROUTINE: SIZE.h
C    !INTERFACE:
C    include SIZE.h
C    !DESCRIPTION: \bv
C     *==========================================================*
C     | SIZE.h Declare size of underlying computational grid.
C     *==========================================================*
C     | The design here support a three-dimensional model grid
C     | with indices I,J and K. The three-dimensional domain
C     | is comprised of nPx*nSx blocks of size sNx along one axis
C     | nPy*nSy blocks of size sNy along another axis and one
C     | block of size Nz along the final axis.
C     | Blocks have overlap regions of size OLx and OLy along the
C     | dimensions that are subdivided.
C     *==========================================================*
C     \ev
CEOP
C     Voodoo numbers controlling data layout.
C     sNx :: No. X points in sub-grid.
C     sNy :: No. Y points in sub-grid.
C     OLx :: Overlap extent in X.
C     OLy :: Overlat extent in Y.
C     nSx :: No. sub-grids in X.
C     nSy :: No. sub-grids in Y.
C     nPx :: No. of processes to use in X.
C     nPy :: No. of processes to use in Y.
C     Nx  :: No. points in X for the total domain.
C     Ny  :: No. points in Y for the total domain.
C     Nr  :: No. points in Z for full process domain.
       INTEGER sNx
       INTEGER sNy
       INTEGER OLx
       INTEGER OLy
       INTEGER nSx
       INTEGER nSy
       INTEGER nPx
       INTEGER nPy
       INTEGER Nx
       INTEGER Ny
       INTEGER Nr
       PARAMETER (
      &           sNx =  224,
      &           sNy =  56,
      &           OLx =   3,
      &           OLy =   3,
      &           nSx =   1,
      &           nSy =   1,
      &           nPx =   2,
      &           nPy =   32,
      &           Nx  = sNx*nSx*nPx,
      &           Ny  = sNy*nSy*nPy,
      &           Nr  =  20)

C     MAX_OLX  - Set to the maximum overlap region size of any array
C     MAX_OLY    that will be exchanged. Controls the sizing of each
C                routine buffers.
       INTEGER MAX_OLX
       INTEGER MAX_OLY
       PARAMETER ( MAX_OLX = OLx,
      &            MAX_OLY = OLy )





More information about the MITgcm-support mailing list