/[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.4 - (show annotations) (download)
Sun Jul 18 01:04:23 2004 UTC (19 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint54c_post
Changes since 1.3: +6 -6 lines
replace surfaceTendency U,V,S,T,Tice by surfaceForcing U,V,S,T,Tice

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

  ViewVC Help
Powered by ViewVC 1.1.22