/[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.27 by heimbach, Wed Mar 8 06:36:39 2006 UTC revision 1.132 by jmc, Tue Jul 9 16:01:01 2013 UTC
# Line 11  C $Name$ Line 11  C $Name$
11  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
12  #  include "KPP_OPTIONS.h"  #  include "KPP_OPTIONS.h"
13  # endif  # endif
14    # ifdef ALLOW_SEAICE
15    #  include "SEAICE_OPTIONS.h"
16    # 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 19  C     !INTERFACE: Line 25  C     !INTERFACE:
25        SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)        SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
26  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
27  C     *==========================================================*  C     *==========================================================*
28  C     | SUBROUTINE DO_OCEANIC_PHYS                                  C     | SUBROUTINE DO_OCEANIC_PHYS
29  C     | o Controlling routine for oceanic physics and  C     | o Controlling routine for oceanic physics and
30  C     |   parameterization  C     |   parameterization
31  C     *==========================================================*  C     *==========================================================*
32  C     | o originally, part of S/R thermodynamics  C     | o originally, part of S/R thermodynamics
# Line 33  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
70  # ifdef EXACT_CONSERV  # ifdef ALLOW_EXF
71  #  include "SURFACE.h"  #  include "ctrl.h"
72    #  include "EXF_FIELDS.h"
73    #  ifdef ALLOW_BULKFORMULAE
74    #   include "EXF_CONSTANTS.h"
75    #  endif
76    # endif
77    # ifdef ALLOW_SEAICE
78    #  include "SEAICE_SIZE.h"
79    #  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 72  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 rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhoKp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
107        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhoKm1  (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 87  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 100  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  #ifdef ALLOW_THSICE        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        IF ( useThSIce .AND. fluidIsWater ) THEN
199  #ifdef ALLOW_DEBUG  # ifdef ALLOW_AUTODIFF_TAMC
200          IF ( debugLevel .GE. debLevB )  CADJ STORE uice,vice         = comlev1, key = ikey_dynamics,
201       &    CALL DEBUG_CALL('THSICE_MAIN',myThid)  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
270          IF ( useSEAICE ) THEN
271    # ifdef ALLOW_AUTODIFF_TAMC
272    cph-adj-test(
273    CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
274    CADJ STORE hsnow  = comlev1, key=ikey_dynamics, kind=isbyte
275    CADJ STORE heff   = comlev1, key=ikey_dynamics, kind=isbyte
276    CADJ STORE tices  = comlev1, key=ikey_dynamics, kind=isbyte
277    CADJ STORE empmr, qnet  = comlev1, key=ikey_dynamics, kind=isbyte
278    CADJ STORE qsw,saltflux = comlev1, key=ikey_dynamics, kind=isbyte
279    CADJ STORE fu, fv = comlev1, key=ikey_dynamics, kind=isbyte
280    cCADJ STORE theta  = comlev1, key=ikey_dynamics, kind=isbyte
281    cCADJ STORE salt   = comlev1, key=ikey_dynamics, kind=isbyte
282    cph-adj-test)
283    c#ifdef ALLOW_EXF
284    CADJ STORE atemp,aqh,precip    = comlev1, key = ikey_dynamics,
285    CADJ &     kind = isbyte
286    CADJ STORE swdown,lwdown       = comlev1, key = ikey_dynamics,
287    CADJ &     kind = isbyte
288    CADJ STORE evap                = comlev1, key = ikey_dynamics,
289    CADJ &     kind = isbyte
290    CADJ STORE uwind,vwind         = comlev1, key = ikey_dynamics,
291    CADJ &     kind = isbyte
292    c#endif
293    CADJ STORE uvel,vvel           = comlev1, key = ikey_dynamics,
294    CADJ &     kind = isbyte
295    #  ifdef SEAICE_CGRID
296    CADJ STORE stressdivergencex   = comlev1, key = ikey_dynamics,
297    CADJ &     kind = isbyte
298    CADJ STORE stressdivergencey   = comlev1, key = ikey_dynamics,
299    CADJ &     kind = isbyte
300    #  endif
301    #  ifdef SEAICE_ALLOW_DYNAMICS
302    CADJ STORE uice                = comlev1, key = ikey_dynamics,
303    CADJ &     kind = isbyte
304    CADJ STORE vice                = comlev1, key = ikey_dynamics,
305    CADJ &     kind = isbyte
306    CADJ STORE dwatn               = comlev1, key = ikey_dynamics,
307    CADJ &     kind = isbyte
308    #   ifdef SEAICE_ALLOW_EVP
309    CADJ STORE seaice_sigma1       = comlev1, key = ikey_dynamics,
310    CADJ &     kind = isbyte
311    CADJ STORE seaice_sigma2       = comlev1, key = ikey_dynamics,
312    CADJ &     kind = isbyte
313    CADJ STORE seaice_sigma12      = comlev1, key = ikey_dynamics,
314    CADJ &     kind = isbyte
315    #   endif
316    #  endif
317    #  ifdef SEAICE_VARIABLE_SALINITY
318    CADJ STORE hsalt               = comlev1, key = ikey_dynamics,
319    CADJ &     kind = isbyte
320    #  endif
321    #  ifdef ATMOSPHERIC_LOADING
322    CADJ STORE pload               = comlev1, key = ikey_dynamics,
323    CADJ &     kind = isbyte
324    CADJ STORE siceload            = comlev1, key = ikey_dynamics,
325    CADJ &     kind = isbyte
326    #  endif
327    #  ifdef NONLIN_FRSURF
328    CADJ STORE recip_hfacc         = comlev1, key = ikey_dynamics,
329    CADJ &     kind = isbyte
330    #  endif
331    #  ifdef ANNUAL_BALANCE
332    CADJ STORE balance_itcount     = comlev1, key = ikey_dynamics,
333    CADJ &     kind = isbyte
334    #  endif /* ANNUAL_BALANCE */
335    #  ifdef ALLOW_THSICE
336    C-- store thSIce vars before advection (called from SEAICE_MODEL) update them:
337    CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
338    CADJ &     kind = isbyte
339    CADJ STORE snowHeight, hOceMxL = comlev1, key = ikey_dynamics,
340    CADJ &     kind = isbyte
341    CADJ STORE Qice1, Qice2  = comlev1, key = ikey_dynamics,
342    CADJ &     kind = isbyte
343    #  endif /* ALLOW_THSICE */
344    # endif /* ALLOW_AUTODIFF_TAMC */
345    # ifdef ALLOW_DEBUG
346            IF (debugMode) CALL DEBUG_CALL('SEAICE_MODEL',myThid)
347    # endif
348            CALL TIMER_START('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
349            CALL SEAICE_MODEL( myTime, myIter, myThid )
350            CALL TIMER_STOP ('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
351    # ifdef ALLOW_COST
352            CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
353    # endif
354          ENDIF
355    #endif /* ALLOW_SEAICE */
356    
357    #ifdef ALLOW_AUTODIFF_TAMC
358    CADJ STORE sst, sss           = comlev1, key = ikey_dynamics,
359    CADJ &     kind = isbyte
360    CADJ STORE qsw                = comlev1, key = ikey_dynamics,
361    CADJ &     kind = isbyte
362    # ifdef ALLOW_SEAICE
363    CADJ STORE area               = comlev1, key = ikey_dynamics,
364    CADJ &     kind = isbyte
365    # endif
366  #endif  #endif
367    
368    #ifdef OLD_THSICE_CALL_SEQUENCE
369    #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
370          IF ( useThSIce .AND. fluidIsWater ) THEN
371    # ifdef ALLOW_AUTODIFF_TAMC
372    cph(
373    #  ifdef NONLIN_FRSURF
374    CADJ STORE uice,vice        = comlev1, key = ikey_dynamics,
375    CADJ &     kind = isbyte
376    CADJ STORE salt,theta       = comlev1, key = ikey_dynamics,
377    CADJ &     kind = isbyte
378    CADJ STORE qnet,qsw, empmr  = comlev1, key = ikey_dynamics,
379    CADJ &     kind = isbyte
380    CADJ STORE hFac_surfC       = comlev1, key = ikey_dynamics,
381    CADJ &     kind = isbyte
382    #  endif
383    # endif
384    # ifdef ALLOW_DEBUG
385            IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
386    # endif
387  C--     Step forward Therm.Sea-Ice variables  C--     Step forward Therm.Sea-Ice variables
388  C       and modify forcing terms including effects from ice  C       and modify forcing terms including effects from ice
389          CALL TIMER_START('THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)          CALL TIMER_START('THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
# Line 129  C       and modify forcing terms includi Line 391  C       and modify forcing terms includi
391          CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)          CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
392        ENDIF        ENDIF
393  #endif /* ALLOW_THSICE */  #endif /* ALLOW_THSICE */
394    #endif /* OLD_THSICE_CALL_SEQUENCE */
395    
396  #ifdef ALLOW_SHELFICE  #ifdef ALLOW_SHELFICE
397    # ifdef ALLOW_AUTODIFF_TAMC
398    CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
399    CADJ &     kind = isbyte
400    # endif
401        IF ( useShelfIce .AND. fluidIsWater ) THEN        IF ( useShelfIce .AND. fluidIsWater ) THEN
402  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
403          IF ( debugLevel .GE. debLevB )         IF (debugMode) CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
      &    CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)  
404  #endif  #endif
405  C     compute temperature and (virtual) salt flux at the  C     compute temperature and (virtual) salt flux at the
406  C     shelf-ice ocean interface  C     shelf-ice ocean interface
407         CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',         CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
408       &       myThid)       &       myThid)
# Line 146  C     shelf-ice ocean interface Line 412  C     shelf-ice ocean interface
412        ENDIF        ENDIF
413  #endif /* ALLOW_SHELFICE */  #endif /* ALLOW_SHELFICE */
414    
415    #ifdef ALLOW_ICEFRONT
416          IF ( useICEFRONT .AND. fluidIsWater ) THEN
417    #ifdef ALLOW_DEBUG
418           IF (debugMode) CALL DEBUG_CALL('ICEFRONT_THERMODYNAMICS',myThid)
419    #endif
420    C     compute temperature and (virtual) salt flux at the
421    C     ice-front ocean interface
422           CALL TIMER_START('ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
423         &       myThid)
424           CALL ICEFRONT_THERMODYNAMICS( myTime, myIter, myThid )
425           CALL TIMER_STOP( 'ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
426         &      myThid)
427          ENDIF
428    #endif /* ALLOW_ICEFRONT */
429    
430    #ifdef ALLOW_SALT_PLUME
431          IF ( useSALT_PLUME ) THEN
432              CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
433          ENDIF
434    #endif /* ALLOW_SALT_PLUME */
435    
436  C--   Freeze water at the surface  C--   Freeze water at the surface
437          IF ( allowFreezing ) THEN
438  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
439  CADJ STORE theta = comlev1, key = ikey_dynamics  CADJ STORE theta = comlev1, key = ikey_dynamics,
440    CADJ &     kind = isbyte
441  #endif  #endif
442        IF ( allowFreezing          CALL FREEZE_SURFACE( myTime, myIter, myThid )
      &                   .AND. .NOT. useSEAICE  
      &                   .AND. .NOT. useThSIce ) THEN  
         CALL FREEZE_SURFACE(  myTime, myIter, myThid )  
443        ENDIF        ENDIF
444    
445  #ifdef COMPONENT_MODULE  #ifdef ALLOW_OCN_COMPON_INTERF
 # ifndef ALLOW_AIM  
446  C--    Apply imported data (from coupled interface) to forcing fields  C--    Apply imported data (from coupled interface) to forcing fields
447  C jmc: do not know precisely where to put this call (bf or af thSIce ?)  C jmc: do not know precisely where to put this call (bf or af thSIce ?)
448         IF ( useCoupler ) THEN        IF ( useCoupler ) THEN
449           CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )           CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
450         ENDIF        ENDIF
451  # endif  #endif /* ALLOW_OCN_COMPON_INTERF */
 #endif /* COMPONENT_MODULE */  
