1 |
jmc |
1.4 |
C$Header: /u/gcmpack/MITgcm/pkg/opps/opps_calc.F,v 1.3 2004/12/04 00:17:18 jmc Exp $ |
2 |
edhill |
1.2 |
C$Name: $ |
3 |
mlosch |
1.1 |
|
4 |
|
|
#include "OPPS_OPTIONS.h" |
5 |
|
|
|
6 |
|
|
CBOP |
7 |
|
|
C !ROUTINE: OPPS_CALC |
8 |
|
|
|
9 |
|
|
C !INTERFACE: ====================================================== |
10 |
|
|
subroutine OPPS_CALC( |
11 |
|
|
U tracerEnv, |
12 |
|
|
I wVel, |
13 |
|
|
I kMax, nTracer, nTracerInuse, |
14 |
|
|
I I, J, bi, bj, myTime, myIter, myThid ) |
15 |
|
|
|
16 |
|
|
C !DESCRIPTION: \bv |
17 |
|
|
C /=====================================================================\ |
18 |
|
|
C | SUBROUTINE OPPS_CALC | |
19 |
|
|
C | o Compute all OPPS fields defined in OPPS.h | |
20 |
|
|
C |=====================================================================| |
21 |
|
|
C | This subroutine is based on the routine 3dconvection.F | |
22 |
|
|
C | by E. Skyllingstad (?) | |
23 |
|
|
C | plenty of modifications to make it work: | |
24 |
|
|
C | - removed many unused parameters and variables | |
25 |
|
|
C | - turned everything (back) into 1D code | |
26 |
|
|
C | - pass variables, that are orginially in common blocks: | |
27 |
|
|
C | maxDepth | |
28 |
|
|
C | - pass vertical velocity, set in OPPS_INTERFACE | |
29 |
|
|
C | - do not use convadj for now (whatever that is) | |
30 |
|
|
C | - changed two .LT. 0 to .LE. 0 statements (because of possible | |
31 |
|
|
C | division) | |
32 |
|
|
C | - replaced statement function state1 by call to a real function | |
33 |
|
|
C | - removed range check, actually moved it up to OPPS_INTERFACE | |
34 |
|
|
C | - avoid division by zero: if (Wd.EQ.0) dt = ...1/Wd | |
35 |
|
|
C | - cleaned-up debugging | |
36 |
|
|
C | - replaced local dz and GridThickness by global drF | |
37 |
|
|
C | - replaced 1/dz by 1*recip_drF | |
38 |
|
|
C | - replaced 9.81 with gravity (=9.81) | |
39 |
|
|
C | - added a lot of comments that relate code to equation in paper | |
40 |
|
|
C | (Paluszkiewicz+Romea, 1997, Dynamics of Atmospheres and Oceans, | |
41 |
|
|
C | 26, pp. 95-130) | |
42 |
|
|
C | - included passive tracer support. This is the main change and may | |
43 |
|
|
C | not improve the readability of the code because of the joint | |
44 |
|
|
C | treatment of active (theta, salt) and passive tracers. The array | |
45 |
|
|
C | tracerEnv(Nr,2+PTRACERS_num) contains | |
46 |
|
|
C | theta = tracerEnv(:,1), | |
47 |
|
|
C | salt = tracerEnv(:,2), and | |
48 |
|
|
C | ptracers = tracerEnv(:,3:PTRACERS_num+2). | |
49 |
|
|
C | All related array names have been changed accordingly, so that | |
50 |
|
|
C | instead of Sd(Nr) and Td(Nr) (plume salinity and temperature), we | |
51 |
|
|
C | have Pd(Nr,nTracer) (tracer in plume), with Sd(:) = Pd(:,2), | |
52 |
|
|
C | Td(:) = Pd(:,1), etc. | |
53 |
|
|
C | o TODO: | |
54 |
|
|
C | clean up the logic of the vertical loops and get rid off the | |
55 |
|
|
C | GOTO statements | |
56 |
|
|
C \=====================================================================/ |
57 |
|
|
IMPLICIT NONE |
58 |
|
|
C |
59 |
|
|
C-------------------------------------------------------------------- |
60 |
|
|
|
61 |
|
|
C \ev |
62 |
|
|
|
63 |
|
|
C !USES: ============================================================ |
64 |
|
|
#include "SIZE.h" |
65 |
|
|
#include "EEPARAMS.h" |
66 |
|
|
#include "PARAMS.h" |
67 |
|
|
#include "OPPS.h" |
68 |
|
|
#include "FFIELDS.h" |
69 |
|
|
#include "GRID.h" |
70 |
|
|
|
71 |
jmc |
1.4 |
EXTERNAL DIFF_BASE_MULTIPLE |
72 |
|
|
LOGICAL DIFF_BASE_MULTIPLE |
73 |
mlosch |
1.1 |
|
74 |
|
|
C !INPUT PARAMETERS: =================================================== |
75 |
|
|
c Routine arguments |
76 |
|
|
c bi, bj - array indices on which to apply calculations |
77 |
|
|
c myTime - Current time in simulation |
78 |
|
|
|
79 |
|
|
INTEGER I, J, bi, bj, KMax, nTracer, nTracerInUse |
80 |
|
|
INTEGER myThid, myIter |
81 |
|
|
_RL myTime |
82 |
|
|
_RL tracerEnv(Nr,nTracer),wVel(Nr) |
83 |
|
|
|
84 |
|
|
#ifdef ALLOW_OPPS |
85 |
|
|
C !LOCAL VARIABLES: ==================================================== |
86 |
|
|
c Local constants |
87 |
|
|
C imin, imax, jmin, jmax - array computation indices |
88 |
|
|
C msgBuf - Informational/error meesage buffer |
89 |
|
|
CHARACTER*(MAX_LEN_MBUF) msgBuf |
90 |
|
|
INTEGER K, K2, K2m1, K2p1, ktr |
91 |
|
|
INTEGER ntime,nn,kmx,ic |
92 |
|
|
INTEGER maxDepth |
93 |
|
|
|
94 |
|
|
_RL wsqr,oldflux,newflux,entrainrate |
95 |
|
|
_RL pmix |
96 |
|
|
_RL D1,D2,state1 |
97 |
|
|
_RL dz1,dz2 |
98 |
|
|
_RL radius,StartingFlux |
99 |
|
|
_RL dtts,dt |
100 |
|
|
C Arrays |
101 |
|
|
_RL Paa(Nr,nTracer) |
102 |
|
|
_RL wda(Nr), mda(Nr), pda(Nr,nTracer) |
103 |
|
|
C |
104 |
|
|
C Pd, Wd - tracers, vertical velocity in plume |
105 |
|
|
C Md - plume mass flux (?) |
106 |
|
|
C Ad - fractional area covered by plume |
107 |
|
|
C Dd - density in plume |
108 |
|
|
C De - density of environment |
109 |
|
|
C PlumeEntrainment - |
110 |
|
|
_RL Ad(Nr),Wd(Nr),Dd(Nr),Md(Nr) |
111 |
|
|
_RL De(Nr) |
112 |
|
|
_RL PlumeEntrainment(Nr) |
113 |
|
|
_RL Pd(Nr,nTracer) |
114 |
|
|
CEOP |
115 |
|
|
|
116 |
|
|
|
117 |
|
|
C-- Check to see if should convect now |
118 |
jmc |
1.4 |
C IF ( DIFF_BASE_MULTIPLE(baseTime,cAdjFreq,myTime,deltaTClock) ) THEN |
119 |
mlosch |
1.1 |
IF ( .true. ) THEN |
120 |
|
|
C local initialization |
121 |
|
|
|
122 |
|
|
C Copy some arrays |
123 |
jmc |
1.3 |
dtts = dTtracerLev(1) |
124 |
mlosch |
1.1 |
C |
125 |
|
|
C start k-loop |
126 |
|
|
C |
127 |
|
|
|
128 |
|
|
DO k=1,KMax-1 |
129 |
|
|
c |
130 |
|
|
c initialize the plume T,S,density, and w velocity |
131 |
|
|
c |
132 |
|
|
DO ktr=1,nTracerInUse |
133 |
|
|
Pd(k,ktr) = tracerEnv(k,ktr) |
134 |
|
|
ENDDO |
135 |
|
|
Dd(k)=state1(Pd(k,2),Pd(k,1),i,j,k,bi,bj,myThid) |
136 |
|
|
De(k)=Dd(k) |
137 |
|
|
CML print *, 'ml-opps:', i,j,k,tracerEnv(k,2),tracerEnv(k,1), |
138 |
|
|
CML & Dd(k),Pd(k,1),Pd(k,2) |
139 |
|
|
CML compute vertical velocity at cell centers from GCM velocity |
140 |
|
|
Wd(k)= - .5*(wVel(K)+wVel(K+1)) |
141 |
|
|
CML( |
142 |
|
|
CML avoid division by zero |
143 |
|
|
CML IF (Wd(K) .EQ. 0.D0) Wd(K) = 2.23e-16 |
144 |
|
|
CML) |
145 |
|
|
c |
146 |
|
|
c guess at initial top grid cell vertical velocity |
147 |
|
|
c |
148 |
|
|
CML Wd(k) = 0.03 |
149 |
|
|
c |
150 |
|
|
c these estimates of initial plume velocity based on plume size and |
151 |
|
|
c top grid cell water mass |
152 |
|
|
c |
153 |
|
|
c Wd(k) = 0.5*drF(k)/(dtts*FRACTIONAL_AREA) |
154 |
|
|
c Wd(k) = 0.5*drF(k)/dtts |
155 |
|
|
c |
156 |
|
|
wsqr=Wd(k)*Wd(k) |
157 |
|
|
PlumeEntrainment(k) = 0.0 |
158 |
|
|
c |
159 |
|
|
c |
160 |
|
|
c |
161 |
|
|
#ifdef ALLOW_OPPS_DEBUG |
162 |
|
|
IF ( OPPSdebugLevel.GE.debLevB ) THEN |
163 |
|
|
WRITE(msgBuf,'(A,I3)') |
164 |
|
|
& 'S/R OPPS_CALC: doing old lowerparcel', k |
165 |
|
|
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
166 |
|
|
& SQUEEZE_RIGHT , 1) |
167 |
|
|
ENDIF |
168 |
|
|
#endif /* ALLOW_OPPS_DEBUG */ |
169 |
|
|
radius=PlumeRadius |
170 |
|
|
StartingFlux=radius*radius*Wd(k)*Dd(k) |
171 |
|
|
oldflux=StartingFlux |
172 |
|
|
|
173 |
|
|
dz2=DrF(k) |
174 |
|
|
DO k2=k,KMax-1 |
175 |
|
|
D1=state1( Pd(k2,2), Pd(k2,1),i,j,k2+1,bi,bj,myThid) |
176 |
|
|
D2=state1( tracerEnv(k2+1,2), tracerEnv(k2+1,1), |
177 |
|
|
& i,j,k2+1,bi,bj,myThid) |
178 |
|
|
De(k2+1)=D2 |
179 |
|
|
c |
180 |
|
|
c To start downward, parcel has to initially be heavier than environment |
181 |
|
|
c but after it has started moving, we continue plume until plume tke or |
182 |
|
|
c flux goes negative |
183 |
|
|
c |
184 |
|
|
CML & _hFacC(i,j,k-1,bi,bj) |
185 |
|
|
CML & *_hFacC(i,j,k,bi,bj) .GT. 0. |
186 |
|
|
CML & .AND. |
187 |
|
|
IF (D2-D1 .LT. STABILITY_THRESHOLD.or.k2.ne.k) THEN |
188 |
|
|
dz1=dz2 |
189 |
|
|
dz2=DrF(k2+1) |
190 |
|
|
c |
191 |
|
|
C find mass flux according to eq.(3) from paper by vertical integration |
192 |
|
|
c |
193 |
|
|
newflux=oldflux+e2*radius*Wd(k2)*Dd(k2)* |
194 |
|
|
& .5*(dz1+dz2) |
195 |
|
|
CML print *, 'ml-opps:', i,j,k,oldflux,newflux,e2,radius, |
196 |
|
|
CML & Wd(k2),Dd(k2),Pd(k2,1),Pd(k2,2),dz1,dz2 |
197 |
|
|
c |
198 |
|
|
PlumeEntrainment(k2+1) = newflux/StartingFlux |
199 |
|
|
c |
200 |
|
|
IF(newflux.LE.0.0) then |
201 |
|
|
#ifdef ALLOW_OPPS_DEBUG |
202 |
|
|
IF ( OPPSdebugLevel.GE.debLevA ) THEN |
203 |
|
|
WRITE(msgBuf,'(A,I3)') |
204 |
|
|
& 'S/R OPPS_CALC: Plume entrained to zero at level ', k2 |
205 |
|
|
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
206 |
|
|
& SQUEEZE_RIGHT , 1) |
207 |
|
|
ENDIF |
208 |
|
|
#endif /* ALLOW_OPPS_DEBUG */ |
209 |
|
|
maxdepth = k2 |
210 |
|
|
if(maxdepth.eq.k) goto 1000 |
211 |
|
|
goto 1 |
212 |
|
|
endif |
213 |
|
|
c |
214 |
|
|
c entrainment rate is basically a scaled mass flux dM/M |
215 |
|
|
c |
216 |
|
|
entrainrate = (newflux - oldflux)/newflux |
217 |
|
|
oldflux = newflux |
218 |
|
|
c |
219 |
|
|
c |
220 |
|
|
c mix var's are the average environmental values over the two grid levels |
221 |
|
|
c |
222 |
|
|
DO ktr=1,nTracerInUse |
223 |
|
|
pmix=(dz1*tracerEnv(k2,ktr)+dz2*tracerEnv(k2+1,ktr)) |
224 |
|
|
& /(dz1+dz2) |
225 |
|
|
Pd(k2+1,ktr)=Pd(k2,ktr) |
226 |
|
|
& - entrainrate*(pmix - Pd(k2,ktr)) |
227 |
|
|
ENDDO |
228 |
|
|
c |
229 |
|
|
c compute the density at this level for the buoyancy term in the |
230 |
|
|
c vertical k.e. equation |
231 |
|
|
c |
232 |
|
|
Dd(k2+1)=state1(Pd(k2+1,2),Pd(k2+1,1),i,j,k2+1,bi,bj,myThid) |
233 |
|
|
c |
234 |
|
|
c next, solve for the vertical velocity k.e. using combined eq. (4) |
235 |
|
|
c and eq (5) from the paper |
236 |
|
|
c |
237 |
|
|
#ifdef ALLOW_OPPS_DEBUG |
238 |
|
|
IF ( OPPSdebugLevel.GE.debLevA ) THEN |
239 |
|
|
WRITE(msgBuf,'(A,3E12.4,I3)') |
240 |
|
|
& 'S/R OPPS_CALC: Dd,De,entr,k ',Dd(k2),De(k2),entrainrate,k2 |
241 |
|
|
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
242 |
|
|
& SQUEEZE_RIGHT , 1) |
243 |
|
|
ENDIF |
244 |
|
|
#endif /* ALLOW_OPPS_DEBUG */ |
245 |
|
|
CML insert Eq. (4) into Eq. (5) to get something like this for wp^2 |
246 |
|
|
wsqr = wsqr - wsqr*abs(entrainrate)+ gravity* |
247 |
|
|
& (dz1*(Dd(k2)-De(k2))/De(k2) |
248 |
|
|
& +dz2*(Dd(k2+1)-De(k2+1))/De(k2+1)) |
249 |
|
|
c |
250 |
|
|
c if negative k.e. then plume has reached max depth, get out of loop |
251 |
|
|
c |
252 |
|
|
IF(wsqr.LE.0.0)then |
253 |
|
|
maxdepth = k2 |
254 |
|
|
#ifdef ALLOW_OPPS_DEBUG |
255 |
|
|
IF ( OPPSdebugLevel.GE.debLevA ) THEN |
256 |
|
|
WRITE(msgBuf,'(A,I3)') |
257 |
|
|
& 'S/R OPPS_CALC: Plume velocity went to zero at level ', k2 |
258 |
|
|
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
259 |
|
|
& SQUEEZE_RIGHT , 1) |
260 |
|
|
WRITE(msgBuf,'(A,4A14)') |
261 |
|
|
& 'S/R OPPS_CALC: ', 'wsqr', 'entrainrate', |
262 |
|
|
& '(Dd-De)/De up', '(Dd-De)/De do' |
263 |
|
|
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
264 |
|
|
& SQUEEZE_RIGHT , 1) |
265 |
|
|
WRITE(msgBuf,'(A,4E14.6)') |
266 |
|
|
& 'S/R OPPS_CALC: ', wsqr, entrainrate, |
267 |
|
|
& (Dd(k2)-De(k2))/De(k2), (Dd(k2+1)-De(k2+1))/De(k2+1) |
268 |
|
|
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
269 |
|
|
& SQUEEZE_RIGHT , 1) |
270 |
|
|
ENDIF |
271 |
|
|
#endif /* ALLOW_OPPS_DEBUG */ |
272 |
|
|
if(maxdepth.eq.k) goto 1000 |
273 |
|
|
goto 1 |
274 |
|
|
endif |
275 |
|
|
Wd(k2+1)=sqrt(wsqr) |
276 |
|
|
C |
277 |
|
|
C compute a new radius based on the new mass flux at this grid level |
278 |
|
|
C from Eq. (4) |
279 |
|
|
C |
280 |
|
|
radius=sqrt(newflux/(Wd(k2)*Dd(k2))) |
281 |
|
|
ELSE |
282 |
|
|
maxdepth=k2 |
283 |
|
|
if(maxdepth.eq.k) goto 1000 |
284 |
|
|
GOTO 1 |
285 |
|
|
ENDIF |
286 |
|
|
ENDDO |
287 |
|
|
c |
288 |
|
|
c plume has reached the bottom |
289 |
|
|
c |
290 |
|
|
MaxDepth=kMax |
291 |
|
|
c |
292 |
|
|
1 CONTINUE |
293 |
|
|
c |
294 |
|
|
Ad(k)=FRACTIONAL_AREA |
295 |
|
|
IC=0 |
296 |
|
|
c |
297 |
|
|
c start iteration on fractional area, not used in OGCM implementation |
298 |
|
|
c |
299 |
|
|
c |
300 |
|
|
DO IC=1,Max_ABE_Iterations |
301 |
|
|
c |
302 |
|
|
c |
303 |
|
|
c next compute the mass flux beteen each grid box using the entrainment |
304 |
|
|
c |
305 |
|
|
Md(k)=Wd(k)*Ad(k) |
306 |
|
|
c |
307 |
|
|
DO k2=k+1,maxDepth |
308 |
|
|
Md(k2)=Md(k)*PlumeEntrainment(k2) |
309 |
|
|
#ifdef ALLOW_OPPS_DEBUG |
310 |
|
|
IF ( OPPSdebugLevel.GE.debLevA ) THEN |
311 |
|
|
WRITE(msgBuf,'(A,2E12.4,I3)') |
312 |
|
|
& 'S/R OPPS_CALC: Md, Wd, and k are ',Md(k2),Wd(k2),k2 |
313 |
|
|
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
314 |
|
|
& SQUEEZE_RIGHT , 1) |
315 |
|
|
ENDIF |
316 |
|
|
#endif /* ALLOW_OPPS_DEBUG */ |
317 |
|
|
ENDDO |
318 |
|
|
c |
319 |
|
|
c Now move on to calculate new temperature using flux from |
320 |
|
|
c Td, Sd, Wd, ta, sa, and we. Values for these variables are at |
321 |
|
|
c center of grid cell, use weighted average to get boundary values |
322 |
|
|
c |
323 |
|
|
c use a timestep limited by the GCM model timestep and the maximum plume |
324 |
|
|
c velocity (CFL criteria) |
325 |
|
|
c |
326 |
|
|
c |
327 |
|
|
c calculate the weighted wd, td, and sd |
328 |
|
|
c |
329 |
|
|
dt = dtts |
330 |
|
|
do k2=k,maxDepth-1 |
331 |
|
|
IF ( Wd(K2) .NE. 0. _d 0 ) dt = min(dt,drF(k2)/Wd(k2)) |
332 |
|
|
c |
333 |
|
|
c time integration will be integer number of steps to get one |
334 |
|
|
c gcm time step |
335 |
|
|
c |
336 |
|
|
ntime = nint(0.5*int(dtts/dt)) |
337 |
|
|
if(ntime.eq.0) then |
338 |
|
|
ntime = 1 |
339 |
|
|
endif |
340 |
|
|
c |
341 |
|
|
c make sure area weighted vertical velocities match; in other words |
342 |
|
|
c make sure mass in equals mass out at the intersection of each grid |
343 |
|
|
c cell. Eq. (20) |
344 |
|
|
c |
345 |
|
|
mda(k2) = (md(k2)*drF(k2)+md(k2+1)*drF(k2+1))/ |
346 |
|
|
* (drF(k2)+drF(k2+1)) |
347 |
|
|
c |
348 |
|
|
wda(k2) = (wd(k2)*drF(k2)+wd(k2+1)*drF(k2+1))/ |
349 |
|
|
* (drF(k2)+drF(k2+1)) |
350 |
|
|
c |
351 |
|
|
DO ktr = 1, nTracerInUse |
352 |
|
|
Pda(k2,ktr) = Pd(k2,ktr) |
353 |
|
|
Paa(k2,ktr) = tracerEnv(k2+1,ktr) |
354 |
|
|
ENDDO |
355 |
|
|
c |
356 |
|
|
enddo |
357 |
|
|
dt = min(dt,dtts) |
358 |
|
|
#ifdef ALLOW_OPPS_DEBUG |
359 |
|
|
IF ( OPPSdebugLevel.GE.debLevA ) THEN |
360 |
|
|
WRITE(msgBuf,'(A,F14.4)') |
361 |
|
|
& 'S/R OPPS_CALC: time step = ', dt |
362 |
|
|
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
363 |
|
|
& SQUEEZE_RIGHT , 1) |
364 |
|
|
ENDIF |
365 |
|
|
#endif /* ALLOW_OPPS_DEBUG */ |
366 |
|
|
DO ktr=1,nTracerInUse |
367 |
|
|
Pda(maxdepth,ktr) = Pd(maxdepth,ktr) |
368 |
|
|
ENDDO |
369 |
|
|
C |
370 |
|
|
kmx = maxdepth-1 |
371 |
|
|
do nn=1,ntime |
372 |
|
|
C |
373 |
|
|
C top point |
374 |
|
|
C |
375 |
|
|
DO ktr = 1,nTracerInUse |
376 |
|
|
tracerEnv(k,ktr) = tracerEnv(k,ktr)- |
377 |
|
|
& (mda(k)*(Pda(k,ktr)-Paa(k,ktr)))*dt*recip_drF(k) |
378 |
|
|
ENDDO |
379 |
|
|
c |
380 |
|
|
c now do inner points if there are any |
381 |
|
|
c |
382 |
|
|
CML if(Maxdepth-k.gt.1) then |
383 |
|
|
CML This if statement is superfluous |
384 |
|
|
CML IF ( k .LT. Maxdepth-1 ) THEN |
385 |
|
|
CML DO k2=k+1,Maxdepth-1 |
386 |
|
|
CML mda(maxDepth) = 0. |
387 |
|
|
DO k2=k+1,kmx |
388 |
|
|
k2m1 = max(k,k2-1) |
389 |
|
|
k2p1 = max(k2+1,maxDepth) |
390 |
|
|
c |
391 |
|
|
DO ktr = 1,nTracerInUse |
392 |
|
|
tracerEnv(k2,ktr) = tracerEnv(k2,ktr) + |
393 |
|
|
& (mda(k2m1)*(Pda(k2m1,ktr)-Paa(k2m1,ktr)) |
394 |
|
|
& -mda(k2) *(Pda(k2,ktr) -Paa(k2,ktr)) ) |
395 |
|
|
& *dt*recip_drF(k2) |
396 |
|
|
ENDDO |
397 |
|
|
ENDDO |
398 |
|
|
CML This if statement is superfluous |
399 |
|
|
CML ENDIF |
400 |
|
|
C |
401 |
|
|
C bottom point |
402 |
|
|
C |
403 |
|
|
DO ktr=1,nTracerInUse |
404 |
|
|
tracerEnv(kmx+1,ktr) = tracerEnv(kmx+1,ktr)+ |
405 |
|
|
& mda(kmx)*(Pda(kmx,ktr)-Paa(kmx,ktr))*dt*recip_drF(kmx+1) |
406 |
|
|
ENDDO |
407 |
|
|
c |
408 |
|
|
c set the environmental temp and salinity to equal new fields |
409 |
|
|
c |
410 |
|
|
DO ktr=1,nTracerInUse |
411 |
|
|
DO k2=1,kmx |
412 |
|
|
paa(k2,ktr) = tracerEnv(k2+1,ktr) |
413 |
|
|
ENDDO |
414 |
|
|
ENDDO |
415 |
|
|
c |
416 |
|
|
c end loop on number of time integration steps |
417 |
|
|
c |
418 |
|
|
enddo |
419 |
|
|
ENDDO |
420 |
|
|
999 continue |
421 |
|
|
C |
422 |
|
|
C count convection event in this grid cell |
423 |
|
|
C |
424 |
|
|
OPPSconvectCount(I,J,K,bi,bj) = |
425 |
|
|
& OPPSconvectCount(I,J,K,bi,bj) + 1. _d 0 |
426 |
|
|
C |
427 |
|
|
C jump here if k = maxdepth or if level not unstable, go to next |
428 |
|
|
C profile point |
429 |
|
|
C |
430 |
|
|
1000 continue |
431 |
|
|
c |
432 |
|
|
C |
433 |
|
|
C end of k-loop |
434 |
|
|
C |
435 |
|
|
ENDDO |
436 |
|
|
|
437 |
jmc |
1.4 |
C-- End IF (DIFF_BASE_MULTIPLE) |
438 |
mlosch |
1.1 |
ENDIF |
439 |
|
|
|
440 |
|
|
RETURN |
441 |
|
|
END |
442 |
|
|
_RL FUNCTION STATE1(sLoc,tLoc,I,J,KREF,bi,bj,mythid) |
443 |
|
|
C !DESCRIPTION: \bv |
444 |
|
|
C *===============================================================* |
445 |
|
|
C | o SUBROUTINE STATE1 |
446 |
|
|
C | Calculates rho(S,T,p) |
447 |
|
|
C | It is absolutely necessary to compute |
448 |
|
|
C | the full rho and not sigma=rho-rhoConst, because |
449 |
|
|
C | density is used as a scale factor for fluxes and velocities |
450 |
|
|
C *===============================================================* |
451 |
|
|
C \ev |
452 |
|
|
|
453 |
|
|
C !USES: |
454 |
|
|
IMPLICIT NONE |
455 |
|
|
C == Global variables == |
456 |
|
|
#include "SIZE.h" |
457 |
|
|
#include "EEPARAMS.h" |
458 |
|
|
#include "PARAMS.h" |
459 |
|
|
#include "EOS.h" |
460 |
|
|
#include "GRID.h" |
461 |
|
|
#include "DYNVARS.h" |
462 |
|
|
|
463 |
|
|
C !INPUT/OUTPUT PARAMETERS: |
464 |
|
|
C == Routine arguments == |
465 |
|
|
INTEGER I,J,kRef,bi,bj,myThid |
466 |
|
|
_RL tLoc,sLoc |
467 |
|
|
|
468 |
|
|
C !LOCAL VARIABLES: |
469 |
|
|
C == Local variables == |
470 |
|
|
_RL rhoLoc, dRho |
471 |
|
|
_RL pLoc |
472 |
|
|
_RL t1, t2, t3, t4, s1, s3o2, p1, p2, sp5, p1t1 |
473 |
|
|
_RL rfresh, rsalt, rhoP0 |
474 |
|
|
_RL bMfresh, bMsalt, bMpres, BulkMod |
475 |
|
|
_RL rhoNum, rhoDen, den, epsln |
476 |
|
|
PARAMETER ( epsln = 0.D0 ) |
477 |
|
|
|
478 |
|
|
character*(max_len_mbuf) msgbuf |
479 |
|
|
|
480 |
|
|
CMLC estimate pressure from depth at cell centers |
481 |
|
|
CML mtoSI = gravity*rhoConst |
482 |
|
|
CML pLoc = ABS(rC(kRef))*mtoSI |
483 |
|
|
|
484 |
|
|
IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN |
485 |
|
|
C in Z coordinates the pressure is rho0 * (hydrostatic) Potential |
486 |
|
|
IF ( useDynP_inEos_Zc ) THEN |
487 |
|
|
C---------- |
488 |
|
|
C NOTE: For now, totPhiHyd only contains the Potential anomaly |
489 |
|
|
C since PhiRef is not available for Atmos and has not (yet) |
490 |
|
|
C been added in S/R DIAGS_PHI_HYD |
491 |
|
|
C---------- |
492 |
|
|
pLoc = rhoConst*( totPhiHyd(i,j,kRef,bi,bj) |
493 |
|
|
& -rC(kRef)*gravity |
494 |
|
|
& )*maskC(i,j,kRef,bi,bj) |
495 |
|
|
ELSE |
496 |
|
|
pLoc = -rhoConst*rC(kRef)*gravity*maskC(i,j,kRef,bi,bj) |
497 |
|
|
ENDIF |
498 |
|
|
ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN |
499 |
|
|
C in P coordinates the pressure is just the coordinate of |
500 |
|
|
C the tracer point |
501 |
|
|
pLoc = rC(kRef)* maskC(i,j,kRef,bi,bj) |
502 |
|
|
ENDIF |
503 |
|
|
|
504 |
|
|
rhoLoc = 0. _d 0 |
505 |
|
|
rhoP0 = 0. _d 0 |
506 |
|
|
bulkMod = 0. _d 0 |
507 |
|
|
rfresh = 0. _d 0 |
508 |
|
|
rsalt = 0. _d 0 |
509 |
|
|
bMfresh = 0. _d 0 |
510 |
|
|
bMsalt = 0. _d 0 |
511 |
|
|
bMpres = 0. _d 0 |
512 |
|
|
rhoNum = 0. _d 0 |
513 |
|
|
rhoDen = 0. _d 0 |
514 |
|
|
den = 0. _d 0 |
515 |
|
|
|
516 |
|
|
t1 = tLoc |
517 |
|
|
t2 = t1*t1 |
518 |
|
|
t3 = t2*t1 |
519 |
|
|
t4 = t3*t1 |
520 |
|
|
|
521 |
|
|
s1 = sLoc |
522 |
|
|
|
523 |
|
|
IF ( equationOfState .EQ. 'LINEAR' ) THEN |
524 |
|
|
|
525 |
|
|
dRho = rhoNil-rhoConst |
526 |
|
|
rhoLoc=rhoNil* ( |
527 |
|
|
& sBeta *(sLoc-sRef(kRef)) |
528 |
|
|
& - tAlpha*(tLoc-tRef(KREF)) ) + dRho |
529 |
|
|
|
530 |
|
|
ELSEIF (equationOfState.EQ.'POLY3') THEN |
531 |
|
|
|
532 |
|
|
C this is not correct, there is a field eosSig0 which should be use here |
533 |
|
|
C but I DO not intent to include the reference level in this routine |
534 |
|
|
WRITE(*,'(a)') |
535 |
|
|
& ' FIND_RHO_SCALAR: for POLY3, the density is not' |
536 |
|
|
WRITE(*,'(a)') |
537 |
|
|
& ' computed correctly in this routine' |
538 |
|
|
rhoLoc = 0. _d 0 |
539 |
|
|
|
540 |
|
|
ELSEIF ( equationOfState(1:5).EQ.'JMD95' |
541 |
|
|
& .OR. equationOfState.EQ.'UNESCO' ) THEN |
542 |
|
|
C nonlinear equation of state in pressure coordinates |
543 |
|
|
|
544 |
|
|
s3o2 = s1*SQRT(s1) |
545 |
|
|
|
546 |
|
|
p1 = pLoc*SItoBar |
547 |
|
|
p2 = p1*p1 |
548 |
|
|
|
549 |
|
|
C density of freshwater at the surface |
550 |
|
|
rfresh = |
551 |
|
|
& eosJMDCFw(1) |
552 |
|
|
& + eosJMDCFw(2)*t1 |
553 |
|
|
& + eosJMDCFw(3)*t2 |
554 |
|
|
& + eosJMDCFw(4)*t3 |
555 |
|
|
& + eosJMDCFw(5)*t4 |
556 |
|
|
& + eosJMDCFw(6)*t4*t1 |
557 |
|
|
C density of sea water at the surface |
558 |
|
|
rsalt = |
559 |
|
|
& s1*( |
560 |
|
|
& eosJMDCSw(1) |
561 |
|
|
& + eosJMDCSw(2)*t1 |
562 |
|
|
& + eosJMDCSw(3)*t2 |
563 |
|
|
& + eosJMDCSw(4)*t3 |
564 |
|
|
& + eosJMDCSw(5)*t4 |
565 |
|
|
& ) |
566 |
|
|
& + s3o2*( |
567 |
|
|
& eosJMDCSw(6) |
568 |
|
|
& + eosJMDCSw(7)*t1 |
569 |
|
|
& + eosJMDCSw(8)*t2 |
570 |
|
|
& ) |
571 |
|
|
& + eosJMDCSw(9)*s1*s1 |
572 |
|
|
|
573 |
|
|
rhoP0 = rfresh + rsalt |
574 |
|
|
|
575 |
|
|
C secant bulk modulus of fresh water at the surface |
576 |
|
|
bMfresh = |
577 |
|
|
& eosJMDCKFw(1) |
578 |
|
|
& + eosJMDCKFw(2)*t1 |
579 |
|
|
& + eosJMDCKFw(3)*t2 |
580 |
|
|
& + eosJMDCKFw(4)*t3 |
581 |
|
|
& + eosJMDCKFw(5)*t4 |
582 |
|
|
C secant bulk modulus of sea water at the surface |
583 |
|
|
bMsalt = |
584 |
|
|
& s1*( eosJMDCKSw(1) |
585 |
|
|
& + eosJMDCKSw(2)*t1 |
586 |
|
|
& + eosJMDCKSw(3)*t2 |
587 |
|
|
& + eosJMDCKSw(4)*t3 |
588 |
|
|
& ) |
589 |
|
|
& + s3o2*( eosJMDCKSw(5) |
590 |
|
|
& + eosJMDCKSw(6)*t1 |
591 |
|
|
& + eosJMDCKSw(7)*t2 |
592 |
|
|
& ) |
593 |
|
|
C secant bulk modulus of sea water at pressure p |
594 |
|
|
bMpres = |
595 |
|
|
& p1*( eosJMDCKP(1) |
596 |
|
|
& + eosJMDCKP(2)*t1 |
597 |
|
|
& + eosJMDCKP(3)*t2 |
598 |
|
|
& + eosJMDCKP(4)*t3 |
599 |
|
|
& ) |
600 |
|
|
& + p1*s1*( eosJMDCKP(5) |
601 |
|
|
& + eosJMDCKP(6)*t1 |
602 |
|
|
& + eosJMDCKP(7)*t2 |
603 |
|
|
& ) |
604 |
|
|
& + p1*s3o2*eosJMDCKP(8) |
605 |
|
|
& + p2*( eosJMDCKP(9) |
606 |
|
|
& + eosJMDCKP(10)*t1 |
607 |
|
|
& + eosJMDCKP(11)*t2 |
608 |
|
|
& ) |
609 |
|
|
& + p2*s1*( eosJMDCKP(12) |
610 |
|
|
& + eosJMDCKP(13)*t1 |
611 |
|
|
& + eosJMDCKP(14)*t2 |
612 |
|
|
& ) |
613 |
|
|
|
614 |
|
|
bulkMod = bMfresh + bMsalt + bMpres |
615 |
|
|
|
616 |
|
|
C density of sea water at pressure p |
617 |
|
|
rhoLoc = rhoP0/(1. _d 0 - p1/bulkMod) - rhoConst |
618 |
|
|
|
619 |
|
|
ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN |
620 |
|
|
|
621 |
|
|
sp5 = SQRT(s1) |
622 |
|
|
|
623 |
|
|
p1 = pLoc*SItodBar |
624 |
|
|
p1t1 = p1*t1 |
625 |
|
|
|
626 |
|
|
rhoNum = eosMDJWFnum(0) |
627 |
|
|
& + t1*(eosMDJWFnum(1) |
628 |
|
|
& + t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) ) |
629 |
|
|
& + s1*(eosMDJWFnum(4) |
630 |
|
|
& + eosMDJWFnum(5)*t1 + eosMDJWFnum(6)*s1) |
631 |
|
|
& + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2 |
632 |
|
|
& + eosMDJWFnum(9)*s1 |
633 |
|
|
& + p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) ) |
634 |
|
|
|
635 |
|
|
|
636 |
|
|
den = eosMDJWFden(0) |
637 |
|
|
& + t1*(eosMDJWFden(1) |
638 |
|
|
& + t1*(eosMDJWFden(2) |
639 |
|
|
& + t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) ) |
640 |
|
|
& + s1*(eosMDJWFden(5) |
641 |
|
|
& + t1*(eosMDJWFden(6) |
642 |
|
|
& + eosMDJWFden(7)*t2) |
643 |
|
|
& + sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) ) |
644 |
|
|
& + p1*(eosMDJWFden(10) |
645 |
|
|
& + p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) ) |
646 |
|
|
|
647 |
|
|
rhoDen = 1.0/(epsln+den) |
648 |
|
|
|
649 |
|
|
rhoLoc = rhoNum*rhoDen - rhoConst |
650 |
|
|
|
651 |
|
|
ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN |
652 |
|
|
C |
653 |
|
|
ELSE |
654 |
|
|
WRITE(msgbuf,'(3A)') |
655 |
|
|
& ' STATE1 : equationOfState = "', |
656 |
|
|
& equationOfState,'"' |
657 |
|
|
CALL PRINT_ERROR( msgbuf, mythid ) |
658 |
|
|
STOP 'ABNORMAL END: S/R STATE1 in OPPS_CALC' |
659 |
|
|
ENDIF |
660 |
|
|
|
661 |
|
|
state1 = rhoLoc + rhoConst |
662 |
|
|
|
663 |
|
|
#endif /* ALLOW_OPPS */ |
664 |
|
|
RETURN |
665 |
|
|
END |
666 |
|
|
|
667 |
|
|
|
668 |
|
|
#undef OPPS_ORGCODE |
669 |
|
|
#ifdef OPPS_ORGCODE |
670 |
|
|
c Listed below is the subroutine for use in parallel 3-d circulation code. |
671 |
|
|
c It has been used in the parallel semtner-chervin code and is now being used |
672 |
|
|
c In the POP code. The subroutine is called nlopps (long story to explain why). |
673 |
|
|
|
674 |
|
|
c I've attached the version of lopps that we've been using in the simulations. |
675 |
|
|
c There is one common block that is different from the standard model commons |
676 |
|
|
c (countc) and it is not needed if the array convadj is not used. The routine |
677 |
|
|
c does need "kmp" which is why the boundc common is included. For massively |
678 |
|
|
c parallel codes (like POP) we think this will work well when converted from a |
679 |
|
|
c "slab" (i=is,ie) to a column, which just means removing the "do i=is,ie" loop. c There are differences between this |
680 |
|
|
c code and the 1-d code and the earlier scheme implemented in 3-d models. These c differences are described below. |
681 |
|
|
|
682 |
|
|
|
683 |
|
|
subroutine nlopps(j,is,ie,ta,sa,gcmdz) |
684 |
|
|
c |
685 |
|
|
parameter (imt = 361 , jmt = 301 , km = 30 ) |
686 |
|
|
c |
687 |
|
|
c Nlopps: E. Skyllingstad and T. Paluszkiewicz |
688 |
|
|
c |
689 |
|
|
c Version: December 11, 1996 |
690 |
|
|
c |
691 |
|
|
c Nlopps: This version of lopps is significantly different from |
692 |
|
|
c the original code developed by R. Romea and T. Paluskiewicz. The |
693 |
|
|
c code uses a flux constraint to control the change in T and S at |
694 |
|
|
c each grid level. First, a plume profile of T,S, and W are |
695 |
|
|
c determined using the standard plume model, but with a detraining |
696 |
|
|
c mass instead of entraining. Thus, the T and S plume |
697 |
|
|
c characteristics still change, but the plume contracts in size |
698 |
|
|
c rather than expanding ala classical entraining plumes. This |
699 |
|
|
c is heuristically more in line with large eddy simulation results. |
700 |
|
|
c At each grid level, the convergence of plume velocity determines |
701 |
|
|
c the flux of T and S, which is conserved by using an upstream |
702 |
|
|
c advection. The vertical velocity is balanced so that the area |
703 |
|
|
c weighted upward velocity equals the area weighted downdraft |
704 |
|
|
c velocity, ensuring mass conservation. The present implementation |
705 |
|
|
c adjusts the plume for a time period equal to the time for 1/2 of |
706 |
|
|
c the mass of the fastest moving level to move downward. As a |
707 |
|
|
c consequence, the model does not completely adjust the profile at |
708 |
|
|
c each model time step, but provides a smooth adjustment over time. |
709 |
|
|
c |
710 |
|
|
c |
711 |
|
|
|
712 |
|
|
|
713 |
|
|
c |
714 |
|
|
|
715 |
|
|
c include "params.h" |
716 |
|
|
c include "plume_fast_inc.h" |
717 |
|
|
c include "plume_fast.h" |
718 |
|
|
c #include "loppsd.h" |
719 |
|
|
|
720 |
|
|
real ta(imt,km),sa(imt,km),gcmdz(km),dz(km) |
721 |
|
|
real pdensity,wsqr,oldflux,newflux,entrainrate,adtemp |
722 |
|
|
REAL Del,D,dza1,dza2,kd,kd1,Smix,Thmix,PlumeS,PlumeT,PlumeD |
723 |
|
|
c |
724 |
|
|
c |
725 |
|
|
|
726 |
|
|
INTEGER i,j,k |
727 |
|
|
clfh |
728 |
|
|
integer is,ie,k2 |
729 |
|
|
clfh |
730 |
|
|
REAL D1,D2,state1,Density |
731 |
|
|
REAL dz1,dz2 |
732 |
|
|
REAL radius,StartingFlux |
733 |
|
|
real ttemp(km),stemp(km),taa(km),saa(km) |
734 |
|
|
real wda(km),tda(km),sda(km),mda(km) |
735 |
|
|
real dtts,dt,sumo,sumn |
736 |
|
|
integer ntime,nn,kmx,ic |
737 |
|
|
c |
738 |
|
|
c |
739 |
|
|
LOGICAL debug,done |
740 |
|
|
INTEGER MAX_ABE_ITERATIONS |
741 |
|
|
PARAMETER(MAX_ABE_ITERATIONS=1) |
742 |
|
|
REAL PlumeRadius |
743 |
|
|
REAL STABILITY_THRESHOLD |
744 |
|
|
REAL FRACTIONAL_AREA |
745 |
|
|
REAL MAX_FRACTIONAL_AREA |
746 |
|
|
REAL VERTICAL_VELOCITY |
747 |
|
|
REAL ENTRAINMENT_RATE |
748 |
|
|
REAL e2 |
749 |
|
|
PARAMETER ( PlumeRadius = 100.D0 ) |
750 |
|
|
PARAMETER ( STABILITY_THRESHOLD = -1.E-4 ) |
751 |
|
|
PARAMETER ( FRACTIONAL_AREA = .1E0 ) |
752 |
|
|
PARAMETER ( MAX_FRACTIONAL_AREA = .8E0 ) |
753 |
|
|
PARAMETER ( VERTICAL_VELOCITY = .02E0 ) |
754 |
|
|
PARAMETER ( ENTRAINMENT_RATE = -.05E0 ) |
755 |
|
|
PARAMETER ( e2 = 2.E0*ENTRAINMENT_RATE ) |
756 |
|
|
! Arrays. |
757 |
|
|
REAL Ad(km),Sd(km),Td(km),Wd(km),Dd(km),Md(km) |
758 |
|
|
REAL Se(km),Te(km),We(km),De(km) |
759 |
|
|
REAL PlumeEntrainment(km) |
760 |
|
|
REAL GridThickness(km) |
761 |
|
|
|
762 |
|
|
c |
763 |
|
|
c input kmp through a common block |
764 |
|
|
c |
765 |
|
|
common / boundc / wsx(imt,jmt),wsy(imt,jmt),hfs(imt,jmt), |
766 |
|
|
1 ple(imt,jmt),kmp(imt,jmt),kmq(imt,jmt) |
767 |
|
|
cwmseas |
768 |
|
|
& ,wsx1(imt,jmt),wsy1(imt,jmt) |
769 |
|
|
1 ,wsx2(imt,jmt),wsy2(imt,jmt) |
770 |
|
|
c |
771 |
|
|
c input the variables through a common |
772 |
|
|
c |
773 |
|
|
logical problem |
774 |
|
|
common /countc/ convadj(imt,jmt,km),ics,depth(km),problem |
775 |
|
|
|
776 |
|
|
|
777 |
|
|
c-----may want to setup an option to get this only on first call |
778 |
|
|
c otherwise it is repetive |
779 |
|
|
c griddz is initialize by call to setupgrid |
780 |
|
|
c |
781 |
|
|
c |
782 |
|
|
|
783 |
|
|
dtts = 2400 |
784 |
|
|
c |
785 |
|
|
do k=1,km |
786 |
|
|
dz(k) = 0.01*gcmdz(k) |
787 |
|
|
enddo |
788 |
|
|
c |
789 |
|
|
do k=1,km |
790 |
|
|
GridThickness(k) = dz(k) |
791 |
|
|
enddo |
792 |
|
|
c |
793 |
|
|
c modified to loop over slab |
794 |
|
|
c |
795 |
|
|
DO i=is,ie |
796 |
|
|
|
797 |
|
|
numgridpoints=kmp(i,j) |
798 |
|
|
c |
799 |
|
|
c go to next column if only 1 grid point or on land |
800 |
|
|
c |
801 |
|
|
if(numgridpoints.le.1) goto 1100 |
802 |
|
|
c |
803 |
|
|
c loop over depth |
804 |
|
|
c |
805 |
|
|
c |
806 |
|
|
debug = .false. |
807 |
|
|
c |
808 |
|
|
c first save copy of initial profile |
809 |
|
|
c |
810 |
|
|
DO k=1,NumGridPoints |
811 |
|
|
stemp(k)=sa(i,k) |
812 |
|
|
ttemp(k)=ta(i,k) |
813 |
|
|
c |
814 |
|
|
c do a check of t and s range, if out of bounds set flag |
815 |
|
|
c |
816 |
|
|
if(problem) then |
817 |
|
|
write(0,*)"Code in trouble before this nlopps call" |
818 |
|
|
return |
819 |
|
|
endif |
820 |
|
|
c |
821 |
|
|
if(sa(i,k).gt.40..or.ta(i,k).lt.-4.0) then |
822 |
|
|
problem = .true. |
823 |
|
|
write(0,*)"t out of range at j ",j |
824 |
|
|
debug = .true. |
825 |
|
|
return |
826 |
|
|
endif |
827 |
|
|
ENDDO |
828 |
|
|
|
829 |
|
|
if(debug) then |
830 |
|
|
write(*,*)"T and S Profile at ",i,j |
831 |
|
|
write(*,*)(ta(i,k),sa(i,k),k=1,NumGridPoints) |
832 |
|
|
endif |
833 |
|
|
|
834 |
|
|
DO k=1,NumGridPoints-1 |
835 |
|
|
c |
836 |
|
|
c initialize the plume T,S,density, and w velocity |
837 |
|
|
c |
838 |
|
|
Sd(k)=stemp(k) |
839 |
|
|
Td(k)=ttemp(k) |
840 |
|
|
Dd(k)=state1(stemp(k),ttemp(k),k) |
841 |
|
|
De(k)=Dd(k) |
842 |
|
|
c Wd(k)=VERTICAL_VELOCITY |
843 |
|
|
c |
844 |
|
|
c guess at initial top grid cell vertical velocity |
845 |
|
|
c |
846 |
|
|
Wd(k) = 0.03 |
847 |
|
|
c |
848 |
|
|
c these estimates of initial plume velocity based on plume size and |
849 |
|
|
c top grid cell water mass |
850 |
|
|
c |
851 |
|
|
c Wd(k) = 0.5*dz(k)/(dtts*FRACTIONAL_AREA) |
852 |
|
|
c Wd(k) = 0.5*dz(k)/dtts |
853 |
|
|
c |
854 |
|
|
wsqr=Wd(k)*Wd(k) |
855 |
|
|
PlumeEntrainment(k) = 0.0 |
856 |
|
|
c |
857 |
|
|
c |
858 |
|
|
c |
859 |
|
|
if(debug) write(0,*) 'Doing old lowerparcel' |
860 |
|
|
radius=PlumeRadius |
861 |
|
|
StartingFlux=radius*radius*Wd(k)*Dd(k) |
862 |
|
|
oldflux=StartingFlux |
863 |
|
|
|
864 |
|
|
dz2=GridThickness(k) |
865 |
|
|
DO k2=k,NumGridPoints-1 |
866 |
|
|
D1=state1(Sd(k2),Td(k2),k2+1) |
867 |
|
|
D2=state1(stemp(k2+1),ttemp(k2+1),k2+1) |
868 |
|
|
De(k2+1)=D2 |
869 |
|
|
c |
870 |
|
|
c To start downward, parcel has to initially be heavier than environment |
871 |
|
|
c but after it has started moving, we continue plume until plume tke or |
872 |
|
|
c flux goes negative |
873 |
|
|
c |
874 |
|
|
IF (D2-D1 .LT. STABILITY_THRESHOLD.or.k2.ne.k) THEN |
875 |
|
|
dz1=dz2 |
876 |
|
|
dz2=GridThickness(k2+1) |
877 |
|
|
c |
878 |
|
|
c define mass flux according to eq. 4 from paper |
879 |
|
|
c |
880 |
|
|
newflux=oldflux+e2*radius*Wd(k2)*Dd(k2)*0.50* |
881 |
|
|
. (dz1+dz2) |
882 |
|
|
c |
883 |
|
|
PlumeEntrainment(k2+1) = newflux/StartingFlux |
884 |
|
|
c |
885 |
|
|
IF(newflux.LT.0.0) then |
886 |
|
|
if(debug) then |
887 |
|
|
write(0,*)"Plume entrained to zero at ",k2 |
888 |
|
|
endif |
889 |
|
|
maxdepth = k2 |
890 |
|
|
if(maxdepth.eq.k) goto 1000 |
891 |
|
|
goto 1 |
892 |
|
|
endif |
893 |
|
|
c |
894 |
|
|
c entrainment rate is basically a scaled mass flux dM/M |
895 |
|
|
c |
896 |
|
|
entrainrate = (newflux - oldflux)/newflux |
897 |
|
|
oldflux = newflux |
898 |
|
|
c |
899 |
|
|
c |
900 |
|
|
c mix var's are the average environmental values over the two grid levels |
901 |
|
|
c |
902 |
|
|
smix=(dz1*stemp(k2)+dz2*stemp(k2+1))/(dz1+dz2) |
903 |
|
|
thmix=(dz1*ttemp(k2)+dz2*ttemp(k2+1))/(dz1+dz2) |
904 |
|
|
c |
905 |
|
|
c first compute the new salinity and temperature for this level |
906 |
|
|
c using equations 3.6 and 3.7 from the paper |
907 |
|
|
c |
908 |
|
|
c |
909 |
|
|
c |
910 |
|
|
sd(k2+1)=sd(k2) - entrainrate*(smix - sd(k2)) |
911 |
|
|
td(k2+1)=td(k2) - entrainrate*(thmix - td(k2)) |
912 |
|
|
c |
913 |
|
|
c |
914 |
|
|
c compute the density at this level for the buoyancy term in the |
915 |
|
|
c vertical k.e. equation |
916 |
|
|
c |
917 |
|
|
Dd(k2+1)=state1(Sd(k2+1),Td(k2+1),k2+1) |
918 |
|
|
c |
919 |
|
|
c next, solve for the vertical velocity k.e. using combined eq. 4 |
920 |
|
|
c and eq 5 from the paper |
921 |
|
|
c |
922 |
|
|
if(debug) then |
923 |
|
|
write(0,*)"Dd,De,entr,k ",Dd(k2),De(k2),entrainrate,k2 |
924 |
|
|
endif |
925 |
|
|
c |
926 |
|
|
wsqr = wsqr - wsqr*abs(entrainrate)+ 9.81* |
927 |
|
|
. (dz1*(Dd(k2)-De(k2))/De(k2) |
928 |
|
|
. +dz2*(Dd(k2+1)-De(k2+1))/De(k2+1)) |
929 |
|
|
c |
930 |
|
|
c if negative k.e. then plume has reached max depth, get out of loop |
931 |
|
|
c |
932 |
|
|
IF(wsqr.LT.0.0)then |
933 |
|
|
maxdepth = k2 |
934 |
|
|
if(debug) then |
935 |
|
|
write(0,*)"Plume velocity went to zero at ",k2 |
936 |
|
|
endif |
937 |
|
|
if(maxdepth.eq.k) goto 1000 |
938 |
|
|
goto 1 |
939 |
|
|
endif |
940 |
|
|
Wd(k2+1)=sqrt(wsqr) |
941 |
|
|
c |
942 |
|
|
c compute a new radius based on the new mass flux at this grid level |
943 |
|
|
c |
944 |
|
|
radius=sqrt(newflux/(Wd(k2)*Dd(k2))) |
945 |
|
|
ELSE |
946 |
|
|
maxdepth=k2 |
947 |
|
|
if(maxdepth.eq.k) goto 1000 |
948 |
|
|
GOTO 1 |
949 |
|
|
ENDIF |
950 |
|
|
ENDDO |
951 |
|
|
c |
952 |
|
|
c plume has reached the bottom |
953 |
|
|
c |
954 |
|
|
MaxDepth=NumGridPoints |
955 |
|
|
c |
956 |
|
|
1 continue |
957 |
|
|
c |
958 |
|
|
Ad(k)=FRACTIONAL_AREA |
959 |
|
|
IC=0 |
960 |
|
|
c |
961 |
|
|
c start iteration on fractional area, not used in OGCM implementation |
962 |
|
|
c |
963 |
|
|
c |
964 |
|
|
DO IC=1,Max_ABE_Iterations |
965 |
|
|
c |
966 |
|
|
c |
967 |
|
|
c next compute the mass flux beteen each grid box using the entrainment |
968 |
|
|
c |
969 |
|
|
92 continue |
970 |
|
|
Md(k)=Wd(k)*Ad(k) |
971 |
|
|
c |
972 |
|
|
DO k2=k+1,MaxDepth |
973 |
|
|
Md(k2)=Md(k)*PlumeEntrainment(k2) |
974 |
|
|
if(debug) then |
975 |
|
|
write(0,*)"Md, Wd, and k are ",Md(k2),Wd(k2),k2 |
976 |
|
|
endif |
977 |
|
|
ENDDO |
978 |
|
|
c |
979 |
|
|
c Now move on to calculate new temperature using flux from |
980 |
|
|
c Td, Sd, Wd, ta, sa, and we. Values for these variables are at |
981 |
|
|
c center of grid cell, use weighted average to get boundary values |
982 |
|
|
c |
983 |
|
|
c use a timestep limited by the GCM model timestep and the maximum plume |
984 |
|
|
c velocity (CFL criteria) |
985 |
|
|
c |
986 |
|
|
c |
987 |
|
|
c calculate the weighted wd, td, and sd |
988 |
|
|
c |
989 |
|
|
dt = dtts |
990 |
|
|
do k2=k,maxdepth-1 |
991 |
|
|
dt = min(dt,dz(k2)/wd(k2)) |
992 |
|
|
c |
993 |
|
|
c time integration will be integer number of steps to get one |
994 |
|
|
c gcm time step |
995 |
|
|
c |
996 |
|
|
ntime = nint(0.5*int(dtts/dt)) |
997 |
|
|
if(ntime.eq.0) then |
998 |
|
|
ntime = 1 |
999 |
|
|
endif |
1000 |
|
|
c |
1001 |
|
|
c make sure area weighted vertical velocities match; in other words |
1002 |
|
|
c make sure mass in equals mass out at the intersection of each grid |
1003 |
|
|
c cell. |
1004 |
|
|
c |
1005 |
|
|
mda(k2) = (md(k2)*dz(k2)+md(k2+1)*dz(k2+1))/ |
1006 |
|
|
* (dz(k2)+dz(k2+1)) |
1007 |
|
|
c |
1008 |
|
|
wda(k2) = (wd(k2)*dz(k2)+wd(k2+1)*dz(k2+1))/ |
1009 |
|
|
* (dz(k2)+dz(k2+1)) |
1010 |
|
|
c |
1011 |
|
|
tda(k2) = td(k2) |
1012 |
|
|
sda(k2) = sd(k2) |
1013 |
|
|
c |
1014 |
|
|
taa(k2) = ttemp(k2+1) |
1015 |
|
|
saa(k2) = stemp(k2+1) |
1016 |
|
|
c |
1017 |
|
|
enddo |
1018 |
|
|
dt = min(dt,dtts) |
1019 |
|
|
if(debug) then |
1020 |
|
|
write(0,*)"Time step is ", dt |
1021 |
|
|
endif |
1022 |
|
|
tda(maxdepth) = td(maxdepth) |
1023 |
|
|
sda(maxdepth) = sd(maxdepth) |
1024 |
|
|
c |
1025 |
|
|
c do top and bottom points first |
1026 |
|
|
c |
1027 |
|
|
kmx = maxdepth-1 |
1028 |
|
|
do nn=1,ntime |
1029 |
|
|
|
1030 |
|
|
ttemp(k) = ttemp(k)- |
1031 |
|
|
* (mda(k)*(tda(k)-taa(k)))*dt*recip_drF(k) |
1032 |
|
|
c |
1033 |
|
|
stemp(k) = stemp(k)- |
1034 |
|
|
* (mda(k)*(sda(k)-saa(k)))*dt*recip_drF(k) |
1035 |
|
|
c |
1036 |
|
|
c |
1037 |
|
|
c now do inner points if there are any |
1038 |
|
|
c |
1039 |
|
|
if(Maxdepth-k.gt.1) then |
1040 |
|
|
do k2=k+1,Maxdepth-1 |
1041 |
|
|
c |
1042 |
|
|
ttemp(k2) = ttemp(k2) + |
1043 |
|
|
* (mda(k2-1)*(tda(k2-1)-taa(k2-1))- |
1044 |
|
|
* mda(k2)*(tda(k2)-taa(k2))) |
1045 |
|
|
* *dt*recip_drF(k2) |
1046 |
|
|
|
1047 |
|
|
c |
1048 |
|
|
stemp(k2) = stemp(k2) + |
1049 |
|
|
* (mda(k2-1)*(sda(k2-1)-saa(k2-1))- |
1050 |
|
|
* mda(k2)*(sda(k2)-saa(k2))) |
1051 |
|
|
* *dt*recip_drF(k2) |
1052 |
|
|
|
1053 |
|
|
c |
1054 |
|
|
enddo |
1055 |
|
|
endif |
1056 |
|
|
ttemp(kmx+1) = ttemp(kmx+1)+ |
1057 |
|
|
* (mda(kmx)*(tda(kmx)-taa(kmx)))* |
1058 |
|
|
* dt*recip_drF(kmx+1) |
1059 |
|
|
c |
1060 |
|
|
stemp(kmx+1) = stemp(kmx+1)+ |
1061 |
|
|
* (mda(kmx)*(sda(kmx)-saa(kmx)))* |
1062 |
|
|
* dt*recip_drF(kmx+1) |
1063 |
|
|
c |
1064 |
|
|
c set the environmental temp and salinity to equal new fields |
1065 |
|
|
c |
1066 |
|
|
do k2=1,maxdepth-1 |
1067 |
|
|
taa(k2) = ttemp(k2+1) |
1068 |
|
|
saa(k2) = stemp(k2+1) |
1069 |
|
|
enddo |
1070 |
|
|
c |
1071 |
|
|
c end loop on number of time integration steps |
1072 |
|
|
c |
1073 |
|
|
enddo |
1074 |
|
|
ENDDO |
1075 |
|
|
999 continue |
1076 |
|
|
c |
1077 |
|
|
c assume that it converged, so update the ta and sa with new fields |
1078 |
|
|
c |
1079 |
|
|
c if(i.gt.180.and.j.gt.200.and.i.lt.240) then |
1080 |
|
|
c write(*,*)"Converged ",i,j,k,maxdepth,ttemp(k+1),ta(i,k+1) |
1081 |
|
|
c endif |
1082 |
|
|
do k2=k,maxdepth |
1083 |
|
|
convadj(i,j,k2) = convadj(i,j,k2) + (ttemp(k2)- |
1084 |
|
|
* ta(i,k2)) |
1085 |
|
|
sa(i,k2) = stemp(k2) |
1086 |
|
|
ta(i,k2) = ttemp(k2) |
1087 |
|
|
c if(i.gt.180.and.j.gt.200.and.i.lt.240) then |
1088 |
|
|
c write(*,*)"convadj ",convadj(i,j,k2) |
1089 |
|
|
c endif |
1090 |
|
|
c |
1091 |
|
|
c see if nlopps messed up |
1092 |
|
|
c |
1093 |
|
|
if(sa(i,k).gt.40..or.ta(i,k).lt.-4.0) then |
1094 |
|
|
problem = .true. |
1095 |
|
|
write(0,*)"t out of range at j after adjust",j |
1096 |
|
|
debug = .true. |
1097 |
|
|
endif |
1098 |
|
|
c |
1099 |
|
|
enddo |
1100 |
|
|
c |
1101 |
|
|
c jump here if k = maxdepth or if level not unstable, go to next |
1102 |
|
|
c profile point |
1103 |
|
|
c |
1104 |
|
|
1000 continue |
1105 |
|
|
c |
1106 |
|
|
c |
1107 |
|
|
c end loop on k, move on to next possible plume |
1108 |
|
|
c |
1109 |
|
|
ENDDO |
1110 |
|
|
1100 continue |
1111 |
|
|
c |
1112 |
|
|
c i loop |
1113 |
|
|
c |
1114 |
|
|
ENDDO |
1115 |
|
|
END |
1116 |
|
|
|
1117 |
edhill |
1.2 |
#endif /* OPPS_ORGCODE */ |