/[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.41 by heimbach, Wed Apr 18 19:54:06 2007 UTC revision 1.131 by jmc, Sun Jul 7 22:27:30 2013 UTC
# Line 14  C $Name$ Line 14  C $Name$
14  # ifdef ALLOW_SEAICE  # ifdef ALLOW_SEAICE
15  #  include "SEAICE_OPTIONS.h"  #  include "SEAICE_OPTIONS.h"
16  # endif  # endif
17    # ifdef ALLOW_EXF
18    #  include "EXF_OPTIONS.h"
19    # endif
20  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
21    
22  CBOP  CBOP
# Line 36  C     == Global variables === Line 39  C     == Global variables ===
39  #include "SIZE.h"  #include "SIZE.h"
40  #include "EEPARAMS.h"  #include "EEPARAMS.h"
41  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
42  #include "GRID.h"  #include "GRID.h"
43    #include "DYNVARS.h"
44  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
45  #include "TIMEAVE_STATV.h"  # include "TIMEAVE_STATV.h"
46  #endif  #endif
47  #if defined (ALLOW_BALANCE_FLUXES) && !(defined ALLOW_AUTODIFF_TAMC)  #ifdef ALLOW_OFFLINE
48  #include "FFIELDS.h"  # include "OFFLINE_SWITCH.h"
49  #endif  #endif
50    
51  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
52    # include "AUTODIFF_MYFIELDS.h"
53  # include "tamc.h"  # include "tamc.h"
54  # include "tamc_keys.h"  # include "tamc_keys.h"
55  # include "FFIELDS.h"  # include "FFIELDS.h"
56    # include "SURFACE.h"
57  # include "EOS.h"  # include "EOS.h"
58  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
59  #  include "KPP.h"  #  include "KPP.h"
60  # endif  # endif
61    # ifdef ALLOW_GGL90
62    #  include "GGL90.h"
63    # endif
64  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
65  #  include "GMREDI.h"  #  include "GMREDI.h"
66  # endif  # endif
67  # ifdef ALLOW_EBM  # ifdef ALLOW_EBM
68  #  include "EBM.h"  #  include "EBM.h"
69  # endif  # endif
 # ifdef EXACT_CONSERV  
 #  include "SURFACE.h"  
 # endif  
70  # ifdef ALLOW_EXF  # ifdef ALLOW_EXF
71  #  include "ctrl.h"  #  include "ctrl.h"
72  #  include "EXF_FIELDS.h"  #  include "EXF_FIELDS.h"
# Line 70  C     == Global variables === Line 75  C     == Global variables ===
75  #  endif  #  endif
76  # endif  # endif
77  # ifdef ALLOW_SEAICE  # ifdef ALLOW_SEAICE
78    #  include "SEAICE_SIZE.h"
79  #  include "SEAICE.h"  #  include "SEAICE.h"
80    #  include "SEAICE_PARAMS.h"
81    # endif
82    # ifdef ALLOW_THSICE
83    #  include "THSICE_VARS.h"
84    # endif
85    # ifdef ALLOW_SALT_PLUME
86    #  include "SALT_PLUME.h"
87  # endif  # endif
88  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
89    
# Line 85  C     myThid :: Thread number for this i Line 98  C     myThid :: Thread number for this i
98    
99  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
100  C     == Local variables  C     == Local variables
101  C     rhoK, rhoKM1  :: Density at current level, and level above  C     rhoK, rhoKm1  :: Density at current level, and level above
102  C     iMin, iMax    :: Ranges and sub-block indices on which calculations  C     iMin, iMax    :: Ranges and sub-block indices on which calculations
103  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
104  C     bi, bj        :: tile indices  C     bi, bj        :: tile indices
105  C     i,j,k         :: loop indices  C     i,j,k         :: loop indices
106        _RL rhokp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhoKp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107        _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)  
108        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
109        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
110        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 101  C     i,j,k         :: loop indices Line 113  C     i,j,k         :: loop indices
113        INTEGER bi, bj        INTEGER bi, bj
114        INTEGER i, j, k        INTEGER i, j, k
115        INTEGER doDiagsRho        INTEGER doDiagsRho
116          LOGICAL calcGMRedi, calcKPP, calcConvect
117  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
118        LOGICAL  DIAGNOSTICS_IS_ON        LOGICAL  DIAGNOSTICS_IS_ON
119        EXTERNAL DIAGNOSTICS_IS_ON        EXTERNAL DIAGNOSTICS_IS_ON
# Line 114  C--   dummy statement to end declaration Line 127  C--   dummy statement to end declaration
127  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
128    
129  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
130        IF ( debugLevel .GE. debLevB )        IF (debugMode) CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
      &     CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)  