452    
453  #ifdef ALLOW_BALANCE_FLUXES        iMin = 1-OLx
454  C     balance fluxes        iMax = sNx+OLx
455         IF ( balanceEmPmR )        jMin = 1-OLy
456       &        CALL REMOVE_MEAN_RS( 1, EmPmR, maskH, maskH, rA, drF,        jMax = sNy+OLy
457       &        'EmPmR', myTime, myThid )  
458         IF ( balanceQnet )  C---  Determines forcing terms based on external fields
459       &        CALL REMOVE_MEAN_RS( 1, Qnet,  maskH, maskH, rA, drF,  C     relaxation terms, etc.
460       &        'Qnet ', myTime, myThid )  #ifdef ALLOW_DEBUG
461  #endif /* ALLOW_BALANCE_FLUXES */        IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
462    #endif
463    #ifdef ALLOW_AUTODIFF_TAMC
464    CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
465    CADJ &     kind = isbyte
466    #else  /* ALLOW_AUTODIFF_TAMC */
467    C--   if fluid is not water, by-pass surfaceForcing, find_rho, gmredi
468    C     and all vertical mixing schemes, but keep OBCS_CALC
469          IF ( fluidIsWater ) THEN
470    #endif /* ALLOW_AUTODIFF_TAMC */
471            CALL EXTERNAL_FORCING_SURF(
472         I             iMin, iMax, jMin, jMax,
473         I             myTime, myIter, myThid )
474    
475  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
476  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 205  C     These inital values do not alter t Line 501  C     These inital values do not alter t
501  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
502  C     point numbers. This prevents spurious hardware signals due to  C     point numbers. This prevents spurious hardware signals due to
503  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  
          ENDDO  
         ENDDO  
   
