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