131  #endif  #endif
132    
133        doDiagsRho = 0        doDiagsRho = 0
134  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
135        IF ( useDiagnostics .AND. fluidIsWater ) THEN        IF ( useDiagnostics .AND. fluidIsWater ) THEN
136          IF ( DIAGNOSTICS_IS_ON('DRHODR  ',myThid) ) doDiagsRho = 1          IF ( DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) )
137          IF ( DIAGNOSTICS_IS_ON('RHOANOSQ',myThid) .OR.       &       doDiagsRho = doDiagsRho + 1
138       &       DIAGNOSTICS_IS_ON('URHOMASS',myThid) .OR.          IF ( DIAGNOSTICS_IS_ON('DRHODR  ',myThid) )
139       &       DIAGNOSTICS_IS_ON('VRHOMASS',myThid) .OR.       &       doDiagsRho = doDiagsRho + 2
140       &       DIAGNOSTICS_IS_ON('WRHOMASS',myThid) .OR.          IF ( DIAGNOSTICS_IS_ON('WdRHO_P ',myThid) )
141       &       DIAGNOSTICS_IS_ON('WRHOMASS',myThid) ) doDiagsRho = 2       &       doDiagsRho = doDiagsRho + 4
142            IF ( DIAGNOSTICS_IS_ON('WdRHOdP ',myThid) )
143         &       doDiagsRho = doDiagsRho + 8
144        ENDIF        ENDIF
145  #endif /* ALLOW_DIAGNOSTICS */  #endif /* ALLOW_DIAGNOSTICS */
146    
147          calcGMRedi  = useGMRedi
148          calcKPP     = useKPP
149          calcConvect = ivdc_kappa.NE.0.
150    #ifdef ALLOW_OFFLINE
151          IF ( useOffLine ) THEN
152            calcGMRedi = useGMRedi .AND. .NOT.offlineLoadGMRedi
153            calcKPP    = useKPP    .AND. .NOT.offlineLoadKPP
154            calcConvect=calcConvect.AND. .NOT.offlineLoadConvec
155          ENDIF
156    #endif /* ALLOW_OFFLINE */
157    
158    #ifdef  ALLOW_OBCS
159          IF (useOBCS) THEN
160    C--   Calculate future values on open boundaries
161    C--   moved before SEAICE_MODEL call since SEAICE_MODEL needs seaice-obcs fields
162    # ifdef ALLOW_AUTODIFF_TAMC
163    CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
164    CADJ STORE salt  = comlev1, key=ikey_dynamics, kind=isbyte
165    # endif
166    # ifdef ALLOW_DEBUG
167           IF (debugMode) CALL DEBUG_CALL('OBCS_CALC',myThid)
168    # endif
169           CALL OBCS_CALC( myTime+deltaTClock, myIter+1,
170         I                 uVel, vVel, wVel, theta, salt, myThid )
171          ENDIF
172    #endif  /* ALLOW_OBCS */
173    
174    #ifdef ALLOW_AUTODIFF_TAMC
175    # ifdef ALLOW_SALT_PLUME
176          DO bj=myByLo(myThid),myByHi(myThid)
177           DO bi=myBxLo(myThid),myBxHi(myThid)
178            DO j=1-OLy,sNy+OLy
179             DO i=1-OLx,sNx+OLx
180              saltPlumeDepth(i,j,bi,bj) = 0. _d 0
181              saltPlumeFlux(i,j,bi,bj)  = 0. _d 0
182             ENDDO
183            ENDDO
184           ENDDO
185          ENDDO
186    # endif
187    #endif /* ALLOW_AUTODIFF_TAMC */
188    
189    #ifdef ALLOW_FRAZIL
190          IF ( useFRAZIL ) THEN
191    C--   Freeze water in the ocean interior and let it rise to the surface
192           CALL FRAZIL_CALC_RHS( myTime, myIter, myThid )
193          ENDIF
194    #endif /* ALLOW_FRAZIL */
195    
196    #ifndef OLD_THSICE_CALL_SEQUENCE
197    #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
198          IF ( useThSIce .AND. fluidIsWater ) THEN
199    # ifdef ALLOW_AUTODIFF_TAMC
200    CADJ STORE uice,vice         = comlev1, key = ikey_dynamics,
201    CADJ &     kind = isbyte
202    CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
203    CADJ &     kind = isbyte
204    CADJ STORE snowHeight, Tsrf  = comlev1, key = ikey_dynamics,
205    CADJ &     kind = isbyte
206    CADJ STORE Qice1, Qice2      = comlev1, key = ikey_dynamics,
207    CADJ &     kind = isbyte
208    CADJ STORE sHeating, snowAge = comlev1, key = ikey_dynamics,
209    CADJ &     kind = isbyte
210    CADJ STORE hocemxl = comlev1, key = ikey_dynamics,
211    CADJ &     kind = isbyte
212    CADJ STORE icflxsw = comlev1, key = ikey_dynamics,
213    CADJ &     kind = isbyte
214    CADJ STORE salt,theta        = comlev1, key = ikey_dynamics,
215    CADJ &     kind = isbyte
216    CADJ STORE uvel,vvel         = comlev1, key = ikey_dynamics,
217    CADJ &     kind = isbyte
218    CADJ STORE qnet,qsw, empmr   = comlev1, key = ikey_dynamics,
219    CADJ &     kind = isbyte
220    CADJ STORE atemp,aqh,precip  = comlev1, key = ikey_dynamics,
221    CADJ &     kind = isbyte
222    CADJ STORE swdown,lwdown     = comlev1, key = ikey_dynamics,
223    CADJ &     kind = isbyte
224    #  ifdef NONLIN_FRSURF
225    CADJ STORE hFac_surfC       = comlev1, key = ikey_dynamics,
226    CADJ &     kind = isbyte
227    #  endif
228    # endif /* ALLOW_AUTODIFF_TAMC */
229    # ifdef ALLOW_DEBUG
230            IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
231    # endif
232    C--     Step forward Therm.Sea-Ice variables
233    C       and modify forcing terms including effects from ice
234            CALL TIMER_START('THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
235            CALL THSICE_MAIN( myTime, myIter, myThid )
236            CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
237          ENDIF
238    #endif /* ALLOW_THSICE */
239    #endif /* ndef OLD_THSICE_CALL_SEQUENCE */
240    
241    #ifdef ALLOW_SEAICE
242    # ifdef ALLOW_AUTODIFF_TAMC
243    CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
244    CADJ STORE fu,fv  = comlev1, key=ikey_dynamics, kind=isbyte
245    CADJ STORE qnet   = comlev1, key=ikey_dynamics, kind=isbyte
246    CADJ STORE qsw    = comlev1, key=ikey_dynamics, kind=isbyte
247    CADJ STORE theta  = comlev1, key=ikey_dynamics, kind=isbyte
248    CADJ STORE salt   = comlev1, key=ikey_dynamics, kind=isbyte
249    #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
250    CADJ STORE evap   = comlev1, key=ikey_dynamics, kind=isbyte
251    #endif
252          IF ( .NOT.useSEAICE ) THEN
253          IF ( SEAICEadjMODE .EQ. -1 ) THEN
254            CALL SEAICE_FAKE( myTime, myIter, myThid )
255          ENDIF
256          ENDIF
257    CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
258    CADJ STORE fu,fv  = comlev1, key=ikey_dynamics, kind=isbyte
259    CADJ STORE qnet   = comlev1, key=ikey_dynamics, kind=isbyte
260    CADJ STORE qsw    = comlev1, key=ikey_dynamics, kind=isbyte
261    CADJ STORE theta  = comlev1, key=ikey_dynamics, kind=isbyte
262    CADJ STORE salt   = comlev1, key=ikey_dynamics, kind=isbyte
263    #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
264    CADJ STORE evap   = comlev1, key=ikey_dynamics, kind=isbyte
265    #endif
266    # endif /* ALLOW_AUTODIFF_TAMC */
267    #endif /* ALLOW_SEAICE */
268    
269  #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.  
270        IF ( useSEAICE ) THEN        IF ( useSEAICE ) THEN
271  #ifdef ALLOW_AUTODIFF_TAMC  # ifdef ALLOW_AUTODIFF_TAMC
272  CADJ STORE atemp,aqh,precip    = comlev1, key = ikey_dynamics  cph-adj-test(
273  CADJ STORE swdown,lwdown       = comlev1, key = ikey_dynamics  CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
274  CADJ STORE theta               = comlev1, key = ikey_dynamics  CADJ STORE hsnow  = comlev1, key=ikey_dynamics, kind=isbyte
275    CADJ STORE heff   = comlev1, key=ikey_dynamics, kind=isbyte
276    CADJ STORE empmr,qsw,theta   = comlev1, key = ikey_dynamics,
277    CADJ &     kind = isbyte
278    cph-adj-test)
279    CADJ STORE atemp,aqh,precip    = comlev1, key = ikey_dynamics,
280    CADJ &     kind = isbyte
281    CADJ STORE swdown,lwdown       = comlev1, key = ikey_dynamics,
282    CADJ &     kind = isbyte
283  cph# ifdef EXF_READ_EVAP  cph# ifdef EXF_READ_EVAP
284  CADJ STORE evap                = comlev1, key = ikey_dynamics  CADJ STORE evap                = comlev1, key = ikey_dynamics,
285  cph# endif  CADJ &     kind = isbyte
 cph# ifdef SEAICE_ALLOW_DYNAMICS  
 CADJ STORE uvel,vvel           = comlev1, key = ikey_dynamics  
286  cph# endif  cph# endif
287  # ifdef SEAICE_CGRID  CADJ STORE uvel,vvel           = comlev1, key = ikey_dynamics,
288  CADJ STORE fu, fv              = comlev1, key = ikey_dynamics  CADJ &     kind = isbyte
289  CADJ STORE seaicemasku         = comlev1, key = ikey_dynamics  #  ifdef SEAICE_CGRID
290  CADJ STORE seaicemaskv         = comlev1, key = ikey_dynamics  CADJ STORE stressdivergencex   = comlev1, key = ikey_dynamics,
291    CADJ &     kind = isbyte
292    CADJ STORE stressdivergencey   = comlev1, key = ikey_dynamics,
293    CADJ &     kind = isbyte
294    #  endif
295    #  ifdef SEAICE_ALLOW_DYNAMICS
296    CADJ STORE uice                = comlev1, key = ikey_dynamics,
297    CADJ &     kind = isbyte
298    CADJ STORE vice                = comlev1, key = ikey_dynamics,
299    CADJ &     kind = isbyte
300    #   ifdef SEAICE_ALLOW_EVP
301    CADJ STORE seaice_sigma1       = comlev1, key = ikey_dynamics,
302    CADJ &     kind = isbyte
303    CADJ STORE seaice_sigma2       = comlev1, key = ikey_dynamics,
304    CADJ &     kind = isbyte
305    CADJ STORE seaice_sigma12      = comlev1, key = ikey_dynamics,
306    CADJ &     kind = isbyte
307    #   endif
308    #  endif
309    cph#  ifdef SEAICE_SALINITY
310    CADJ STORE salt                = comlev1, key = ikey_dynamics,
311    CADJ &     kind = isbyte
312    cph#  endif
313  #  ifdef ATMOSPHERIC_LOADING  #  ifdef ATMOSPHERIC_LOADING
314  CADJ STORE siceload            = comlev1, key = ikey_dynamics  CADJ STORE pload               = comlev1, key = ikey_dynamics,
315    CADJ &     kind = isbyte
316    CADJ STORE siceload            = comlev1, key = ikey_dynamics,
317    CADJ &     kind = isbyte
318    #  endif
319    #  ifdef NONLIN_FRSURF
320    CADJ STORE recip_hfacc         = comlev1, key = ikey_dynamics,
321    CADJ &     kind = isbyte
322    #  endif
323    #  ifdef ANNUAL_BALANCE
324    CADJ STORE balance_itcount     = comlev1, key = ikey_dynamics,
325    CADJ &     kind = isbyte
326    #  endif /* ANNUAL_BALANCE */
327    #  ifdef ALLOW_THSICE
328    CADJ STORE fu, fv              = comlev1, key = ikey_dynamics,
329    CADJ &     kind = isbyte
330  #  endif  #  endif
331    # endif /* ALLOW_AUTODIFF_TAMC */
332    # ifdef ALLOW_DEBUG
333            IF (debugMode) CALL DEBUG_CALL('SEAICE_MODEL',myThid)
334  # endif  # endif
 #endif  
 #ifdef ALLOW_DEBUG  
         IF ( debugLevel .GE. debLevB )  
      &    CALL DEBUG_CALL('SEAICE_MODEL',myThid)  
 #endif  
335          CALL TIMER_START('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)          CALL TIMER_START('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
336          CALL SEAICE_MODEL( myTime, myIter, myThid )          CALL SEAICE_MODEL( myTime, myIter, myThid )
337          CALL TIMER_STOP ('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)          CALL TIMER_STOP ('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
338  #ifdef ALLOW_COST_ICE  # ifdef ALLOW_COST
339          CALL COST_ICE_TEST ( myTime, myIter, myThid )          CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
340  #endif  # endif
341        ENDIF        ENDIF
342  #endif /* ALLOW_SEAICE */  #endif /* ALLOW_SEAICE */
343    
344    #ifdef ALLOW_AUTODIFF_TAMC
345    CADJ STORE sst, sss           = comlev1, key = ikey_dynamics,
346    CADJ &     kind = isbyte
347    CADJ STORE qsw                = comlev1, key = ikey_dynamics,
348    CADJ &     kind = isbyte
349    # ifdef ALLOW_SEAICE
350    CADJ STORE area               = comlev1, key = ikey_dynamics,
351    CADJ &     kind = isbyte
352    # endif
353    #endif
354    
355    #ifdef OLD_THSICE_CALL_SEQUENCE
356  #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)  #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
357        IF ( useThSIce .AND. fluidIsWater ) THEN        IF ( useThSIce .AND. fluidIsWater ) THEN
358  #ifdef ALLOW_DEBUG  # ifdef ALLOW_AUTODIFF_TAMC
359          IF ( debugLevel .GE. debLevB )  cph(
360       &    CALL DEBUG_CALL('THSICE_MAIN',myThid)  #  ifdef NONLIN_FRSURF
361  #endif  CADJ STORE uice,vice        = comlev1, key = ikey_dynamics,
362    CADJ &     kind = isbyte
363    CADJ STORE salt,theta       = comlev1, key = ikey_dynamics,
364    CADJ &     kind = isbyte
365    CADJ STORE qnet,qsw, empmr  = comlev1, key = ikey_dynamics,
366    CADJ &     kind = isbyte
367    CADJ STORE hFac_surfC       = comlev1, key = ikey_dynamics,
368    CADJ &     kind = isbyte
369    #  endif
370    # endif
371    # ifdef ALLOW_DEBUG
372            IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
373    # endif
374  C--     Step forward Therm.Sea-Ice variables  C--     Step forward Therm.Sea-Ice variables
375  C       and modify forcing terms including effects from ice  C       and modify forcing terms including effects from ice
376          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 378  C       and modify forcing terms includi
378          CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)          CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
379        ENDIF        ENDIF
380  #endif /* ALLOW_THSICE */  #endif /* ALLOW_THSICE */
381    #endif /* OLD_THSICE_CALL_SEQUENCE */
382    
383  #ifdef ALLOW_SHELFICE  #ifdef ALLOW_SHELFICE
384    # ifdef ALLOW_AUTODIFF_TAMC
385    CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
386    CADJ &     kind = isbyte
387    # endif
388        IF ( useShelfIce .AND. fluidIsWater ) THEN        IF ( useShelfIce .AND. fluidIsWater ) THEN
389  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
390          IF ( debugLevel .GE. debLevB )         IF (debugMode) CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
      &    CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)  
391  #endif  #endif
392  C     compute temperature and (virtual) salt flux at the  C     compute temperature and (virtual) salt flux at the
393  C     shelf-ice ocean interface  C     shelf-ice ocean interface
394         CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',         CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
395       &       myThid)       &       myThid)
# Line 204  C     shelf-ice ocean interface Line 399  C     shelf-ice ocean interface
399        ENDIF        ENDIF
400  #endif /* ALLOW_SHELFICE */  #endif /* ALLOW_SHELFICE */
401    
402    #ifdef ALLOW_ICEFRONT
403          IF ( useICEFRONT .AND. fluidIsWater ) THEN
404    #ifdef ALLOW_DEBUG
405           IF (debugMode) CALL DEBUG_CALL('ICEFRONT_THERMODYNAMICS',myThid)
406    #endif
407    C     compute temperature and (virtual) salt flux at the
408    C     ice-front ocean interface
409           CALL TIMER_START('ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
410         &       myThid)
411           CALL ICEFRONT_THERMODYNAMICS( myTime, myIter, myThid )
412           CALL TIMER_STOP( 'ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
413         &      myThid)
414          ENDIF
415    #endif /* ALLOW_ICEFRONT */
416    
417    #ifdef ALLOW_SALT_PLUME
418          IF ( useSALT_PLUME ) THEN
419              CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
420          ENDIF
421    #endif /* ALLOW_SALT_PLUME */
422    
423  C--   Freeze water at the surface  C--   Freeze water at the surface
424          IF ( allowFreezing ) THEN
425  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
426  CADJ STORE theta = comlev1, key = ikey_dynamics  CADJ STORE theta = comlev1, key = ikey_dynamics,
427    CADJ &     kind = isbyte
428  #endif  #endif
429        IF ( allowFreezing          CALL FREEZE_SURFACE( myTime, myIter, myThid )
      &                   .AND. .NOT. useSEAICE  
      &                   .AND. .NOT. useThSIce ) THEN  
         CALL FREEZE_SURFACE(  myTime, myIter, myThid )  
