/[MITgcm]/MITgcm/model/src/do_oceanic_phys.F
ViewVC logotype

Diff of /MITgcm/model/src/do_oceanic_phys.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.42 by heimbach, Tue Apr 24 18:45:11 2007 UTC revision 1.148 by heimbach, Fri Aug 12 14:15:58 2016 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "PACKAGES_CONFIG.h"  #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    #ifdef ALLOW_AUTODIFF
7    # include "AUTODIFF_OPTIONS.h"
8    #endif
9    #ifdef ALLOW_CTRL
10    # include "CTRL_OPTIONS.h"
11    #endif
12    #ifdef ALLOW_SALT_PLUME
13    # include "SALT_PLUME_OPTIONS.h"
14    #endif
15    #ifdef ALLOW_ECCO
16    # include "ECCO_OPTIONS.h"
17    #endif
18    
19  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
20    # ifdef ALLOW_GGL90
21    #  include "GGL90_OPTIONS.h"
22    # endif
23  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
24  #  include "GMREDI_OPTIONS.h"  #  include "GMREDI_OPTIONS.h"
25  # endif  # endif
# Line 14  C $Name$ Line 29  C $Name$
29  # ifdef ALLOW_SEAICE  # ifdef ALLOW_SEAICE
30  #  include "SEAICE_OPTIONS.h"  #  include "SEAICE_OPTIONS.h"
31  # endif  # endif
32  #endif /* ALLOW_AUTODIFF_TAMC */  # ifdef ALLOW_EXF
33    #  include "EXF_OPTIONS.h"
34    # endif
35    #endif /* ALLOW_AUTODIFF */
36    
37  CBOP  CBOP
38  C     !ROUTINE: DO_OCEANIC_PHYS  C     !ROUTINE: DO_OCEANIC_PHYS
# Line 30  C     | o originally, part of S/R thermo Line 48  C     | o originally, part of S/R thermo
48  C     *==========================================================*  C     *==========================================================*
49  C     \ev  C     \ev
50    
51    C     !CALLING SEQUENCE:
52    C     DO_OCEANIC_PHYS
53    C       |
54    C       |-- OBCS_CALC
55    C       |
56    C       |-- OCN_APPLY_IMPORT
57    C       |
58    C       |-- FRAZIL_CALC_RHS
59    C       |
60    C       |-- THSICE_MAIN
61    C       |
62    C       |-- SEAICE_FAKE
63    C       |-- SEAICE_MODEL
64    C       |-- SEAICE_COST_SENSI
65    C       |
66    C       |-- OCN_EXPORT_DATA
67    C       |
68    C       |-- SHELFICE_THERMODYNAMICS
69    C       |
70    C       |-- ICEFRONT_THERMODYNAMICS
71    C       |
72    C       |-- SALT_PLUME_DO_EXCH
73    C       |
74    C       |-- FREEZE_SURFACE
75    C       |
76    C       |-- EXTERNAL_FORCING_SURF
77    C       |
78    C       |- k loop (Nr:1):
79    C       | - DWNSLP_CALC_RHO
80    C       | - BBL_CALC_RHO
81    C       | - FIND_RHO_2D @ p(k)
82    C       | - FIND_RHO_2D @ p(k-1)
83    C       | - GRAD_SIGMA
84    C       | - CALC_IVDC
85    C       | - DIAGS_RHO_L
86    C       |- end k loop.
87    C       |
88    C       |-- CALC_OCE_MXLAYER
89    C       |
90    C       |-- SALT_PLUME_CALC_DEPTH
91    C       |-- SALT_PLUME_VOLFRAC
92    C       |-- SALT_PLUME_APPLY
93    C       |-- SALT_PLUME_APPLY
94    C       |-- SALT_PLUME_FORCING_SURF
95    C       |
96    C       |-- KPP_CALC
97    C       |-- KPP_CALC_DUMMY
98    C       |
99    C       |-- PP81_CALC
100    C       |
101    C       |-- KL10_CALC
102    C       |
103    C       |-- MY82_CALC
104    C       |
105    C       |-- GGL90_CALC
106    C       |
107    C       |-- TIMEAVE_SURF_FLUX
108    C       |
109    C       |-- GMREDI_CALC_TENSOR
110    C       |-- GMREDI_CALC_TENSOR_DUMMY
111    C       |
112    C       |-- DWNSLP_CALC_FLOW
113    C       |-- DWNSLP_CALC_FLOW
114    C       |
115    C       |-- OFFLINE_GET_DIFFUS
116    C       |
117    C       |-- BBL_CALC_RHS
118    C       |
119    C       |-- MYPACKAGE_CALC_RHS
120    C       |
121    C       |-- GMREDI_DO_EXCH
122    C       |
123    C       |-- KPP_DO_EXCH
124    C       |
125    C       |-- DIAGS_RHO_G
126    C       |-- DIAGS_OCEANIC_SURF_FLUX
127    C       |-- SALT_PLUME_DIAGNOSTICS_FILL
128    C       |
129    C       |-- ECCO_PHYS
130    
131  C     !USES:  C     !USES:
132        IMPLICIT NONE        IMPLICIT NONE
133  C     == Global variables ===  C     == Global variables ===
134  #include "SIZE.h"  #include "SIZE.h"
135  #include "EEPARAMS.h"  #include "EEPARAMS.h"
136  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
137  #include "GRID.h"  #include "GRID.h"
138    #include "DYNVARS.h"
139  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
140  #include "TIMEAVE_STATV.h"  # include "TIMEAVE_STATV.h"
141  #endif  #endif
142  #if defined (ALLOW_BALANCE_FLUXES) && !(defined ALLOW_AUTODIFF_TAMC)  #ifdef ALLOW_OFFLINE
143  #include "FFIELDS.h"  # include "OFFLINE_SWITCH.h"
144  #endif  #endif
145    
146  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
147    # include "AUTODIFF_MYFIELDS.h"
148  # include "tamc.h"  # include "tamc.h"
149  # include "tamc_keys.h"  # include "tamc_keys.h"
150  # include "FFIELDS.h"  # include "FFIELDS.h"
151    # include "SURFACE.h"
152  # include "EOS.h"  # include "EOS.h"
153    # ifdef ALLOW_GMREDI
154    #  include "GMREDI.h"
155    # endif
156  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
157  #  include "KPP.h"  #  include "KPP.h"
158  # endif  # endif
159  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GGL90
160  #  include "GMREDI.h"  #  include "GGL90.h"
161  # endif  # endif
162  # ifdef ALLOW_EBM  # ifdef ALLOW_EBM
163  #  include "EBM.h"  #  include "EBM.h"
164  # endif  # endif
 # ifdef EXACT_CONSERV  
 #  include "SURFACE.h"  
 # endif  