504          DO k=1,Nr          DO k=1,Nr
505           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
506            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
507  C This is currently also used by IVDC and Diagnostics  C This is currently used by GMRedi, IVDC, MXL-depth  and Diagnostics
508             sigmaX(i,j,k) = 0. _d 0             sigmaX(i,j,k) = 0. _d 0
509             sigmaY(i,j,k) = 0. _d 0             sigmaY(i,j,k) = 0. _d 0
510             sigmaR(i,j,k) = 0. _d 0             sigmaR(i,j,k) = 0. _d 0
511              ENDDO
512             ENDDO
513            ENDDO
514    
515  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
516            DO j=1-OLy,sNy+OLy
517             DO i=1-OLx,sNx+OLx
518              rhoKm1 (i,j)   = 0. _d 0
519              rhoKp1 (i,j)   = 0. _d 0
520             ENDDO
521            ENDDO
522  cph all the following init. are necessary for TAF  cph all the following init. are necessary for TAF
523  cph although some of these are re-initialised later.  cph although some of these are re-initialised later.
524            DO k=1,Nr
525             DO j=1-OLy,sNy+OLy
526              DO i=1-OLx,sNx+OLx
527               rhoInSitu(i,j,k,bi,bj) = 0.
528    # ifdef ALLOW_GGL90
529               GGL90viscArU(i,j,k,bi,bj)  = 0. _d 0
530               GGL90viscArV(i,j,k,bi,bj)  = 0. _d 0
531               GGL90diffKr(i,j,k,bi,bj)  = 0. _d 0
532    # endif /* ALLOW_GGL90 */
533              ENDDO
534             ENDDO
535            ENDDO
536    #ifdef ALLOW_OFFLINE
537           IF ( calcConvect ) THEN
538    #endif
539            DO k=1,Nr
540             DO j=1-OLy,sNy+OLy
541              DO i=1-OLx,sNx+OLx
542             IVDConvCount(i,j,k,bi,bj) = 0.             IVDConvCount(i,j,k,bi,bj) = 0.
543              ENDDO
544             ENDDO
545            ENDDO
546    #ifdef ALLOW_OFFLINE
547           ENDIF
548           IF ( calcGMRedi ) THEN
549    #endif
550  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
551            DO k=1,Nr
552             DO j=1-OLy,sNy+OLy
553              DO i=1-OLx,sNx+OLx
554             Kwx(i,j,k,bi,bj)  = 0. _d 0             Kwx(i,j,k,bi,bj)  = 0. _d 0
555             Kwy(i,j,k,bi,bj)  = 0. _d 0             Kwy(i,j,k,bi,bj)  = 0. _d 0
556             Kwz(i,j,k,bi,bj)  = 0. _d 0             Kwz(i,j,k,bi,bj)  = 0. _d 0
# Line 243  cph although some of these are re-initia Line 569  cph although some of these are re-initia
569  #  ifdef GM_VISBECK_VARIABLE_K  #  ifdef GM_VISBECK_VARIABLE_K
570             VisbeckK(i,j,bi,bj)   = 0. _d 0             VisbeckK(i,j,bi,bj)   = 0. _d 0
571  #  endif  #  endif
572              ENDDO
573             ENDDO
574            ENDDO
575  # endif /* ALLOW_GMREDI */  # endif /* ALLOW_GMREDI */
576  #endif /* ALLOW_AUTODIFF_TAMC */  #ifdef ALLOW_OFFLINE
577           ENDIF
578           IF ( calcKPP ) THEN
579    #endif
580    # ifdef ALLOW_KPP
581            DO k=1,Nr
582             DO j=1-OLy,sNy+OLy
583              DO i=1-OLx,sNx+OLx
584               KPPdiffKzS(i,j,k,bi,bj)  = 0. _d 0
585               KPPdiffKzT(i,j,k,bi,bj)  = 0. _d 0
586            ENDDO            ENDDO
587           ENDDO           ENDDO
588          ENDDO          ENDDO
589    # endif /* ALLOW_KPP */
590          iMin = 1-OLx  #ifdef ALLOW_OFFLINE
591          iMax = sNx+OLx         ENDIF
592          jMin = 1-OLy  #endif
593          jMax = sNy+OLy  #endif /* ALLOW_AUTODIFF_TAMC */
594    
595  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
596  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
597  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     kind = isbyte
598    CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
599    CADJ &     kind = isbyte
600  CADJ STORE totphihyd(:,:,:,bi,bj)  CADJ STORE totphihyd(:,:,:,bi,bj)
601  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
602    CADJ &     kind = isbyte
603  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
604  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
605  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     kind = isbyte
606    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
607    CADJ &     kind = isbyte
608    # endif
609    # ifdef ALLOW_SALT_PLUME
610    CADJ STORE saltplumedepth(:,:,bi,bj) = comlev1_bibj, key=itdkey,
611    CADJ &     kind = isbyte
612  # endif  # endif
613  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
614    
615    C--   Always compute density (stored in common block) here; even when it is not
616    C     needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
617    #ifdef ALLOW_DEBUG
618            IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D (xNr)',myThid)
619    #endif
620    #ifdef ALLOW_AUTODIFF_TAMC
621            IF ( fluidIsWater ) THEN
622    #endif /* ALLOW_AUTODIFF_TAMC */
623    #ifdef ALLOW_DOWN_SLOPE
624             IF ( useDOWN_SLOPE ) THEN
625               DO k=1,Nr
626                CALL DWNSLP_CALC_RHO(
627         I                  theta, salt,
628         O                  rhoInSitu(1-OLx,1-OLy,k,bi,bj),
629         I                  k, bi, bj, myTime, myIter, myThid )
630               ENDDO
631             ENDIF
632    #endif /* ALLOW_DOWN_SLOPE */
633    #ifdef ALLOW_BBL
634             IF ( useBBL ) THEN
635    C     pkg/bbl requires in-situ bbl density for depths equal to and deeper than the bbl.
636    C     To reduce computation and storage requirement, these densities are stored in the
637    C     dry grid boxes of rhoInSitu.  See BBL_CALC_RHO for details.
638               DO k=Nr,1,-1
639                CALL BBL_CALC_RHO(
640         I                  theta, salt,
641         O                  rhoInSitu,
642         I                  k, bi, bj, myTime, myIter, myThid )
643    
644               ENDDO
645             ENDIF
646    #endif /* ALLOW_BBL */
647             IF ( .NOT. ( useDOWN_SLOPE .OR. useBBL ) ) THEN
648               DO k=1,Nr
649                CALL FIND_RHO_2D(
650         I                iMin, iMax, jMin, jMax, k,
651         I                theta(1-OLx,1-OLy,k,bi,bj),
652         I                salt (1-OLx,1-OLy,k,bi,bj),
653         O                rhoInSitu(1-OLx,1-OLy,k,bi,bj),
654         I                k, bi, bj, myThid )
655               ENDDO
656             ENDIF
657    #ifdef ALLOW_AUTODIFF_TAMC
658            ELSE
659    C-        fluid is not water:
660              DO k=1,Nr
661               DO j=1-OLy,sNy+OLy
662                DO i=1-OLx,sNx+OLx
663                  rhoInSitu(i,j,k,bi,bj) = 0.
664                ENDDO
665               ENDDO
666              ENDDO
667            ENDIF
668    #endif /* ALLOW_AUTODIFF_TAMC */
669    
670  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
671          IF ( debugLevel .GE. debLevB )          IF (debugMode) CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
      &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)  