430        ENDIF        ENDIF
431    
432  #ifdef ALLOW_OCN_COMPON_INTERF  #ifdef ALLOW_OCN_COMPON_INTERF
# Line 222  C jmc: do not know precisely where to pu Line 437  C jmc: do not know precisely where to pu
437        ENDIF        ENDIF
438  #endif /* ALLOW_OCN_COMPON_INTERF */  #endif /* ALLOW_OCN_COMPON_INTERF */
439    
440  #ifdef ALLOW_BALANCE_FLUXES        iMin = 1-OLx
441  C     balance fluxes        iMax = sNx+OLx
442        IF ( balanceEmPmR )        jMin = 1-OLy
443       &        CALL REMOVE_MEAN_RS( 1, EmPmR, maskH, maskH, rA, drF,        jMax = sNy+OLy
444       &        'EmPmR', myTime, myThid )  
445        IF ( balanceQnet )  C---  Determines forcing terms based on external fields
446       &        CALL REMOVE_MEAN_RS( 1, Qnet,  maskH, maskH, rA, drF,  C     relaxation terms, etc.
447       &        'Qnet ', myTime, myThid )  #ifdef ALLOW_DEBUG
448  #endif /* ALLOW_BALANCE_FLUXES */        IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
449    #endif
450    #ifdef ALLOW_AUTODIFF_TAMC
451    CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
452    CADJ &     kind = isbyte
453    #else  /* ALLOW_AUTODIFF_TAMC */
454    C--   if fluid is not water, by-pass surfaceForcing, find_rho, gmredi
455    C     and all vertical mixing schemes, but keep OBCS_CALC
456          IF ( fluidIsWater ) THEN
457    #endif /* ALLOW_AUTODIFF_TAMC */
458            CALL EXTERNAL_FORCING_SURF(
459         I             iMin, iMax, jMin, jMax,
460         I             myTime, myIter, myThid )
461    
462  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
463  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 254  CHPF$ INDEPENDENT Line 481  CHPF$ INDEPENDENT
481            itdkey = (act1 + 1) + act2*max1            itdkey = (act1 + 1) + act2*max1
482       &                      + act3*max1*max2       &                      + act3*max1*max2
483       &                      + 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  
484  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
485    
486  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 488  C     These inital values do not alter t
488  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
489  C     point numbers. This prevents spurious hardware signals due to  C     point numbers. This prevents spurious hardware signals due to
490  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  
   
