[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