165  # ifdef ALLOW_EXF  # ifdef ALLOW_EXF
166  #  include "ctrl.h"  #  include "ctrl.h"
167  #  include "EXF_FIELDS.h"  #  include "EXF_FIELDS.h"
# Line 70  C     == Global variables === Line 170  C     == Global variables ===
170  #  endif  #  endif
171  # endif  # endif
172  # ifdef ALLOW_SEAICE  # ifdef ALLOW_SEAICE
173    #  include "SEAICE_SIZE.h"
174  #  include "SEAICE.h"  #  include "SEAICE.h"
175    #  include "SEAICE_PARAMS.h"
176  # endif  # endif
177  #endif /* ALLOW_AUTODIFF_TAMC */  # ifdef ALLOW_THSICE
178    #  include "THSICE_VARS.h"
179    # endif
180    # ifdef ALLOW_SALT_PLUME
181    #  include "SALT_PLUME.h"
182    # endif
183    # ifdef ALLOW_ECCO
184    #  ifdef ALLOW_SIGMAR_COST_CONTRIBUTION
185    #   include "ecco_cost.h"
186    #  endif
187    # endif
188    #endif /* ALLOW_AUTODIFF */
189    
190  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
191  C     == Routine arguments ==  C     == Routine arguments ==
# Line 85  C     myThid :: Thread number for this i Line 198  C     myThid :: Thread number for this i
198    
199  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
200  C     == Local variables  C     == Local variables
201  C     rhoK, rhoKM1  :: Density at current level, and level above  C     rhoK, rhoKm1  :: Density at current level, and level above
202  C     iMin, iMax    :: Ranges and sub-block indices on which calculations  C     iMin, iMax    :: Ranges and sub-block indices on which calculations
203  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
204  C     bi, bj        :: tile indices  C     bi, bj        :: tile indices
205    C     msgBuf        :: Temp. for building output string
206  C     i,j,k         :: loop indices  C     i,j,k         :: loop indices
207        _RL rhokp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhoKp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
208        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhoKm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
209        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
210        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
211        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 100  C     i,j,k         :: loop indices Line 213  C     i,j,k         :: loop indices
213        INTEGER jMin, jMax        INTEGER jMin, jMax
214        INTEGER bi, bj        INTEGER bi, bj
215        INTEGER i, j, k        INTEGER i, j, k
216          CHARACTER*(MAX_LEN_MBUF) msgBuf
217        INTEGER doDiagsRho        INTEGER doDiagsRho
218          LOGICAL calcGMRedi, calcKPP, calcConvect
219  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
220        LOGICAL  DIAGNOSTICS_IS_ON        LOGICAL  DIAGNOSTICS_IS_ON
221        EXTERNAL DIAGNOSTICS_IS_ON        EXTERNAL DIAGNOSTICS_IS_ON
222  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
223    #ifdef ALLOW_AUTODIFF
224          _RL thetaRef
225    #endif /* ALLOW_AUTODIFF */
226  CEOP  CEOP
227    
228  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 114  C--   dummy statement to end declaration Line 231  C--   dummy statement to end declaration
231  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
232    
233  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
234        IF ( debugLevel .GE. debLevB )        IF (debugMode) CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
      &     CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)  