491          DO k=1,Nr          DO k=1,Nr
492           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
493            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
494  C This is currently also used by IVDC and Diagnostics  C This is currently used by GMRedi, IVDC, MXL-depth  and Diagnostics
495             sigmaX(i,j,k) = 0. _d 0             sigmaX(i,j,k) = 0. _d 0
496             sigmaY(i,j,k) = 0. _d 0             sigmaY(i,j,k) = 0. _d 0
497             sigmaR(i,j,k) = 0. _d 0             sigmaR(i,j,k) = 0. _d 0
498              ENDDO
499             ENDDO
500            ENDDO
501    
502  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
503            DO j=1-OLy,sNy+OLy
504             DO i=1-OLx,sNx+OLx
505              rhoKm1 (i,j)   = 0. _d 0
506              rhoKp1 (i,j)   = 0. _d 0
507             ENDDO
508            ENDDO
509  cph all the following init. are necessary for TAF  cph all the following init. are necessary for TAF
510  cph although some of these are re-initialised later.  cph although some of these are re-initialised later.
511            DO k=1,Nr
512             DO j=1-OLy,sNy+OLy
513              DO i=1-OLx,sNx+OLx
514               rhoInSitu(i,j,k,bi,bj) = 0.
515    # ifdef ALLOW_GGL90
516               GGL90viscArU(i,j,k,bi,bj)  = 0. _d 0
517               GGL90viscArV(i,j,k,bi,bj)  = 0. _d 0
518               GGL90diffKr(i,j,k,bi,bj)  = 0. _d 0
519    # endif /* ALLOW_GGL90 */
520              ENDDO
521             ENDDO
522            ENDDO
523    #ifdef ALLOW_OFFLINE
524           IF ( calcConvect ) THEN
525    #endif
526            DO k=1,Nr
527             DO j=1-OLy,sNy+OLy
528              DO i=1-OLx,sNx+OLx
529             IVDConvCount(i,j,k,bi,bj) = 0.             IVDConvCount(i,j,k,bi,bj) = 0.
530              ENDDO
531             ENDDO
532            ENDDO
533    #ifdef ALLOW_OFFLINE
534           ENDIF
535           IF ( calcGMRedi ) THEN
536    #endif
537  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
538            DO k=1,Nr
539             DO j=1-OLy,sNy+OLy
540              DO i=1-OLx,sNx+OLx
541             Kwx(i,j,k,bi,bj)  = 0. _d 0             Kwx(i,j,k,bi,bj)  = 0. _d 0
542             Kwy(i,j,k,bi,bj)  = 0. _d 0             Kwy(i,j,k,bi,bj)  = 0. _d 0
543             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 556  cph although some of these are re-initia
556  #  ifdef GM_VISBECK_VARIABLE_K  #  ifdef GM_VISBECK_VARIABLE_K
557             VisbeckK(i,j,bi,bj)   = 0. _d 0             VisbeckK(i,j,bi,bj)   = 0. _d 0
558  #  endif  #  endif
559              ENDDO
560             ENDDO
561            ENDDO
562  # endif /* ALLOW_GMREDI */  # endif /* ALLOW_GMREDI */
563  #endif /* ALLOW_AUTODIFF_TAMC */  #ifdef ALLOW_OFFLINE
564           ENDIF
565           IF ( calcKPP ) THEN
566    #endif
567    # ifdef ALLOW_KPP
568            DO k=1,Nr
569             DO j=1-OLy,sNy+OLy
570              DO i=1-OLx,sNx+OLx
571               KPPdiffKzS(i,j,k,bi,bj)  = 0. _d 0
572               KPPdiffKzT(i,j,k,bi,bj)  = 0. _d 0
573            ENDDO            ENDDO
574           ENDDO           ENDDO
575          ENDDO          ENDDO
576    # endif /* ALLOW_KPP */
577          iMin = 1-OLx  #ifdef ALLOW_OFFLINE
578          iMax = sNx+OLx         ENDIF
579          jMin = 1-OLy  #endif
580          jMax = sNy+OLy  #endif /* ALLOW_AUTODIFF_TAMC */
581    
582  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
583  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
584  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     kind = isbyte
585    CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
586    CADJ &     kind = isbyte
587  CADJ STORE totphihyd(:,:,:,bi,bj)  CADJ STORE totphihyd(:,:,:,bi,bj)
588  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
589    CADJ &     kind = isbyte
590  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
591  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
592  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     kind = isbyte
593    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
594    CADJ &     kind = isbyte
595    # endif
596    # ifdef ALLOW_SALT_PLUME
597    CADJ STORE saltplumedepth(:,:,bi,bj) = comlev1_bibj, key=itdkey,
598    CADJ &     kind = isbyte
599  # endif  # endif
600  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
601    
602    C--   Always compute density (stored in common block) here; even when it is not
603    C     needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
604  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
605          IF ( debugLevel .GE. debLevB )          IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D (xNr)',myThid)
606       &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)  #endif
607    #ifdef ALLOW_AUTODIFF_TAMC
608            IF ( fluidIsWater ) THEN
609    #endif /* ALLOW_AUTODIFF_TAMC */
610    #ifdef ALLOW_DOWN_SLOPE
611             IF ( useDOWN_SLOPE ) THEN
612               DO k=1,Nr
613                CALL DWNSLP_CALC_RHO(
614         I                  theta, salt,
615         O                  rhoInSitu(1-OLx,1-OLy,k,bi,bj),
616         I                  k, bi, bj, myTime, myIter, myThid )
617               ENDDO
618             ENDIF
619    #endif /* ALLOW_DOWN_SLOPE */
620    #ifdef ALLOW_BBL
621             IF ( useBBL ) THEN
622    C     pkg/bbl requires in-situ bbl density for depths equal to and deeper than the bbl.
623    C     To reduce computation and storage requirement, these densities are stored in the
624    C     dry grid boxes of rhoInSitu.  See BBL_CALC_RHO for details.
625               DO k=Nr,1,-1
626                CALL BBL_CALC_RHO(
627         I                  theta, salt,
628         O                  rhoInSitu,
629         I                  k, bi, bj, myTime, myIter, myThid )
630    
631               ENDDO
632             ENDIF
633    #endif /* ALLOW_BBL */
634             IF ( .NOT. ( useDOWN_SLOPE .OR. useBBL ) ) THEN
635               DO k=1,Nr
636                CALL FIND_RHO_2D(
637         I                iMin, iMax, jMin, jMax, k,
638         I                theta(1-OLx,1-OLy,k,bi,bj),
639         I                salt (1-OLx,1-OLy,k,bi,bj),
640         O                rhoInSitu(1-OLx,1-OLy,k,bi,bj),
641         I                k, bi, bj, myThid )
642               ENDDO
643             ENDIF
644    #ifdef ALLOW_AUTODIFF_TAMC
645            ELSE
646    C-        fluid is not water:
647              DO k=1,Nr
648               DO j=1-OLy,sNy+OLy
649                DO i=1-OLx,sNx+OLx
650                  rhoInSitu(i,j,k,bi,bj) = 0.
651                ENDDO
652               ENDDO
653              ENDDO
654            ENDIF
655    #endif /* ALLOW_AUTODIFF_TAMC */
656    
657    #ifdef ALLOW_DEBUG
658            IF (debugMode) CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
659  #endif  #endif
660    
661  C--     Start of diagnostic loop  C--     Start of diagnostic loop
# Line 337  C--     Start of diagnostic loop Line 664  C--     Start of diagnostic loop
664  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
665  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?
666  C? Do we still need this?  C? Do we still need this?
667  cph kkey formula corrected.  cph kkey formula corrected.
668  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhoK, rhoKm1, in the case useGMREDI.
669           kkey = (itdkey-1)*Nr + k            kkey = (itdkey-1)*Nr + k
670  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
671    
672    c#ifdef ALLOW_AUTODIFF_TAMC
673    cCADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
674    cCADJ &     kind = isbyte
675    cCADJ STORE salt(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey,
676    cCADJ &     kind = isbyte
677    c#endif /* ALLOW_AUTODIFF_TAMC */
678    
679  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
680  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
681            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)            IF ( calcGMRedi .OR. (k.GT.1 .AND. calcConvect)
682       &                   .OR. doDiagsRho.GE.1 ) THEN       &         .OR. usePP81 .OR. useMY82 .OR. useGGL90
683  #ifdef ALLOW_DEBUG       &         .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
             IF ( debugLevel .GE. debLevB )  
      &       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 )  
   
