[MITgcm-support] (no subject)
Abolfazl Shouryabi
ashouryabi at yahoo.com
Mon Dec 5 04:47:23 EST 2016
HelloMy 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 = 0.04842, 0.0433231, 0.0382263, 0.0331295, 0.0280326, 0.0229358, 0.0178389, 0.0127421, 0.00764526, 0.00254842, -0.00254842, -0.00764526, -0.0127421, -0.0178389, -0.0229358, -0.0280326, -0.0331295, -0.0382263, -0.0433231, -0.04842, sRef= 20*35., viscAz=1.E-3, viscAh=1.E-2, no_slip_sides=.FALSE., no_slip_bottom=.FALSE., 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=.FALSE., hFacMin=0.2, implicSurfPress=0.5, implicDiv2DFlow=0.5, nonlinFreeSurf=3, hFacInf=0.2, hFacSup=1.8, saltStepping=.FALSE.,#- 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=100, deltaT=500., abEps=0.1, pChkptFreq=0., chkptFreq=0., dumpFreq=50000., monitorFreq=2500., monitorSelect=2, &
# Gridding parameters &PARM04 usingCartesianGrid=.TRUE., delXfile='delXvar', delY=5.E3, delZ=20*10., &
# Input datasets &PARM05 hydrogThetaFile='T.init', bathyFile='topog.slope', &data.obcs is :# Open-boundaries &OBCS_PARM01# OB_Ieast=-1, OB_Iwest=1, &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 $CC /==========================================================\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 = 30, & sNy = 1, & OLx = 2, & OLy = 2, & nSx = 2, & nSy = 1, & nPx = 1, & nPy = 1, & Nx = sNx*nSx*nPx, & Ny = sNy*nSy*nPy, & Nr = 20)
C MAX_OLX - Set to the maximum overlap region size of any arrayC MAX_OLY that will be exchanged. Controls the sizing of exchC routine buufers. INTEGER MAX_OLX INTEGER MAX_OLY PARAMETER ( MAX_OLX = OLx, & MAX_OLY = OLy )
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_CALCC | o Calculate future boundary data at open boundariesC | at time = futureTimeC *==========================================================* 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 0C Stratification strat = 1.0 _d -6 / (gravity*tAlpha)
C Create a vertical structure function with zero mean tmpsum=0. do K=1,Nr vertStructWst(K)=cos(mz*PI* (rC(K)/rF(Nr+1)) ) tmpsum=tmpsum+vertStructWst(K)*drF(K) enddo tmpsum=tmpsum/rF(Nr+1) do K=1,Nr vertStructWst(K)=vertStructWst(K)-tmpsum enddoc obTimeScale = 44567.0 _d 0 kx=mz*2. _d 0*pi/400.0 _d 0 & *sqrt((2.0 _d 0*pi*2.0 _d 0*pi/(obTimeScale*obTimeScale) & - f0*f0)/(1.0 _d -6 & - 2.0 _d 0*pi*2.0 _d 0*pi/(obTimeScale*obTimeScale))) Uinflow = 0.024 _d 0C *NOTE* I have commented out the ramp function belowC just to speed things up. You will probably want to use itC for smoother looking solutions. rampTime2 = 4. _d 0*44567.0 _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=1,Nr 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=1,Nr DO J=1-Oly,sNy+Oly OBWu(J,K,bi,bj)=0. _d 0 & +Uinflow & *vertStructWst(K) & *sin(2. _d 0*PI*futureTime/obTimeScale)c & *(exp(futureTime/rampTime2)c & - exp(-futureTime/rampTime2))c & /(exp(futureTime/rampTime2)c & + exp(-futureTime/rampTime2)) & *cos(kx*(3. _d 0-2. _d 0-0.5 _d 0)*delX(1)) OBWv(J,K,bi,bj)=0. _d 0 & +Uinflow & *f0/(2.0 _d 0*PI/obTimeScale) & *vertStructWst(K) & *cos(2. _d 0*PI*futureTime/obTimeScale ) & * (exp(futureTime/rampTime2) & - exp(-futureTime/rampTime2)) & /(exp(futureTime/rampTime2) & + exp(-futureTime/rampTime2)) OBWt(J,K,bi,bj)=tRef(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 - f0*f0) & /(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)=-Uinflow & *sqrt(2.0 _d 0*PI/obTimeScale*2.0 _d 0*PI/obTimeScale - f0*f0) & /sqrt(strat*strat - & 2.0 _d 0*PI/obTimeScale*2.0 _d 0*PI/obTimeScale) & *sin(mz*PI*(float(k)-0.5 _d 0)/float(Nr)) & *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
gendata.m is :
% This is a matlab script that generates the input data% variable x resolutionprec='real*8';ieee='b';
% Dimensions of gridnx=60;ny=1;nz=20;% Nominal depth of model (meters)H=200.0;% Size of domainLx=13.3e3;
% Horizontal resolution (m)dx=zeros(nx,1);for i=1:nxdx(i) = Lx/(nx+1);end
dy = Lx/nx% Stratificationgravity=9.81;talpha=2.0e-4;N2=1e-6;
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:nxx(i)=x(i-1) + dx(i);endz=-dz/2:-dz:-H;
%[Y,X]=meshgrid(y,x);
% Temperature profileTref=Tz*z-mean(Tz*z);[sprintf('Tref =') sprintf(' %8.6g,',Tref)]t=0.0*rand([nx,ny,nz]);for k=1:nzt(:,:,k) = t(:,:,k) + Tref(k);endfid=fopen('T.init','w',ieee); fwrite(fid,t,prec); fclose(fid);
% Sloping channel slope=0.03offset=2.5e3;dmax=-40.0;d=0.0*rand([nx,ny]);for i=1:nxfor j=1:nyd(i,j) = -H;d(i,j) = d(i,j) + slope*(x(i) - offset);if d(i,j) < -H ;d(i,j) = -H;endif d(i,j) > dmax;d(i,j) = dmax;endendendd(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 slopenxslope=(dmax + H)/(slope)d1=zeros(nx,ny);hamp=(H-dmax)/5.0pi=4.0*atan(1.0)for i=1:nxfor j=1:nyif x(i) < (offset + nxslope)if x(i) < offsetd1(i,j) = d(i,j);elsed1(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;endif d1(i,j) > dmax;d1(i,j) = dmax;endendelsed1(i,j) = d(i,j);endendend%d1(end-1:end,:)=d1(1:2,:); % debug by ajafid=fopen('topog.convex','w',ieee); fwrite(fid,d1,prec); fclose(fid);plot(x,d1(:,1),'g')
%convex sloped2=zeros(nx,ny);for i=1:nxfor j=1:nyif x(i) < (offset + nxslope)if x(i) < offsetd2(i,j) = d(i,j);elsed2(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;endif d2(i,j) > dmax;d2(i,j) = dmax;endendelsed2(i,j) = d(i,j);endendend%d2(end-1:end,:)=d2(1:2,:); % debug by ajafid=fopen('topog.concave','w',ieee); fwrite(fid,d2,prec); fclose(fid);hold onplot(x,d2(:,1),'r')hold off
fid=fopen('delXvar','w',ieee); fwrite(fid,dx,prec); fclose(fid);
Abolfazl Shouryabi
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20161205/99aa6494/attachment-0001.htm>
More information about the MITgcm-support
mailing list