/[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.104 - (hide annotations) (download)
Tue May 3 18:05:03 2011 UTC (13 years, 1 month ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint63, checkpoint62z, checkpoint62y, checkpoint62x
Changes since 1.103: +8 -6 lines
Remove unpleasant #ifndef ALLOW_AUTODIFF_TAMC
(hard to track down later)

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

  ViewVC Help
Powered by ViewVC 1.1.22