684              IF (k.GT.1) THEN              IF (k.GT.1) THEN
685  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
686  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,
687  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ &     kind = isbyte
688  #endif /* ALLOW_AUTODIFF_TAMC */  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
689               CALL FIND_RHO(  CADJ &     kind = isbyte
690       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,  CADJ STORE rhokm1 (bi,bj)       = comlev1_bibj_k, key=kkey,
691       I        theta, salt,  CADJ &     kind = isbyte
692       O        rhoKm1,  #endif /* ALLOW_AUTODIFF_TAMC */
693       I        myThid )               CALL FIND_RHO_2D(
694         I                 iMin, iMax, jMin, jMax, k,
695         I                 theta(1-OLx,1-OLy,k-1,bi,bj),
696         I                 salt (1-OLx,1-OLy,k-1,bi,bj),
697         O                 rhoKm1,
698         I                 k-1, bi, bj, myThid )
699              ENDIF              ENDIF
700  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
701              IF ( debugLevel .GE. debLevB )              IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
      &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)  
702  #endif  #endif
703  cph Avoid variable aliasing for adjoint !!!  cph Avoid variable aliasing for adjoint !!!
704              DO j=jMin,jMax              DO j=jMin,jMax
705               DO i=iMin,iMax               DO i=iMin,iMax
706                rhoKP1(i,j) = rhoK(i,j)                rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
707               ENDDO               ENDDO
708              ENDDO              ENDDO
709              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
710       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
711       I             rhoK, rhoKm1, rhoKp1,       I             rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
712       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
713       I             myThid )       I             myThid )
           ENDIF  
   
