/[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.74 - (hide annotations) (download)
Fri Oct 3 02:26:49 2008 UTC (15 years, 8 months ago) by jmc
Branch: MAIN
Changes since 1.73: +18 -8 lines 
AD: back to version 1.72 ;
  + compute rhoInSitu if fluidIswater, otherwise set it to zero.
1 jmc 1.74 C $Header: /u/gcmpack/MITgcm/model/src/do_oceanic_phys.F,v 1.73 2008/09/25 14:52:25 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 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.74 IF ( DIAGNOSTICS_IS_ON('WRHOMASS',myThid) )
122 jmc 1.71 & 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.74 #else /* ALLOW_AUTODIFF_TAMC */
263 jmc 1.47 C if fluid is not water, by-pass find_rho, gmredi, surfaceForcing
264 jmc 1.36 C and all vertical mixing schemes, but keep OBCS_CALC
265     IF ( fluidIsWater ) THEN
266 jmc 1.74 #endif /* ALLOW_AUTODIFF_TAMC */
267 jmc 1.1
268     C-- Set up work arrays with valid (i.e. not NaN) values
269     C These inital values do not alter the numerical results. They
270     C just ensure that all memory references are to valid floating
271     C point numbers. This prevents spurious hardware signals due to
272     C uninitialised but inert locations.
273    
274 jmc 1.69 #ifdef ALLOW_AUTODIFF_TAMC
275 jmc 1.1 DO j=1-OLy,sNy+OLy
276     DO i=1-OLx,sNx+OLx
277 jmc 1.69 rhoKm1 (i,j) = 0. _d 0
278 jmc 1.47 rhoKp1 (i,j) = 0. _d 0
279 jmc 1.1 ENDDO
280     ENDDO
281 jmc 1.69 #endif /* ALLOW_AUTODIFF_TAMC */
282 jmc 1.1
283     DO k=1,Nr
284     DO j=1-OLy,sNy+OLy
285     DO i=1-OLx,sNx+OLx
286 jmc 1.69 C This is currently used by GMRedi, IVDC, MXL-depth and Diagnostics
287 jmc 1.1 sigmaX(i,j,k) = 0. _d 0
288     sigmaY(i,j,k) = 0. _d 0
289     sigmaR(i,j,k) = 0. _d 0
290     #ifdef ALLOW_AUTODIFF_TAMC
291     cph all the following init. are necessary for TAF
292     cph although some of these are re-initialised later.
293 jmc 1.71 c rhoInSitu(i,j,k,bi,bj) = 0.
294 jmc 1.1 IVDConvCount(i,j,k,bi,bj) = 0.
295     # ifdef ALLOW_GMREDI
296     Kwx(i,j,k,bi,bj) = 0. _d 0
297     Kwy(i,j,k,bi,bj) = 0. _d 0
298     Kwz(i,j,k,bi,bj) = 0. _d 0
299     # ifdef GM_NON_UNITY_DIAGONAL
300     Kux(i,j,k,bi,bj) = 0. _d 0
301     Kvy(i,j,k,bi,bj) = 0. _d 0
302     # endif
303     # ifdef GM_EXTRA_DIAGONAL
304     Kuz(i,j,k,bi,bj) = 0. _d 0
305     Kvz(i,j,k,bi,bj) = 0. _d 0
306     # endif
307     # ifdef GM_BOLUS_ADVEC
308     GM_PsiX(i,j,k,bi,bj) = 0. _d 0
309     GM_PsiY(i,j,k,bi,bj) = 0. _d 0
310     # endif
311     # ifdef GM_VISBECK_VARIABLE_K
312     VisbeckK(i,j,bi,bj) = 0. _d 0
313     # endif
314     # endif /* ALLOW_GMREDI */
315 heimbach 1.42 # ifdef ALLOW_KPP
316     KPPdiffKzS(i,j,k,bi,bj) = 0. _d 0
317     KPPdiffKzT(i,j,k,bi,bj) = 0. _d 0
318     # endif /* ALLOW_KPP */
319 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
320     ENDDO
321     ENDDO
322     ENDDO
323    
324     iMin = 1-OLx
325     iMax = sNx+OLx
326     jMin = 1-OLy
327     jMax = sNy+OLy
328    
329     #ifdef ALLOW_AUTODIFF_TAMC
330     CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
331     CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
332 heimbach 1.12 CADJ STORE totphihyd(:,:,:,bi,bj)
333 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
334 heimbach 1.10 # ifdef ALLOW_KPP
335 jmc 1.1 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
336     CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
337 heimbach 1.10 # endif
338 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
339    
340     #ifdef ALLOW_DEBUG
341 jmc 1.36 IF ( debugLevel .GE. debLevB )
342 jmc 1.1 & CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
343     #endif
344    
345     C-- Start of diagnostic loop
346     DO k=Nr,1,-1
347    
348     #ifdef ALLOW_AUTODIFF_TAMC
349     C? Patrick, is this formula correct now that we change the loop range?
350     C? Do we still need this?
351 jmc 1.69 cph kkey formula corrected.
352 jmc 1.47 cph Needed for rhoK, rhoKm1, in the case useGMREDI.
353 jmc 1.71 kkey = (itdkey-1)*Nr + k
354 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
355    
356 jmc 1.71 C-- Always compute density (stored in common block) here; even when it is not
357     C needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
358     #ifdef ALLOW_DEBUG
359     IF ( debugLevel .GE. debLevB )
360     & CALL DEBUG_CALL('FIND_RHO_2D',myThid)
361     #endif
362     #ifdef ALLOW_AUTODIFF_TAMC
363 jmc 1.74 IF ( fluidIsWater ) THEN
364 jmc 1.71 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 jmc 1.74 #ifdef ALLOW_AUTODIFF_TAMC
385     ELSE
386     C- fluid is not water:
387     DO j=1-OLy,sNy+OLy
388     DO i=1-OLx,sNx+OLx
389     rhoInSitu(i,j,k,bi,bj) = 0.
390     ENDDO
391     ENDDO
392     ENDIF
393     #endif /* ALLOW_AUTODIFF_TAMC */
394 jmc 1.69
395 jmc 1.1 C-- Calculate gradients of potential density for isoneutral
396     C slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
397 jmc 1.17 IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)
398 dimitri 1.61 & .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
399 jmc 1.1 IF (k.GT.1) THEN
400     #ifdef ALLOW_AUTODIFF_TAMC
401     CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
402     CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
403 heimbach 1.72 CADJ STORE rhokm1 (bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
404 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
405 jmc 1.68 CALL FIND_RHO_2D(
406     I iMin, iMax, jMin, jMax, k,
407     I theta(1-OLx,1-OLy,k-1,bi,bj),
408     I salt (1-OLx,1-OLy,k-1,bi,bj),
409     O rhoKm1,
410     I k-1, bi, bj, myThid )
411 jmc 1.1 ENDIF
412     #ifdef ALLOW_DEBUG
413 jmc 1.36 IF ( debugLevel .GE. debLevB )
414 jmc 1.1 & CALL DEBUG_CALL('GRAD_SIGMA',myThid)
415     #endif
416 heimbach 1.31 cph Avoid variable aliasing for adjoint !!!
417     DO j=jMin,jMax
418     DO i=iMin,iMax
419 jmc 1.71 rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
420 heimbach 1.31 ENDDO
421     ENDDO
422 jmc 1.1 CALL GRAD_SIGMA(
423     I bi, bj, iMin, iMax, jMin, jMax, k,
424 jmc 1.71 I rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
425 jmc 1.1 O sigmaX, sigmaY, sigmaR,
426     I myThid )
427 gforget 1.66 #ifdef ALLOW_AUTODIFF_TAMC
428 jmc 1.69 #ifdef GMREDI_WITH_STABLE_ADJOINT
429 gforget 1.66 cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
430     cgf -> cuts adjoint dependency from slope to state
431 jmc 1.69 CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
432     CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
433     CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
434 gforget 1.66 #endif
435     #endif /* ALLOW_AUTODIFF_TAMC */
436 jmc 1.1 ENDIF
437    
438     C-- Implicit Vertical Diffusion for Convection
439     c ==> should use sigmaR !!!
440     IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
441     #ifdef ALLOW_DEBUG
442 jmc 1.36 IF ( debugLevel .GE. debLevB )
443 jmc 1.1 & CALL DEBUG_CALL('CALC_IVDC',myThid)
444     #endif
445     CALL CALC_IVDC(
446     I bi, bj, iMin, iMax, jMin, jMax, k,
447 jmc 1.71 I rhoKm1, rhoInSitu(1-OLx,1-OLy,k,bi,bj),
448 jmc 1.1 I myTime, myIter, myThid)
449     ENDIF
450    
451 jmc 1.17 #ifdef ALLOW_DIAGNOSTICS
452 jmc 1.71 IF ( MOD(doDiagsRho,2).EQ.1 ) THEN
453     CALL DIAGS_RHO_L( k, bi, bj,
454     I rhoInSitu(1-OLx,1-OLy,k,bi,bj),
455 jmc 1.74 I rhoKm1, wVel,
456 jmc 1.71 I myTime, myIter, myThid )
457 jmc 1.17 ENDIF
458     #endif
459    
460 jmc 1.1 C-- end of diagnostic k loop (Nr:1)
461     ENDDO
462    
463 heimbach 1.57 #ifdef ALLOW_AUTODIFF_TAMC
464 jmc 1.69 CADJ STORE IVDConvCount(:,:,:,bi,bj)
465 heimbach 1.57 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
466     #endif
467    
468 jmc 1.47 C-- Diagnose Mixed Layer Depth:
469 jmc 1.71 IF ( useGMRedi .OR. doDiagsRho.GE.4 ) THEN
470     CALL CALC_OCE_MXLAYER(
471     I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
472     I bi, bj, myTime, myIter, myThid )
473 jmc 1.47 ENDIF
474 heimbach 1.53
475 dimitri 1.52 #ifdef ALLOW_SALT_PLUME
476 dimitri 1.61 IF ( useSALT_PLUME ) THEN
477 jmc 1.71 CALL SALT_PLUME_CALC_DEPTH(
478     I rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
479     I bi, bj, myTime, myIter, myThid )
480 dimitri 1.60 ENDIF
481 dimitri 1.61 #endif /* ALLOW_SALT_PLUME */
482    
483 jmc 1.8 #ifdef ALLOW_DIAGNOSTICS
484 jmc 1.71 IF ( MOD(doDiagsRho,4).GE.2 ) THEN
485 jmc 1.16 CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR ', 0, Nr,
486     & 2, bi, bj, myThid)
487 jmc 1.8 ENDIF
488 dimitri 1.61 #endif /* ALLOW_DIAGNOSTICS */
489 jmc 1.8
490 jmc 1.1 C-- Determines forcing terms based on external fields
491     C relaxation terms, etc.
492     #ifdef ALLOW_DEBUG
493 jmc 1.36 IF ( debugLevel .GE. debLevB )
494 jmc 1.1 & CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
495     #endif
496 heimbach 1.23 #ifdef ALLOW_AUTODIFF_TAMC
497     CADJ STORE EmPmR(:,:,bi,bj)
498     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
499 heimbach 1.26 # ifdef EXACT_CONSERV
500 heimbach 1.23 CADJ STORE PmEpR(:,:,bi,bj)
501     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
502 heimbach 1.26 # endif
503 heimbach 1.27 # ifdef NONLIN_FRSURF
504     CADJ STORE hFac_surfC(:,:,bi,bj)
505     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
506     CADJ STORE recip_hFacC(:,:,:,bi,bj)
507     CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
508     # endif
509 heimbach 1.23 #endif
510 jmc 1.36 CALL EXTERNAL_FORCING_SURF(
511 jmc 1.1 I bi, bj, iMin, iMax, jMin, jMax,
512     I myTime, myIter, myThid )
513 heimbach 1.27 #ifdef ALLOW_AUTODIFF_TAMC
514     # ifdef EXACT_CONSERV
515     cph-test
516     cphCADJ STORE PmEpR(:,:,bi,bj)
517     cphCADJ & = comlev1_bibj, key=itdkey, byte=isbyte
518     # endif
519     #endif
520 jmc 1.1
521     #ifdef ALLOW_AUTODIFF_TAMC
522     cph needed for KPP
523 jmc 1.4 CADJ STORE surfaceForcingU(:,:,bi,bj)
524 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
525 jmc 1.4 CADJ STORE surfaceForcingV(:,:,bi,bj)
526 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
527 jmc 1.4 CADJ STORE surfaceForcingS(:,:,bi,bj)
528 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
529 jmc 1.4 CADJ STORE surfaceForcingT(:,:,bi,bj)
530 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
531 jmc 1.4 CADJ STORE surfaceForcingTice(:,:,bi,bj)
532 jmc 1.1 CADJ & = comlev1_bibj, key=itdkey, byte=isbyte
533     #endif /* ALLOW_AUTODIFF_TAMC */
534    
535     #ifdef ALLOW_KPP
536     C-- Compute KPP mixing coefficients
537     IF (useKPP) THEN
538     #ifdef ALLOW_DEBUG
539 jmc 1.36 IF ( debugLevel .GE. debLevB )
540 jmc 1.1 & CALL DEBUG_CALL('KPP_CALC',myThid)
541     #endif
542     CALL KPP_CALC(
543 jmc 1.44 I bi, bj, myTime, myIter, myThid )
544 jmc 1.1 #ifdef ALLOW_AUTODIFF_TAMC
545     ELSE
546     CALL KPP_CALC_DUMMY(
547 jmc 1.44 I bi, bj, myTime, myIter, myThid )
548 jmc 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
549     ENDIF
550    
551     #endif /* ALLOW_KPP */
552    
553 mlosch 1.6 #ifdef ALLOW_PP81
554     C-- Compute PP81 mixing coefficients
555     IF (usePP81) THEN
556     #ifdef ALLOW_DEBUG
557 jmc 1.36 IF ( debugLevel .GE. debLevB )
558 mlosch 1.6 & CALL DEBUG_CALL('PP81_CALC',myThid)
559     #endif
560     CALL PP81_CALC(
561     I bi, bj, myTime, myThid )
562     ENDIF
563     #endif /* ALLOW_PP81 */
564    
565     #ifdef ALLOW_MY82
566     C-- Compute MY82 mixing coefficients
567     IF (useMY82) THEN
568     #ifdef ALLOW_DEBUG
569 jmc 1.36 IF ( debugLevel .GE. debLevB )
570 mlosch 1.6 & CALL DEBUG_CALL('MY82_CALC',myThid)
571     #endif
572     CALL MY82_CALC(
573     I bi, bj, myTime, myThid )
574     ENDIF
575     #endif /* ALLOW_MY82 */
576    
577 mlosch 1.9 #ifdef ALLOW_GGL90
578     C-- Compute GGL90 mixing coefficients
579     IF (useGGL90) THEN
580     #ifdef ALLOW_DEBUG
581 jmc 1.36 IF ( debugLevel .GE. debLevB )
582 mlosch 1.9 & CALL DEBUG_CALL('GGL90_CALC',myThid)
583     #endif
584     CALL GGL90_CALC(
585     I bi, bj, myTime, myThid )
586     ENDIF
587     #endif /* ALLOW_GGL90 */
588    
589 jmc 1.20 #ifdef ALLOW_TIMEAVE
590 jmc 1.36 IF ( taveFreq.GT. 0. _d 0 ) THEN
591 jmc 1.20 CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
592     ENDIF
593     IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
594     CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
595     I Nr, deltaTclock, bi, bj, myThid)
596     ENDIF
597     #endif /* ALLOW_TIMEAVE */
598    
599 jmc 1.69 #ifdef ALLOW_GMREDI
600 jmc 1.47 #ifdef ALLOW_AUTODIFF_TAMC
601     # ifndef GM_EXCLUDE_CLIPPING
602     cph storing here is needed only for one GMREDI_OPTIONS:
603     cph define GM_BOLUS_ADVEC
604     cph keep it although TAF says you dont need to.
605     cph but I've avoided the #ifdef for now, in case more things change
606     CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
607     CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
608     CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key=itdkey, byte=isbyte
609     # endif
610     #endif /* ALLOW_AUTODIFF_TAMC */
611    
612     C-- Calculate iso-neutral slopes for the GM/Redi parameterisation
613     IF (useGMRedi) THEN
614     #ifdef ALLOW_DEBUG
615     IF ( debugLevel .GE. debLevB )
616     & CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
617     #endif
618     CALL GMREDI_CALC_TENSOR(
619 jmc 1.51 I iMin, iMax, jMin, jMax,
620 jmc 1.47 I sigmaX, sigmaY, sigmaR,
621 jmc 1.51 I bi, bj, myTime, myIter, myThid )
622 jmc 1.47 #ifdef ALLOW_AUTODIFF_TAMC
623     ELSE
624     CALL GMREDI_CALC_TENSOR_DUMMY(
625 jmc 1.51 I iMin, iMax, jMin, jMax,
626 jmc 1.47 I sigmaX, sigmaY, sigmaR,
627 jmc 1.51 I bi, bj, myTime, myIter, myThid )
628 jmc 1.47 #endif /* ALLOW_AUTODIFF_TAMC */
629     ENDIF
630 jmc 1.69 #endif /* ALLOW_GMREDI */
631    
632     #ifdef ALLOW_DOWN_SLOPE
633     IF ( useDOWN_SLOPE ) THEN
634     C-- Calculate Downsloping Flow for Down_Slope parameterization
635     IF ( usingPCoords ) THEN
636     CALL DWNSLP_CALC_FLOW(
637 jmc 1.71 I bi, bj, kSurfC, rhoInSitu,
638 jmc 1.69 I myTime, myIter, myThid )
639     ELSE
640     CALL DWNSLP_CALC_FLOW(
641 jmc 1.71 I bi, bj, kLowC, rhoInSitu,
642 jmc 1.69 I myTime, myIter, myThid )
643     ENDIF
644     ENDIF
645     #endif /* ALLOW_DOWN_SLOPE */
646 jmc 1.47
647 jmc 1.74 #ifndef ALLOW_AUTODIFF_TAMC
648 jmc 1.36 C--- if fluid Is Water: end
649     ENDIF
650     #endif
651    
652     #ifdef ALLOW_OBCS
653     C-- Calculate future values on open boundaries
654     IF (useOBCS) THEN
655     #ifdef ALLOW_DEBUG
656     IF ( debugLevel .GE. debLevB )
657     & CALL DEBUG_CALL('OBCS_CALC',myThid)
658     #endif
659     CALL OBCS_CALC( bi, bj, myTime+deltaTclock, myIter+1,
660     I uVel, vVel, wVel, theta, salt,
661     I myThid )
662     ENDIF
663     #endif /* ALLOW_OBCS */
664    
665 jmc 1.1 C-- end bi,bj loops.
666     ENDDO
667     ENDDO
668    
669 jmc 1.45 #ifdef ALLOW_KPP
670     IF (useKPP) THEN
671     CALL KPP_DO_EXCH( myThid )
672     ENDIF
673     #endif /* ALLOW_KPP */
674    
675 jmc 1.18 #ifdef ALLOW_DIAGNOSTICS
676     IF ( fluidIsWater .AND. useDiagnostics ) THEN
677 jmc 1.74 CALL DIAGS_RHO_G(
678 jmc 1.71 I rhoInSitu, uVel, vVel,
679     I myTime, myIter, myThid )
680 jmc 1.18 CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
681     ENDIF
682 jmc 1.19 IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN
683 jmc 1.71 CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
684     & 0, Nr, 0, 1, 1, myThid )
685 jmc 1.19 ENDIF
686 jmc 1.18 #endif
687    
688 jmc 1.1 #ifdef ALLOW_DEBUG
689 jmc 1.29 IF ( debugLevel .GE. debLevB )
690     & CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
691 jmc 1.1 #endif
692    
693     RETURN
694     END
  ViewVC Help
Powered by ViewVC 1.1.22