/[MITgcm]/MITgcm/verification/internal_wave/code/obcs_calc.F
ViewVC logotype

Contents of /MITgcm/verification/internal_wave/code/obcs_calc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.9 - (show annotations) (download)
Mon Dec 12 19:04:25 2011 UTC (8 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint64, checkpoint65, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63g, checkpoint65o, checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, HEAD
Changes since 1.8: +5 -1 lines
move delX,delY to new header file (SET_GRID.h) and adjust length to 1rst
face dimensions.

1 C $Header: /u/gcmpack/MITgcm/verification/internal_wave/code/obcs_calc.F,v 1.8 2011/05/24 20:34:01 jmc Exp $
2 C $Name: $
3
4 #include "OBCS_OPTIONS.h"
5
6 SUBROUTINE OBCS_CALC( futureTime, futureIter,
7 & uVel, vVel, wVel, theta, salt,
8 & myThid )
9 C *==========================================================*
10 C | SUBROUTINE OBCS_CALC
11 C | o Calculate future boundary data at open boundaries
12 C | at time = futureTime
13 C *==========================================================*
14 IMPLICIT NONE
15
16 C === Global variables ===
17 #include "SIZE.h"
18 #include "EEPARAMS.h"
19 #include "PARAMS.h"
20 #ifdef ALLOW_EXCH2
21 # include "W2_EXCH2_SIZE.h"
22 #endif /* ALLOW_EXCH2 */
23 #include "SET_GRID.h"
24 #include "GRID.h"
25 #include "OBCS_PARAMS.h"
26 #include "OBCS_GRID.h"
27 #include "OBCS_FIELDS.h"
28 #include "EOS.h"
29
30 C == Routine arguments ==
31 INTEGER futureIter
32 _RL futureTime
33 _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
34 _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
35 _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
36 _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
37 _RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
38 INTEGER myThid
39
40 #ifdef ALLOW_OBCS
41
42 C == Local variables ==
43 INTEGER bi, bj
44 INTEGER I, J ,K
45 _RL obTimeScale,Uinflow,rampTime2
46 _RL vertStructWst(Nr)
47 _RL mz,strat,kx
48 _RL tmpsum
49
50 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
51
52 #ifdef ALLOW_DEBUG
53 IF (debugMode) CALL DEBUG_ENTER('OBCS_CALC',myThid)
54 #endif
55
56 C Vertical mode number
57 mz=1.0 _d 0
58 C Stratification
59 strat = 1.0 _d -6 / (gravity*tAlpha)
60
61 C Create a vertical structure function with zero mean
62 tmpsum=0.
63 do K=1,Nr
64 vertStructWst(K)=cos(mz*PI* (rC(K)/rF(Nr+1)) )
65 tmpsum=tmpsum+vertStructWst(K)*drF(K)
66 enddo
67 tmpsum=tmpsum/rF(Nr+1)
68 do K=1,Nr
69 vertStructWst(K)=vertStructWst(K)-tmpsum
70 enddo
71 c
72 obTimeScale = 44567.0 _d 0
73 kx=mz*2. _d 0*pi/400.0 _d 0
74 & *sqrt((2.0 _d 0*pi*2.0 _d 0*pi/(obTimeScale*obTimeScale)
75 & - f0*f0)/(1.0 _d -6
76 & - 2.0 _d 0*pi*2.0 _d 0*pi/(obTimeScale*obTimeScale)))
77 Uinflow = 0.024 _d 0
78 C *NOTE* I have commented out the ramp function below
79 C just to speed things up. You will probably want to use it
80 C for smoother looking solutions.
81 rampTime2 = 4. _d 0*44567.0 _d 0
82
83 DO bj=myByLo(myThid),myByHi(myThid)
84 DO bi=myBxLo(myThid),myBxHi(myThid)
85
86 C Eastern OB
87 IF (useOrlanskiEast) THEN
88 CALL ORLANSKI_EAST(
89 & bi, bj, futureTime,
90 & uVel, vVel, wVel, theta, salt,
91 & myThid )
92 ELSE
93 DO K=1,Nr
94 DO J=1-Oly,sNy+Oly
95 OBEu(J,K,bi,bj)=0.
96 OBEv(J,K,bi,bj)=0.
97 OBEt(J,K,bi,bj)=tRef(K)
98 OBEs(J,K,bi,bj)=sRef(K)
99 #ifdef ALLOW_NONHYDROSTATIC
100 OBEw(J,K,bi,bj)=0.
101 #endif
102 ENDDO
103 ENDDO
104 ENDIF
105
106 C Western OB
107 IF (useOrlanskiWest) THEN
108 CALL ORLANSKI_WEST(
109 & bi, bj, futureTime,
110 & uVel, vVel, wVel, theta, salt,
111 & myThid )
112 ELSE
113 DO K=1,Nr
114 DO J=1-Oly,sNy+Oly
115 OBWu(J,K,bi,bj)=0. _d 0
116 & +Uinflow
117 & *vertStructWst(K)
118 & *sin(2. _d 0*PI*futureTime/obTimeScale)
119 c & *(exp(futureTime/rampTime2)
120 c & - exp(-futureTime/rampTime2))
121 c & /(exp(futureTime/rampTime2)
122 c & + exp(-futureTime/rampTime2))
123 & *cos(kx*(3. _d 0-2. _d 0-0.5 _d 0)*delX(1))
124 OBWv(J,K,bi,bj)=0. _d 0
125 & +Uinflow
126 & *f0/(2.0 _d 0*PI/obTimeScale)
127 & *vertStructWst(K)
128 & *cos(2. _d 0*PI*futureTime/obTimeScale )
129 & * (exp(futureTime/rampTime2)
130 & - exp(-futureTime/rampTime2))
131 & /(exp(futureTime/rampTime2)
132 & + exp(-futureTime/rampTime2))
133 OBWt(J,K,bi,bj)=tRef(K)
134 & + Uinflow*sin(mz*PI*(float(k)-0.5 _d 0)/float(Nr))
135 & * sin(2.0 _d 0*PI*futureTime/obTimeScale)
136 & *sqrt(strat/(tAlpha*gravity))
137 & *sqrt(2.0 _d 0*PI/obTimeScale*2.0*PI/obTimeScale - f0*f0)
138 & /(2.0 _d 0*PI/obTimeScale)
139 c & * (exp(futureTime/rampTime2)
140 c & - exp(-futureTime/rampTime2))
141 c & /(exp(futureTime/rampTime2)
142 c & + exp(-futureTime/rampTime2))
143 #ifdef ALLOW_NONHYDROSTATIC
144 OBWw(J,K,bi,bj)=-Uinflow
145 & *sqrt(2.0 _d 0*PI/obTimeScale*2.0 _d 0*PI/obTimeScale - f0*f0)
146 & /sqrt(strat*strat -
147 & 2.0 _d 0*PI/obTimeScale*2.0 _d 0*PI/obTimeScale)
148 & *sin(mz*PI*(float(k)-0.5 _d 0)/float(Nr))
149 & *cos(2. _d 0*PI*futureTime/obTimeScale)
150 c & *(exp(futureTime/rampTime2)
151 c & - exp(-futureTime/rampTime2))
152 c & /(exp(futureTime/rampTime2)
153 c & + exp(-futureTime/rampTime2))
154 #endif
155 ENDDO
156 ENDDO
157 ENDIF
158
159 C Northern OB, template for forcing
160 IF (useOrlanskiNorth) THEN
161 CALL ORLANSKI_NORTH(
162 & bi, bj, futureTime,
163 & uVel, vVel, wVel, theta, salt,
164 & myThid )
165 ELSE
166 DO K=1,Nr
167 DO I=1-Olx,sNx+Olx
168 OBNv(I,K,bi,bj)=0.
169 OBNu(I,K,bi,bj)=0.
170 OBNt(I,K,bi,bj)=tRef(K)
171 OBNs(I,K,bi,bj)=sRef(K)
172 #ifdef ALLOW_NONHYDROSTATIC
173 OBNw(I,K,bi,bj)=0.
174 #endif
175 ENDDO
176 ENDDO
177 ENDIF
178
179 C Southern OB, template for forcing
180 IF (useOrlanskiSouth) THEN
181 CALL ORLANSKI_SOUTH(
182 & bi, bj, futureTime,
183 & uVel, vVel, wVel, theta, salt,
184 & myThid )
185 ELSE
186 DO K=1,Nr
187 DO I=1-Olx,sNx+Olx
188 OBSu(I,K,bi,bj)=0.
189 OBSv(I,K,bi,bj)=0.
190 OBSt(I,K,bi,bj)=tRef(K)
191 OBSs(I,K,bi,bj)=sRef(K)
192 #ifdef ALLOW_NONHYDROSTATIC
193 OBSw(I,K,bi,bj)=0.
194 #endif
195 ENDDO
196 ENDDO
197 ENDIF
198
199 C-- end bi,bj loops.
200 ENDDO
201 ENDDO
202
203 #ifdef ALLOW_OBCS_BALANCE
204 IF ( useOBCSbalance ) THEN
205 CALL OBCS_BALANCE_FLOW( futureTime, futureIter, myThid )
206 ENDIF
207 #endif /* ALLOW_OBCS_BALANCE */
208
209 #ifdef ALLOW_DEBUG
210 IF (debugMode) CALL DEBUG_LEAVE('OBCS_CALC',myThid)
211 #endif
212 #endif /* ALLOW_OBCS */
213
214 RETURN
215 END

  ViewVC Help
Powered by ViewVC 1.1.22