714  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
715  ctest# ifndef GM_EXCLUDE_CLIPPING  #ifdef GMREDI_WITH_STABLE_ADJOINT
716  CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
717  ctest# endif  cgf -> cuts adjoint dependency from slope to state
718  CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte              CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
719                CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
720                CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
721    #endif
722  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
723              ENDIF
724    
725  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
726  c ==> should use sigmaR !!!            IF (k.GT.1 .AND. calcConvect) THEN
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
727  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
728              IF ( debugLevel .GE. debLevB )              IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid)
      &       CALL DEBUG_CALL('CALC_IVDC',myThid)  
729  #endif  #endif
730              CALL CALC_IVDC(              CALL CALC_IVDC(
731       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
732       I        rhoKm1, rhoK,       I        sigmaR,
733       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
734            ENDIF            ENDIF
735    
736  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
737            IF ( doDiagsRho.GE.2 ) THEN            IF ( doDiagsRho.GE.4 ) THEN
738              CALL DIAGS_RHO( k, bi, bj,              CALL DIAGS_RHO_L( doDiagsRho, k, bi, bj,
739       I                      rhoK, rhoKm1,       I                        rhoInSitu(1-OLx,1-OLy,1,bi,bj),
740       I                      myTime, myIter, myThid)       I                        rhoKm1, wVel,
741         I                        myTime, myIter, myThid )
742            ENDIF            ENDIF
743  #endif  #endif
744    
745  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
746          ENDDO          ENDDO
747    
748    #ifdef ALLOW_AUTODIFF_TAMC
749    CADJ STORE IVDConvCount(:,:,:,bi,bj)
750    CADJ &     = comlev1_bibj, key=itdkey,
751    CADJ &     kind = isbyte
752    #endif
753    
754    C--     Diagnose Mixed Layer Depth:
755            IF ( calcGMRedi .OR. MOD(doDiagsRho,2).EQ.1 ) THEN
756              CALL CALC_OCE_MXLAYER(
757         I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
758         I              bi, bj, myTime, myIter, myThid )
759            ENDIF
760    
761    #ifdef ALLOW_SALT_PLUME
762            IF ( useSALT_PLUME ) THEN
763              CALL SALT_PLUME_CALC_DEPTH(
764         I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
765         I              bi, bj, myTime, myIter, myThid )
766            ENDIF
767    #endif /* ALLOW_SALT_PLUME */
768    
769  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
770          IF ( doDiagsRho.GE.1 ) THEN          IF ( MOD(doDiagsRho,4).GE.2 ) THEN
771            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,
772       &         2, bi, bj, myThid)       &         2, bi, bj, myThid)
773          ENDIF          ENDIF
774  #endif  #endif /* ALLOW_DIAGNOSTICS */
775    
776  C--     Determines forcing terms based on external fields  C--    This is where EXTERNAL_FORCING_SURF(bi,bj) used to be called;
777  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  
778    
779  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
780  cph needed for KPP  cph needed for KPP
781  CADJ STORE surfaceForcingU(:,:,bi,bj)  CADJ STORE surfaceForcingU(:,:,bi,bj)
782  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
783    CADJ &     kind = isbyte
784  CADJ STORE surfaceForcingV(:,:,bi,bj)  CADJ STORE surfaceForcingV(:,:,bi,bj)
785  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
786    CADJ &     kind = isbyte
787  CADJ STORE surfaceForcingS(:,:,bi,bj)  CADJ STORE surfaceForcingS(:,:,bi,bj)
788  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
789    CADJ &     kind = isbyte
790  CADJ STORE surfaceForcingT(:,:,bi,bj)  CADJ STORE surfaceForcingT(:,:,bi,bj)
791  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
792    CADJ &     kind = isbyte
793  CADJ STORE surfaceForcingTice(:,:,bi,bj)  CADJ STORE surfaceForcingTice(:,:,bi,bj)
794  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
795    CADJ &     kind = isbyte
796  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
797    
 #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 */  
   