672  #endif  #endif
673    
674  C--     Start of diagnostic loop  C--     Start of diagnostic loop
# Line 276  C--     Start of diagnostic loop Line 677  C--     Start of diagnostic loop
677  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
678  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?
679  C? Do we still need this?  C? Do we still need this?
680  cph kkey formula corrected.  cph kkey formula corrected.
681  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhoK, rhoKm1, in the case useGMREDI.
682           kkey = (itdkey-1)*Nr + k            kkey = (itdkey-1)*Nr + k
683  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
684    
685    c#ifdef ALLOW_AUTODIFF_TAMC
686    cCADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
687    cCADJ &     kind = isbyte
688    cCADJ STORE salt(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey,
689    cCADJ &     kind = isbyte
690    c#endif /* ALLOW_AUTODIFF_TAMC */
691    
692  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
693  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
694  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN            IF ( calcGMRedi .OR. (k.GT.1 .AND. calcConvect)
695            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)       &         .OR. usePP81 .OR. useMY82 .OR. useGGL90
696       &                   .OR. doDiagsRho.GE.1 ) THEN       &         .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
 #ifdef ALLOW_DEBUG  
             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 )  
   
697              IF (k.GT.1) THEN              IF (k.GT.1) THEN
698  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
699  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,
700  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ &     kind = isbyte
701  #endif /* ALLOW_AUTODIFF_TAMC */  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
702               CALL FIND_RHO(  CADJ &     kind = isbyte
703       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,  CADJ STORE rhokm1 (bi,bj)       = comlev1_bibj_k, key=kkey,
704       I        theta, salt,  CADJ &     kind = isbyte
705       O        rhoKm1,  #endif /* ALLOW_AUTODIFF_TAMC */
706       I        myThid )               CALL FIND_RHO_2D(
707         I                 iMin, iMax, jMin, jMax, k,
708         I                 theta(1-OLx,1-OLy,k-1,bi,bj),
709         I                 salt (1-OLx,1-OLy,k-1,bi,bj),
710         O                 rhoKm1,
711         I                 k-1, bi, bj, myThid )
712              ENDIF              ENDIF
713  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
714              IF ( debugLevel .GE. debLevB )              IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
      &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)  