235  #endif  #endif
236    
237        doDiagsRho = 0        doDiagsRho = 0
238  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
239        IF ( useDiagnostics .AND. fluidIsWater ) THEN        IF ( useDiagnostics .AND. fluidIsWater ) THEN
240          IF ( DIAGNOSTICS_IS_ON('DRHODR  ',myThid) ) doDiagsRho = 1          IF ( DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) )
241          IF ( DIAGNOSTICS_IS_ON('RHOANOSQ',myThid) .OR.       &       doDiagsRho = doDiagsRho + 1
242       &       DIAGNOSTICS_IS_ON('URHOMASS',myThid) .OR.          IF ( DIAGNOSTICS_IS_ON('DRHODR  ',myThid) )
243       &       DIAGNOSTICS_IS_ON('VRHOMASS',myThid) .OR.       &       doDiagsRho = doDiagsRho + 2
244       &       DIAGNOSTICS_IS_ON('WRHOMASS',myThid) .OR.          IF ( DIAGNOSTICS_IS_ON('WdRHO_P ',myThid) )
245       &       DIAGNOSTICS_IS_ON('WRHOMASS',myThid) ) doDiagsRho = 2       &       doDiagsRho = doDiagsRho + 4
246            IF ( DIAGNOSTICS_IS_ON('WdRHOdP ',myThid) )
247         &       doDiagsRho = doDiagsRho + 8
248        ENDIF        ENDIF
249  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
250    
251          calcGMRedi  = useGMRedi
252          calcKPP     = useKPP
253          calcConvect = ivdc_kappa.NE.0.
254    #ifdef ALLOW_OFFLINE
255          IF ( useOffLine ) THEN
256            calcGMRedi = useGMRedi .AND. .NOT.offlineLoadGMRedi
257            calcKPP    = useKPP    .AND. .NOT.offlineLoadKPP
258            calcConvect=calcConvect.AND. .NOT.offlineLoadConvec
259          ENDIF
260    #endif /* ALLOW_OFFLINE */
261    
262    #ifdef  ALLOW_OBCS
263          IF (useOBCS) THEN
264    C--   Calculate future values on open boundaries
265    C--   moved before SEAICE_MODEL call since SEAICE_MODEL needs seaice-obcs fields
266    # ifdef ALLOW_AUTODIFF_TAMC
267    CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
268    CADJ STORE salt  = comlev1, key=ikey_dynamics, kind=isbyte
269    # endif
270    # ifdef ALLOW_DEBUG
271           IF (debugMode) CALL DEBUG_CALL('OBCS_CALC',myThid)
272    # endif
273           CALL OBCS_CALC( myTime+deltaTClock, myIter+1,
274         I                 uVel, vVel, wVel, theta, salt, myThid )
275          ENDIF
276    #endif  /* ALLOW_OBCS */
277    
278    #ifdef ALLOW_OCN_COMPON_INTERF
279    C--    Apply imported data (from coupled interface) to forcing fields
280    C jmc: moved here before any freezing/seaice pkg adjustment of surf-fluxes
281          IF ( useCoupler ) THEN
282             CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
283          ENDIF
284    #endif /* ALLOW_OCN_COMPON_INTERF */
285    
286    #ifdef ALLOW_AUTODIFF
287    # ifdef ALLOW_SALT_PLUME
288          DO bj=myByLo(myThid),myByHi(myThid)
289           DO bi=myBxLo(myThid),myBxHi(myThid)
290            DO j=1-OLy,sNy+OLy
291             DO i=1-OLx,sNx+OLx
292              saltPlumeDepth(i,j,bi,bj) = 0. _d 0
293              saltPlumeFlux(i,j,bi,bj)  = 0. _d 0
294             ENDDO
295            ENDDO
296           ENDDO
297          ENDDO
298    # endif
299    # ifdef ALLOW_ECCO
300    #  ifdef ALLOW_SIGMAR_COST_CONTRIBUTION
301          DO bj=myByLo(myThid),myByHi(myThid)
302           DO bi=myBxLo(myThid),myBxHi(myThid)
303            DO k=1,Nr
304             DO j=1-OLy,sNy+OLy
305              DO i=1-OLx,sNx+OLx
306               sigmaRfield(i,j,k,bi,bj) = 0. _d 0
307              ENDDO
308             ENDDO
309            ENDDO
310           ENDDO
311          ENDDO
312    #  endif
313    # endif
314    #endif /* ALLOW_AUTODIFF */
315    
316    #ifdef ALLOW_FRAZIL
317          IF ( useFRAZIL ) THEN
318    C--   Freeze water in the ocean interior and let it rise to the surface
319    CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
320    CADJ STORE salt  = comlev1, key=ikey_dynamics, kind=isbyte
321           CALL FRAZIL_CALC_RHS( myTime, myIter, myThid )
322          ENDIF
323    #endif /* ALLOW_FRAZIL */
324    
325    #ifndef OLD_THSICE_CALL_SEQUENCE
326    #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
327          IF ( useThSIce .AND. fluidIsWater ) THEN
328    # ifdef ALLOW_AUTODIFF_TAMC
329    CADJ STORE uice,vice         = comlev1, key = ikey_dynamics,
330    CADJ &     kind = isbyte
331    CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
332    CADJ &     kind = isbyte
333    CADJ STORE snowHeight, Tsrf  = comlev1, key = ikey_dynamics,
334    CADJ &     kind = isbyte
335    CADJ STORE Qice1, Qice2      = comlev1, key = ikey_dynamics,
336    CADJ &     kind = isbyte
337    CADJ STORE sHeating, snowAge = comlev1, key = ikey_dynamics,
338    CADJ &     kind = isbyte
339    CADJ STORE hocemxl = comlev1, key = ikey_dynamics,
340    CADJ &     kind = isbyte
341    CADJ STORE icflxsw = comlev1, key = ikey_dynamics,
342    CADJ &     kind = isbyte
343    CADJ STORE salt,theta        = comlev1, key = ikey_dynamics,
344    CADJ &     kind = isbyte
345    CADJ STORE uvel,vvel         = comlev1, key = ikey_dynamics,
346    CADJ &     kind = isbyte
347    CADJ STORE qnet,qsw, empmr   = comlev1, key = ikey_dynamics,
348    CADJ &     kind = isbyte
349    CADJ STORE atemp,aqh,precip  = comlev1, key = ikey_dynamics,
350    CADJ &     kind = isbyte
351    CADJ STORE swdown,lwdown     = comlev1, key = ikey_dynamics,
352    CADJ &     kind = isbyte
353    #  ifdef NONLIN_FRSURF
354    CADJ STORE hFac_surfC       = comlev1, key = ikey_dynamics,
355    CADJ &     kind = isbyte
356    #  endif
357    # endif /* ALLOW_AUTODIFF_TAMC */
358    # ifdef ALLOW_DEBUG
359            IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
360    # endif
361    C--     Step forward Therm.Sea-Ice variables
362    C       and modify forcing terms including effects from ice
363            CALL TIMER_START('THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
364            CALL THSICE_MAIN( myTime, myIter, myThid )
365            CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
366          ENDIF
367    #endif /* ALLOW_THSICE */
368    #endif /* ndef OLD_THSICE_CALL_SEQUENCE */
369    
370    #ifdef ALLOW_SEAICE
371    # ifdef ALLOW_AUTODIFF
372    CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
373    CADJ STORE fu,fv  = comlev1, key=ikey_dynamics, kind=isbyte
374    CADJ STORE qnet   = comlev1, key=ikey_dynamics, kind=isbyte
375    CADJ STORE qsw    = comlev1, key=ikey_dynamics, kind=isbyte
376    CADJ STORE theta  = comlev1, key=ikey_dynamics, kind=isbyte
377    CADJ STORE salt   = comlev1, key=ikey_dynamics, kind=isbyte
378    #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
379    CADJ STORE evap   = comlev1, key=ikey_dynamics, kind=isbyte
380    #endif
381          IF ( .NOT.useSEAICE .AND. SEAICEadjMODE .EQ. -1 ) THEN
382            CALL SEAICE_FAKE( myTime, myIter, myThid )
383          ENDIF
384    CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
385    CADJ STORE fu,fv  = comlev1, key=ikey_dynamics, kind=isbyte
386    CADJ STORE qnet   = comlev1, key=ikey_dynamics, kind=isbyte
387    CADJ STORE qsw    = comlev1, key=ikey_dynamics, kind=isbyte
388    CADJ STORE theta  = comlev1, key=ikey_dynamics, kind=isbyte
389    CADJ STORE salt   = comlev1, key=ikey_dynamics, kind=isbyte
390    #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
391    CADJ STORE evap   = comlev1, key=ikey_dynamics, kind=isbyte
392    #endif
393    # endif /* ALLOW_AUTODIFF */
394    #endif /* ALLOW_SEAICE */
395    
396  #ifdef ALLOW_SEAICE  #ifdef ALLOW_SEAICE
 C--   Call sea ice model to compute forcing/external data fields.  In  
 C     addition to computing prognostic sea-ice variables and diagnosing the  
 C     forcing/external data fields that drive the ocean model, SEAICE_MODEL  
 C     also sets theta to the freezing point under sea-ice.  The implied  
 C     surface heat flux is then stored in variable surfaceTendencyTice,  
 C     which is needed by KPP package (kpp_calc.F and kpp_transport_t.F)  
 C     to diagnose surface buoyancy fluxes and for the non-local transport  
 C     term.  Because this call precedes model thermodynamics, temperature  
 C     under sea-ice may not be "exactly" at the freezing point by the time  
 C     theta is dumped or time-averaged.  
397        IF ( useSEAICE ) THEN        IF ( useSEAICE ) THEN
398  #ifdef ALLOW_AUTODIFF_TAMC  # ifdef ALLOW_AUTODIFF_TAMC
399  CADJ STORE atemp,aqh,precip    = comlev1, key = ikey_dynamics  cph-adj-test(
400  CADJ STORE swdown,lwdown       = comlev1, key = ikey_dynamics  CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
401  CADJ STORE theta               = comlev1, key = ikey_dynamics  CADJ STORE hsnow  = comlev1, key=ikey_dynamics, kind=isbyte
402  cph# ifdef EXF_READ_EVAP  CADJ STORE heff   = comlev1, key=ikey_dynamics, kind=isbyte
403  CADJ STORE evap                = comlev1, key = ikey_dynamics  CADJ STORE tices  = comlev1, key=ikey_dynamics, kind=isbyte
404  cph# endif  CADJ STORE empmr, qnet  = comlev1, key=ikey_dynamics, kind=isbyte
405  cph# ifdef SEAICE_ALLOW_DYNAMICS  CADJ STORE qsw,saltflux = comlev1, key=ikey_dynamics, kind=isbyte
406  CADJ STORE uvel,vvel           = comlev1, key = ikey_dynamics  CADJ STORE fu, fv = comlev1, key=ikey_dynamics, kind=isbyte
407  cph# endif  cCADJ STORE theta  = comlev1, key=ikey_dynamics, kind=isbyte
408  # ifdef SEAICE_CGRID  cCADJ STORE salt   = comlev1, key=ikey_dynamics, kind=isbyte
409  CADJ STORE fu, fv              = comlev1, key = ikey_dynamics  cph-adj-test)
410  CADJ STORE seaicemasku         = comlev1, key = ikey_dynamics  c#ifdef ALLOW_EXF
411  CADJ STORE seaicemaskv         = comlev1, key = ikey_dynamics  CADJ STORE atemp,aqh,precip    = comlev1, key = ikey_dynamics,
412    CADJ &     kind = isbyte
413    CADJ STORE swdown,lwdown       = comlev1, key = ikey_dynamics,
414    CADJ &     kind = isbyte
415    CADJ STORE evap                = comlev1, key = ikey_dynamics,
416    CADJ &     kind = isbyte
417    CADJ STORE uwind,vwind         = comlev1, key = ikey_dynamics,
418    CADJ &     kind = isbyte
419    c#endif
420    CADJ STORE uvel,vvel           = comlev1, key = ikey_dynamics,
421    CADJ &     kind = isbyte
422    #  ifdef SEAICE_CGRID
423    CADJ STORE stressdivergencex   = comlev1, key = ikey_dynamics,
424    CADJ &     kind = isbyte
425    CADJ STORE stressdivergencey   = comlev1, key = ikey_dynamics,
426    CADJ &     kind = isbyte
427    #  endif
428    #  ifdef SEAICE_ALLOW_DYNAMICS
429    CADJ STORE uice                = comlev1, key = ikey_dynamics,
430    CADJ &     kind = isbyte
431    CADJ STORE vice                = comlev1, key = ikey_dynamics,
432    CADJ &     kind = isbyte
433    CADJ STORE dwatn               = comlev1, key = ikey_dynamics,
434    CADJ &     kind = isbyte
435    #   ifdef SEAICE_ALLOW_EVP
436    CADJ STORE seaice_sigma1       = comlev1, key = ikey_dynamics,
437    CADJ &     kind = isbyte
438    CADJ STORE seaice_sigma2       = comlev1, key = ikey_dynamics,
439    CADJ &     kind = isbyte
440    CADJ STORE seaice_sigma12      = comlev1, key = ikey_dynamics,
441    CADJ &     kind = isbyte
442    #   endif
443    #  endif
444    #  ifdef SEAICE_VARIABLE_SALINITY
445    CADJ STORE hsalt               = comlev1, key = ikey_dynamics,
446    CADJ &     kind = isbyte
447    #  endif
448  #  ifdef ATMOSPHERIC_LOADING  #  ifdef ATMOSPHERIC_LOADING
449  CADJ STORE siceload            = comlev1, key = ikey_dynamics  CADJ STORE pload               = comlev1, key = ikey_dynamics,
450    CADJ &     kind = isbyte
451    CADJ STORE siceload            = comlev1, key = ikey_dynamics,
452    CADJ &     kind = isbyte
453    #  endif
454    #  ifdef NONLIN_FRSURF
455    CADJ STORE recip_hfacc         = comlev1, key = ikey_dynamics,
456    CADJ &     kind = isbyte
457  #  endif  #  endif
458    #  ifdef ANNUAL_BALANCE
459    CADJ STORE balance_itcount     = comlev1, key = ikey_dynamics,
460    CADJ &     kind = isbyte
461    #  endif /* ANNUAL_BALANCE */
462    #  ifdef ALLOW_THSICE
463    C-- store thSIce vars before advection (called from SEAICE_MODEL) update them:
464    CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
465    CADJ &     kind = isbyte
466    CADJ STORE snowHeight, hOceMxL = comlev1, key = ikey_dynamics,
467    CADJ &     kind = isbyte
468    CADJ STORE Qice1, Qice2  = comlev1, key = ikey_dynamics,
469    CADJ &     kind = isbyte
470    #  endif /* ALLOW_THSICE */
471    # endif /* ALLOW_AUTODIFF_TAMC */
472    # ifdef ALLOW_DEBUG
473            IF (debugMode) CALL DEBUG_CALL('SEAICE_MODEL',myThid)
474  # endif  # endif
 #endif  
 #ifdef ALLOW_DEBUG  
         IF ( debugLevel .GE. debLevB )  
      &    CALL DEBUG_CALL('SEAICE_MODEL',myThid)  
 #endif  
475          CALL TIMER_START('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)          CALL TIMER_START('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
476          CALL SEAICE_MODEL( myTime, myIter, myThid )          CALL SEAICE_MODEL( myTime, myIter, myThid )
477          CALL TIMER_STOP ('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)          CALL TIMER_STOP ('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
478  #ifdef ALLOW_COST_ICE  # ifdef ALLOW_COST
479          CALL COST_ICE_TEST ( myTime, myIter, myThid )          CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
480  #endif  # endif
481        ENDIF        ENDIF
482  #endif /* ALLOW_SEAICE */  #endif /* ALLOW_SEAICE */
483    
484    #if (defined ALLOW_OCN_COMPON_INTERF) && (defined ALLOW_THSICE)
485    C--   After seaice-dyn and advection of pkg/thsice fields,
486    C     Export ocean coupling fields to coupled interface (only with pkg/thsice)
487          IF ( useCoupler ) THEN
488    # ifdef ALLOW_DEBUG
489            IF (debugMode) CALL DEBUG_CALL('OCN_EXPORT_DATA',myThid)
490    # endif
491             CALL TIMER_START('OCN_EXPORT_DATA [DO_OCEANIC_PHYS]', myThid)
492             CALL OCN_EXPORT_DATA( myTime, myIter, myThid )
493             CALL TIMER_STOP ('OCN_EXPORT_DATA [DO_OCEANIC_PHYS]', myThid)
494          ENDIF
495    #endif /* ALLOW_OCN_COMPON_INTERF & ALLOW_THSICE */
496    
497    #ifdef ALLOW_AUTODIFF_TAMC
498    CADJ STORE sst, sss           = comlev1, key = ikey_dynamics,
499    CADJ &     kind = isbyte
500    CADJ STORE qsw                = comlev1, key = ikey_dynamics,
501    CADJ &     kind = isbyte
502    # ifdef ALLOW_SEAICE
503    CADJ STORE area               = comlev1, key = ikey_dynamics,
504    CADJ &     kind = isbyte
505    # endif
506    #endif
507    
508    #ifdef OLD_THSICE_CALL_SEQUENCE
509  #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)  #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
510        IF ( useThSIce .AND. fluidIsWater ) THEN        IF ( useThSIce .AND. fluidIsWater ) THEN
511  #ifdef ALLOW_DEBUG  # ifdef ALLOW_AUTODIFF_TAMC
512          IF ( debugLevel .GE. debLevB )  cph(
513       &    CALL DEBUG_CALL('THSICE_MAIN',myThid)  #  ifdef NONLIN_FRSURF
514  #endif  CADJ STORE uice,vice        = comlev1, key = ikey_dynamics,
515    CADJ &     kind = isbyte
516    CADJ STORE salt,theta       = comlev1, key = ikey_dynamics,
517    CADJ &     kind = isbyte
518    CADJ STORE qnet,qsw, empmr  = comlev1, key = ikey_dynamics,
519    CADJ &     kind = isbyte
520    CADJ STORE hFac_surfC       = comlev1, key = ikey_dynamics,
521    CADJ &     kind = isbyte
522    #  endif
523    # endif
524    # ifdef ALLOW_DEBUG
525            IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
526    # endif
527  C--     Step forward Therm.Sea-Ice variables  C--     Step forward Therm.Sea-Ice variables
528  C       and modify forcing terms including effects from ice  C       and modify forcing terms including effects from ice
529          CALL TIMER_START('THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)          CALL TIMER_START('THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
# Line 187  C       and modify forcing terms includi Line 531  C       and modify forcing terms includi
531          CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)          CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
532        ENDIF        ENDIF
533  #endif /* ALLOW_THSICE */  #endif /* ALLOW_THSICE */
534    #endif /* OLD_THSICE_CALL_SEQUENCE */
535    
536  #ifdef ALLOW_SHELFICE  #ifdef ALLOW_SHELFICE
537        IF ( useShelfIce .AND. fluidIsWater ) THEN        IF ( useShelfIce .AND. fluidIsWater ) THEN
538  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
539          IF ( debugLevel .GE. debLevB )         IF (debugMode) CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
540       &    CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)  #endif
541    #ifdef ALLOW_AUTODIFF_TAMC
542    CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
543    CADJ &     kind = isbyte
544    CADJ STORE uvel, vvel = comlev1, key = ikey_dynamics,
545    CADJ &     kind = isbyte
546  #endif  #endif
547  C     compute temperature and (virtual) salt flux at the  C     compute temperature and (virtual) salt flux at the
548  C     shelf-ice ocean interface  C     shelf-ice ocean interface
549         CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',         CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
550       &       myThid)       &       myThid)
# Line 204  C     shelf-ice ocean interface Line 554  C     shelf-ice ocean interface
554        ENDIF        ENDIF
555  #endif /* ALLOW_SHELFICE */  #endif /* ALLOW_SHELFICE */
556    
557  C--   Freeze water at the surface  #ifdef ALLOW_ICEFRONT
558  #ifdef ALLOW_AUTODIFF_TAMC        IF ( useICEFRONT .AND. fluidIsWater ) THEN
559  CADJ STORE theta = comlev1, key = ikey_dynamics  #ifdef ALLOW_DEBUG
560           IF (debugMode) CALL DEBUG_CALL('ICEFRONT_THERMODYNAMICS',myThid)
561  #endif  #endif
562        IF ( allowFreezing  C     compute temperature and (virtual) salt flux at the
563       &                   .AND. .NOT. useSEAICE  C     ice-front ocean interface
564       &                   .AND. .NOT. useThSIce ) THEN         CALL TIMER_START('ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
565          CALL FREEZE_SURFACE(  myTime, myIter, myThid )       &       myThid)
566           CALL ICEFRONT_THERMODYNAMICS( myTime, myIter, myThid )
567           CALL TIMER_STOP( 'ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
568         &      myThid)
569        ENDIF        ENDIF
570    #endif /* ALLOW_ICEFRONT */
571    
572  #ifdef ALLOW_OCN_COMPON_INTERF  #ifdef ALLOW_SALT_PLUME
573  C--    Apply imported data (from coupled interface) to forcing fields        IF ( useSALT_PLUME ) THEN
574  C jmc: do not know precisely where to put this call (bf or af thSIce ?)  Catn: exchanging saltPlumeFlux:
575        IF ( useCoupler ) THEN            CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
          CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )  
