[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