715  #endif  #endif
716    cph Avoid variable aliasing for adjoint !!!
717                DO j=jMin,jMax
718                 DO i=iMin,iMax
719                  rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
720                 ENDDO
721                ENDDO
722              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
723       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
724       I             rhoK, rhoKm1, rhoK,       I             rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
725       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
726       I             myThid )       I             myThid )
           ENDIF  
   
727  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
728  ctest# ifndef GM_EXCLUDE_CLIPPING  #ifdef GMREDI_WITH_STABLE_ADJOINT
729  CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
730  ctest# endif  cgf -> cuts adjoint dependency from slope to state
731  CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte              CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
732                CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
733                CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
734    #endif
735  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
736              ENDIF
737    
738  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
739  c ==> should use sigmaR !!!            IF (k.GT.1 .AND. calcConvect) THEN
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
740  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
741              IF ( debugLevel .GE. debLevB )              IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid)
      &       CALL DEBUG_CALL('CALC_IVDC',myThid)  
742  #endif  #endif
743              CALL CALC_IVDC(              CALL CALC_IVDC(
744       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
745       I        rhoKm1, rhoK,       I        sigmaR,
746       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
747            ENDIF            ENDIF
748    
749  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
750            IF ( doDiagsRho.GE.2 ) THEN            IF ( doDiagsRho.GE.4 ) THEN
751              CALL DIAGS_RHO( k, bi, bj,              CALL DIAGS_RHO_L( doDiagsRho, k, bi, bj,
752       I                      rhoK, rhoKm1,       I                        rhoInSitu(1-OLx,1-OLy,1,bi,bj),
753       I                      myTime, myIter, myThid)       I                        rhoKm1, wVel,
754         I                        myTime, myIter, myThid )
755            ENDIF            ENDIF
756  #endif  #endif
757    
758  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
759          ENDDO          ENDDO
760    
761  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_AUTODIFF_TAMC
762  c       IF ( useDiagnostics .AND.  CADJ STORE IVDConvCount(:,:,:,bi,bj)
763  c    &       (useGMRedi .OR. ivdc_kappa.NE.0.) ) THEN  CADJ &     = comlev1_bibj, key=itdkey,
764          IF ( doDiagsRho.GE.1 ) THEN  CADJ &     kind = isbyte
           CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,  
      &         2, bi, bj, myThid)  
         ENDIF  
765  #endif  #endif
766    
767  #ifdef  ALLOW_OBCS  C--     Diagnose Mixed Layer Depth:
768  C--     Calculate future values on open boundaries          IF ( calcGMRedi .OR. MOD(doDiagsRho,2).EQ.1 ) THEN
769          IF (useOBCS) THEN            CALL CALC_OCE_MXLAYER(
770  #ifdef ALLOW_DEBUG       I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
771            IF ( debugLevel .GE. debLevB )       I              bi, bj, myTime, myIter, myThid )
      &     CALL DEBUG_CALL('OBCS_CALC',myThid)  
 #endif  
           CALL OBCS_CALC( bi, bj, myTime+deltaTclock, myIter+1,  
      I            uVel, vVel, wVel, theta, salt,  
      I            myThid )  
772          ENDIF          ENDIF
 #endif  /* ALLOW_OBCS */  
773    
774  #ifndef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_SALT_PLUME
775          IF ( fluidIsWater ) THEN          IF ( useSALT_PLUME ) THEN
776  #endif            CALL SALT_PLUME_CALC_DEPTH(
777  C--     Determines forcing terms based on external fields       I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
778  C       relaxation terms, etc.       I              bi, bj, myTime, myIter, myThid )
 #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 )  
 #ifndef ALLOW_AUTODIFF_TAMC  