576        ENDIF        ENDIF
577  #endif /* ALLOW_OCN_COMPON_INTERF */  #endif /* ALLOW_SALT_PLUME */
578    
579  #ifdef ALLOW_BALANCE_FLUXES  C--   Freeze water at the surface
580  C     balance fluxes        IF ( allowFreezing ) THEN
581        IF ( balanceEmPmR )  #ifdef ALLOW_AUTODIFF_TAMC
582       &        CALL REMOVE_MEAN_RS( 1, EmPmR, maskH, maskH, rA, drF,  CADJ STORE theta = comlev1, key = ikey_dynamics,
583       &        'EmPmR', myTime, myThid )  CADJ &     kind = isbyte
584        IF ( balanceQnet )  #endif
585       &        CALL REMOVE_MEAN_RS( 1, Qnet,  maskH, maskH, rA, drF,          CALL FREEZE_SURFACE( myTime, myIter, myThid )
586       &        'Qnet ', myTime, myThid )        ENDIF
587  #endif /* ALLOW_BALANCE_FLUXES */  
588          iMin = 1-OLx
589          iMax = sNx+OLx
590          jMin = 1-OLy
591          jMax = sNy+OLy
592    
593    C---  Determines forcing terms based on external fields
594    C     relaxation terms, etc.
595    #ifdef ALLOW_AUTODIFF
596    CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
597    CADJ &     kind = isbyte
598    #else  /* ALLOW_AUTODIFF */
599    C--   if fluid is not water, by-pass surfaceForcing, find_rho, gmredi
600    C     and all vertical mixing schemes, but keep OBCS_CALC
601          IF ( fluidIsWater ) THEN
602    #endif /* ALLOW_AUTODIFF */
603    #ifdef ALLOW_DEBUG
604          IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
605    #endif
606            CALL EXTERNAL_FORCING_SURF(
607         I             iMin, iMax, jMin, jMax,
608         I             myTime, myIter, myThid )
609    
610  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
611  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 254  CHPF$ INDEPENDENT Line 629  CHPF$ INDEPENDENT
629            itdkey = (act1 + 1) + act2*max1            itdkey = (act1 + 1) + act2*max1
630       &                      + act3*max1*max2       &                      + act3*max1*max2
631       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
 #else  /* ALLOW_AUTODIFF_TAMC */  
 C     if fluid is not water, by-pass find_rho, gmredi, surfaceForcing  
 C     and all vertical mixing schemes, but keep OBCS_CALC  
         IF ( fluidIsWater ) THEN  
