[MITgcm-support] (no subject)
Martin Losch
martin.losch at awi.de
Mon Dec 19 03:30:25 EST 2016
Hi Mohammad,
I don’t think that you can expect any reply to your question, because it is totally unclear what you want to know. Please describe your configuration (what you are trying to do), and what happens. I assume that “next time” means “next time step”, but it’s not clear how many time steps you can run with your configuration, etc.
There’s no buoyancy gradient in your configuration and no external forcing.
I don't think anyone will want to spend time on trying to understand what you are trying to do just from a subset of your configuration files, please be more specific and someone may be able to help.
Martin
> On 19 Dec 2016, at 07:21, Mohammad Akbarinasab <akbarinasabmohamad at gmail.com> wrote:
>
>
>
>
>
> 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
>
>
> <123.rar>_______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
More information about the MITgcm-support
mailing list