779          ENDIF          ENDIF
780  #endif  #endif /* ALLOW_SALT_PLUME */
781  #ifdef ALLOW_AUTODIFF_TAMC  
782  # ifdef EXACT_CONSERV  #ifdef ALLOW_DIAGNOSTICS
783  cph-test          IF ( MOD(doDiagsRho,4).GE.2 ) THEN
784  cphCADJ STORE PmEpR(:,:,bi,bj)            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,
785  cphCADJ &     = comlev1_bibj, key=itdkey, byte=isbyte       &         2, bi, bj, myThid)
786  # endif          ENDIF
787  #endif  #endif /* ALLOW_DIAGNOSTICS */
788    
789    C--    This is where EXTERNAL_FORCING_SURF(bi,bj) used to be called;
790    C      now called earlier, before bi,bj loop.
791    
792  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
793  cph needed for KPP  cph needed for KPP
794  CADJ STORE surfaceForcingU(:,:,bi,bj)  CADJ STORE surfaceForcingU(:,:,bi,bj)
795  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
796    CADJ &     kind = isbyte
797  CADJ STORE surfaceForcingV(:,:,bi,bj)  CADJ STORE surfaceForcingV(:,:,bi,bj)
798  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
799    CADJ &     kind = isbyte
800  CADJ STORE surfaceForcingS(:,:,bi,bj)  CADJ STORE surfaceForcingS(:,:,bi,bj)
801  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
802    CADJ &     kind = isbyte
803  CADJ STORE surfaceForcingT(:,:,bi,bj)  CADJ STORE surfaceForcingT(:,:,bi,bj)
804  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
805    CADJ &     kind = isbyte
806  CADJ STORE surfaceForcingTice(:,:,bi,bj)  CADJ STORE surfaceForcingTice(:,:,bi,bj)
807  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
808    CADJ &     kind = isbyte
809  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
810    
 #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 */  
   
