[MITgcm-support] Fwd: (no subject)
Mohammad Akbarinasab
akbarinasabmohamad at gmail.com
Mon Dec 19 01:21:36 EST 2016
Hello
My problem is that the next time there is an increase in salinity
Please advise where's the problem?
data is :
# ====================
# | Model parameters |
# ====================
#
# Continuous equation parameters
&PARM01
Tref = 30*25.,
sRef= 10.54707989, 11.37348106, 12.39817393, 13.11977195, 14.4565467,
15.66976952, 16.92519145, 18.25814781, 19.56011301, 20.57347013,
21.52479808, 22.54991112, 23.44316005, 24.20244213, 25.52131167,
26.28910276, 27.62253252, 28.75248372, 29.46196077, 31.17486667,
32.32473282, 33.11892958, 33.48091282, 34.13403135, 34.20672289,
34.86205116, 35.37312434, 35.51936584, 35.73891183, 36.39887237,
viscAz=1.E-3,
viscAh=1.E-2,
no_slip_sides=.TRUE.,
no_slip_bottom=.TRUE.,
diffKhT=1.E-2
diffKzT=1.E-3,
f0=0.0,
beta=0.E-11,
eosType='LINEAR',
tAlpha=2.E-4,
sBeta =0.E-4,
gravity=9.81,
implicitFreeSurface=.TRUE.,
exactConserv=.TRUE.
nonHydrostatic=.TRUE.,
hFacMin=0.2,
implicSurfPress=0.5,
implicDiv2DFlow=0.5,
# nonlinFreeSurf=3,
hFacInf=0.2,
hFacSup=1.8,
saltStepping=.TRUE.,
#- not safe to use globalFiles in multi-processors runs
#globalFiles=.TRUE.,
readBinaryPrec=64,
writeBinaryPrec=64,
writeStatePrec=64,
&
# Elliptic solver parameters
&PARM02
cg2dMaxIters=1000,
cg2dTargetResidual=1.E-13,
cg3dMaxIters=400,
cg3dTargetResidual=1.E-13,
&
# Time stepping parameters
&PARM03
nIter0=0,
nTimeSteps=10000,
deltaT=0.01,
abEps=0.1,
pChkptFreq=2.5,
chkptFreq=2.5,
dumpFreq=0.4,
monitorFreq=2.5,
monitorSelect=2,
&
# Gridding parameters
&PARM04
usingCartesianGrid=.TRUE.,
delX=120*0.025,
delY=5,
delZ=30*0.03,
&
data.obcs is:
# Input datasets
&PARM05
# hydrogThetaFile='T.init',
bathyFile='topog.slope',
&
# Open-boundaries
&OBCS_PARM01
# OB_Ieast=-1,
OB_Iwest=1,
&
obcs_calc.F is:
C $Header: /u/gcmpack/MITgcm/verification/internal_wave/code/obcs_calc.F,v
1.9 2011/12/12 19:04:25 jmc Exp $
C $Name: checkpoint64m $
#include "OBCS_OPTIONS.h"
SUBROUTINE OBCS_CALC( futureTime, futureIter,
& uVel, vVel, wVel, theta, salt,
& myThid )
C *==========================================================*
C | SUBROUTINE OBCS_CALC
C | o Calculate future boundary data at open boundaries
C | at time = futureTime
C *==========================================================*
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#ifdef ALLOW_EXCH2
# include "W2_EXCH2_SIZE.h"
#endif /* ALLOW_EXCH2 */
#include "SET_GRID.h"
#include "GRID.h"
#include "OBCS_PARAMS.h"
#include "OBCS_GRID.h"
#include "OBCS_FIELDS.h"
#include "EOS.h"
C == Routine arguments ==
INTEGER futureIter
_RL futureTime
_RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
INTEGER myThid
#ifdef ALLOW_OBCS
C == Local variables ==
INTEGER bi, bj
INTEGER I, J ,K
_RL obTimeScale,Uinflow,rampTime2
_RL vertStructWst(Nr)
_RL mz,strat,kx
_RL tmpsum
C---+----1----+----2----+----3----+----4----+----5----+----
6----+----7-|--+----|
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_ENTER('OBCS_CALC',myThid)
#endif
C Vertical mode number
mz=1.0 _d 0
C Stratification
strat = 28.0 _d -2 / (gravity*tAlpha)
C Create a vertical structure function with zero mean
tmpsum=0.
do K=15,25
vertStructWst(K)=cos(mz*PI* (rC(K)/rF(Nr+1)) )
tmpsum=tmpsum+vertStructWst(K)*drF(K)
enddo
tmpsum=tmpsum/rF(Nr+1)
do K=15,25
vertStructWst(K)=vertStructWst(K)-tmpsum
enddo
c
obTimeScale = 127.57 _d 0
kx=mz*2. _d 0*pi/2.0 _d 0
& *sqrt((24232859 _d -10
& - 0)/(259727 _d -8
& - 64931.75 _d -8))
Uinflow = 0.0 _d 0
C *NOTE* I have commented out the ramp function below
C just to speed things up. You will probably want to use it
C for smoother looking solutions.
rampTime2 = 4. _d 0*142.29 _d 0
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
C Eastern OB
IF (useOrlanskiEast) THEN
CALL ORLANSKI_EAST(
& bi, bj, futureTime,
& uVel, vVel, wVel, theta, salt,
& myThid )
ELSE
DO K=15,25
DO J=1-Oly,sNy+Oly
OBEu(J,K,bi,bj)=0.
OBEv(J,K,bi,bj)=0.
OBEt(J,K,bi,bj)=tRef(K)
OBEs(J,K,bi,bj)=sRef(K)
#ifdef ALLOW_NONHYDROSTATIC
OBEw(J,K,bi,bj)=0.
#endif
ENDDO
ENDDO
ENDIF
C Western OB
IF (useOrlanskiWest) THEN
CALL ORLANSKI_WEST(
& bi, bj, futureTime,
& uVel, vVel, wVel, theta, salt,
& myThid )
ELSE
DO K=15,25
DO J=1-Oly,sNy+Oly
OBWu(J,K,bi,bj)=0. _d 0
c & +Uinflow
c & *vertStructWst(K)
c & *sin(2. _d 0*PI*futureTime/obTimeScale)
c & *(exp(futureTime/rampTime2)
c & - exp(-futureTime/rampTime2))
c & /(exp(futureTime/rampTime2)
c & + exp(-futureTime/rampTime2))
c & *cos(kx*(3. _d 0-2. _d 0-0.5 _d 0)*delX(1))
OBWv(J,K,bi,bj)=0. _d 0
c & +Uinflow
c & *f0/(2.0 _d 0*PI/obTimeScale)
c & *vertStructWst(K)
c & *cos(2. _d 0*PI*futureTime/obTimeScale )
c & * (exp(futureTime/rampTime2)
c & - exp(-futureTime/rampTime2))
c & /(exp(futureTime/rampTime2)
c & + exp(-futureTime/rampTime2))
OBWt(J,K,bi,bj)=tRef(K)
c & + Uinflow*sin(mz*PI*(float(k)-0.5 _d 0)/float(Nr))
c & * sin(2.0 _d 0*PI*futureTime/obTimeScale)
c & *sqrt(strat/(tAlpha*gravity))
c & *sqrt(2.0 _d 0*PI/obTimeScale*2.0*PI/obTimeScale - f0*f0)
c & /(2.0 _d 0*PI/obTimeScale)
c & * (exp(futureTime/rampTime2)
c & - exp(-futureTime/rampTime2))
c & /(exp(futureTime/rampTime2)
c & + exp(-futureTime/rampTime2))
OBWs(J,K,bi,bj)=sRef(K)
& + Uinflow*sin(mz*PI*(float(k)-0.5 _d 0)/float(Nr))
& * sin(2.0 _d 0*PI*futureTime/obTimeScale)
& *sqrt(strat/(tAlpha*gravity))
& *sqrt(2.0 _d 0*PI/obTimeScale*2.0*PI/obTimeScale - 0)
& /(2.0 _d 0*PI/obTimeScale)
c & * (exp(futureTime/rampTime2)
c & - exp(-futureTime/rampTime2))
c & /(exp(futureTime/rampTime2)
c & + exp(-futureTime/rampTime2))
#ifdef ALLOW_NONHYDROSTATIC
OBWw(J,K,bi,bj)=0. _d 0
c -Uinflow
c & *sqrt(2.0 _d 0*PI/obTimeScale*2.0 _d 0*PI/obTimeScale - 0)
c & /sqrt(strat*strat -
c & 2.0 _d 0*PI/obTimeScale*2.0 _d 0*PI/obTimeScale)
c & *sin(mz*PI*(float(k)-0.5 _d 0)/float(Nr))
c & *cos(2. _d 0*PI*futureTime/obTimeScale)
c & *(exp(futureTime/rampTime2)
c & - exp(-futureTime/rampTime2))
c & /(exp(futureTime/rampTime2)
c & + exp(-futureTime/rampTime2))
#endif
ENDDO
ENDDO
ENDIF
C Northern OB, template for forcing
IF (useOrlanskiNorth) THEN
CALL ORLANSKI_NORTH(
& bi, bj, futureTime,
& uVel, vVel, wVel, theta, salt,
& myThid )
ELSE
DO K=1,Nr
DO I=1-Olx,sNx+Olx
OBNv(I,K,bi,bj)=0.
OBNu(I,K,bi,bj)=0.
OBNt(I,K,bi,bj)=tRef(K)
OBNs(I,K,bi,bj)=sRef(K)
#ifdef ALLOW_NONHYDROSTATIC
OBNw(I,K,bi,bj)=0.
#endif
ENDDO
ENDDO
ENDIF
C Southern OB, template for forcing
IF (useOrlanskiSouth) THEN
CALL ORLANSKI_SOUTH(
& bi, bj, futureTime,
& uVel, vVel, wVel, theta, salt,
& myThid )
ELSE
DO K=1,Nr
DO I=1-Olx,sNx+Olx
OBSu(I,K,bi,bj)=0.
OBSv(I,K,bi,bj)=0.
OBSt(I,K,bi,bj)=tRef(K)
OBSs(I,K,bi,bj)=sRef(K)
#ifdef ALLOW_NONHYDROSTATIC
OBSw(I,K,bi,bj)=0.
#endif
ENDDO
ENDDO
ENDIF
C-- end bi,bj loops.
ENDDO
ENDDO
#ifdef ALLOW_OBCS_BALANCE
IF ( useOBCSbalance ) THEN
CALL OBCS_BALANCE_FLOW( futureTime, futureIter, myThid )
ENDIF
#endif /* ALLOW_OBCS_BALANCE */
#ifdef ALLOW_DEBUG
IF (debugMode) CALL DEBUG_LEAVE('OBCS_CALC',myThid)
#endif
#endif /* ALLOW_OBCS */
RETURN
END
size.h is:
C $Header: /u/gcmpack/MITgcm/verification/internal_wave/code/SIZE.h,v 1.6
2008/12/18 19:48:20 jmc Exp $
C $Name: checkpoint64m $
C
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 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 = 120,
& sNy = 1,
& OLx = 2,
& OLy = 2,
& nSx = 1,
& nSy = 1,
& nPx = 1,
& nPy = 1,
& Nx = sNx*nSx*nPx,
& Ny = sNy*nSy*nPy,
& Nr = 30)
C MAX_OLX - Set to the maximum overlap region size of any array
C MAX_OLY that will be exchanged. Controls the sizing of exch
C routine buufers.
INTEGER MAX_OLX
INTEGER MAX_OLY
PARAMETER ( MAX_OLX = OLx,
& MAX_OLY = OLy )
gendata.m is :
% This is a matlab script that generates the input data
% variable x resolution
prec='real*8';
ieee='b';
% Dimensions of grid
nx=120;
ny=5;
nz=30;
% Nominal depth of model (meters)
H=1.0;
% Size of domain
Lx=3;
% Horizontal resolution (m)
dx=zeros(nx,1);
for i=1:nx
dx(i) = Lx/(nx+1);
end
dy = Lx/nx
% Stratification
gravity=9.81;
talpha=2.0e-4;
N2=28e-2;
Tz=N2/(gravity*talpha);
dz=H/nz;
sprintf('delZ = %d * %7.6g,',nz,dz)
x=zeros(nx,1);
x(1) = dx(1);
for i=2:nx
x(i)=x(i-1) + dx(i);
end
z=-dz/2:-dz:-H;
%[Y,X]=meshgrid(y,x);
% Temperature profile
Tref=Tz*z-mean(Tz*z);
[sprintf('Tref =') sprintf(' %8.6g,',Tref)]
t=0.0*rand([nx,ny,nz]);
for k=1:nz
t(:,:,k) = t(:,:,k) + Tref(k);
end
fid=fopen('T.init','w',ieee); fwrite(fid,t,prec); fclose(fid);
% Sloping channel
slope=0.82;
offset=1.78;
dmax=0;
d=0.0*rand([nx,ny]);
for i=1:nx
for j=1:ny
d(i,j) = -H;
d(i,j) = d(i,j) + slope*(x(i) - offset);
if d(i,j) < -H ;
d(i,j) = -H;
end
if d(i,j) > dmax;
d(i,j) = dmax;
end
end
end
d(nx,:)=0.0;
fid=fopen('topog.slope','w',ieee); fwrite(fid,d,prec); fclose(fid);
plot(x,d(:,1))
fid=fopen('delXvar','w',ieee); fwrite(fid,dx,prec); fclose(fid);
% %convex slope
% nxslope=(dmax + H)/(slope)
% d1=zeros(nx,ny);
% hamp=(H-dmax)/5.0
% pi=4.0*atan(1.0)
% for i=1:nx
% for j=1:ny
% if x(i) < (offset + nxslope)
% if x(i) < offset
% d1(i,j) = d(i,j);
% else
% d1(i,j) = -H;
% d1(i,j) = d(i,j) + hamp*sin(pi*(x(i)-offset)/nxslope);
% if d1(i,j) < -H ;
% d1(i,j) = -H;
% end
% if d1(i,j) > dmax;
% d1(i,j) = dmax;
% end
% end
% else
% d1(i,j) = d(i,j);
% end
% end
% end
% %d1(end-1:end,:)=d1(1:2,:); % debug by aja
% fid=fopen('topog.convex','w',ieee); fwrite(fid,d1,prec); fclose(fid);
% plot(x,d1(:,1),'g')
%
% %convex slope
% d2=zeros(nx,ny);
% for i=1:nx
% for j=1:ny
% if x(i) < (offset + nxslope)
% if x(i) < offset
% d2(i,j) = d(i,j);
% else
% d2(i,j) = -H;
% d2(i,j) = d(i,j) - hamp*sin(pi*(x(i)-offset)/nxslope);
% if d2(i,j) < -H ;
% d2(i,j) = -H;
% end
% if d2(i,j) > dmax;
% d2(i,j) = dmax;
% end
% end
% else
% d2(i,j) = d(i,j);
% end
% end
% end
% %d2(end-1:end,:)=d2(1:2,:); % debug by aja
% fid=fopen('topog.concave','w',ieee); fwrite(fid,d2,prec); fclose(fid);
% hold on
% plot(x,d2(:,1),'r')
% hold off
%
%
% fid=fopen('delXvar','w',ieee); fwrite(fid,dx,prec); fclose(fid);
Thank you
_______________________________________________
MITgcm-support mailing list
MITgcm-support at mitgcm.org
http://mitgcm.org/mailman/listinfo/mitgcm-support
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20161219/e80257fc/attachment-0001.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 123.rar
Type: application/octet-stream
Size: 4578 bytes
Desc: not available
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20161219/e80257fc/attachment-0001.obj>
More information about the MITgcm-support
mailing list