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

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

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


Revision 1.99 - (hide annotations) (download)
Thu Jan 27 20:34:53 2011 UTC (13 years, 4 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62s
Changes since 1.98: +7 -1 lines
add call to new S/R GMREDI_DO_EXCH

1 jmc 1.99 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.98 2011/01/26 03:46:01 jmc Exp $
2 jmc 1.1 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 jmc 1.29 # ifdef ALLOW_SEAICE
15     # include "SEAICE_OPTIONS.h"
16     # endif
17 jmc 1.1 #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 jmc 1.28 C | SUBROUTINE DO_OCEANIC_PHYS
26     C | o Controlling routine for oceanic physics and
27 jmc 1.1 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 jmc 1.69 #include "GRID.h"
40 jmc 1.1 #include "DYNVARS.h"
41 jmc 1.20 #ifdef ALLOW_TIMEAVE
42     #include "TIMEAVE_STATV.h"
43     #endif
44 mlosch 1.22 #if defined (ALLOW_BALANCE_FLUXES) && !(defined ALLOW_AUTODIFF_TAMC)
45     #include "FFIELDS.h"
46     #endif
47 jmc 1.1
48     #ifdef ALLOW_AUTODIFF_TAMC
49 heimbach 1.80 # include "AUTODIFF_MYFIELDS.h"
50 jmc 1.1 # include "tamc.h"
51     # include "tamc_keys.h"
52     # include "FFIELDS.h"
53 heimbach 1.54 # include "SURFACE.h"
54 jmc 1.1 # include "EOS.h"
55     # ifdef ALLOW_KPP
56     # include "KPP.h"
57     # endif
58 gforget 1.91 # ifdef ALLOW_GGL90
59     # include "GGL90.h"
60     # endif
61 jmc 1.1 # ifdef ALLOW_GMREDI
62     # include "GMREDI.h"
63     # endif
64     # ifdef ALLOW_EBM
65     # include "EBM.h"
66     # endif
67 jmc 1.29 # ifdef ALLOW_EXF
68     # include "ctrl.h"
69 jmc 1.40 # include "EXF_FIELDS.h"
70 jmc 1.29 # ifdef ALLOW_BULKFORMULAE
71 jmc 1.40 # include "EXF_CONSTANTS.h"
72 jmc 1.29 # endif
73     # endif
74     # ifdef ALLOW_SEAICE
75     # include "SEAICE.h"
76     # endif
77 heimbach 1.75 # ifdef ALLOW_SALT_PLUME
78     # include "SALT_PLUME.h"
79     # endif
80 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
81    
82     C !INPUT/OUTPUT PARAMETERS:
83     C == Routine arguments ==
84 jmc 1.18 C myTime :: Current time in simulation
85     C myIter :: Current iteration number in simulation
86     C myThid :: Thread number for this instance of the routine.
87 jmc 1.1 _RL myTime
88     INTEGER myIter
89     INTEGER myThid
90    
91     C !LOCAL VARIABLES:
92     C == Local variables
93 jmc 1.47 C rhoK, rhoKm1 :: Density at current level, and level above
94 jmc 1.18 C iMin, iMax :: Ranges and sub-block indices on which calculations
95 jmc 1.1 C jMin, jMax are applied.
96 jmc 1.18 C bi, bj :: tile indices
97     C i,j,k :: loop indices
98 jmc 1.47 _RL rhoKp1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99     _RL rhoKm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100 jmc 1.1 _RL sigmaX (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101     _RL sigmaY (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
102     _RL sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
103     INTEGER iMin, iMax
104     INTEGER jMin, jMax
105     INTEGER bi, bj
106 jmc 1.18 INTEGER i, j, k
107 jmc 1.17 INTEGER doDiagsRho
108     #ifdef ALLOW_DIAGNOSTICS
109     LOGICAL DIAGNOSTICS_IS_ON
110     EXTERNAL DIAGNOSTICS_IS_ON
111     #endif /* ALLOW_DIAGNOSTICS */
112 jmc 1.1
113     CEOP
114    
115 heimbach 1.12 #ifdef ALLOW_AUTODIFF_TAMC
116     C-- dummy statement to end declaration part
117     itdkey = 1
118     #endif /* ALLOW_AUTODIFF_TAMC */
119    
120 jmc 1.1 #ifdef ALLOW_DEBUG
121 jmc 1.96 IF (debugMode) CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
122 jmc 1.1 #endif
123 jmc 1.36
124 jmc 1.17 doDiagsRho = 0
125     #ifdef ALLOW_DIAGNOSTICS
126     IF ( useDiagnostics .AND. fluidIsWater ) THEN
127 jmc 1.74 IF ( DIAGNOSTICS_IS_ON('WRHOMASS',myThid) )
128 jmc 1.71 & doDiagsRho = doDiagsRho + 1
129     IF ( DIAGNOSTICS_IS_ON('DRHODR ',myThid) )
130     & doDiagsRho = doDiagsRho + 2
131     IF ( DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) )
132     & doDiagsRho = doDiagsRho + 4
133 jmc 1.17 ENDIF
134     #endif /* ALLOW_DIAGNOSTICS */
135    
136 jmc 1.82 #ifdef ALLOW_OBCS
137     IF (useOBCS) THEN
138     C-- Calculate future values on open boundaries
139     C-- moved before SEAICE_MODEL call since SEAICE_MODEL needs seaice-obcs fields
140     #ifdef ALLOW_DEBUG
141 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('OBCS_CALC',myThid)
142 jmc 1.82 #endif
143 jahn 1.83 CALL OBCS_CALC( myTime+deltaTclock, myIter+1,
144     I uVel, vVel, wVel, theta, salt, myThid )
145 jmc 1.82 ENDIF
146     #endif /* ALLOW_OBCS */
147 jmc 1.69
148 jmc 1.98 #ifdef ALLOW_ADDFLUID
149     IF ( fluidIsWater ) THEN
150     DO bj=myByLo(myThid),myByHi(myThid)
151     DO bi=myBxLo(myThid),myBxHi(myThid)
152     DO k=1,Nr
153     DO j=1-OLy,sNy+OLy
154     DO i=1-OLx,sNx+OLx
155     addMass(i,j,k,bi,bj) = 0. _d 0
156     ENDDO
157     ENDDO
158     ENDDO
159     ENDDO
160     ENDDO
161     ENDIF
162     #endif /* ALLOW_ADDFLUID */
163    
164 gforget 1.87 #ifdef ALLOW_AUTODIFF_TAMC
165     # ifdef ALLOW_SALT_PLUME
166     DO bj=myByLo(myThid),myByHi(myThid)
167     DO bi=myBxLo(myThid),myBxHi(myThid)
168     DO j=1-OLy,sNy+OLy
169     DO i=1-OLx,sNx+OLx
170     saltPlumeDepth(i,j,bi,bj) = 0. _d 0
171     saltPlumeFlux(i,j,bi,bj) = 0. _d 0
172     ENDDO
173     ENDDO
174     ENDDO
175     ENDDO
176     # endif
177     #endif /* ALLOW_AUTODIFF_TAMC */
178    
179 jmc 1.29 #ifdef ALLOW_SEAICE
180     IF ( useSEAICE ) THEN
181 heimbach 1.62 # ifdef ALLOW_AUTODIFF_TAMC
182 heimbach 1.65 cph-adj-test(
183 heimbach 1.81 CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
184     CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte
185 heimbach 1.88 CADJ STORE heff = comlev1, key=ikey_dynamics, kind=isbyte
186 heimbach 1.81 CADJ STORE empmr,qsw,theta = comlev1, key = ikey_dynamics,
187 heimbach 1.77 CADJ & kind = isbyte
188 heimbach 1.65 cph-adj-test)
189 heimbach 1.77 CADJ STORE atemp,aqh,precip = comlev1, key = ikey_dynamics,
190     CADJ & kind = isbyte
191     CADJ STORE swdown,lwdown = comlev1, key = ikey_dynamics,
192     CADJ & kind = isbyte
193 heimbach 1.34 cph# ifdef EXF_READ_EVAP
194 heimbach 1.77 CADJ STORE evap = comlev1, key = ikey_dynamics,
195     CADJ & kind = isbyte
196 heimbach 1.34 cph# endif
197 heimbach 1.77 CADJ STORE uvel,vvel = comlev1, key = ikey_dynamics,
198     CADJ & kind = isbyte
199 heimbach 1.95 # ifdef SEAICE_CGRID
200 heimbach 1.93 CADJ STORE stressdivergencex = comlev1, key = ikey_dynamics,
201     CADJ & kind = isbyte
202     CADJ STORE stressdivergencey = comlev1, key = ikey_dynamics,
203     CADJ & kind = isbyte
204     # endif
205 heimbach 1.62 # ifdef SEAICE_ALLOW_DYNAMICS
206 heimbach 1.77 CADJ STORE uice = comlev1, key = ikey_dynamics,
207     CADJ & kind = isbyte
208     CADJ STORE vice = comlev1, key = ikey_dynamics,
209     CADJ & kind = isbyte
210 heimbach 1.62 # ifdef SEAICE_ALLOW_EVP
211 heimbach 1.77 CADJ STORE seaice_sigma1 = comlev1, key = ikey_dynamics,
212     CADJ & kind = isbyte
213     CADJ STORE seaice_sigma2 = comlev1, key = ikey_dynamics,
214     CADJ & kind = isbyte
215     CADJ STORE seaice_sigma12 = comlev1, key = ikey_dynamics,
216     CADJ & kind = isbyte
217 heimbach 1.62 # endif
218     # endif
219     # ifdef SEAICE_SALINITY
220 heimbach 1.77 CADJ STORE salt = comlev1, key = ikey_dynamics,
221     CADJ & kind = isbyte
222 heimbach 1.62 # endif
223     # ifdef ATMOSPHERIC_LOADING
224 heimbach 1.77 CADJ STORE pload = comlev1, key = ikey_dynamics,
225     CADJ & kind = isbyte
226     CADJ STORE siceload = comlev1, key = ikey_dynamics,
227     CADJ & kind = isbyte
228 heimbach 1.62 # endif
229     # ifdef NONLIN_FRSURF
230 heimbach 1.77 CADJ STORE recip_hfacc = comlev1, key = ikey_dynamics,
231     CADJ & kind = isbyte
232 heimbach 1.62 # endif
233 heimbach 1.78 # ifdef ANNUAL_BALANCE
234     CADJ STORE balance_itcount = comlev1, key = ikey_dynamics,
235     CADJ & kind = isbyte
236     # endif /* ANNUAL_BALANCE */
237 heimbach 1.55 # endif
238 heimbach 1.62 # ifdef ALLOW_DEBUG
239 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('SEAICE_MODEL',myThid)
240 heimbach 1.62 # endif
241 jmc 1.29 CALL TIMER_START('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid)
242     CALL SEAICE_MODEL( myTime, myIter, myThid )
243     CALL TIMER_STOP ('SEAICE_MODEL [DO_OCEANIC_PHYS]', myThid)
244 heimbach 1.62 # ifdef ALLOW_COST
245 heimbach 1.57 CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
246 heimbach 1.62 # endif
247 heimbach 1.35 ENDIF
248 jmc 1.29 #endif /* ALLOW_SEAICE */
249    
250 heimbach 1.64 #ifdef ALLOW_AUTODIFF_TAMC
251 heimbach 1.77 CADJ STORE sst, sss = comlev1, key = ikey_dynamics,
252     CADJ & kind = isbyte
253     CADJ STORE qsw = comlev1, key = ikey_dynamics,
254     CADJ & kind = isbyte
255 heimbach 1.64 # ifdef ALLOW_SEAICE
256 heimbach 1.77 CADJ STORE area = comlev1, key = ikey_dynamics,
257     CADJ & kind = isbyte
258 heimbach 1.64 # endif
259     #endif
260    
261 jscott 1.30 #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
262 jmc 1.14 IF ( useThSIce .AND. fluidIsWater ) THEN
263 jmc 1.5 #ifdef ALLOW_DEBUG
264 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
265 jmc 1.5 #endif
266     C-- Step forward Therm.Sea-Ice variables
267     C and modify forcing terms including effects from ice
268     CALL TIMER_START('THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
269     CALL THSICE_MAIN( myTime, myIter, myThid )
270     CALL TIMER_STOP( 'THSICE_MAIN [DO_OCEANIC_PHYS]', myThid)
271     ENDIF
272     #endif /* ALLOW_THSICE */
273    
274 mlosch 1.21 #ifdef ALLOW_SHELFICE
275 heimbach 1.92 # ifdef ALLOW_AUTODIFF_TAMC
276     CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
277     CADJ & kind = isbyte
278     # endif
279 mlosch 1.21 IF ( useShelfIce .AND. fluidIsWater ) THEN
280     #ifdef ALLOW_DEBUG
281 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
282 mlosch 1.21 #endif
283 jmc 1.47 C compute temperature and (virtual) salt flux at the
284 mlosch 1.21 C shelf-ice ocean interface
285     CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
286     & myThid)
287     CALL SHELFICE_THERMODYNAMICS( myTime, myIter, myThid )
288     CALL TIMER_STOP( 'SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
289     & myThid)
290     ENDIF
291     #endif /* ALLOW_SHELFICE */
292    
293 dimitri 1.85 #ifdef ALLOW_ICEFRONT
294     IF ( useICEFRONT .AND. fluidIsWater ) THEN
295     #ifdef ALLOW_DEBUG
296 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('ICEFRONT_THERMODYNAMICS',myThid)
297 dimitri 1.85 #endif
298     C compute temperature and (virtual) salt flux at the
299     C ice-front ocean interface
300     CALL TIMER_START('ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
301     & myThid)
302     CALL ICEFRONT_THERMODYNAMICS( myTime, myIter, myThid )
303     CALL TIMER_STOP( 'ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
304     & myThid)
305     ENDIF
306     #endif /* ALLOW_ICEFRONT */
307    
308 jmc 1.5 C-- Freeze water at the surface
309     #ifdef ALLOW_AUTODIFF_TAMC
310 heimbach 1.77 CADJ STORE theta = comlev1, key = ikey_dynamics,
311     CADJ & kind = isbyte
312 jmc 1.5 #endif
313 jmc 1.70 IF ( allowFreezing ) THEN
314 jmc 1.5 CALL FREEZE_SURFACE( myTime, myIter, myThid )
315     ENDIF
316    
317 jmc 1.28 #ifdef ALLOW_OCN_COMPON_INTERF
318 jmc 1.5 C-- Apply imported data (from coupled interface) to forcing fields
319 jmc 1.28 C jmc: do not know precisely where to put this call (bf or af thSIce ?)
320 jmc 1.36 IF ( useCoupler ) THEN
321 jmc 1.5 CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
322 jmc 1.36 ENDIF
323 jmc 1.28 #endif /* ALLOW_OCN_COMPON_INTERF */
324 jmc 1.5
325 jmc 1.25 #ifdef ALLOW_BALANCE_FLUXES
326 jmc 1.36 C balance fluxes
327     IF ( balanceEmPmR )
328 jmc 1.84 & CALL REMOVE_MEAN_RS( 1, EmPmR, maskInC, maskInC, rA, drF,
329 jmc 1.25 & 'EmPmR', myTime, myThid )
330 jmc 1.36 IF ( balanceQnet )
331 jmc 1.84 & CALL REMOVE_MEAN_RS( 1, Qnet, maskInC, maskInC, rA, drF,
332 jmc 1.25 & 'Qnet ', myTime, myThid )
333     #endif /* ALLOW_BALANCE_FLUXES */
334    
335 jmc 1.1 #ifdef ALLOW_AUTODIFF_TAMC
336     C-- HPF directive to help TAMC
337     CHPF$ INDEPENDENT
338     #endif /* ALLOW_AUTODIFF_TAMC */
339     DO bj=myByLo(myThid),myByHi(myThid)
340     #ifdef ALLOW_AUTODIFF_TAMC
341 heimbach 1.15 C-- HPF directive to help TAMC
342     CHPF$ INDEPENDENT
343 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
344     DO bi=myBxLo(myThid),myBxHi(myThid)
345    
346     #ifdef ALLOW_AUTODIFF_TAMC
347     act1 = bi - myBxLo(myThid)
348     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
349     act2 = bj - myByLo(myThid)
350     max2 = myByHi(myThid) - myByLo(myThid) + 1
351     act3 = myThid - 1
352     max3 = nTx*nTy
353     act4 = ikey_dynamics - 1
354     itdkey = (act1 + 1) + act2*max1
355     & + act3*max1*max2
356     & + act4*max1*max2*max3
357 jmc 1.74 #else /* ALLOW_AUTODIFF_TAMC */
358 jmc 1.47 C if fluid is not water, by-pass find_rho, gmredi, surfaceForcing
359 jmc 1.36 C and all vertical mixing schemes, but keep OBCS_CALC
360     IF ( fluidIsWater ) THEN
361 jmc 1.74 #endif /* ALLOW_AUTODIFF_TAMC */
362 jmc 1.1
363     C-- Set up work arrays with valid (i.e. not NaN) values
364     C These inital values do not alter the numerical results. They
365     C just ensure that all memory references are to valid floating
366     C point numbers. This prevents spurious hardware signals due to
367     C uninitialised but inert locations.
368    
369 jmc 1.69 #ifdef ALLOW_AUTODIFF_TAMC
370 jmc 1.1 DO j=1-OLy,sNy+OLy
371     DO i=1-OLx,sNx+OLx
372 jmc 1.69 rhoKm1 (i,j) = 0. _d 0
373 jmc 1.47 rhoKp1 (i,j) = 0. _d 0
374 jmc 1.1 ENDDO
375     ENDDO
376 jmc 1.69 #endif /* ALLOW_AUTODIFF_TAMC */
377 jmc 1.1
378     DO k=1,Nr
379     DO j=1-OLy,sNy+OLy
380     DO i=1-OLx,sNx+OLx
381 jmc 1.69 C This is currently used by GMRedi, IVDC, MXL-depth and Diagnostics
382 jmc 1.1 sigmaX(i,j,k) = 0. _d 0
383     sigmaY(i,j,k) = 0. _d 0
384     sigmaR(i,j,k) = 0. _d 0
385     #ifdef ALLOW_AUTODIFF_TAMC
386     cph all the following init. are necessary for TAF
387     cph although some of these are re-initialised later.
388 jmc 1.71 c rhoInSitu(i,j,k,bi,bj) = 0.
389 jmc 1.1 IVDConvCount(i,j,k,bi,bj) = 0.
390     # ifdef ALLOW_GMREDI
391     Kwx(i,j,k,bi,bj) = 0. _d 0
392     Kwy(i,j,k,bi,bj) = 0. _d 0
393     Kwz(i,j,k,bi,bj) = 0. _d 0
394     # ifdef GM_NON_UNITY_DIAGONAL
395     Kux(i,j,k,bi,bj) = 0. _d 0
396     Kvy(i,j,k,bi,bj) = 0. _d 0
397     # endif
398     # ifdef GM_EXTRA_DIAGONAL
399     Kuz(i,j,k,bi,bj) = 0. _d 0
400     Kvz(i,j,k,bi,bj) = 0. _d 0
401     # endif
402     # ifdef GM_BOLUS_ADVEC
403     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
404     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
405     # endif
406     # ifdef GM_VISBECK_VARIABLE_K
407     VisbeckK(i,j,bi,bj) = 0. _d 0
408     # endif
409     # endif /* ALLOW_GMREDI */
410 heimbach 1.42 # ifdef ALLOW_KPP
411     KPPdiffKzS(i,j,k,bi,bj) = 0. _d 0
412     KPPdiffKzT(i,j,k,bi,bj) = 0. _d 0
413     # endif /* ALLOW_KPP */
414 gforget 1.91 # ifdef ALLOW_GGL90
415     GGL90viscArU(i,j,k,bi,bj) = 0. _d 0
416     GGL90viscArV(i,j,k,bi,bj) = 0. _d 0
417     GGL90diffKr(i,j,k,bi,bj) = 0. _d 0
418     # endif /* ALLOW_GGL90 */
419 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
420     ENDDO
421     ENDDO
422     ENDDO
423    
424     iMin = 1-OLx
425     iMax = sNx+OLx
426     jMin = 1-OLy
427     jMax = sNy+OLy
428    
429     #ifdef ALLOW_AUTODIFF_TAMC
430 jmc 1.96 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
431 heimbach 1.77 CADJ & kind = isbyte
432 jmc 1.96 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
433 heimbach 1.77 CADJ & kind = isbyte
434 heimbach 1.12 CADJ STORE totphihyd(:,:,:,bi,bj)
435 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
436 heimbach 1.77 CADJ & kind = isbyte
437 heimbach 1.10 # ifdef ALLOW_KPP
438 jmc 1.96 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
439 heimbach 1.77 CADJ & kind = isbyte
440 jmc 1.96 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
441 heimbach 1.77 CADJ & kind = isbyte
442 heimbach 1.10 # endif
443 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
444    
445     #ifdef ALLOW_DEBUG
446 jmc 1.96 IF (debugMode) CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
447 jmc 1.1 #endif
448    
449     C-- Start of diagnostic loop
450     DO k=Nr,1,-1
451    
452     #ifdef ALLOW_AUTODIFF_TAMC
453     C? Patrick, is this formula correct now that we change the loop range?
454     C? Do we still need this?
455 jmc 1.69 cph kkey formula corrected.
456 jmc 1.47 cph Needed for rhoK, rhoKm1, in the case useGMREDI.
457 jmc 1.71 kkey = (itdkey-1)*Nr + k
458 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
459    
460 jmc 1.71 C-- Always compute density (stored in common block) here; even when it is not
461     C needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
462     #ifdef ALLOW_DEBUG
463 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D',myThid)
464 jmc 1.71 #endif
465     #ifdef ALLOW_AUTODIFF_TAMC
466 jmc 1.74 IF ( fluidIsWater ) THEN
467 jmc 1.96 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
468 heimbach 1.77 CADJ & kind = isbyte
469 jmc 1.96 CADJ STORE salt(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
470 heimbach 1.77 CADJ & kind = isbyte
471 jmc 1.71 #endif /* ALLOW_AUTODIFF_TAMC */
472 jmc 1.69 #ifdef ALLOW_DOWN_SLOPE
473     IF ( useDOWN_SLOPE ) THEN
474     CALL DWNSLP_CALC_RHO(
475     I theta, salt,
476 jmc 1.71 O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
477 jmc 1.69 I k, bi, bj, myTime, myIter, myThid )
478 jmc 1.71 ELSE
479     #endif /* ALLOW_DOWN_SLOPE */
480     CALL FIND_RHO_2D(
481     I iMin, iMax, jMin, jMax, k,
482     I theta(1-OLx,1-OLy,k,bi,bj),
483     I salt (1-OLx,1-OLy,k,bi,bj),
484     O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
485     I k, bi, bj, myThid )
486     #ifdef ALLOW_DOWN_SLOPE
487 jmc 1.69 ENDIF
488     #endif /* ALLOW_DOWN_SLOPE */
489 jmc 1.74 #ifdef ALLOW_AUTODIFF_TAMC
490     ELSE
491     C- fluid is not water:
492     DO j=1-OLy,sNy+OLy
493     DO i=1-OLx,sNx+OLx
494     rhoInSitu(i,j,k,bi,bj) = 0.
495     ENDDO
496     ENDDO
497     ENDIF
498     #endif /* ALLOW_AUTODIFF_TAMC */
499 jmc 1.69
500 jmc 1.1 C-- Calculate gradients of potential density for isoneutral
501     C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
502 jmc 1.17 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)
503 dimitri 1.61 & .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
504 jmc 1.1 IF (k.GT.1) THEN
505     #ifdef ALLOW_AUTODIFF_TAMC
506 jmc 1.96 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
507 heimbach 1.77 CADJ & kind = isbyte
508 jmc 1.96 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
509 heimbach 1.77 CADJ & kind = isbyte
510 jmc 1.96 CADJ STORE rhokm1 (bi,bj) = comlev1_bibj_k, key=kkey,
511 heimbach 1.77 CADJ & kind = isbyte
512 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
513 jmc 1.68 CALL FIND_RHO_2D(
514     I iMin, iMax, jMin, jMax, k,
515     I theta(1-OLx,1-OLy,k-1,bi,bj),
516     I salt (1-OLx,1-OLy,k-1,bi,bj),
517     O rhoKm1,
518     I k-1, bi, bj, myThid )
519 jmc 1.1 ENDIF
520     #ifdef ALLOW_DEBUG
521 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
522 jmc 1.1 #endif
523 heimbach 1.31 cph Avoid variable aliasing for adjoint !!!
524     DO j=jMin,jMax
525     DO i=iMin,iMax
526 jmc 1.71 rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
527 heimbach 1.31 ENDDO
528     ENDDO
529 jmc 1.1 CALL GRAD_SIGMA(
530     I bi, bj, iMin, iMax, jMin, jMax, k,
531 jmc 1.71 I rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
532 jmc 1.1 O sigmaX, sigmaY, sigmaR,
533     I myThid )
534 gforget 1.66 #ifdef ALLOW_AUTODIFF_TAMC
535 jmc 1.69 #ifdef GMREDI_WITH_STABLE_ADJOINT
536 gforget 1.66 cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
537     cgf -> cuts adjoint dependency from slope to state
538 jmc 1.69 CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
539     CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
540     CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
541 gforget 1.66 #endif
542     #endif /* ALLOW_AUTODIFF_TAMC */
543 jmc 1.1 ENDIF
544    
545     C-- Implicit Vertical Diffusion for Convection
546     c ==> should use sigmaR !!!
547     IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
548     #ifdef ALLOW_DEBUG
549 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid)
550 jmc 1.1 #endif
551     CALL CALC_IVDC(
552     I bi, bj, iMin, iMax, jMin, jMax, k,
553 jmc 1.71 I rhoKm1, rhoInSitu(1-OLx,1-OLy,k,bi,bj),
554 jmc 1.1 I myTime, myIter, myThid)
555     ENDIF
556    
557 jmc 1.17 #ifdef ALLOW_DIAGNOSTICS
558 jmc 1.71 IF ( MOD(doDiagsRho,2).EQ.1 ) THEN
559     CALL DIAGS_RHO_L( k, bi, bj,
560     I rhoInSitu(1-OLx,1-OLy,k,bi,bj),
561 jmc 1.74 I rhoKm1, wVel,
562 jmc 1.71 I myTime, myIter, myThid )
563 jmc 1.17 ENDIF
564     #endif
565    
566 jmc 1.1 C-- end of diagnostic k loop (Nr:1)
567     ENDDO
568    
569 heimbach 1.57 #ifdef ALLOW_AUTODIFF_TAMC
570 jmc 1.69 CADJ STORE IVDConvCount(:,:,:,bi,bj)
571 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
572 heimbach 1.77 CADJ & kind = isbyte
573 heimbach 1.57 #endif
574    
575 jmc 1.47 C-- Diagnose Mixed Layer Depth:
576 jmc 1.71 IF ( useGMRedi .OR. doDiagsRho.GE.4 ) THEN
577     CALL CALC_OCE_MXLAYER(
578     I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
579     I bi, bj, myTime, myIter, myThid )
580 jmc 1.47 ENDIF
581 heimbach 1.53
582 dimitri 1.52 #ifdef ALLOW_SALT_PLUME
583 dimitri 1.61 IF ( useSALT_PLUME ) THEN
584 jmc 1.71 CALL SALT_PLUME_CALC_DEPTH(
585     I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
586     I bi, bj, myTime, myIter, myThid )
587 dimitri 1.60 ENDIF
588 dimitri 1.61 #endif /* ALLOW_SALT_PLUME */
589    
590 jmc 1.8 #ifdef ALLOW_DIAGNOSTICS
591 jmc 1.71 IF ( MOD(doDiagsRho,4).GE.2 ) THEN
592 jmc 1.16 CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR ', 0, Nr,
593     & 2, bi, bj, myThid)
594 jmc 1.8 ENDIF
595 dimitri 1.61 #endif /* ALLOW_DIAGNOSTICS */
596 jmc 1.8
597 jmc 1.1 C-- Determines forcing terms based on external fields
598     C relaxation terms, etc.
599     #ifdef ALLOW_DEBUG
600 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
601 jmc 1.1 #endif
602 heimbach 1.23 #ifdef ALLOW_AUTODIFF_TAMC
603     CADJ STORE EmPmR(:,:,bi,bj)
604 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
605 heimbach 1.77 CADJ & kind = isbyte
606 heimbach 1.26 # ifdef EXACT_CONSERV
607 heimbach 1.23 CADJ STORE PmEpR(:,:,bi,bj)
608 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
609 heimbach 1.77 CADJ & kind = isbyte
610 heimbach 1.26 # endif
611 heimbach 1.27 # ifdef NONLIN_FRSURF
612     CADJ STORE hFac_surfC(:,:,bi,bj)
613 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
614 heimbach 1.77 CADJ & kind = isbyte
615 heimbach 1.27 CADJ STORE recip_hFacC(:,:,:,bi,bj)
616 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
617 heimbach 1.77 CADJ & kind = isbyte
618 heimbach 1.97 # if (defined (ALLOW_PTRACERS))
619     CADJ STORE surfaceForcingS(:,:,bi,bj) = comlev1_bibj, key=itdkey,
620     CADJ & kind = isbyte
621     CADJ STORE surfaceForcingT(:,:,bi,bj) = comlev1_bibj, key=itdkey,
622     CADJ & kind = isbyte
623     # endif /* ALLOW_PTRACERS */
624     # endif /* NONLIN_FRSURF */
625 heimbach 1.23 #endif
626 jmc 1.36 CALL EXTERNAL_FORCING_SURF(
627 jmc 1.1 I bi, bj, iMin, iMax, jMin, jMax,
628     I myTime, myIter, myThid )
629 heimbach 1.27 #ifdef ALLOW_AUTODIFF_TAMC
630     # ifdef EXACT_CONSERV
631     cph-test
632     cphCADJ STORE PmEpR(:,:,bi,bj)
633 jmc 1.96 cphCADJ & = comlev1_bibj, key=itdkey,
634 heimbach 1.77 cphCADJ & kind = isbyte
635 heimbach 1.27 # endif
636     #endif
637 jmc 1.1
638     #ifdef ALLOW_AUTODIFF_TAMC
639     cph needed for KPP
640 jmc 1.4 CADJ STORE surfaceForcingU(:,:,bi,bj)
641 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
642 heimbach 1.77 CADJ & kind = isbyte
643 jmc 1.4 CADJ STORE surfaceForcingV(:,:,bi,bj)
644 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
645 heimbach 1.77 CADJ & kind = isbyte
646 jmc 1.4 CADJ STORE surfaceForcingS(:,:,bi,bj)
647 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
648 heimbach 1.77 CADJ & kind = isbyte
649 jmc 1.4 CADJ STORE surfaceForcingT(:,:,bi,bj)
650 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
651 heimbach 1.77 CADJ & kind = isbyte
652 jmc 1.4 CADJ STORE surfaceForcingTice(:,:,bi,bj)
653 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
654 heimbach 1.77 CADJ & kind = isbyte
655 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
656    
657     #ifdef ALLOW_KPP
658     C-- Compute KPP mixing coefficients
659     IF (useKPP) THEN
660     #ifdef ALLOW_DEBUG
661 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid)
662 jmc 1.1 #endif
663 dfer 1.76 CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
664 jmc 1.1 CALL KPP_CALC(
665 jmc 1.44 I bi, bj, myTime, myIter, myThid )
666 dfer 1.76 CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
667 jmc 1.1 #ifdef ALLOW_AUTODIFF_TAMC
668     ELSE
669     CALL KPP_CALC_DUMMY(
670 jmc 1.44 I bi, bj, myTime, myIter, myThid )
671 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
672     ENDIF
673    
674     #endif /* ALLOW_KPP */
675    
676 mlosch 1.6 #ifdef ALLOW_PP81
677     C-- Compute PP81 mixing coefficients
678     IF (usePP81) THEN
679     #ifdef ALLOW_DEBUG
680 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid)
681 mlosch 1.6 #endif
682     CALL PP81_CALC(
683     I bi, bj, myTime, myThid )
684     ENDIF
685     #endif /* ALLOW_PP81 */
686    
687     #ifdef ALLOW_MY82
688     C-- Compute MY82 mixing coefficients
689     IF (useMY82) THEN
690     #ifdef ALLOW_DEBUG
691 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid)
692 mlosch 1.6 #endif
693     CALL MY82_CALC(
694     I bi, bj, myTime, myThid )
695     ENDIF
696     #endif /* ALLOW_MY82 */
697    
698 mlosch 1.9 #ifdef ALLOW_GGL90
699 gforget 1.91 #ifdef ALLOW_AUTODIFF_TAMC
700     CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
701     CADJ & kind = isbyte
702     #endif /* ALLOW_AUTODIFF_TAMC */
703 mlosch 1.9 C-- Compute GGL90 mixing coefficients
704     IF (useGGL90) THEN
705     #ifdef ALLOW_DEBUG
706 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
707 mlosch 1.9 #endif
708 dfer 1.76 CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
709 mlosch 1.9 CALL GGL90_CALC(
710 jmc 1.90 I bi, bj, myTime, myIter, myThid )
711 dfer 1.76 CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
712 mlosch 1.9 ENDIF
713     #endif /* ALLOW_GGL90 */
714    
715 jmc 1.20 #ifdef ALLOW_TIMEAVE
716 jmc 1.36 IF ( taveFreq.GT. 0. _d 0 ) THEN
717 jmc 1.20 CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
718     ENDIF
719     IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
720     CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
721     I Nr, deltaTclock, bi, bj, myThid)
722     ENDIF
723     #endif /* ALLOW_TIMEAVE */
724    
725 jmc 1.69 #ifdef ALLOW_GMREDI
726 jmc 1.47 #ifdef ALLOW_AUTODIFF_TAMC
727     # ifndef GM_EXCLUDE_CLIPPING
728     cph storing here is needed only for one GMREDI_OPTIONS:
729     cph define GM_BOLUS_ADVEC
730     cph keep it although TAF says you dont need to.
731 jmc 1.86 cph but I have avoided the #ifdef for now, in case more things change
732 jmc 1.96 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey,
733 heimbach 1.77 CADJ & kind = isbyte
734 jmc 1.96 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey,
735 heimbach 1.77 CADJ & kind = isbyte
736 jmc 1.96 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey,
737 heimbach 1.77 CADJ & kind = isbyte
738 jmc 1.47 # endif
739     #endif /* ALLOW_AUTODIFF_TAMC */
740    
741     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
742     IF (useGMRedi) THEN
743     #ifdef ALLOW_DEBUG
744 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
745 jmc 1.47 #endif
746     CALL GMREDI_CALC_TENSOR(
747 jmc 1.51 I iMin, iMax, jMin, jMax,
748 jmc 1.47 I sigmaX, sigmaY, sigmaR,
749 jmc 1.51 I bi, bj, myTime, myIter, myThid )
750 jmc 1.47 #ifdef ALLOW_AUTODIFF_TAMC
751     ELSE
752     CALL GMREDI_CALC_TENSOR_DUMMY(
753 jmc 1.51 I iMin, iMax, jMin, jMax,
754 jmc 1.47 I sigmaX, sigmaY, sigmaR,
755 jmc 1.51 I bi, bj, myTime, myIter, myThid )
756 jmc 1.47 #endif /* ALLOW_AUTODIFF_TAMC */
757     ENDIF
758 jmc 1.69 #endif /* ALLOW_GMREDI */
759    
760     #ifdef ALLOW_DOWN_SLOPE
761     IF ( useDOWN_SLOPE ) THEN
762     C-- Calculate Downsloping Flow for Down_Slope parameterization
763     IF ( usingPCoords ) THEN
764     CALL DWNSLP_CALC_FLOW(
765 jmc 1.71 I bi, bj, kSurfC, rhoInSitu,
766 jmc 1.69 I myTime, myIter, myThid )
767     ELSE
768     CALL DWNSLP_CALC_FLOW(
769 jmc 1.71 I bi, bj, kLowC, rhoInSitu,
770 jmc 1.69 I myTime, myIter, myThid )
771     ENDIF
772     ENDIF
773     #endif /* ALLOW_DOWN_SLOPE */
774 jmc 1.47
775 jmc 1.98 #ifndef ALLOW_AUTODIFF_TAMC
776     C--- if fluid Is Water: end
777     ENDIF
778     #endif
779    
780 dimitri 1.94 #ifdef ALLOW_MYPACKAGE
781     IF ( useMYPACKAGE ) THEN
782     CALL MYPACKAGE_CALC_RHS(
783     I bi, bj, myTime, myIter, myThid )
784     ENDIF
785     #endif /* ALLOW_MYPACKAGE */
786    
787 jmc 1.1 C-- end bi,bj loops.
788     ENDDO
789     ENDDO
790    
791 jmc 1.99 #ifdef ALLOW_GMREDI
792     IF ( useGMRedi ) THEN
793     CALL GMREDI_DO_EXCH( myTime, myIter, myThid )
794     ENDIF
795     #endif /* ALLOW_GMREDI */
796    
797 jmc 1.98 #ifdef ALLOW_KPP
798 jmc 1.45 IF (useKPP) THEN
799     CALL KPP_DO_EXCH( myThid )
800     ENDIF
801 jmc 1.98 #endif /* ALLOW_KPP */
802 jmc 1.45
803 jmc 1.18 #ifdef ALLOW_DIAGNOSTICS
804     IF ( fluidIsWater .AND. useDiagnostics ) THEN
805 jmc 1.74 CALL DIAGS_RHO_G(
806 jmc 1.71 I rhoInSitu, uVel, vVel,
807     I myTime, myIter, myThid )
808 jmc 1.18 CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
809     ENDIF
810 jmc 1.19 IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN
811 jmc 1.71 CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
812     & 0, Nr, 0, 1, 1, myThid )
813 jmc 1.19 ENDIF
814 jmc 1.18 #endif
815    
816 jmc 1.1 #ifdef ALLOW_DEBUG
817 jmc 1.96 IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
818 jmc 1.1 #endif
819    
820     RETURN
821     END

  ViewVC Help
Powered by ViewVC 1.1.22