798  #ifdef  ALLOW_KPP  #ifdef  ALLOW_KPP
799  C--     Compute KPP mixing coefficients  C--     Compute KPP mixing coefficients
800          IF (useKPP) THEN          IF ( calcKPP ) THEN
801  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
802            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid)
      &     CALL DEBUG_CALL('KPP_CALC',myThid)  
803  #endif  #endif
804              CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
805            CALL KPP_CALC(            CALL KPP_CALC(
806       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myIter, myThid )
807  #ifdef ALLOW_AUTODIFF_TAMC            CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
808    #if (defined ALLOW_AUTODIFF_TAMC) && !(defined ALLOW_OFFLINE)
809          ELSE          ELSE
810            CALL KPP_CALC_DUMMY(            CALL KPP_CALC_DUMMY(
811       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myIter, myThid )
812  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC and not ALLOW_OFFLINE */
813          ENDIF          ENDIF
   
814  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
815    
816  #ifdef  ALLOW_PP81  #ifdef  ALLOW_PP81
817  C--     Compute PP81 mixing coefficients  C--     Compute PP81 mixing coefficients
818          IF (usePP81) THEN          IF (usePP81) THEN
819  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
820            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid)
      &     CALL DEBUG_CALL('PP81_CALC',myThid)  
821  #endif  #endif
822            CALL PP81_CALC(            CALL PP81_CALC(
823       I                  bi, bj, myTime, myThid )       I                     bi, bj, sigmaR, myTime, myIter, myThid )
824          ENDIF          ENDIF
825  #endif /* ALLOW_PP81 */  #endif /* ALLOW_PP81 */
826    
# Line 539  C--     Compute PP81 mixing coefficients Line 828  C--     Compute PP81 mixing coefficients
828  C--     Compute MY82 mixing coefficients  C--     Compute MY82 mixing coefficients
829          IF (useMY82) THEN          IF (useMY82) THEN
830  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
831            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid)
      &     CALL DEBUG_CALL('MY82_CALC',myThid)  
832  #endif  #endif
833            CALL MY82_CALC(            CALL MY82_CALC(
834       I                  bi, bj, myTime, myThid )       I                     bi, bj, sigmaR, myTime, myIter, myThid )
835          ENDIF          ENDIF
836  #endif /* ALLOW_MY82 */  #endif /* ALLOW_MY82 */
837    
838  #ifdef  ALLOW_GGL90  #ifdef  ALLOW_GGL90
839    #ifdef ALLOW_AUTODIFF_TAMC
840    CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
841    CADJ &     kind = isbyte
842    #endif /* ALLOW_AUTODIFF_TAMC */
843  C--     Compute GGL90 mixing coefficients  C--     Compute GGL90 mixing coefficients
844          IF (useGGL90) THEN          IF (useGGL90) THEN
845  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
846            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
      &     CALL DEBUG_CALL('GGL90_CALC',myThid)  