632  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
633    
634  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
# Line 265  C     These inital values do not alter t Line 636  C     These inital values do not alter t
636  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
637  C     point numbers. This prevents spurious hardware signals due to  C     point numbers. This prevents spurious hardware signals due to
638  C     uninitialised but inert locations.  C     uninitialised but inert locations.
   
         DO j=1-OLy,sNy+OLy  
          DO i=1-OLx,sNx+OLx  
           rhok   (i,j)   = 0. _d 0  
           rhoKM1 (i,j)   = 0. _d 0  
           rhoKP1 (i,j)   = 0. _d 0  
          ENDDO  
         ENDDO  
   
639          DO k=1,Nr          DO k=1,Nr
640           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
641            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
642  C This is currently also used by IVDC and Diagnostics  C This is currently used by GMRedi, IVDC, MXL-depth  and Diagnostics
643             sigmaX(i,j,k) = 0. _d 0             sigmaX(i,j,k) = 0. _d 0
644             sigmaY(i,j,k) = 0. _d 0             sigmaY(i,j,k) = 0. _d 0
645             sigmaR(i,j,k) = 0. _d 0             sigmaR(i,j,k) = 0. _d 0
646  #ifdef ALLOW_AUTODIFF_TAMC            ENDDO
647             ENDDO
648            ENDDO
649    
650    #ifdef ALLOW_AUTODIFF
651            DO j=1-OLy,sNy+OLy
652             DO i=1-OLx,sNx+OLx
653              rhoKm1 (i,j)   = 0. _d 0
654              rhoKp1 (i,j)   = 0. _d 0
655             ENDDO
656            ENDDO
657  cph all the following init. are necessary for TAF  cph all the following init. are necessary for TAF
658  cph although some of these are re-initialised later.  cph although some of these are re-initialised later.
659            DO k=1,Nr
660             DO j=1-OLy,sNy+OLy
661              DO i=1-OLx,sNx+OLx
662               rhoInSitu(i,j,k,bi,bj) = 0.
663    # ifdef ALLOW_GGL90
664               GGL90viscArU(i,j,k,bi,bj)  = 0. _d 0
665               GGL90viscArV(i,j,k,bi,bj)  = 0. _d 0
666               GGL90diffKr(i,j,k,bi,bj)  = 0. _d 0
667    # endif /* ALLOW_GGL90 */
668    # ifdef ALLOW_SALT_PLUME
669    #  ifdef SALT_PLUME_VOLUME
670               SPforcingS(i,j,k,bi,bj) = 0. _d 0
671               SPforcingT(i,j,k,bi,bj) = 0. _d 0
672    #  endif
673    # endif /* ALLOW_SALT_PLUME */
674              ENDDO
675             ENDDO
676            ENDDO
677    #ifdef ALLOW_OFFLINE
678           IF ( calcConvect ) THEN
679    #endif
680            DO k=1,Nr
681             DO j=1-OLy,sNy+OLy
682              DO i=1-OLx,sNx+OLx
683             IVDConvCount(i,j,k,bi,bj) = 0.             IVDConvCount(i,j,k,bi,bj) = 0.
684              ENDDO
685             ENDDO
686            ENDDO
687    #ifdef ALLOW_OFFLINE
688           ENDIF
689           IF ( calcGMRedi ) THEN
690    #endif
691  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
692            DO k=1,Nr
693             DO j=1-OLy,sNy+OLy
694              DO i=1-OLx,sNx+OLx
695             Kwx(i,j,k,bi,bj)  = 0. _d 0             Kwx(i,j,k,bi,bj)  = 0. _d 0
696             Kwy(i,j,k,bi,bj)  = 0. _d 0             Kwy(i,j,k,bi,bj)  = 0. _d 0
697             Kwz(i,j,k,bi,bj)  = 0. _d 0             Kwz(i,j,k,bi,bj)  = 0. _d 0
# Line 304  cph although some of these are re-initia Line 710  cph although some of these are re-initia
710  #  ifdef GM_VISBECK_VARIABLE_K  #  ifdef GM_VISBECK_VARIABLE_K
711             VisbeckK(i,j,bi,bj)   = 0. _d 0             VisbeckK(i,j,bi,bj)   = 0. _d 0
712  #  endif  #  endif
713              ENDDO
714             ENDDO
715            ENDDO
716  # endif /* ALLOW_GMREDI */  # endif /* ALLOW_GMREDI */
717    #ifdef ALLOW_OFFLINE
718           ENDIF
719           IF ( calcKPP ) THEN
720    #endif
721  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
722            DO k=1,Nr
723             DO j=1-OLy,sNy+OLy
724              DO i=1-OLx,sNx+OLx
725             KPPdiffKzS(i,j,k,bi,bj)  = 0. _d 0             KPPdiffKzS(i,j,k,bi,bj)  = 0. _d 0
726             KPPdiffKzT(i,j,k,bi,bj)  = 0. _d 0             KPPdiffKzT(i,j,k,bi,bj)  = 0. _d 0
 # endif /* ALLOW_KPP */  
 #endif /* ALLOW_AUTODIFF_TAMC */  