811  #ifdef  ALLOW_KPP  #ifdef  ALLOW_KPP
812  C--     Compute KPP mixing coefficients  C--     Compute KPP mixing coefficients
813          IF (useKPP) THEN          IF ( calcKPP ) THEN
814  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
815            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid)
      &     CALL DEBUG_CALL('KPP_CALC',myThid)  
816  #endif  #endif
817              CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
818            CALL KPP_CALC(            CALL KPP_CALC(
819       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myIter, myThid )
820  #ifdef ALLOW_AUTODIFF_TAMC            CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
821    #if (defined ALLOW_AUTODIFF_TAMC) && !(defined ALLOW_OFFLINE)
822          ELSE          ELSE
823            CALL KPP_CALC_DUMMY(            CALL KPP_CALC_DUMMY(
824       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myIter, myThid )
825  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC and not ALLOW_OFFLINE */
826          ENDIF          ENDIF
   
827  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
828    
829  #ifdef  ALLOW_PP81  #ifdef  ALLOW_PP81
830  C--     Compute PP81 mixing coefficients  C--     Compute PP81 mixing coefficients
831          IF (usePP81) THEN          IF (usePP81) THEN
832  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
833            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid)
      &     CALL DEBUG_CALL('PP81_CALC',myThid)  
834  #endif  #endif
835            CALL PP81_CALC(            CALL PP81_CALC(
836       I                  bi, bj, myTime, myThid )       I                     bi, bj, sigmaR, myTime, myIter, myThid )
837          ENDIF          ENDIF
838  #endif /* ALLOW_PP81 */  #endif /* ALLOW_PP81 */
839    
# Line 494  C--     Compute PP81 mixing coefficients Line 841  C--     Compute PP81 mixing coefficients
841  C--     Compute MY82 mixing coefficients  C--     Compute MY82 mixing coefficients
842          IF (useMY82) THEN          IF (useMY82) THEN
843  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
844            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid)
      &     CALL DEBUG_CALL('MY82_CALC',myThid)  
845  #endif  #endif
846            CALL MY82_CALC(            CALL MY82_CALC(
847       I                  bi, bj, myTime, myThid )       I                     bi, bj, sigmaR, myTime, myIter, myThid )
848          ENDIF          ENDIF
849  #endif /* ALLOW_MY82 */  #endif /* ALLOW_MY82 */
850    
851  #ifdef  ALLOW_GGL90  #ifdef  ALLOW_GGL90
852    #ifdef ALLOW_AUTODIFF_TAMC
853    CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
854    CADJ &     kind = isbyte
855    #endif /* ALLOW_AUTODIFF_TAMC */
856  C--     Compute GGL90 mixing coefficients  C--     Compute GGL90 mixing coefficients
857          IF (useGGL90) THEN          IF (useGGL90) THEN
858  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
859            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
      &     CALL DEBUG_CALL('GGL90_CALC',myThid)  
