[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