/[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.73 - (hide annotations) (download)
Thu Sep 25 14:52:25 2008 UTC (15 years, 7 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint61d
Changes since 1.72: +5 -4 lines
try to put back "if fluidIsWater" (needed for hs94 AD test)

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

  ViewVC Help
Powered by ViewVC 1.1.22