727            ENDDO            ENDDO
728           ENDDO           ENDDO
729          ENDDO          ENDDO
730    # endif /* ALLOW_KPP */
731          iMin = 1-OLx  #ifdef ALLOW_OFFLINE
732          iMax = sNx+OLx         ENDIF
733          jMin = 1-OLy  #endif
734          jMax = sNy+OLy  #endif /* ALLOW_AUTODIFF */
735    
736  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
737  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
738  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     kind = isbyte
739    CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
740    CADJ &     kind = isbyte
741  CADJ STORE totphihyd(:,:,:,bi,bj)  CADJ STORE totphihyd(:,:,:,bi,bj)
742  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
743    CADJ &     kind = isbyte
744  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
745  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
746  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     kind = isbyte
747    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
748    CADJ &     kind = isbyte
749    # endif
750    # ifdef ALLOW_SALT_PLUME
751    CADJ STORE saltplumedepth(:,:,bi,bj) = comlev1_bibj, key=itdkey,
752    CADJ &     kind = isbyte
753    CADJ STORE saltplumeflux(:,:,bi,bj) = comlev1_bibj, key=itdkey,
754    CADJ &     kind = isbyte
755  # endif  # endif
756  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
757    
758    C--   Always compute density (stored in common block) here; even when it is not
759    C     needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
760  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
761          IF ( debugLevel .GE. debLevB )          IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D (xNr)',myThid)
762       &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)  #endif
763    #ifdef ALLOW_AUTODIFF
764            IF ( fluidIsWater ) THEN
765    #endif /* ALLOW_AUTODIFF */
766    #ifdef ALLOW_DOWN_SLOPE
767             IF ( useDOWN_SLOPE ) THEN
768               DO k=1,Nr
769                CALL DWNSLP_CALC_RHO(
770         I                  theta, salt,
771         O                  rhoInSitu(1-OLx,1-OLy,k,bi,bj),
772         I                  k, bi, bj, myTime, myIter, myThid )
773               ENDDO
774             ENDIF
775    #endif /* ALLOW_DOWN_SLOPE */
776    #ifdef ALLOW_BBL
777             IF ( useBBL ) THEN
778    C     pkg/bbl requires in-situ bbl density for depths equal to and deeper than the bbl.
779    C     To reduce computation and storage requirement, these densities are stored in the
780    C     dry grid boxes of rhoInSitu.  See BBL_CALC_RHO for details.
781               DO k=Nr,1,-1
782                CALL BBL_CALC_RHO(
783         I                  theta, salt,
784         O                  rhoInSitu,
785         I                  k, bi, bj, myTime, myIter, myThid )
786    
787               ENDDO
788             ENDIF
789    #endif /* ALLOW_BBL */
790             IF ( .NOT. ( useDOWN_SLOPE .OR. useBBL ) ) THEN
791               DO k=1,Nr
792                CALL FIND_RHO_2D(
793         I                iMin, iMax, jMin, jMax, k,
794         I                theta(1-OLx,1-OLy,k,bi,bj),
795         I                salt (1-OLx,1-OLy,k,bi,bj),
796         O                rhoInSitu(1-OLx,1-OLy,k,bi,bj),
797         I                k, bi, bj, myThid )
798               ENDDO
799             ENDIF
800    #ifdef ALLOW_AUTODIFF
801            ELSE
802    C-        fluid is not water:
803              DO k=1,Nr
804               IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN
805    C-    isothermal (theta=const) reference state
806                 thetaRef = thetaConst
807               ELSE
808    C-    horizontally uniform (tRef) reference state
809                 thetaRef = tRef(k)
810               ENDIF
811               DO j=1-OLy,sNy+OLy
812                DO i=1-OLx,sNx+OLx
813                 rhoInSitu(i,j,k,bi,bj) =
814         &         ( theta(i,j,k,bi,bj)
815         &              *( salt(i,j,k,bi,bj)*atm_Rq + oneRL )
816         &         - thetaRef )*maskC(i,j,k,bi,bj)
817                ENDDO
818               ENDDO
819              ENDDO
820            ENDIF
821    #endif /* ALLOW_AUTODIFF */
822    
823    #ifdef ALLOW_DEBUG
824            IF (debugMode) THEN
825              WRITE(msgBuf,'(A,2(I4,A))')
826         &         'ENTERING UPWARD K LOOP (bi=', bi, ', bj=', bj,')'
827              CALL DEBUG_MSG(msgBuf(1:43),myThid)
828            ENDIF
829  #endif  #endif
830    
831  C--     Start of diagnostic loop  C--     Start of diagnostic loop
# Line 341  C--     Start of diagnostic loop Line 834  C--     Start of diagnostic loop
834  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
835  C? Patrick, is this formula correct now that we change the loop range?  C? Patrick, is this formula correct now that we change the loop range?
836  C? Do we still need this?  C? Do we still need this?
837  cph kkey formula corrected.  cph kkey formula corrected.
838  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhoK, rhoKm1, in the case useGMREDI.
839           kkey = (itdkey-1)*Nr + k            kkey = (itdkey-1)*Nr + k
840  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
841    
842    c#ifdef ALLOW_AUTODIFF_TAMC
843    cCADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
844    cCADJ &     kind = isbyte
845    cCADJ STORE salt(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey,
846    cCADJ &     kind = isbyte
847    c#endif /* ALLOW_AUTODIFF_TAMC */
848    
849  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
850  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
851            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)            IF ( calcGMRedi .OR. (k.GT.1 .AND. calcConvect)
852       &                   .OR. doDiagsRho.GE.1 ) THEN       &         .OR. usePP81 .OR. useKL10
853  #ifdef ALLOW_DEBUG       &         .OR. useMY82 .OR. useGGL90
854              IF ( debugLevel .GE. debLevB )       &         .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
      &       CALL DEBUG_CALL('FIND_RHO',myThid)  
 #endif  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
   
855              IF (k.GT.1) THEN              IF (k.GT.1) THEN
856  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
857  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
858  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ &     kind = isbyte
859    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
860    CADJ &     kind = isbyte
861    CADJ STORE rhokm1 (bi,bj)       = comlev1_bibj_k, key=kkey,
862    CADJ &     kind = isbyte
863  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
864               CALL FIND_RHO(               CALL FIND_RHO_2D(
865       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,       I                 iMin, iMax, jMin, jMax, k,
866       I        theta, salt,       I                 theta(1-OLx,1-OLy,k-1,bi,bj),
867       O        rhoKm1,       I                 salt (1-OLx,1-OLy,k-1,bi,bj),
868       I        myThid )       O                 rhoKm1,
869         I                 k-1, bi, bj, myThid )
870              ENDIF              ENDIF
871  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
872              IF ( debugLevel .GE. debLevB )              IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
      &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)  
873  #endif  #endif
874  cph Avoid variable aliasing for adjoint !!!  cph Avoid variable aliasing for adjoint !!!
875              DO j=jMin,jMax              DO j=jMin,jMax
876               DO i=iMin,iMax               DO i=iMin,iMax
877                rhoKP1(i,j) = rhoK(i,j)                rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
878               ENDDO               ENDDO
879              ENDDO              ENDDO
880              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
881       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
882       I             rhoK, rhoKm1, rhoKp1,       I             rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
883       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
884       I             myThid )       I             myThid )
885    #ifdef ALLOW_ECCO
886    # ifdef ALLOW_SIGMAR_COST_CONTRIBUTION
887                DO j=jMin,jMax
888                 DO i=iMin,iMax
889                  sigmaRfield(i,j,k,bi,bj)=sigmaR(i,j,k)
890                 ENDDO
891                ENDDO
892    # endif
893    #endif /* ALLOW_ECCO */
894    #ifdef ALLOW_AUTODIFF
895    #ifdef GMREDI_WITH_STABLE_ADJOINT
896    cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
897    cgf -> cuts adjoint dependency from slope to state
898                CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
899                CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
900                CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
901    #endif
902    #endif /* ALLOW_AUTODIFF */
903            ENDIF            ENDIF
904    
 #ifdef ALLOW_AUTODIFF_TAMC  
 ctest# ifndef GM_EXCLUDE_CLIPPING  
 CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 ctest# endif  
 CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
905  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
906  c ==> should use sigmaR !!!            IF (k.GT.1 .AND. calcConvect) THEN
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
907  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
908              IF ( debugLevel .GE. debLevB )              IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid)
      &       CALL DEBUG_CALL('CALC_IVDC',myThid)  
