/[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.93 - (hide annotations) (download)
Thu Sep 16 20:10:22 2010 UTC (13 years, 8 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint62k, checkpoint62l
Changes since 1.92: +7 -5 lines
Change store dir. logic following common block changes (05/2009)

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

  ViewVC Help
Powered by ViewVC 1.1.22