[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