847  #endif  #endif
848              CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
849            CALL GGL90_CALC(            CALL GGL90_CALC(
850       I                  bi, bj, myTime, myThid )       I                     bi, bj, sigmaR, myTime, myIter, myThid )
851              CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
852          ENDIF          ENDIF
853  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
854    
# Line 563  C--     Compute GGL90 mixing coefficient Line 856  C--     Compute GGL90 mixing coefficient
856          IF ( taveFreq.GT. 0. _d 0 ) THEN          IF ( taveFreq.GT. 0. _d 0 ) THEN
857            CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)            CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
858          ENDIF          ENDIF
859          IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN          IF ( taveFreq.GT.0. .AND. calcConvect ) THEN
860            CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,            CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
861       I                           Nr, deltaTclock, bi, bj, myThid)       I                           Nr, deltaTClock, bi, bj, myThid)
862          ENDIF          ENDIF
863  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
864    
865  #ifndef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_GMREDI
866  C---  if fluid Is Water: end  #ifdef ALLOW_AUTODIFF_TAMC
867          ENDIF  # ifndef GM_EXCLUDE_CLIPPING
868  #endif  cph storing here is needed only for one GMREDI_OPTIONS:
869    cph define GM_BOLUS_ADVEC
870    cph keep it although TAF says you dont need to.
871    cph but I have avoided the #ifdef for now, in case more things change
872    CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey,
873    CADJ &     kind = isbyte
874    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey,
875    CADJ &     kind = isbyte
876    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey,
877    CADJ &     kind = isbyte
878    # endif
879    #endif /* ALLOW_AUTODIFF_TAMC */
880    
881  #ifdef  ALLOW_OBCS  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
882  C--     Calculate future values on open boundaries          IF ( calcGMRedi ) THEN
         IF (useOBCS) THEN  
883  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
884            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
      &     CALL DEBUG_CALL('OBCS_CALC',myThid)  
885  #endif  #endif
886            CALL OBCS_CALC( bi, bj, myTime+deltaTclock, myIter+1,            CALL GMREDI_CALC_TENSOR(
887       I            uVel, vVel, wVel, theta, salt,       I             iMin, iMax, jMin, jMax,
888       I            myThid )       I             sigmaX, sigmaY, sigmaR,
889         I             bi, bj, myTime, myIter, myThid )
890    #if (defined ALLOW_AUTODIFF_TAMC) && !(defined ALLOW_OFFLINE)
891            ELSE
892              CALL GMREDI_CALC_TENSOR_DUMMY(
893         I             iMin, iMax, jMin, jMax,
894         I             sigmaX, sigmaY, sigmaR,
895         I             bi, bj, myTime, myIter, myThid )
896    #endif /* ALLOW_AUTODIFF_TAMC and not ALLOW_OFFLINE */
897          ENDIF          ENDIF
898  #endif  /* ALLOW_OBCS */  #endif /* ALLOW_GMREDI */
899    
900    #ifdef ALLOW_DOWN_SLOPE
901            IF ( useDOWN_SLOPE ) THEN
902    C--     Calculate Downsloping Flow for Down_Slope parameterization
903             IF ( usingPCoords ) THEN
904              CALL DWNSLP_CALC_FLOW(
905         I                bi, bj, kSurfC, rhoInSitu,
906         I                myTime, myIter, myThid )
907             ELSE
908              CALL DWNSLP_CALC_FLOW(
909         I                bi, bj, kLowC, rhoInSitu,
910         I                myTime, myIter, myThid )
911             ENDIF
912            ENDIF
913    #endif /* ALLOW_DOWN_SLOPE */
914    
915  C--   end bi,bj loops.  C--   end bi,bj loops.
916         ENDDO         ENDDO
917        ENDDO        ENDDO
918    
919    #ifndef ALLOW_AUTODIFF_TAMC
920    C---  if fluid Is Water: end
921          ENDIF
922    #endif
923    
924    #ifdef ALLOW_BBL
925          IF ( useBBL ) THEN
926           CALL BBL_CALC_RHS(
927         I                          myTime, myIter, myThid )
928          ENDIF
929    #endif /* ALLOW_BBL */
930    
931    #ifdef ALLOW_MYPACKAGE
932          IF ( useMYPACKAGE ) THEN
933           CALL MYPACKAGE_CALC_RHS(
934         I                          myTime, myIter, myThid )
935          ENDIF
936    #endif /* ALLOW_MYPACKAGE */
937    
938    #ifdef ALLOW_GMREDI
939          IF ( calcGMRedi ) THEN
940            CALL GMREDI_DO_EXCH( myTime, myIter, myThid )
941          ENDIF
942    #endif /* ALLOW_GMREDI */
943    
944    #ifdef ALLOW_KPP
945          IF ( calcKPP ) THEN
946            CALL KPP_DO_EXCH( myThid )
947          ENDIF
948    #endif /* ALLOW_KPP */
949    
950  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
951        IF ( fluidIsWater .AND. useDiagnostics ) THEN        IF ( fluidIsWater .AND. useDiagnostics ) THEN
952            CALL DIAGS_RHO_G(
953         I                    rhoInSitu, uVel, vVel, wVel,
954         I                    myTime, myIter, myThid )
955          CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )          CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
956        ENDIF        ENDIF
957        IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN        IF ( calcConvect .AND. useDiagnostics ) THEN
958          CALL DIAGNOSTICS_FILL( IVDConvCount,'CONVADJ ',          CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
959       &                         0, Nr, 0, 1, 1, myThid )       &                               0, Nr, 0, 1, 1, myThid )
960        ENDIF        ENDIF
961  #endif  #endif
962    
963    #ifdef ALLOW_ECCO
964          CALL ECCO_PHYS( myThid )
965    #endif
966    
967  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
968        IF ( debugLevel .GE. debLevB )        IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
      &     CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)  
969  #endif  #endif
970    
971        RETURN        RETURN
Legend:
Removed from v.1.41  
changed lines
  Added in v.1.131
  ViewVC Help
Powered by ViewVC 1.1.22