[MITgcm-support] (no subject)

Abolfazl Shouryabi ashouryabi at yahoo.com
Mon Dec 5 05:18:24 EST 2016


HelloMy problem is that the next time there is an increase in salinityPlease 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_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 = 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      enddoc      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 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*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 0c     &       +Uinflowc     &       *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 0c     &       +Uinflowc     &       *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 0c       -Uinflowc     & *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      ENDsize.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 =  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 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 )
gendata.m is :

% This is a matlab script that generates the input data% variable x resolutionprec='real*8';ieee='b';

% Dimensions of gridnx=120;ny=5;nz=30;% Nominal depth of model (meters)H=1.0;% Size of domainLx=3;
% 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=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: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.82;offset=1.78;dmax=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 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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20161205/546d12e0/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/20161205/546d12e0/attachment-0001.obj>


More information about the MITgcm-support mailing list