/[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.6 - (show annotations) (download)
Tue Dec 15 17:03:29 2009 UTC (14 years, 4 months ago) by jahn
Branch: MAIN
CVS Tags: checkpoint62c, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62, checkpoint62b
Changes since 1.5: +9 -3 lines
move bi,bj loops into obcs_calc, so obcs_prescribe_read is called only once

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

  ViewVC Help
Powered by ViewVC 1.1.22