/[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.78 - (hide annotations) (download)
Sun May 24 18:00:39 2009 UTC (15 years ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint61o
Changes since 1.77: +5 -1 lines
One temp. store.

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

  ViewVC Help
Powered by ViewVC 1.1.22