/[MITgcm]/MITgcm/model/src/do_oceanic_phys.F
ViewVC logotype

Contents of /MITgcm/model/src/do_oceanic_phys.F

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


Revision 1.2 - (show annotations) (download)
Tue Jul 6 21:12:51 2004 UTC (19 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint54b_post
Changes since 1.1: +1 -20 lines
put atmospheric physics & state-vars diagnostics calls in 2 dedicated S/R.

1 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.1 2004/07/06 00:58:40 jmc Exp $
2 C $Name: $
3
4 #include "PACKAGES_CONFIG.h"
5 #include "CPP_OPTIONS.h"
6 c #ifdef ALLOW_PTRACERS
7 c # include "PTRACERS_OPTIONS.h"
8 c #endif
9
10 #ifdef ALLOW_AUTODIFF_TAMC
11 # ifdef ALLOW_GMREDI
12 # include "GMREDI_OPTIONS.h"
13 # endif
14 # ifdef ALLOW_KPP
15 # include "KPP_OPTIONS.h"
16 # endif
17 #endif /* ALLOW_AUTODIFF_TAMC */
18
19 CBOP
20 C !ROUTINE: DO_OCEANIC_PHYS
21 C !INTERFACE:
22 SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
23 C !DESCRIPTION: \bv
24 C *==========================================================*
25 C | SUBROUTINE DO_OCEANIC_PHYS
26 C | o Controlling routine for oceanic physics and
27 C | parameterization
28 C *==========================================================*
29 C | o originally, part of S/R thermodynamics
30 C *==========================================================*
31 C \ev
32
33 C !USES:
34 IMPLICIT NONE
35 C == Global variables ===
36 #include "SIZE.h"
37 #include "EEPARAMS.h"
38 #include "PARAMS.h"
39 #include "DYNVARS.h"
40 #include "GRID.h"
41 c #include "GAD.h"
42 c #ifdef ALLOW_PASSIVE_TRACER
43 c #include "TR1.h"
44 c #endif
45 c #ifdef ALLOW_PTRACERS
46 c #include "PTRACERS.h"
47 c #endif
48 c #ifdef ALLOW_TIMEAVE
49 c #include "TIMEAVE_STATV.h"
50 c #endif
51
52 #ifdef ALLOW_AUTODIFF_TAMC
53 # include "tamc.h"
54 # include "tamc_keys.h"
55 # include "FFIELDS.h"
56 # include "EOS.h"
57 # ifdef ALLOW_KPP
58 # include "KPP.h"
59 # endif
60 # ifdef ALLOW_GMREDI
61 # include "GMREDI.h"
62 # endif
63 # ifdef ALLOW_EBM
64 # include "EBM.h"
65 # endif
66 #endif /* ALLOW_AUTODIFF_TAMC */
67
68 C !INPUT/OUTPUT PARAMETERS:
69 C == Routine arguments ==
70 C myTime - Current time in simulation
71 C myIter - Current iteration number in simulation
72 C myThid - Thread number for this instance of the routine.
73 _RL myTime
74 INTEGER myIter
75 INTEGER myThid
76
77 C !LOCAL VARIABLES:
78 C == Local variables
79 C rhoK, rhoKM1 - Density at current level, and level above
80 C useVariableK = T when vertical diffusion is not constant
81 C iMin, iMax - Ranges and sub-block indices on which calculations
82 C jMin, jMax are applied.
83 C bi, bj
84 C k, kup, - Index for layer above and below. kup and kDown
85 C kDown, km1 are switched with layer to be the appropriate
86 C index into fVerTerm.
87 _RL rhokm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88 _RL rhok (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89 _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
90 _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
91 _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
92 _RL kp1Msk
93 LOGICAL useVariableK
94 INTEGER iMin, iMax
95 INTEGER jMin, jMax
96 INTEGER bi, bj
97 INTEGER i, j
98 INTEGER k, km1, kup, kDown
99 INTEGER iTracer, ip
100
101 CEOP
102
103 #ifdef ALLOW_DEBUG
104 IF ( debugLevel .GE. debLevB )
105 & CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
106 #endif
107
108 #ifdef ALLOW_AUTODIFF_TAMC
109 C-- dummy statement to end declaration part
110 ikey = 1
111 itdkey = 1
112 #endif /* ALLOW_AUTODIFF_TAMC */
113
114 #ifdef ALLOW_AUTODIFF_TAMC
115 C-- HPF directive to help TAMC
116 CHPF$ INDEPENDENT
117 #endif /* ALLOW_AUTODIFF_TAMC */
118
119 DO bj=myByLo(myThid),myByHi(myThid)
120
121 #ifdef ALLOW_AUTODIFF_TAMC
122 C-- HPF directive to help TAMC
123 CHPF$ INDEPENDENT, NEW (rTrans,fVerT,fVerS
124 CHPF$& ,utrans,vtrans,xA,yA
125 CHPF$& )
126 #endif /* ALLOW_AUTODIFF_TAMC */
127
128 DO bi=myBxLo(myThid),myBxHi(myThid)
129
130 #ifdef ALLOW_AUTODIFF_TAMC
131 act1 = bi - myBxLo(myThid)
132 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
133 act2 = bj - myByLo(myThid)
134 max2 = myByHi(myThid) - myByLo(myThid) + 1
135 act3 = myThid - 1
136 max3 = nTx*nTy
137 act4 = ikey_dynamics - 1
138 itdkey = (act1 + 1) + act2*max1
139 & + act3*max1*max2
140 & + act4*max1*max2*max3
141 #endif /* ALLOW_AUTODIFF_TAMC */
142
143 C-- Set up work arrays with valid (i.e. not NaN) values
144 C These inital values do not alter the numerical results. They
145 C just ensure that all memory references are to valid floating
146 C point numbers. This prevents spurious hardware signals due to
147 C uninitialised but inert locations.
148
149 DO j=1-OLy,sNy+OLy
150 DO i=1-OLx,sNx+OLx
151 rhok (i,j) = 0. _d 0
152 rhoKM1 (i,j) = 0. _d 0
153 ENDDO
154 ENDDO
155
156 DO k=1,Nr
157 DO j=1-OLy,sNy+OLy
158 DO i=1-OLx,sNx+OLx
159 C This is currently also used by IVDC and Diagnostics
160 sigmaX(i,j,k) = 0. _d 0
161 sigmaY(i,j,k) = 0. _d 0
162 sigmaR(i,j,k) = 0. _d 0
163 #ifdef ALLOW_AUTODIFF_TAMC
164 cph all the following init. are necessary for TAF
165 cph although some of these are re-initialised later.
166 IVDConvCount(i,j,k,bi,bj) = 0.
167 # ifdef ALLOW_GMREDI
168 Kwx(i,j,k,bi,bj) = 0. _d 0
169 Kwy(i,j,k,bi,bj) = 0. _d 0
170 Kwz(i,j,k,bi,bj) = 0. _d 0
171 # ifdef GM_NON_UNITY_DIAGONAL
172 Kux(i,j,k,bi,bj) = 0. _d 0
173 Kvy(i,j,k,bi,bj) = 0. _d 0
174 # endif
175 # ifdef GM_EXTRA_DIAGONAL
176 Kuz(i,j,k,bi,bj) = 0. _d 0
177 Kvz(i,j,k,bi,bj) = 0. _d 0
178 # endif
179 # ifdef GM_BOLUS_ADVEC
180 GM_PsiX(i,j,k,bi,bj) = 0. _d 0
181 GM_PsiY(i,j,k,bi,bj) = 0. _d 0
182 # endif
183 # ifdef GM_VISBECK_VARIABLE_K
184 VisbeckK(i,j,bi,bj) = 0. _d 0
185 # endif
186 # endif /* ALLOW_GMREDI */
187 #endif /* ALLOW_AUTODIFF_TAMC */
188 ENDDO
189 ENDDO
190 ENDDO
191
192 iMin = 1-OLx
193 iMax = sNx+OLx
194 jMin = 1-OLy
195 jMax = sNy+OLy
196
197 #ifdef ALLOW_AUTODIFF_TAMC
198 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
199 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
200 CADJ STORE totphihyd
201 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
202 #ifdef ALLOW_KPP
203 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
204 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
205 #endif
206 #endif /* ALLOW_AUTODIFF_TAMC */
207
208 #ifdef ALLOW_DEBUG
209 IF ( debugLevel .GE. debLevB )
210 & CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
211 #endif
212
213 C-- Start of diagnostic loop
214 DO k=Nr,1,-1
215
216 #ifdef ALLOW_AUTODIFF_TAMC
217 C? Patrick, is this formula correct now that we change the loop range?
218 C? Do we still need this?
219 cph kkey formula corrected.
220 cph Needed for rhok, rhokm1, in the case useGMREDI.
221 kkey = (itdkey-1)*Nr + k
222 #endif /* ALLOW_AUTODIFF_TAMC */
223
224 C-- Calculate gradients of potential density for isoneutral
225 C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
226 c IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
227 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
228 #ifdef ALLOW_DEBUG
229 IF ( debugLevel .GE. debLevB )
230 & CALL DEBUG_CALL('FIND_RHO',myThid)
231 #endif
232 #ifdef ALLOW_AUTODIFF_TAMC
233 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
234 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
235 #endif /* ALLOW_AUTODIFF_TAMC */
236 CALL FIND_RHO(
237 I bi, bj, iMin, iMax, jMin, jMax, k, k,
238 I theta, salt,
239 O rhoK,
240 I myThid )
241
242 IF (k.GT.1) THEN
243 #ifdef ALLOW_AUTODIFF_TAMC
244 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
245 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
246 #endif /* ALLOW_AUTODIFF_TAMC */
247 CALL FIND_RHO(
248 I bi, bj, iMin, iMax, jMin, jMax, k-1, k,
249 I theta, salt,
250 O rhoKm1,
251 I myThid )
252 ENDIF
253 #ifdef ALLOW_DEBUG
254 IF ( debugLevel .GE. debLevB )
255 & CALL DEBUG_CALL('GRAD_SIGMA',myThid)
256 #endif
257 CALL GRAD_SIGMA(
258 I bi, bj, iMin, iMax, jMin, jMax, k,
259 I rhoK, rhoKm1, rhoK,
260 O sigmaX, sigmaY, sigmaR,
261 I myThid )
262 ENDIF
263
264 #ifdef ALLOW_AUTODIFF_TAMC
265 CADJ STORE rhok (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
266 CADJ STORE rhokm1 (:,:) = comlev1_bibj_k , key=kkey, byte=isbyte
267 #endif /* ALLOW_AUTODIFF_TAMC */
268 C-- Implicit Vertical Diffusion for Convection
269 c ==> should use sigmaR !!!
270 IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
271 #ifdef ALLOW_DEBUG
272 IF ( debugLevel .GE. debLevB )
273 & CALL DEBUG_CALL('CALC_IVDC',myThid)
274 #endif
275 CALL CALC_IVDC(
276 I bi, bj, iMin, iMax, jMin, jMax, k,
277 I rhoKm1, rhoK,
278 I myTime, myIter, myThid)
279 ENDIF
280
281 C-- end of diagnostic k loop (Nr:1)
282 ENDDO
283
284 #ifdef ALLOW_AUTODIFF_TAMC
285 cph avoids recomputation of integrate_for_w
286 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
287 #endif /* ALLOW_AUTODIFF_TAMC */
288
289 #ifdef ALLOW_OBCS
290 C-- Calculate future values on open boundaries
291 IF (useOBCS) THEN
292 #ifdef ALLOW_DEBUG
293 IF ( debugLevel .GE. debLevB )
294 & CALL DEBUG_CALL('OBCS_CALC',myThid)
295 #endif
296 CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
297 I uVel, vVel, wVel, theta, salt,
298 I myThid )
299 ENDIF
300 #endif /* ALLOW_OBCS */
301
302 #ifndef ALLOW_AUTODIFF_TAMC
303 IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN
304 #endif
305 C-- Determines forcing terms based on external fields
306 C relaxation terms, etc.
307 #ifdef ALLOW_DEBUG
308 IF ( debugLevel .GE. debLevB )
309 & CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
310 #endif
311 CALL EXTERNAL_FORCING_SURF(
312 I bi, bj, iMin, iMax, jMin, jMax,
313 I myTime, myIter, myThid )
314 #ifndef ALLOW_AUTODIFF_TAMC
315 ENDIF
316 #endif
317
318 #ifdef ALLOW_AUTODIFF_TAMC
319 cph needed for KPP
320 CADJ STORE surfacetendencyU(:,:,bi,bj)
321 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
322 CADJ STORE surfacetendencyV(:,:,bi,bj)
323 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
324 CADJ STORE surfacetendencyS(:,:,bi,bj)
325 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
326 CADJ STORE surfacetendencyT(:,:,bi,bj)
327 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
328 # ifdef ALLOW_SEAICE
329 CADJ STORE surfacetendencyTice(:,:,bi,bj)
330 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
331 # endif
332 #endif /* ALLOW_AUTODIFF_TAMC */
333
334 #ifdef ALLOW_GMREDI
335
336 #ifdef ALLOW_AUTODIFF_TAMC
337 cph storing here is needed only for one GMREDI_OPTIONS:
338 cph define GM_BOLUS_ADVEC
339 cph but I've avoided the #ifdef for now, in case more things change
340 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
341 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
342 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
343 #endif /* ALLOW_AUTODIFF_TAMC */
344
345 C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
346 IF (useGMRedi) THEN
347 #ifdef ALLOW_DEBUG
348 IF ( debugLevel .GE. debLevB )
349 & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
350 #endif
351 CALL GMREDI_CALC_TENSOR(
352 I bi, bj, iMin, iMax, jMin, jMax,
353 I sigmaX, sigmaY, sigmaR,
354 I myThid )
355 #ifdef ALLOW_AUTODIFF_TAMC
356 ELSE
357 CALL GMREDI_CALC_TENSOR_DUMMY(
358 I bi, bj, iMin, iMax, jMin, jMax,
359 I sigmaX, sigmaY, sigmaR,
360 I myThid )
361 #endif /* ALLOW_AUTODIFF_TAMC */
362 ENDIF
363
364 #ifdef ALLOW_AUTODIFF_TAMC
365 CADJ STORE Kwx(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
366 CADJ STORE Kwy(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
367 CADJ STORE Kwz(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
368 #endif /* ALLOW_AUTODIFF_TAMC */
369
370 #endif /* ALLOW_GMREDI */
371
372 #ifdef ALLOW_KPP
373 C-- Compute KPP mixing coefficients
374 IF (useKPP) THEN
375 #ifdef ALLOW_DEBUG
376 IF ( debugLevel .GE. debLevB )
377 & CALL DEBUG_CALL('KPP_CALC',myThid)
378 #endif
379 CALL KPP_CALC(
380 I bi, bj, myTime, myThid )
381 #ifdef ALLOW_AUTODIFF_TAMC
382 ELSE
383 CALL KPP_CALC_DUMMY(
384 I bi, bj, myTime, myThid )
385 #endif /* ALLOW_AUTODIFF_TAMC */
386 ENDIF
387
388 #ifdef ALLOW_AUTODIFF_TAMC
389 CADJ STORE KPPghat (:,:,:,bi,bj)
390 CADJ & , KPPfrac (:,: ,bi,bj)
391 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
392 #endif /* ALLOW_AUTODIFF_TAMC */
393
394 #endif /* ALLOW_KPP */
395
396 #ifdef ALLOW_AUTODIFF_TAMC
397 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
398 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
399 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
400 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
401 #ifdef ALLOW_PASSIVE_TRACER
402 CADJ STORE tr1 (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
403 #endif
404 #ifdef ALLOW_PTRACERS
405 cph-- moved to forward_step to avoid key computation
406 cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
407 cphCADJ & key=itdkey, byte=isbyte
408 #endif
409 #endif /* ALLOW_AUTODIFF_TAMC */
410
411 C-- end bi,bj loops.
412 ENDDO
413 ENDDO
414
415 #ifdef ALLOW_DEBUG
416 IF ( debugLevel .GE. debLevB )
417 & CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
418 #endif
419
420 RETURN
421 END

  ViewVC Help
Powered by ViewVC 1.1.22