/[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.103 - (hide annotations) (download)
Wed Apr 20 16:20:37 2011 UTC (13 years, 1 month ago) by dimitri
Branch: MAIN
CVS Tags: checkpoint62w
Changes since 1.102: +5 -1 lines
temporarily exclude from adjoint computations until
impact has been vetted in forward integrations

1 dimitri 1.103 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.102 2011/04/20 13:45:14 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     # ifdef SEAICE_SALINITY
224 heimbach 1.77 CADJ STORE salt = comlev1, key = ikey_dynamics,
225     CADJ & kind = isbyte
226 heimbach 1.62 # endif
227     # 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 #ifndef ALLOW_AUTODIFF_TAMC
327     C temporarily exclude from adjoint computations until
328     C impact has been vetted in forward integrations
329 dimitri 1.102 IF ( allowInteriorFreezing ) THEN
330     CALL FREEZE_INTERIOR( myTime, myIter, myThid )
331     ENDIF
332 dimitri 1.103 #endif
333 dimitri 1.102
334 jmc 1.5 C-- Freeze water at the surface
335     #ifdef ALLOW_AUTODIFF_TAMC
336 heimbach 1.77 CADJ STORE theta = comlev1, key = ikey_dynamics,
337     CADJ & kind = isbyte
338 jmc 1.5 #endif
339 jmc 1.70 IF ( allowFreezing ) THEN
340 jmc 1.5 CALL FREEZE_SURFACE( myTime, myIter, myThid )
341     ENDIF
342    
343 jmc 1.28 #ifdef ALLOW_OCN_COMPON_INTERF
344 jmc 1.5 C-- Apply imported data (from coupled interface) to forcing fields
345 jmc 1.28 C jmc: do not know precisely where to put this call (bf or af thSIce ?)
346 jmc 1.36 IF ( useCoupler ) THEN
347 jmc 1.5 CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
348 jmc 1.36 ENDIF
349 jmc 1.28 #endif /* ALLOW_OCN_COMPON_INTERF */
350 jmc 1.5
351 jmc 1.25 #ifdef ALLOW_BALANCE_FLUXES
352 jmc 1.36 C balance fluxes
353     IF ( balanceEmPmR )
354 jmc 1.84 & CALL REMOVE_MEAN_RS( 1, EmPmR, maskInC, maskInC, rA, drF,
355 jmc 1.25 & 'EmPmR', myTime, myThid )
356 jmc 1.36 IF ( balanceQnet )
357 jmc 1.84 & CALL REMOVE_MEAN_RS( 1, Qnet, maskInC, maskInC, rA, drF,
358 jmc 1.25 & 'Qnet ', myTime, myThid )
359     #endif /* ALLOW_BALANCE_FLUXES */
360    
361 jmc 1.1 #ifdef ALLOW_AUTODIFF_TAMC
362     C-- HPF directive to help TAMC
363     CHPF$ INDEPENDENT
364     #endif /* ALLOW_AUTODIFF_TAMC */
365     DO bj=myByLo(myThid),myByHi(myThid)
366     #ifdef ALLOW_AUTODIFF_TAMC
367 heimbach 1.15 C-- HPF directive to help TAMC
368     CHPF$ INDEPENDENT
369 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
370     DO bi=myBxLo(myThid),myBxHi(myThid)
371    
372     #ifdef ALLOW_AUTODIFF_TAMC
373     act1 = bi - myBxLo(myThid)
374     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
375     act2 = bj - myByLo(myThid)
376     max2 = myByHi(myThid) - myByLo(myThid) + 1
377     act3 = myThid - 1
378     max3 = nTx*nTy
379     act4 = ikey_dynamics - 1
380     itdkey = (act1 + 1) + act2*max1
381     & + act3*max1*max2
382     & + act4*max1*max2*max3
383 jmc 1.74 #else /* ALLOW_AUTODIFF_TAMC */
384 jmc 1.47 C if fluid is not water, by-pass find_rho, gmredi, surfaceForcing
385 jmc 1.36 C and all vertical mixing schemes, but keep OBCS_CALC
386     IF ( fluidIsWater ) THEN
387 jmc 1.74 #endif /* ALLOW_AUTODIFF_TAMC */
388 jmc 1.1
389     C-- Set up work arrays with valid (i.e. not NaN) values
390     C These inital values do not alter the numerical results. They
391     C just ensure that all memory references are to valid floating
392     C point numbers. This prevents spurious hardware signals due to
393     C uninitialised but inert locations.
394    
395 jmc 1.69 #ifdef ALLOW_AUTODIFF_TAMC
396 jmc 1.1 DO j=1-OLy,sNy+OLy
397     DO i=1-OLx,sNx+OLx
398 jmc 1.69 rhoKm1 (i,j) = 0. _d 0
399 jmc 1.47 rhoKp1 (i,j) = 0. _d 0
400 jmc 1.1 ENDDO
401     ENDDO
402 jmc 1.69 #endif /* ALLOW_AUTODIFF_TAMC */
403 jmc 1.1
404     DO k=1,Nr
405     DO j=1-OLy,sNy+OLy
406     DO i=1-OLx,sNx+OLx
407 jmc 1.69 C This is currently used by GMRedi, IVDC, MXL-depth and Diagnostics
408 jmc 1.1 sigmaX(i,j,k) = 0. _d 0
409     sigmaY(i,j,k) = 0. _d 0
410     sigmaR(i,j,k) = 0. _d 0
411     #ifdef ALLOW_AUTODIFF_TAMC
412     cph all the following init. are necessary for TAF
413     cph although some of these are re-initialised later.
414 jmc 1.71 c rhoInSitu(i,j,k,bi,bj) = 0.
415 jmc 1.1 IVDConvCount(i,j,k,bi,bj) = 0.
416     # ifdef ALLOW_GMREDI
417     Kwx(i,j,k,bi,bj) = 0. _d 0
418     Kwy(i,j,k,bi,bj) = 0. _d 0
419     Kwz(i,j,k,bi,bj) = 0. _d 0
420     # ifdef GM_NON_UNITY_DIAGONAL
421     Kux(i,j,k,bi,bj) = 0. _d 0
422     Kvy(i,j,k,bi,bj) = 0. _d 0
423     # endif
424     # ifdef GM_EXTRA_DIAGONAL
425     Kuz(i,j,k,bi,bj) = 0. _d 0
426     Kvz(i,j,k,bi,bj) = 0. _d 0
427     # endif
428     # ifdef GM_BOLUS_ADVEC
429     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
430     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
431     # endif
432     # ifdef GM_VISBECK_VARIABLE_K
433     VisbeckK(i,j,bi,bj) = 0. _d 0
434     # endif
435     # endif /* ALLOW_GMREDI */
436 heimbach 1.42 # ifdef ALLOW_KPP
437     KPPdiffKzS(i,j,k,bi,bj) = 0. _d 0
438     KPPdiffKzT(i,j,k,bi,bj) = 0. _d 0
439     # endif /* ALLOW_KPP */
440 gforget 1.91 # ifdef ALLOW_GGL90
441     GGL90viscArU(i,j,k,bi,bj) = 0. _d 0
442     GGL90viscArV(i,j,k,bi,bj) = 0. _d 0
443     GGL90diffKr(i,j,k,bi,bj) = 0. _d 0
444     # endif /* ALLOW_GGL90 */
445 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
446     ENDDO
447     ENDDO
448     ENDDO
449    
450     iMin = 1-OLx
451     iMax = sNx+OLx
452     jMin = 1-OLy
453     jMax = sNy+OLy
454    
455     #ifdef ALLOW_AUTODIFF_TAMC
456 jmc 1.96 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
457 heimbach 1.77 CADJ & kind = isbyte
458 jmc 1.96 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
459 heimbach 1.77 CADJ & kind = isbyte
460 heimbach 1.12 CADJ STORE totphihyd(:,:,:,bi,bj)
461 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
462 heimbach 1.77 CADJ & kind = isbyte
463 heimbach 1.10 # ifdef ALLOW_KPP
464 jmc 1.96 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
465 heimbach 1.77 CADJ & kind = isbyte
466 jmc 1.96 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
467 heimbach 1.77 CADJ & kind = isbyte
468 heimbach 1.10 # endif
469 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
470    
471     #ifdef ALLOW_DEBUG
472 jmc 1.96 IF (debugMode) CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
473 jmc 1.1 #endif
474    
475     C-- Start of diagnostic loop
476     DO k=Nr,1,-1
477    
478     #ifdef ALLOW_AUTODIFF_TAMC
479     C? Patrick, is this formula correct now that we change the loop range?
480     C? Do we still need this?
481 jmc 1.69 cph kkey formula corrected.
482 jmc 1.47 cph Needed for rhoK, rhoKm1, in the case useGMREDI.
483 jmc 1.71 kkey = (itdkey-1)*Nr + k
484 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
485    
486 jmc 1.71 C-- Always compute density (stored in common block) here; even when it is not
487     C needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
488     #ifdef ALLOW_DEBUG
489 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D',myThid)
490 jmc 1.71 #endif
491     #ifdef ALLOW_AUTODIFF_TAMC
492 jmc 1.74 IF ( fluidIsWater ) THEN
493 jmc 1.96 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
494 heimbach 1.77 CADJ & kind = isbyte
495 jmc 1.96 CADJ STORE salt(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
496 heimbach 1.77 CADJ & kind = isbyte
497 jmc 1.71 #endif /* ALLOW_AUTODIFF_TAMC */
498 jmc 1.69 #ifdef ALLOW_DOWN_SLOPE
499     IF ( useDOWN_SLOPE ) THEN
500     CALL DWNSLP_CALC_RHO(
501     I theta, salt,
502 jmc 1.71 O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
503 jmc 1.69 I k, bi, bj, myTime, myIter, myThid )
504 jmc 1.71 ELSE
505     #endif /* ALLOW_DOWN_SLOPE */
506     CALL FIND_RHO_2D(
507     I iMin, iMax, jMin, jMax, k,
508     I theta(1-OLx,1-OLy,k,bi,bj),
509     I salt (1-OLx,1-OLy,k,bi,bj),
510     O rhoInSitu(1-OLx,1-OLy,k,bi,bj),
511     I k, bi, bj, myThid )
512     #ifdef ALLOW_DOWN_SLOPE
513 jmc 1.69 ENDIF
514     #endif /* ALLOW_DOWN_SLOPE */
515 jmc 1.74 #ifdef ALLOW_AUTODIFF_TAMC
516     ELSE
517     C- fluid is not water:
518     DO j=1-OLy,sNy+OLy
519     DO i=1-OLx,sNx+OLx
520     rhoInSitu(i,j,k,bi,bj) = 0.
521     ENDDO
522     ENDDO
523     ENDIF
524     #endif /* ALLOW_AUTODIFF_TAMC */
525 jmc 1.69
526 jmc 1.1 C-- Calculate gradients of potential density for isoneutral
527     C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
528 jmc 1.17 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)
529 dimitri 1.61 & .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
530 jmc 1.1 IF (k.GT.1) THEN
531     #ifdef ALLOW_AUTODIFF_TAMC
532 jmc 1.96 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
533 heimbach 1.77 CADJ & kind = isbyte
534 jmc 1.96 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
535 heimbach 1.77 CADJ & kind = isbyte
536 jmc 1.96 CADJ STORE rhokm1 (bi,bj) = comlev1_bibj_k, key=kkey,
537 heimbach 1.77 CADJ & kind = isbyte
538 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
539 jmc 1.68 CALL FIND_RHO_2D(
540     I iMin, iMax, jMin, jMax, k,
541     I theta(1-OLx,1-OLy,k-1,bi,bj),
542     I salt (1-OLx,1-OLy,k-1,bi,bj),
543     O rhoKm1,
544     I k-1, bi, bj, myThid )
545 jmc 1.1 ENDIF
546     #ifdef ALLOW_DEBUG
547 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
548 jmc 1.1 #endif
549 heimbach 1.31 cph Avoid variable aliasing for adjoint !!!
550     DO j=jMin,jMax
551     DO i=iMin,iMax
552 jmc 1.71 rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
553 heimbach 1.31 ENDDO
554     ENDDO
555 jmc 1.1 CALL GRAD_SIGMA(
556     I bi, bj, iMin, iMax, jMin, jMax, k,
557 jmc 1.71 I rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
558 jmc 1.1 O sigmaX, sigmaY, sigmaR,
559     I myThid )
560 gforget 1.66 #ifdef ALLOW_AUTODIFF_TAMC
561 jmc 1.69 #ifdef GMREDI_WITH_STABLE_ADJOINT
562 gforget 1.66 cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
563     cgf -> cuts adjoint dependency from slope to state
564 jmc 1.69 CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
565     CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
566     CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
567 gforget 1.66 #endif
568     #endif /* ALLOW_AUTODIFF_TAMC */
569 jmc 1.1 ENDIF
570    
571     C-- Implicit Vertical Diffusion for Convection
572     c ==> should use sigmaR !!!
573     IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
574     #ifdef ALLOW_DEBUG
575 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid)
576 jmc 1.1 #endif
577     CALL CALC_IVDC(
578     I bi, bj, iMin, iMax, jMin, jMax, k,
579 jmc 1.71 I rhoKm1, rhoInSitu(1-OLx,1-OLy,k,bi,bj),
580 jmc 1.1 I myTime, myIter, myThid)
581     ENDIF
582    
583 jmc 1.17 #ifdef ALLOW_DIAGNOSTICS
584 jmc 1.71 IF ( MOD(doDiagsRho,2).EQ.1 ) THEN
585     CALL DIAGS_RHO_L( k, bi, bj,
586     I rhoInSitu(1-OLx,1-OLy,k,bi,bj),
587 jmc 1.74 I rhoKm1, wVel,
588 jmc 1.71 I myTime, myIter, myThid )
589 jmc 1.17 ENDIF
590     #endif
591    
592 jmc 1.1 C-- end of diagnostic k loop (Nr:1)
593     ENDDO
594    
595 heimbach 1.57 #ifdef ALLOW_AUTODIFF_TAMC
596 jmc 1.69 CADJ STORE IVDConvCount(:,:,:,bi,bj)
597 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
598 heimbach 1.77 CADJ & kind = isbyte
599 heimbach 1.57 #endif
600    
601 jmc 1.47 C-- Diagnose Mixed Layer Depth:
602 jmc 1.71 IF ( useGMRedi .OR. doDiagsRho.GE.4 ) THEN
603     CALL CALC_OCE_MXLAYER(
604     I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
605     I bi, bj, myTime, myIter, myThid )
606 jmc 1.47 ENDIF
607 heimbach 1.53
608 dimitri 1.52 #ifdef ALLOW_SALT_PLUME
609 dimitri 1.61 IF ( useSALT_PLUME ) THEN
610 jmc 1.71 CALL SALT_PLUME_CALC_DEPTH(
611     I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
612     I bi, bj, myTime, myIter, myThid )
613 dimitri 1.60 ENDIF
614 dimitri 1.61 #endif /* ALLOW_SALT_PLUME */
615    
616 jmc 1.8 #ifdef ALLOW_DIAGNOSTICS
617 jmc 1.71 IF ( MOD(doDiagsRho,4).GE.2 ) THEN
618 jmc 1.16 CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR ', 0, Nr,
619     & 2, bi, bj, myThid)
620 jmc 1.8 ENDIF
621 dimitri 1.61 #endif /* ALLOW_DIAGNOSTICS */
622 jmc 1.8
623 jmc 1.1 C-- Determines forcing terms based on external fields
624     C relaxation terms, etc.
625     #ifdef ALLOW_DEBUG
626 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
627 jmc 1.1 #endif
628 heimbach 1.23 #ifdef ALLOW_AUTODIFF_TAMC
629     CADJ STORE EmPmR(:,:,bi,bj)
630 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
631 heimbach 1.77 CADJ & kind = isbyte
632 heimbach 1.26 # ifdef EXACT_CONSERV
633 heimbach 1.23 CADJ STORE PmEpR(:,:,bi,bj)
634 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
635 heimbach 1.77 CADJ & kind = isbyte
636 heimbach 1.26 # endif
637 heimbach 1.27 # ifdef NONLIN_FRSURF
638     CADJ STORE hFac_surfC(:,:,bi,bj)
639 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
640 heimbach 1.77 CADJ & kind = isbyte
641 heimbach 1.27 CADJ STORE recip_hFacC(:,:,:,bi,bj)
642 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
643 heimbach 1.77 CADJ & kind = isbyte
644 heimbach 1.97 # if (defined (ALLOW_PTRACERS))
645     CADJ STORE surfaceForcingS(:,:,bi,bj) = comlev1_bibj, key=itdkey,
646     CADJ & kind = isbyte
647     CADJ STORE surfaceForcingT(:,:,bi,bj) = comlev1_bibj, key=itdkey,
648     CADJ & kind = isbyte
649     # endif /* ALLOW_PTRACERS */
650     # endif /* NONLIN_FRSURF */
651 heimbach 1.23 #endif
652 jmc 1.36 CALL EXTERNAL_FORCING_SURF(
653 jmc 1.1 I bi, bj, iMin, iMax, jMin, jMax,
654     I myTime, myIter, myThid )
655 heimbach 1.27 #ifdef ALLOW_AUTODIFF_TAMC
656     # ifdef EXACT_CONSERV
657     cph-test
658     cphCADJ STORE PmEpR(:,:,bi,bj)
659 jmc 1.96 cphCADJ & = comlev1_bibj, key=itdkey,
660 heimbach 1.77 cphCADJ & kind = isbyte
661 heimbach 1.27 # endif
662     #endif
663 jmc 1.1
664     #ifdef ALLOW_AUTODIFF_TAMC
665     cph needed for KPP
666 jmc 1.4 CADJ STORE surfaceForcingU(:,:,bi,bj)
667 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
668 heimbach 1.77 CADJ & kind = isbyte
669 jmc 1.4 CADJ STORE surfaceForcingV(:,:,bi,bj)
670 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
671 heimbach 1.77 CADJ & kind = isbyte
672 jmc 1.4 CADJ STORE surfaceForcingS(:,:,bi,bj)
673 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
674 heimbach 1.77 CADJ & kind = isbyte
675 jmc 1.4 CADJ STORE surfaceForcingT(:,:,bi,bj)
676 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
677 heimbach 1.77 CADJ & kind = isbyte
678 jmc 1.4 CADJ STORE surfaceForcingTice(:,:,bi,bj)
679 jmc 1.96 CADJ & = comlev1_bibj, key=itdkey,
680 heimbach 1.77 CADJ & kind = isbyte
681 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
682    
683     #ifdef ALLOW_KPP
684     C-- Compute KPP mixing coefficients
685     IF (useKPP) THEN
686     #ifdef ALLOW_DEBUG
687 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid)
688 jmc 1.1 #endif
689 dfer 1.76 CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
690 jmc 1.1 CALL KPP_CALC(
691 jmc 1.44 I bi, bj, myTime, myIter, myThid )
692 dfer 1.76 CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
693 jmc 1.1 #ifdef ALLOW_AUTODIFF_TAMC
694     ELSE
695     CALL KPP_CALC_DUMMY(
696 jmc 1.44 I bi, bj, myTime, myIter, myThid )
697 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
698     ENDIF
699    
700     #endif /* ALLOW_KPP */
701    
702 mlosch 1.6 #ifdef ALLOW_PP81
703     C-- Compute PP81 mixing coefficients
704     IF (usePP81) THEN
705     #ifdef ALLOW_DEBUG
706 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid)
707 mlosch 1.6 #endif
708     CALL PP81_CALC(
709     I bi, bj, myTime, myThid )
710     ENDIF
711     #endif /* ALLOW_PP81 */
712    
713     #ifdef ALLOW_MY82
714     C-- Compute MY82 mixing coefficients
715     IF (useMY82) THEN
716     #ifdef ALLOW_DEBUG
717 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid)
718 mlosch 1.6 #endif
719     CALL MY82_CALC(
720     I bi, bj, myTime, myThid )
721     ENDIF
722     #endif /* ALLOW_MY82 */
723    
724 mlosch 1.9 #ifdef ALLOW_GGL90
725 gforget 1.91 #ifdef ALLOW_AUTODIFF_TAMC
726     CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
727     CADJ & kind = isbyte
728     #endif /* ALLOW_AUTODIFF_TAMC */
729 mlosch 1.9 C-- Compute GGL90 mixing coefficients
730     IF (useGGL90) THEN
731     #ifdef ALLOW_DEBUG
732 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
733 mlosch 1.9 #endif
734 dfer 1.76 CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
735 mlosch 1.9 CALL GGL90_CALC(
736 jmc 1.90 I bi, bj, myTime, myIter, myThid )
737 dfer 1.76 CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
738 mlosch 1.9 ENDIF
739     #endif /* ALLOW_GGL90 */
740    
741 jmc 1.20 #ifdef ALLOW_TIMEAVE
742 jmc 1.36 IF ( taveFreq.GT. 0. _d 0 ) THEN
743 jmc 1.20 CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
744     ENDIF
745     IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
746     CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
747     I Nr, deltaTclock, bi, bj, myThid)
748     ENDIF
749     #endif /* ALLOW_TIMEAVE */
750    
751 jmc 1.69 #ifdef ALLOW_GMREDI
752 jmc 1.47 #ifdef ALLOW_AUTODIFF_TAMC
753     # ifndef GM_EXCLUDE_CLIPPING
754     cph storing here is needed only for one GMREDI_OPTIONS:
755     cph define GM_BOLUS_ADVEC
756     cph keep it although TAF says you dont need to.
757 jmc 1.86 cph but I have avoided the #ifdef for now, in case more things change
758 jmc 1.96 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey,
759 heimbach 1.77 CADJ & kind = isbyte
760 jmc 1.96 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey,
761 heimbach 1.77 CADJ & kind = isbyte
762 jmc 1.96 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey,
763 heimbach 1.77 CADJ & kind = isbyte
764 jmc 1.47 # endif
765     #endif /* ALLOW_AUTODIFF_TAMC */
766    
767     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
768     IF (useGMRedi) THEN
769     #ifdef ALLOW_DEBUG
770 jmc 1.96 IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
771 jmc 1.47 #endif
772     CALL GMREDI_CALC_TENSOR(
773 jmc 1.51 I iMin, iMax, jMin, jMax,
774 jmc 1.47 I sigmaX, sigmaY, sigmaR,
775 jmc 1.51 I bi, bj, myTime, myIter, myThid )
776 jmc 1.47 #ifdef ALLOW_AUTODIFF_TAMC
777     ELSE
778     CALL GMREDI_CALC_TENSOR_DUMMY(
779 jmc 1.51 I iMin, iMax, jMin, jMax,
780 jmc 1.47 I sigmaX, sigmaY, sigmaR,
781 jmc 1.51 I bi, bj, myTime, myIter, myThid )
782 jmc 1.47 #endif /* ALLOW_AUTODIFF_TAMC */
783     ENDIF
784 jmc 1.69 #endif /* ALLOW_GMREDI */
785    
786     #ifdef ALLOW_DOWN_SLOPE
787     IF ( useDOWN_SLOPE ) THEN
788     C-- Calculate Downsloping Flow for Down_Slope parameterization
789     IF ( usingPCoords ) THEN
790     CALL DWNSLP_CALC_FLOW(
791 jmc 1.71 I bi, bj, kSurfC, rhoInSitu,
792 jmc 1.69 I myTime, myIter, myThid )
793     ELSE
794     CALL DWNSLP_CALC_FLOW(
795 jmc 1.71 I bi, bj, kLowC, rhoInSitu,
796 jmc 1.69 I myTime, myIter, myThid )
797     ENDIF
798     ENDIF
799     #endif /* ALLOW_DOWN_SLOPE */
800 jmc 1.47
801 jmc 1.98 #ifndef ALLOW_AUTODIFF_TAMC
802     C--- if fluid Is Water: end
803     ENDIF
804     #endif
805    
806 dimitri 1.94 #ifdef ALLOW_MYPACKAGE
807     IF ( useMYPACKAGE ) THEN
808     CALL MYPACKAGE_CALC_RHS(
809     I bi, bj, myTime, myIter, myThid )
810     ENDIF
811     #endif /* ALLOW_MYPACKAGE */
812    
813 jmc 1.1 C-- end bi,bj loops.
814     ENDDO
815     ENDDO
816    
817 jmc 1.99 #ifdef ALLOW_GMREDI
818     IF ( useGMRedi ) THEN
819     CALL GMREDI_DO_EXCH( myTime, myIter, myThid )
820     ENDIF
821     #endif /* ALLOW_GMREDI */
822    
823 jmc 1.98 #ifdef ALLOW_KPP
824 jmc 1.45 IF (useKPP) THEN
825     CALL KPP_DO_EXCH( myThid )
826     ENDIF
827 jmc 1.98 #endif /* ALLOW_KPP */
828 jmc 1.45
829 jmc 1.18 #ifdef ALLOW_DIAGNOSTICS
830     IF ( fluidIsWater .AND. useDiagnostics ) THEN
831 jmc 1.74 CALL DIAGS_RHO_G(
832 jmc 1.71 I rhoInSitu, uVel, vVel,
833     I myTime, myIter, myThid )
834 jmc 1.18 CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
835     ENDIF
836 jmc 1.19 IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN
837 jmc 1.71 CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
838     & 0, Nr, 0, 1, 1, myThid )
839 jmc 1.19 ENDIF
840 jmc 1.18 #endif
841    
842 jmc 1.1 #ifdef ALLOW_DEBUG
843 jmc 1.96 IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
844 jmc 1.1 #endif
845    
846     RETURN
847     END

  ViewVC Help
Powered by ViewVC 1.1.22