860  #endif  #endif
861              CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
862            CALL GGL90_CALC(            CALL GGL90_CALC(
863       I                  bi, bj, myTime, myThid )       I                     bi, bj, sigmaR, myTime, myIter, myThid )
864              CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
865          ENDIF          ENDIF
866  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
867    
868  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
869          IF ( taveFreq.GT. 0. _d 0 .AND. fluidIsWater ) THEN          IF ( taveFreq.GT. 0. _d 0 ) THEN
870            CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)            CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
871          ENDIF          ENDIF
872          IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN          IF ( taveFreq.GT.0. .AND. calcConvect ) THEN
873            CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,            CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
874       I                           Nr, deltaTclock, bi, bj, myThid)       I                           Nr, deltaTClock, bi, bj, myThid)
875          ENDIF          ENDIF
876  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
877    
878    #ifdef ALLOW_GMREDI
879    #ifdef ALLOW_AUTODIFF_TAMC
880    # ifndef GM_EXCLUDE_CLIPPING
881    cph storing here is needed only for one GMREDI_OPTIONS:
882    cph define GM_BOLUS_ADVEC
883    cph keep it although TAF says you dont need to.
884    cph but I have avoided the #ifdef for now, in case more things change
885    CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey,
886    CADJ &     kind = isbyte
887    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey,
888    CADJ &     kind = isbyte
889    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey,
890    CADJ &     kind = isbyte
891    # endif
892    #endif /* ALLOW_AUTODIFF_TAMC */
893    
894    C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
895            IF ( calcGMRedi ) THEN
896    #ifdef ALLOW_DEBUG
897              IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
898    #endif
899              CALL GMREDI_CALC_TENSOR(
900         I             iMin, iMax, jMin, jMax,
901         I             sigmaX, sigmaY, sigmaR,
902         I             bi, bj, myTime, myIter, myThid )
903    #if (defined ALLOW_AUTODIFF_TAMC) && !(defined ALLOW_OFFLINE)
904            ELSE
905              CALL GMREDI_CALC_TENSOR_DUMMY(
906         I             iMin, iMax, jMin, jMax,
907         I             sigmaX, sigmaY, sigmaR,
908         I             bi, bj, myTime, myIter, myThid )
909    #endif /* ALLOW_AUTODIFF_TAMC and not ALLOW_OFFLINE */
910            ENDIF
911    #endif /* ALLOW_GMREDI */
912    
913    #ifdef ALLOW_DOWN_SLOPE
914            IF ( useDOWN_SLOPE ) THEN
915    C--     Calculate Downsloping Flow for Down_Slope parameterization
916             IF ( usingPCoords ) THEN
917              CALL DWNSLP_CALC_FLOW(
918         I                bi, bj, kSurfC, rhoInSitu,
919         I                myTime, myIter, myThid )
920             ELSE
921              CALL DWNSLP_CALC_FLOW(
922         I                bi, bj, kLowC, rhoInSitu,
923         I                myTime, myIter, myThid )
924             ENDIF
925            ENDIF
926    #endif /* ALLOW_DOWN_SLOPE */
927    
928  C--   end bi,bj loops.  C--   end bi,bj loops.
929         ENDDO         ENDDO
930        ENDDO        ENDDO
931    
932    #ifndef ALLOW_AUTODIFF_TAMC
933    C---  if fluid Is Water: end
934          ENDIF
935    #endif
936    
937    #ifdef ALLOW_BBL
938          IF ( useBBL ) THEN
939           CALL BBL_CALC_RHS(
940         I                          myTime, myIter, myThid )
941          ENDIF
942    #endif /* ALLOW_BBL */
943    
944    #ifdef ALLOW_MYPACKAGE
945          IF ( useMYPACKAGE ) THEN
946           CALL MYPACKAGE_CALC_RHS(
947         I                          myTime, myIter, myThid )
948          ENDIF
949    #endif /* ALLOW_MYPACKAGE */
950    
951    #ifdef ALLOW_GMREDI
952          IF ( calcGMRedi ) THEN
953            CALL GMREDI_DO_EXCH( myTime, myIter, myThid )
954          ENDIF
955    #endif /* ALLOW_GMREDI */
956    
957    #ifdef ALLOW_KPP
958          IF ( calcKPP ) THEN
959            CALL KPP_DO_EXCH( myThid )
960          ENDIF
961    #endif /* ALLOW_KPP */
962    
963  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
964        IF ( fluidIsWater .AND. useDiagnostics ) THEN        IF ( fluidIsWater .AND. useDiagnostics ) THEN
965            CALL DIAGS_RHO_G(
966         I                    rhoInSitu, uVel, vVel, wVel,
967         I                    myTime, myIter, myThid )
968          CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )          CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
969        ENDIF        ENDIF
970        IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN        IF ( calcConvect .AND. useDiagnostics ) THEN
971          CALL DIAGNOSTICS_FILL( IVDConvCount,'CONVADJ ',          CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
972       &                         0, Nr, 0, 1, 1, myThid )       &                               0, Nr, 0, 1, 1, myThid )
973        ENDIF        ENDIF
974  #endif  #endif
975    
976    #ifdef ALLOW_ECCO
977          CALL ECCO_PHYS( myThid )
978    #endif
979    
980  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
981           IF ( debugLevel .GE. debLevB )        IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
      &    CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)  
982  #endif  #endif
983    
984        RETURN        RETURN
Legend:
Removed from v.1.27  
changed lines
  Added in v.1.132
  ViewVC Help
Powered by ViewVC 1.1.22