909  #endif  #endif
910              CALL CALC_IVDC(              CALL CALC_IVDC(
911       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
912       I        rhoKm1, rhoK,       I        sigmaR,
913       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
914            ENDIF            ENDIF
915    
916  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
917            IF ( doDiagsRho.GE.2 ) THEN            IF ( doDiagsRho.GE.4 ) THEN
918              CALL DIAGS_RHO( k, bi, bj,              CALL DIAGS_RHO_L( doDiagsRho, k, bi, bj,
919       I                      rhoK, rhoKm1,       I                        rhoInSitu(1-OLx,1-OLy,1,bi,bj),
920       I                      myTime, myIter, myThid)       I                        rhoKm1, wVel,
921         I                        myTime, myIter, myThid )
922            ENDIF            ENDIF
923  #endif  #endif
924    
925  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
926          ENDDO          ENDDO
927    
928    #ifdef ALLOW_AUTODIFF_TAMC
929    CADJ STORE IVDConvCount(:,:,:,bi,bj)
930    CADJ &     = comlev1_bibj, key=itdkey,
931    CADJ &     kind = isbyte
932    #endif
933    
934    C--     Diagnose Mixed Layer Depth:
935            IF ( calcGMRedi .OR. MOD(doDiagsRho,2).EQ.1 ) THEN
936              CALL CALC_OCE_MXLAYER(
937         I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
938         I              bi, bj, myTime, myIter, myThid )
939            ENDIF
940    
941    #ifdef ALLOW_SALT_PLUME
942            IF ( useSALT_PLUME ) THEN
943              CALL SALT_PLUME_CALC_DEPTH(
944         I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
945         I              bi, bj, myTime, myIter, myThid )
946    #ifdef SALT_PLUME_VOLUME
947              CALL SALT_PLUME_VOLFRAC(
948         I              bi, bj, myTime, myIter, myThid )
949    C-- get forcings for kpp
950              CALL SALT_PLUME_APPLY(
951         I              1, bi, bj, recip_hFacC(1-OLx,1-OLy,1,bi,bj),
952         I              theta, 0,
953         I              myTime, myIter, myThid )
954              CALL SALT_PLUME_APPLY(
955         I              2, bi, bj, recip_hFacC(1-OLx,1-OLy,1,bi,bj),
956         I              salt, 0,
957         I              myTime, myIter, myThid )
958    C-- need to call this S/R from here to apply just before kpp
959              CALL SALT_PLUME_FORCING_SURF(
960         I              bi, bj, iMin, iMax, jMin, jMax,
961         I              myTime, myIter, myThid )
962    #endif /* SALT_PLUME_VOLUME */
963            ENDIF
964    #endif /* ALLOW_SALT_PLUME */
965    
966  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
967          IF ( doDiagsRho.GE.1 ) THEN          IF ( MOD(doDiagsRho,4).GE.2 ) THEN
968            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,
969       &         2, bi, bj, myThid)       &         2, bi, bj, myThid)
970          ENDIF          ENDIF
971  #endif  #endif /* ALLOW_DIAGNOSTICS */
972    
973  C--     Determines forcing terms based on external fields  C--    This is where EXTERNAL_FORCING_SURF(bi,bj) used to be called;
974  C       relaxation terms, etc.  C      now called earlier, before bi,bj loop.
 #ifdef ALLOW_DEBUG  
         IF ( debugLevel .GE. debLevB )  
      &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)  
 #endif  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE EmPmR(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 # ifdef EXACT_CONSERV  
 CADJ STORE PmEpR(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 # endif  
 # ifdef NONLIN_FRSURF  
 CADJ STORE hFac_surfC(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE recip_hFacC(:,:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 # endif  
 #endif  
         CALL EXTERNAL_FORCING_SURF(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             myTime, myIter, myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
 # ifdef EXACT_CONSERV  
 cph-test  
 cphCADJ STORE PmEpR(:,:,bi,bj)  
 cphCADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 # endif  
 #endif  
975    
976  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
977  cph needed for KPP  cph needed for KPP
978  CADJ STORE surfaceForcingU(:,:,bi,bj)  CADJ STORE surfaceForcingU(:,:,bi,bj)
979  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
980    CADJ &     kind = isbyte
981  CADJ STORE surfaceForcingV(:,:,bi,bj)  CADJ STORE surfaceForcingV(:,:,bi,bj)
982  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
983    CADJ &     kind = isbyte
984  CADJ STORE surfaceForcingS(:,:,bi,bj)  CADJ STORE surfaceForcingS(:,:,bi,bj)
985  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
986    CADJ &     kind = isbyte
987  CADJ STORE surfaceForcingT(:,:,bi,bj)  CADJ STORE surfaceForcingT(:,:,bi,bj)
988  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
989    CADJ &     kind = isbyte
990  CADJ STORE surfaceForcingTice(:,:,bi,bj)  CADJ STORE surfaceForcingTice(:,:,bi,bj)
991  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
992    CADJ &     kind = isbyte
993  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
994    
 #ifdef  ALLOW_GMREDI  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 # ifndef GM_EXCLUDE_CLIPPING  
 cph storing here is needed only for one GMREDI_OPTIONS:  
 cph define GM_BOLUS_ADVEC  
 cph keep it although TAF says you dont need to.  
 cph but I've avoided the #ifdef for now, in case more things change  
 CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  
 # endif  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  
         IF (useGMRedi) THEN  
 #ifdef ALLOW_DEBUG  
           IF ( debugLevel .GE. debLevB )  
      &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)  
 #endif  
           CALL GMREDI_CALC_TENSOR(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           CALL GMREDI_CALC_TENSOR_DUMMY(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
   
 #endif  /* ALLOW_GMREDI */  
   
995  #ifdef  ALLOW_KPP  #ifdef  ALLOW_KPP
996  C--     Compute KPP mixing coefficients  C--     Compute KPP mixing coefficients
997          IF (useKPP) THEN          IF ( calcKPP ) THEN
998  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
999            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid)
      &     CALL DEBUG_CALL('KPP_CALC',myThid)  
1000  #endif  #endif
1001              CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
1002            CALL KPP_CALC(            CALL KPP_CALC(
1003       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myIter, myThid )
1004  #ifdef ALLOW_AUTODIFF_TAMC            CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
1005    #if (defined ALLOW_AUTODIFF) && !(defined ALLOW_OFFLINE)
1006          ELSE          ELSE
1007            CALL KPP_CALC_DUMMY(            CALL KPP_CALC_DUMMY(
1008       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myIter, myThid )
1009  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF and not ALLOW_OFFLINE */
1010          ENDIF          ENDIF
   
1011  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
1012    
1013  #ifdef  ALLOW_PP81  #ifdef  ALLOW_PP81
1014  C--     Compute PP81 mixing coefficients  C--     Compute PP81 mixing coefficients
1015          IF (usePP81) THEN          IF (usePP81) THEN
1016  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
1017            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid)
      &     CALL DEBUG_CALL('PP81_CALC',myThid)  
1018  #endif  #endif
1019            CALL PP81_CALC(            CALL PP81_CALC(
1020       I                  bi, bj, myTime, myThid )       I                     bi, bj, sigmaR, myTime, myIter, myThid )
1021          ENDIF          ENDIF
1022  #endif /* ALLOW_PP81 */  #endif /* ALLOW_PP81 */
1023    
1024    #ifdef  ALLOW_KL10
1025    C--     Compute KL10 mixing coefficients
1026            IF (useKL10) THEN
1027    #ifdef ALLOW_DEBUG
1028              IF (debugMode) CALL DEBUG_CALL('KL10_CALC',myThid)
1029    #endif
1030              CALL KL10_CALC(
1031         I                     bi, bj, sigmaR, myTime, myIter, myThid )
1032            ENDIF
1033    #endif /* ALLOW_KL10 */
1034    
1035  #ifdef  ALLOW_MY82  #ifdef  ALLOW_MY82
1036  C--     Compute MY82 mixing coefficients  C--     Compute MY82 mixing coefficients
1037          IF (useMY82) THEN          IF (useMY82) THEN
1038  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
1039            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid)
      &     CALL DEBUG_CALL('MY82_CALC',myThid)  
1040  #endif  #endif
1041            CALL MY82_CALC(            CALL MY82_CALC(
1042       I                  bi, bj, myTime, myThid )       I                     bi, bj, sigmaR, myTime, myIter, myThid )
1043          ENDIF          ENDIF
1044  #endif /* ALLOW_MY82 */  #endif /* ALLOW_MY82 */
1045    
1046  #ifdef  ALLOW_GGL90  #ifdef  ALLOW_GGL90
1047    #ifdef ALLOW_AUTODIFF_TAMC
1048    CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
1049    CADJ &     kind = isbyte
1050    #endif /* ALLOW_AUTODIFF_TAMC */
1051  C--     Compute GGL90 mixing coefficients  C--     Compute GGL90 mixing coefficients
1052          IF (useGGL90) THEN          IF (useGGL90) THEN
1053  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
1054            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
      &     CALL DEBUG_CALL('GGL90_CALC',myThid)  
1055  #endif  #endif
1056              CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
1057            CALL GGL90_CALC(            CALL GGL90_CALC(
1058       I                  bi, bj, myTime, myThid )       I                     bi, bj, sigmaR, myTime, myIter, myThid )
1059              CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
1060          ENDIF          ENDIF
1061  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
1062    
# Line 567  C--     Compute GGL90 mixing coefficient Line 1064  C--     Compute GGL90 mixing coefficient
1064          IF ( taveFreq.GT. 0. _d 0 ) THEN          IF ( taveFreq.GT. 0. _d 0 ) THEN
1065            CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)            CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
1066          ENDIF          ENDIF
1067          IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN          IF ( taveFreq.GT.0. .AND. calcConvect ) THEN
1068            CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,            CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
1069       I                           Nr, deltaTclock, bi, bj, myThid)       I                           Nr, deltaTClock, bi, bj, myThid)
1070          ENDIF          ENDIF
1071  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
1072    
1073  #ifndef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_GMREDI
1074  C---  if fluid Is Water: end  #ifdef ALLOW_AUTODIFF_TAMC
1075          ENDIF  # ifndef GM_EXCLUDE_CLIPPING
1076  #endif  cph storing here is needed only for one GMREDI_OPTIONS:
1077    cph define GM_BOLUS_ADVEC
1078    cph keep it although TAF says you dont need to.
1079    cph but I have avoided the #ifdef for now, in case more things change
1080    CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey,
1081    CADJ &     kind = isbyte
1082    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey,
1083    CADJ &     kind = isbyte
1084    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey,
1085    CADJ &     kind = isbyte
1086    # endif
1087    #endif /* ALLOW_AUTODIFF_TAMC */
1088    
1089  #ifdef  ALLOW_OBCS  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
1090  C--     Calculate future values on open boundaries          IF ( calcGMRedi ) THEN
         IF (useOBCS) THEN  
1091  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
1092            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
      &     CALL DEBUG_CALL('OBCS_CALC',myThid)  
1093  #endif  #endif
1094            CALL OBCS_CALC( bi, bj, myTime+deltaTclock, myIter+1,            CALL GMREDI_CALC_TENSOR(
1095       I            uVel, vVel, wVel, theta, salt,       I             iMin, iMax, jMin, jMax,
1096       I            myThid )       I             sigmaX, sigmaY, sigmaR,
1097         I             bi, bj, myTime, myIter, myThid )
1098    #if (defined ALLOW_AUTODIFF) && !(defined ALLOW_OFFLINE)
1099            ELSE
1100              CALL GMREDI_CALC_TENSOR_DUMMY(
1101         I             iMin, iMax, jMin, jMax,
1102         I             sigmaX, sigmaY, sigmaR,
1103         I             bi, bj, myTime, myIter, myThid )
1104    #endif /* ALLOW_AUTODIFF and not ALLOW_OFFLINE */
1105          ENDIF          ENDIF
1106  #endif  /* ALLOW_OBCS */  #endif /* ALLOW_GMREDI */
1107    
1108    #ifdef ALLOW_DOWN_SLOPE
1109            IF ( useDOWN_SLOPE ) THEN
1110    C--     Calculate Downsloping Flow for Down_Slope parameterization
1111             IF ( usingPCoords ) THEN
1112              CALL DWNSLP_CALC_FLOW(
1113         I                bi, bj, kSurfC, rhoInSitu,
1114         I                myTime, myIter, myThid )
1115             ELSE
1116              CALL DWNSLP_CALC_FLOW(
1117         I                bi, bj, kLowC, rhoInSitu,
1118         I                myTime, myIter, myThid )
1119             ENDIF
1120            ENDIF
1121    #endif /* ALLOW_DOWN_SLOPE */
1122    
1123  C--   end bi,bj loops.  C--   end bi,bj loops.
1124         ENDDO         ENDDO
1125        ENDDO        ENDDO
1126    
1127    #ifndef ALLOW_AUTODIFF
1128    C---  if fluid Is Water: end
1129          ENDIF
1130    #endif
1131    
1132    #ifdef ALLOW_OFFLINE
1133          IF ( useOffLine ) THEN
1134    #ifdef ALLOW_DEBUG
1135            IF (debugMode) CALL DEBUG_CALL('OFFLINE_GET_DIFFUS',myThid)
1136    #endif /* ALLOW_DEBUG */
1137            CALL OFFLINE_GET_DIFFUS( myTime, myIter, myThid )
1138          ENDIF
1139    #endif /* ALLOW_OFFLINE */
1140    
1141    #ifdef ALLOW_BBL
1142          IF ( useBBL ) THEN
1143           CALL BBL_CALC_RHS(
1144         I                          myTime, myIter, myThid )
1145          ENDIF
1146    #endif /* ALLOW_BBL */
1147    
1148    #ifdef ALLOW_MYPACKAGE
1149          IF ( useMYPACKAGE ) THEN
1150           CALL MYPACKAGE_CALC_RHS(
1151         I                          myTime, myIter, myThid )
1152          ENDIF
1153    #endif /* ALLOW_MYPACKAGE */
1154    
1155    #ifdef ALLOW_GMREDI
1156          IF ( calcGMRedi ) THEN
1157            CALL GMREDI_DO_EXCH( myTime, myIter, myThid )
1158          ENDIF
1159    #endif /* ALLOW_GMREDI */
1160    
1161    #ifdef ALLOW_KPP
1162          IF ( calcKPP ) THEN
1163            CALL KPP_DO_EXCH( myThid )
1164          ENDIF
1165    #endif /* ALLOW_KPP */
1166    
1167  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
1168        IF ( fluidIsWater .AND. useDiagnostics ) THEN        IF ( fluidIsWater .AND. useDiagnostics ) THEN
1169            CALL DIAGS_RHO_G(
1170         I                    rhoInSitu, uVel, vVel, wVel,
1171         I                    myTime, myIter, myThid )
1172          ENDIF
1173          IF ( useDiagnostics ) THEN
1174          CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )          CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
1175        ENDIF        ENDIF
1176        IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN        IF ( calcConvect .AND. useDiagnostics ) THEN
1177          CALL DIAGNOSTICS_FILL( IVDConvCount,'CONVADJ ',          CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
1178       &                         0, Nr, 0, 1, 1, myThid )       &                               0, Nr, 0, 1, 1, myThid )
1179        ENDIF        ENDIF
1180    #ifdef ALLOW_SALT_PLUME
1181          IF ( useDiagnostics )
1182         &      CALL SALT_PLUME_DIAGNOSTICS_FILL(bi,bj,myThid)
1183    #endif
1184    #endif
1185    
1186    #ifdef ALLOW_ECCO
1187          IF ( useECCO ) CALL ECCO_PHYS( myThid )
1188  #endif  #endif
1189    
1190  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
1191        IF ( debugLevel .GE. debLevB )        IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
      &     CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)  
1192  #endif  #endif
1193    
1194        RETURN        RETURN
Legend:
Removed from v.1.42  
changed lines
  Added in v.1.148
  ViewVC Help
Powered by ViewVC 1.1.22