/[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.3 by jmc, Tue Jul 13 16:48:48 2004 UTC revision 1.136 by heimbach, Wed May 21 10:44:59 2014 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "PACKAGES_CONFIG.h"  #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    #ifdef ALLOW_AUTODIFF
7    # include "AUTODIFF_OPTIONS.h"
8    #endif
9    #ifdef ALLOW_CTRL
10    # include "CTRL_OPTIONS.h"
11    #endif
12    #ifdef ALLOW_SALT_PLUME
13    # include "SALT_PLUME_OPTIONS.h"
14    #endif
15    
16  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
17  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
18  #  include "GMREDI_OPTIONS.h"  #  include "GMREDI_OPTIONS.h"
19  # endif  # endif
20  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
21  #  include "KPP_OPTIONS.h"  #  include "KPP_OPTIONS.h"
22  # endif  # endif
23  #endif /* ALLOW_AUTODIFF_TAMC */  # ifdef ALLOW_SEAICE
24    #  include "SEAICE_OPTIONS.h"
25    # endif
26    # ifdef ALLOW_EXF
27    #  include "EXF_OPTIONS.h"
28    # endif
29    #endif /* ALLOW_AUTODIFF */
30    
31  CBOP  CBOP
32  C     !ROUTINE: DO_OCEANIC_PHYS  C     !ROUTINE: DO_OCEANIC_PHYS
# Line 19  C     !INTERFACE: Line 34  C     !INTERFACE:
34        SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)        SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
35  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
36  C     *==========================================================*  C     *==========================================================*
37  C     | SUBROUTINE DO_OCEANIC_PHYS                                  C     | SUBROUTINE DO_OCEANIC_PHYS
38  C     | o Controlling routine for oceanic physics and  C     | o Controlling routine for oceanic physics and
39  C     |   parameterization  C     |   parameterization
40  C     *==========================================================*  C     *==========================================================*
41  C     | o originally, part of S/R thermodynamics  C     | o originally, part of S/R thermodynamics
# Line 33  C     == Global variables === Line 48  C     == Global variables ===
48  #include "SIZE.h"  #include "SIZE.h"
49  #include "EEPARAMS.h"  #include "EEPARAMS.h"
50  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
51  #include "GRID.h"  #include "GRID.h"
52  c #include "GAD.h"  #include "DYNVARS.h"
53  c #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_TIMEAVE
54  c #include "TR1.h"  # include "TIMEAVE_STATV.h"
55  c #endif  #endif
56  c #ifdef ALLOW_PTRACERS  #ifdef ALLOW_OFFLINE
57  c #include "PTRACERS_SIZE.h"  # include "OFFLINE_SWITCH.h"
58  c #include "PTRACERS.h"  #endif
 c #endif  
 c #ifdef ALLOW_TIMEAVE  
 c #include "TIMEAVE_STATV.h"  
 c #endif  
59    
60  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
61    # include "AUTODIFF_MYFIELDS.h"
62  # include "tamc.h"  # include "tamc.h"
63  # include "tamc_keys.h"  # include "tamc_keys.h"
64  # include "FFIELDS.h"  # include "FFIELDS.h"
65    # include "SURFACE.h"
66  # include "EOS.h"  # include "EOS.h"
67  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
68  #  include "KPP.h"  #  include "KPP.h"
69  # endif  # endif
70    # ifdef ALLOW_GGL90
71    #  include "GGL90.h"
72    # endif
73  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
74  #  include "GMREDI.h"  #  include "GMREDI.h"
75  # endif  # endif
76  # ifdef ALLOW_EBM  # ifdef ALLOW_EBM
77  #  include "EBM.h"  #  include "EBM.h"
78  # endif  # endif
79  #endif /* ALLOW_AUTODIFF_TAMC */  # ifdef ALLOW_EXF
80    #  include "ctrl.h"
81    #  include "EXF_FIELDS.h"
82    #  ifdef ALLOW_BULKFORMULAE
83    #   include "EXF_CONSTANTS.h"
84    #  endif
85    # endif
86    # ifdef ALLOW_SEAICE
87    #  include "SEAICE_SIZE.h"
88    #  include "SEAICE.h"
89    #  include "SEAICE_PARAMS.h"
90    # endif
91    # ifdef ALLOW_THSICE
92    #  include "THSICE_VARS.h"
93    # endif
94    # ifdef ALLOW_SALT_PLUME
95    #  include "SALT_PLUME.h"
96    # endif
97    #endif /* ALLOW_AUTODIFF */
98    
99  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
100  C     == Routine arguments ==  C     == Routine arguments ==
101  C     myTime - Current time in simulation  C     myTime :: Current time in simulation
102  C     myIter - Current iteration number in simulation  C     myIter :: Current iteration number in simulation
103  C     myThid - Thread number for this instance of the routine.  C     myThid :: Thread number for this instance of the routine.
104        _RL myTime        _RL myTime
105        INTEGER myIter        INTEGER myIter
106        INTEGER myThid        INTEGER myThid
107    
108  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
109  C     == Local variables  C     == Local variables
110  C     rhoK, rhoKM1   - Density at current level, and level above  C     rhoK, rhoKm1  :: Density at current level, and level above
111  C     useVariableK   = T when vertical diffusion is not constant  C     iMin, iMax    :: Ranges and sub-block indices on which calculations
 C     iMin, iMax     - Ranges and sub-block indices on which calculations  
112  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
113  C     bi, bj  C     bi, bj        :: tile indices
114  C     k, kup,        - Index for layer above and below. kup and kDown  C     i,j,k         :: loop indices
115  C     kDown, km1       are switched with layer to be the appropriate        _RL rhoKp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116  C                      index into fVerTerm.        _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)  
117        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
118        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
119        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL kp1Msk  
       LOGICAL useVariableK  
120        INTEGER iMin, iMax        INTEGER iMin, iMax
121        INTEGER jMin, jMax        INTEGER jMin, jMax
122        INTEGER bi, bj        INTEGER bi, bj
123        INTEGER i, j        INTEGER i, j, k
124        INTEGER k, km1, kup, kDown        INTEGER doDiagsRho
125        INTEGER iTracer, ip        LOGICAL calcGMRedi, calcKPP, calcConvect
126          CHARACTER*(MAX_LEN_MBUF) msgBuf
127    #ifdef ALLOW_DIAGNOSTICS
128          LOGICAL  DIAGNOSTICS_IS_ON
129          EXTERNAL DIAGNOSTICS_IS_ON
130    #endif /* ALLOW_DIAGNOSTICS */
131    
132    C#ifdef ALLOW_SALT_PLUME
133    C#ifdef SALT_PLUME_VOLUME
134    C      _RL temparray(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
135    C#endif
136    C#endif
137    
138  CEOP  CEOP
139    
 #ifdef ALLOW_DEBUG  
          IF ( debugLevel .GE. debLevB )  
      &    CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)  
 #endif  
   
140  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
141  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
       ikey = 1  
142        itdkey = 1        itdkey = 1
143  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
144    
145    #ifdef ALLOW_DEBUG
146          IF (debugMode) CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
147    #endif
148    
149          doDiagsRho = 0
150    #ifdef ALLOW_DIAGNOSTICS
151          IF ( useDiagnostics .AND. fluidIsWater ) THEN
152            IF ( DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) )
153         &       doDiagsRho = doDiagsRho + 1
154            IF ( DIAGNOSTICS_IS_ON('DRHODR  ',myThid) )
155         &       doDiagsRho = doDiagsRho + 2
156            IF ( DIAGNOSTICS_IS_ON('WdRHO_P ',myThid) )
157         &       doDiagsRho = doDiagsRho + 4
158            IF ( DIAGNOSTICS_IS_ON('WdRHOdP ',myThid) )
159         &       doDiagsRho = doDiagsRho + 8
160          ENDIF
161    #endif /* ALLOW_DIAGNOSTICS */
162    
163          calcGMRedi  = useGMRedi
164          calcKPP     = useKPP
165          calcConvect = ivdc_kappa.NE.0.
166    #ifdef ALLOW_OFFLINE
167          IF ( useOffLine ) THEN
168            calcGMRedi = useGMRedi .AND. .NOT.offlineLoadGMRedi
169            calcKPP    = useKPP    .AND. .NOT.offlineLoadKPP
170            calcConvect=calcConvect.AND. .NOT.offlineLoadConvec
171          ENDIF
172    #endif /* ALLOW_OFFLINE */
173    
174    #ifdef  ALLOW_OBCS
175          IF (useOBCS) THEN
176    C--   Calculate future values on open boundaries
177    C--   moved before SEAICE_MODEL call since SEAICE_MODEL needs seaice-obcs fields
178    # ifdef ALLOW_AUTODIFF_TAMC
179    CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
180    CADJ STORE salt  = comlev1, key=ikey_dynamics, kind=isbyte
181    # endif
182    # ifdef ALLOW_DEBUG
183           IF (debugMode) CALL DEBUG_CALL('OBCS_CALC',myThid)
184    # endif
185           CALL OBCS_CALC( myTime+deltaTClock, myIter+1,
186         I                 uVel, vVel, wVel, theta, salt, myThid )
187          ENDIF
188    #endif  /* ALLOW_OBCS */
189    
190    #ifdef ALLOW_AUTODIFF
191    # ifdef ALLOW_SALT_PLUME
192          DO bj=myByLo(myThid),myByHi(myThid)
193           DO bi=myBxLo(myThid),myBxHi(myThid)
194            DO j=1-OLy,sNy+OLy
195             DO i=1-OLx,sNx+OLx
196              saltPlumeDepth(i,j,bi,bj) = 0. _d 0
197              saltPlumeFlux(i,j,bi,bj)  = 0. _d 0
198             ENDDO
199            ENDDO
200           ENDDO
201          ENDDO
202    # endif
203    #endif /* ALLOW_AUTODIFF */
204    
205    #ifdef ALLOW_FRAZIL
206          IF ( useFRAZIL ) THEN
207    C--   Freeze water in the ocean interior and let it rise to the surface
208           CALL FRAZIL_CALC_RHS( myTime, myIter, myThid )
209          ENDIF
210    #endif /* ALLOW_FRAZIL */
211    
212    #ifndef OLD_THSICE_CALL_SEQUENCE
213    #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
214          IF ( useThSIce .AND. fluidIsWater ) THEN
215    # ifdef ALLOW_AUTODIFF_TAMC
216    CADJ STORE uice,vice         = comlev1, key = ikey_dynamics,
217    CADJ &     kind = isbyte
218    CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
219    CADJ &     kind = isbyte
220    CADJ STORE snowHeight, Tsrf  = comlev1, key = ikey_dynamics,
221    CADJ &     kind = isbyte
222    CADJ STORE Qice1, Qice2      = comlev1, key = ikey_dynamics,
223    CADJ &     kind = isbyte
224    CADJ STORE sHeating, snowAge = comlev1, key = ikey_dynamics,
225    CADJ &     kind = isbyte
226    CADJ STORE hocemxl = comlev1, key = ikey_dynamics,
227    CADJ &     kind = isbyte
228    CADJ STORE icflxsw = comlev1, key = ikey_dynamics,
229    CADJ &     kind = isbyte
230    CADJ STORE salt,theta        = comlev1, key = ikey_dynamics,
231    CADJ &     kind = isbyte
232    CADJ STORE uvel,vvel         = comlev1, key = ikey_dynamics,
233    CADJ &     kind = isbyte
234    CADJ STORE qnet,qsw, empmr   = comlev1, key = ikey_dynamics,
235    CADJ &     kind = isbyte
236    CADJ STORE atemp,aqh,precip  = comlev1, key = ikey_dynamics,
237    CADJ &     kind = isbyte
238    CADJ STORE swdown,lwdown     = comlev1, key = ikey_dynamics,
239    CADJ &     kind = isbyte
240    #  ifdef NONLIN_FRSURF
241    CADJ STORE hFac_surfC       = comlev1, key = ikey_dynamics,
242    CADJ &     kind = isbyte
243    #  endif
244    # endif /* ALLOW_AUTODIFF_TAMC */
245    # ifdef ALLOW_DEBUG
246            IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
247    # endif
248    C--     Step forward Therm.Sea-Ice variables
249    C       and modify forcing terms including effects from ice
250            CALL TIMER_START('THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
251            CALL THSICE_MAIN( myTime, myIter, myThid )
252            CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
253          ENDIF
254    #endif /* ALLOW_THSICE */
255    #endif /* ndef OLD_THSICE_CALL_SEQUENCE */
256    
257    #ifdef ALLOW_SEAICE
258    # ifdef ALLOW_AUTODIFF
259    CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
260    CADJ STORE fu,fv  = comlev1, key=ikey_dynamics, kind=isbyte
261    CADJ STORE qnet   = comlev1, key=ikey_dynamics, kind=isbyte
262    CADJ STORE qsw    = comlev1, key=ikey_dynamics, kind=isbyte
263    CADJ STORE theta  = comlev1, key=ikey_dynamics, kind=isbyte
264    CADJ STORE salt   = comlev1, key=ikey_dynamics, kind=isbyte
265    #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
266    CADJ STORE evap   = comlev1, key=ikey_dynamics, kind=isbyte
267    #endif
268          IF ( .NOT.useSEAICE .AND. SEAICEadjMODE .EQ. -1 ) THEN
269            CALL SEAICE_FAKE( myTime, myIter, myThid )
270          ENDIF
271    CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
272    CADJ STORE fu,fv  = comlev1, key=ikey_dynamics, kind=isbyte
273    CADJ STORE qnet   = comlev1, key=ikey_dynamics, kind=isbyte
274    CADJ STORE qsw    = comlev1, key=ikey_dynamics, kind=isbyte
275    CADJ STORE theta  = comlev1, key=ikey_dynamics, kind=isbyte
276    CADJ STORE salt   = comlev1, key=ikey_dynamics, kind=isbyte
277    #if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
278    CADJ STORE evap   = comlev1, key=ikey_dynamics, kind=isbyte
279    #endif
280    # endif /* ALLOW_AUTODIFF */
281    #endif /* ALLOW_SEAICE */
282    
283    #ifdef ALLOW_SEAICE
284          IF ( useSEAICE ) THEN
285    # ifdef ALLOW_AUTODIFF_TAMC
286    cph-adj-test(
287    CADJ STORE area   = comlev1, key=ikey_dynamics, kind=isbyte
288    CADJ STORE hsnow  = comlev1, key=ikey_dynamics, kind=isbyte
289    CADJ STORE heff   = comlev1, key=ikey_dynamics, kind=isbyte
290    CADJ STORE tices  = comlev1, key=ikey_dynamics, kind=isbyte
291    CADJ STORE empmr, qnet  = comlev1, key=ikey_dynamics, kind=isbyte
292    CADJ STORE qsw,saltflux = comlev1, key=ikey_dynamics, kind=isbyte
293    CADJ STORE fu, fv = comlev1, key=ikey_dynamics, kind=isbyte
294    cCADJ STORE theta  = comlev1, key=ikey_dynamics, kind=isbyte
295    cCADJ STORE salt   = comlev1, key=ikey_dynamics, kind=isbyte
296    cph-adj-test)
297    c#ifdef ALLOW_EXF
298    CADJ STORE atemp,aqh,precip    = comlev1, key = ikey_dynamics,
299    CADJ &     kind = isbyte
300    CADJ STORE swdown,lwdown       = comlev1, key = ikey_dynamics,
301    CADJ &     kind = isbyte
302    CADJ STORE evap                = comlev1, key = ikey_dynamics,
303    CADJ &     kind = isbyte
304    CADJ STORE uwind,vwind         = comlev1, key = ikey_dynamics,
305    CADJ &     kind = isbyte
306    c#endif
307    CADJ STORE uvel,vvel           = comlev1, key = ikey_dynamics,
308    CADJ &     kind = isbyte
309    #  ifdef SEAICE_CGRID
310    CADJ STORE stressdivergencex   = comlev1, key = ikey_dynamics,
311    CADJ &     kind = isbyte
312    CADJ STORE stressdivergencey   = comlev1, key = ikey_dynamics,
313    CADJ &     kind = isbyte
314    #  endif
315    #  ifdef SEAICE_ALLOW_DYNAMICS
316    CADJ STORE uice                = comlev1, key = ikey_dynamics,
317    CADJ &     kind = isbyte
318    CADJ STORE vice                = comlev1, key = ikey_dynamics,
319    CADJ &     kind = isbyte
320    CADJ STORE dwatn               = comlev1, key = ikey_dynamics,
321    CADJ &     kind = isbyte
322    #   ifdef SEAICE_ALLOW_EVP
323    CADJ STORE seaice_sigma1       = comlev1, key = ikey_dynamics,
324    CADJ &     kind = isbyte
325    CADJ STORE seaice_sigma2       = comlev1, key = ikey_dynamics,
326    CADJ &     kind = isbyte
327    CADJ STORE seaice_sigma12      = comlev1, key = ikey_dynamics,
328    CADJ &     kind = isbyte
329    #   endif
330    #  endif
331    #  ifdef SEAICE_VARIABLE_SALINITY
332    CADJ STORE hsalt               = comlev1, key = ikey_dynamics,
333    CADJ &     kind = isbyte
334    #  endif
335    #  ifdef ATMOSPHERIC_LOADING
336    CADJ STORE pload               = comlev1, key = ikey_dynamics,
337    CADJ &     kind = isbyte
338    CADJ STORE siceload            = comlev1, key = ikey_dynamics,
339    CADJ &     kind = isbyte
340    #  endif
341    #  ifdef NONLIN_FRSURF
342    CADJ STORE recip_hfacc         = comlev1, key = ikey_dynamics,
343    CADJ &     kind = isbyte
344    #  endif
345    #  ifdef ANNUAL_BALANCE
346    CADJ STORE balance_itcount     = comlev1, key = ikey_dynamics,
347    CADJ &     kind = isbyte
348    #  endif /* ANNUAL_BALANCE */
349    #  ifdef ALLOW_THSICE
350    C-- store thSIce vars before advection (called from SEAICE_MODEL) update them:
351    CADJ STORE iceMask,iceHeight = comlev1, key = ikey_dynamics,
352    CADJ &     kind = isbyte
353    CADJ STORE snowHeight, hOceMxL = comlev1, key = ikey_dynamics,
354    CADJ &     kind = isbyte
355    CADJ STORE Qice1, Qice2  = comlev1, key = ikey_dynamics,
356    CADJ &     kind = isbyte
357    #  endif /* ALLOW_THSICE */
358    # endif /* ALLOW_AUTODIFF_TAMC */
359    # ifdef ALLOW_DEBUG
360            IF (debugMode) CALL DEBUG_CALL('SEAICE_MODEL',myThid)
361    # endif
362            CALL TIMER_START('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
363            CALL SEAICE_MODEL( myTime, myIter, myThid )
364            CALL TIMER_STOP ('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
365    # ifdef ALLOW_COST
366            CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
367    # endif
368          ENDIF
369    #endif /* ALLOW_SEAICE */
370    
371    #ifdef ALLOW_AUTODIFF_TAMC
372    CADJ STORE sst, sss           = comlev1, key = ikey_dynamics,
373    CADJ &     kind = isbyte
374    CADJ STORE qsw                = comlev1, key = ikey_dynamics,
375    CADJ &     kind = isbyte
376    # ifdef ALLOW_SEAICE
377    CADJ STORE area               = comlev1, key = ikey_dynamics,
378    CADJ &     kind = isbyte
379    # endif
380    #endif
381    
382    #ifdef OLD_THSICE_CALL_SEQUENCE
383    #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
384          IF ( useThSIce .AND. fluidIsWater ) THEN
385    # ifdef ALLOW_AUTODIFF_TAMC
386    cph(
387    #  ifdef NONLIN_FRSURF
388    CADJ STORE uice,vice        = comlev1, key = ikey_dynamics,
389    CADJ &     kind = isbyte
390    CADJ STORE salt,theta       = comlev1, key = ikey_dynamics,
391    CADJ &     kind = isbyte
392    CADJ STORE qnet,qsw, empmr  = comlev1, key = ikey_dynamics,
393    CADJ &     kind = isbyte
394    CADJ STORE hFac_surfC       = comlev1, key = ikey_dynamics,
395    CADJ &     kind = isbyte
396    #  endif
397    # endif
398    # ifdef ALLOW_DEBUG
399            IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
400    # endif
401    C--     Step forward Therm.Sea-Ice variables
402    C       and modify forcing terms including effects from ice
403            CALL TIMER_START('THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
404            CALL THSICE_MAIN( myTime, myIter, myThid )
405            CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
406          ENDIF
407    #endif /* ALLOW_THSICE */
408    #endif /* OLD_THSICE_CALL_SEQUENCE */
409    
410    #ifdef ALLOW_SHELFICE
411    # ifdef ALLOW_AUTODIFF_TAMC
412    CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
413    CADJ &     kind = isbyte
414    # endif
415          IF ( useShelfIce .AND. fluidIsWater ) THEN
416    #ifdef ALLOW_DEBUG
417           IF (debugMode) CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
418    #endif
419    C     compute temperature and (virtual) salt flux at the
420    C     shelf-ice ocean interface
421           CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
422         &       myThid)
423           CALL SHELFICE_THERMODYNAMICS( myTime, myIter, myThid )
424           CALL TIMER_STOP( 'SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
425         &      myThid)
426          ENDIF
427    #endif /* ALLOW_SHELFICE */
428    
429    #ifdef ALLOW_ICEFRONT
430          IF ( useICEFRONT .AND. fluidIsWater ) THEN
431    #ifdef ALLOW_DEBUG
432           IF (debugMode) CALL DEBUG_CALL('ICEFRONT_THERMODYNAMICS',myThid)
433    #endif
434    C     compute temperature and (virtual) salt flux at the
435    C     ice-front ocean interface
436           CALL TIMER_START('ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
437         &       myThid)
438           CALL ICEFRONT_THERMODYNAMICS( myTime, myIter, myThid )
439           CALL TIMER_STOP( 'ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
440         &      myThid)
441          ENDIF
442    #endif /* ALLOW_ICEFRONT */
443    
444    #ifdef ALLOW_SALT_PLUME
445          IF ( useSALT_PLUME ) THEN
446    Catn: exchanging saltPlumeFlux:
447              CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
448          ENDIF
449    #endif /* ALLOW_SALT_PLUME */
450    
451    C--   Freeze water at the surface
452          IF ( allowFreezing ) THEN
453    #ifdef ALLOW_AUTODIFF_TAMC
454    CADJ STORE theta = comlev1, key = ikey_dynamics,
455    CADJ &     kind = isbyte
456    #endif
457            CALL FREEZE_SURFACE( myTime, myIter, myThid )
458          ENDIF
459    
460    #ifdef ALLOW_OCN_COMPON_INTERF
461    C--    Apply imported data (from coupled interface) to forcing fields
462    C jmc: do not know precisely where to put this call (bf or af thSIce ?)
463          IF ( useCoupler ) THEN
464             CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
465          ENDIF
466    #endif /* ALLOW_OCN_COMPON_INTERF */
467    
468          iMin = 1-OLx
469          iMax = sNx+OLx
470          jMin = 1-OLy
471          jMax = sNy+OLy
472    
473    C---  Determines forcing terms based on external fields
474    C     relaxation terms, etc.
475    #ifdef ALLOW_DEBUG
476          IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
477    #endif
478    #ifdef ALLOW_AUTODIFF
479    CADJ STORE salt, theta = comlev1, key = ikey_dynamics,
480    CADJ &     kind = isbyte
481    #else  /* ALLOW_AUTODIFF */
482    C--   if fluid is not water, by-pass surfaceForcing, find_rho, gmredi
483    C     and all vertical mixing schemes, but keep OBCS_CALC
484          IF ( fluidIsWater ) THEN
485    #endif /* ALLOW_AUTODIFF */
486            CALL EXTERNAL_FORCING_SURF(
487         I             iMin, iMax, jMin, jMax,
488         I             myTime, myIter, myThid )
489    
490  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
491  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
492  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
493  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
494        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
   
495  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
496  C--    HPF directive to help TAMC  C--   HPF directive to help TAMC
497  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$ INDEPENDENT
 CHPF$&                  ,utrans,vtrans,xA,yA  
 CHPF$&                  )  
498  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
499         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
500    
501  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 143  C     These inital values do not alter t Line 516  C     These inital values do not alter t
516  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
517  C     point numbers. This prevents spurious hardware signals due to  C     point numbers. This prevents spurious hardware signals due to
518  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  
   
519          DO k=1,Nr          DO k=1,Nr
520           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
521            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
522  C This is currently also used by IVDC and Diagnostics  C This is currently used by GMRedi, IVDC, MXL-depth  and Diagnostics
523             sigmaX(i,j,k) = 0. _d 0             sigmaX(i,j,k) = 0. _d 0
524             sigmaY(i,j,k) = 0. _d 0             sigmaY(i,j,k) = 0. _d 0
525             sigmaR(i,j,k) = 0. _d 0             sigmaR(i,j,k) = 0. _d 0
526  #ifdef ALLOW_AUTODIFF_TAMC            ENDDO
527             ENDDO
528            ENDDO
529    
530    #ifdef ALLOW_AUTODIFF
531            DO j=1-OLy,sNy+OLy
532             DO i=1-OLx,sNx+OLx
533              rhoKm1 (i,j)   = 0. _d 0
534              rhoKp1 (i,j)   = 0. _d 0
535             ENDDO
536            ENDDO
537  cph all the following init. are necessary for TAF  cph all the following init. are necessary for TAF
538  cph although some of these are re-initialised later.  cph although some of these are re-initialised later.
539            DO k=1,Nr
540             DO j=1-OLy,sNy+OLy
541              DO i=1-OLx,sNx+OLx
542               rhoInSitu(i,j,k,bi,bj) = 0.
543    # ifdef ALLOW_GGL90
544               GGL90viscArU(i,j,k,bi,bj)  = 0. _d 0
545               GGL90viscArV(i,j,k,bi,bj)  = 0. _d 0
546               GGL90diffKr(i,j,k,bi,bj)  = 0. _d 0
547    # endif /* ALLOW_GGL90 */
548    # ifdef ALLOW_SALT_PLUME
549    #  ifdef SALT_PLUME_VOLUME
550               SPforcingS(i,j,k,bi,bj) = 0. _d 0
551               SPforcingT(i,j,k,bi,bj) = 0. _d 0
552    #  endif
553    # endif /* ALLOW_SALT_PLUME */
554              ENDDO
555             ENDDO
556            ENDDO
557    #ifdef ALLOW_OFFLINE
558           IF ( calcConvect ) THEN
559    #endif
560            DO k=1,Nr
561             DO j=1-OLy,sNy+OLy
562              DO i=1-OLx,sNx+OLx
563             IVDConvCount(i,j,k,bi,bj) = 0.             IVDConvCount(i,j,k,bi,bj) = 0.
564              ENDDO
565             ENDDO
566            ENDDO
567    #ifdef ALLOW_OFFLINE
568           ENDIF
569           IF ( calcGMRedi ) THEN
570    #endif
571  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
572            DO k=1,Nr
573             DO j=1-OLy,sNy+OLy
574              DO i=1-OLx,sNx+OLx
575             Kwx(i,j,k,bi,bj)  = 0. _d 0             Kwx(i,j,k,bi,bj)  = 0. _d 0
576             Kwy(i,j,k,bi,bj)  = 0. _d 0             Kwy(i,j,k,bi,bj)  = 0. _d 0
577             Kwz(i,j,k,bi,bj)  = 0. _d 0             Kwz(i,j,k,bi,bj)  = 0. _d 0
# Line 181  cph although some of these are re-initia Line 590  cph although some of these are re-initia
590  #  ifdef GM_VISBECK_VARIABLE_K  #  ifdef GM_VISBECK_VARIABLE_K
591             VisbeckK(i,j,bi,bj)   = 0. _d 0             VisbeckK(i,j,bi,bj)   = 0. _d 0
592  #  endif  #  endif
593              ENDDO
594             ENDDO
595            ENDDO
596  # endif /* ALLOW_GMREDI */  # endif /* ALLOW_GMREDI */
597  #endif /* ALLOW_AUTODIFF_TAMC */  #ifdef ALLOW_OFFLINE
598           ENDIF
599           IF ( calcKPP ) THEN
600    #endif
601    # ifdef ALLOW_KPP
602            DO k=1,Nr
603             DO j=1-OLy,sNy+OLy
604              DO i=1-OLx,sNx+OLx
605               KPPdiffKzS(i,j,k,bi,bj)  = 0. _d 0
606               KPPdiffKzT(i,j,k,bi,bj)  = 0. _d 0
607            ENDDO            ENDDO
608           ENDDO           ENDDO
609          ENDDO          ENDDO
610    # endif /* ALLOW_KPP */
611          iMin = 1-OLx  #ifdef ALLOW_OFFLINE
612          iMax = sNx+OLx         ENDIF
613          jMin = 1-OLy  #endif
614          jMax = sNy+OLy  #endif /* ALLOW_AUTODIFF */
615    
616  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
617  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
618  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     kind = isbyte
619  CADJ STORE totphihyd  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
620  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     kind = isbyte
621  #ifdef ALLOW_KPP  CADJ STORE totphihyd(:,:,:,bi,bj)
622  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
623  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     kind = isbyte
624  #endif  # ifdef ALLOW_KPP
625    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
626    CADJ &     kind = isbyte
627    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
628    CADJ &     kind = isbyte
629    # endif
630    # ifdef ALLOW_SALT_PLUME
631    CADJ STORE saltplumedepth(:,:,bi,bj) = comlev1_bibj, key=itdkey,
632    CADJ &     kind = isbyte
633    CADJ STORE saltplumeflux(:,:,bi,bj) = comlev1_bibj, key=itdkey,
634    CADJ &     kind = isbyte
635    # endif
636  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
637    
638    C--   Always compute density (stored in common block) here; even when it is not
639    C     needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
640    #ifdef ALLOW_DEBUG
641            IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D (xNr)',myThid)
642    #endif
643    #ifdef ALLOW_AUTODIFF
644            IF ( fluidIsWater ) THEN
645    #endif /* ALLOW_AUTODIFF */
646    #ifdef ALLOW_DOWN_SLOPE
647             IF ( useDOWN_SLOPE ) THEN
648               DO k=1,Nr
649                CALL DWNSLP_CALC_RHO(
650         I                  theta, salt,
651         O                  rhoInSitu(1-OLx,1-OLy,k,bi,bj),
652         I                  k, bi, bj, myTime, myIter, myThid )
653               ENDDO
654             ENDIF
655    #endif /* ALLOW_DOWN_SLOPE */
656    #ifdef ALLOW_BBL
657             IF ( useBBL ) THEN
658    C     pkg/bbl requires in-situ bbl density for depths equal to and deeper than the bbl.
659    C     To reduce computation and storage requirement, these densities are stored in the
660    C     dry grid boxes of rhoInSitu.  See BBL_CALC_RHO for details.
661               DO k=Nr,1,-1
662                CALL BBL_CALC_RHO(
663         I                  theta, salt,
664         O                  rhoInSitu,
665         I                  k, bi, bj, myTime, myIter, myThid )
666    
667               ENDDO
668             ENDIF
669    #endif /* ALLOW_BBL */
670             IF ( .NOT. ( useDOWN_SLOPE .OR. useBBL ) ) THEN
671               DO k=1,Nr
672                CALL FIND_RHO_2D(
673         I                iMin, iMax, jMin, jMax, k,
674         I                theta(1-OLx,1-OLy,k,bi,bj),
675         I                salt (1-OLx,1-OLy,k,bi,bj),
676         O                rhoInSitu(1-OLx,1-OLy,k,bi,bj),
677         I                k, bi, bj, myThid )
678               ENDDO
679             ENDIF
680    #ifdef ALLOW_AUTODIFF
681            ELSE
682    C-        fluid is not water:
683              DO k=1,Nr
684               DO j=1-OLy,sNy+OLy
685                DO i=1-OLx,sNx+OLx
686                  rhoInSitu(i,j,k,bi,bj) = 0.
687                ENDDO
688               ENDDO
689              ENDDO
690            ENDIF
691    #endif /* ALLOW_AUTODIFF */
692    
693  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
694          IF ( debugLevel .GE. debLevB )          IF (debugMode) CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
      &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)  
695  #endif  #endif
696    
697  C--     Start of diagnostic loop  C--     Start of diagnostic loop
# Line 214  C--     Start of diagnostic loop Line 700  C--     Start of diagnostic loop
700  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
701  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?
702  C? Do we still need this?  C? Do we still need this?
703  cph kkey formula corrected.  cph kkey formula corrected.
704  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhoK, rhoKm1, in the case useGMREDI.
705           kkey = (itdkey-1)*Nr + k            kkey = (itdkey-1)*Nr + k
706  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
707    
708    c#ifdef ALLOW_AUTODIFF_TAMC
709    cCADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey,
710    cCADJ &     kind = isbyte
711    cCADJ STORE salt(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey,
712    cCADJ &     kind = isbyte
713    c#endif /* ALLOW_AUTODIFF_TAMC */
714    
715  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
716  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
717  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN            IF ( calcGMRedi .OR. (k.GT.1 .AND. calcConvect)
718            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN       &         .OR. usePP81 .OR. useMY82 .OR. useGGL90
719  #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 )  
   
720              IF (k.GT.1) THEN              IF (k.GT.1) THEN
721  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
722  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,
723  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ &     kind = isbyte
724    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey,
725    CADJ &     kind = isbyte
726    CADJ STORE rhokm1 (bi,bj)       = comlev1_bibj_k, key=kkey,
727    CADJ &     kind = isbyte
728  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
729               CALL FIND_RHO(               CALL FIND_RHO_2D(
730       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,       I                 iMin, iMax, jMin, jMax, k,
731       I        theta, salt,       I                 theta(1-OLx,1-OLy,k-1,bi,bj),
732       O        rhoKm1,       I                 salt (1-OLx,1-OLy,k-1,bi,bj),
733       I        myThid )       O                 rhoKm1,
734         I                 k-1, bi, bj, myThid )
735              ENDIF              ENDIF
736  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
737              IF ( debugLevel .GE. debLevB )              IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
      &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)  
738  #endif  #endif
739    cph Avoid variable aliasing for adjoint !!!
740                DO j=jMin,jMax
741                 DO i=iMin,iMax
742                  rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
743                 ENDDO
744                ENDDO
745              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
746       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
747       I             rhoK, rhoKm1, rhoK,       I             rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
748       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
749       I             myThid )       I             myThid )
750    #ifdef ALLOW_AUTODIFF
751    #ifdef GMREDI_WITH_STABLE_ADJOINT
752    cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
753    cgf -> cuts adjoint dependency from slope to state
754                CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
755                CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
756                CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
757    #endif
758    #endif /* ALLOW_AUTODIFF */
759            ENDIF            ENDIF
760    
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
761  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
762  c ==> should use sigmaR !!!            IF (k.GT.1 .AND. calcConvect) THEN
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
763  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
764              IF ( debugLevel .GE. debLevB )              IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid)
      &       CALL DEBUG_CALL('CALC_IVDC',myThid)  
765  #endif  #endif
766              CALL CALC_IVDC(              CALL CALC_IVDC(
767       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
768       I        rhoKm1, rhoK,       I        sigmaR,
769       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
770            ENDIF            ENDIF
771    
772    #ifdef ALLOW_DIAGNOSTICS
773              IF ( doDiagsRho.GE.4 ) THEN
774                CALL DIAGS_RHO_L( doDiagsRho, k, bi, bj,
775         I                        rhoInSitu(1-OLx,1-OLy,1,bi,bj),
776         I                        rhoKm1, wVel,
777         I                        myTime, myIter, myThid )
778              ENDIF
779    #endif
780    
781  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
782          ENDDO          ENDDO
783    
784  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
785  cph avoids recomputation of integrate_for_w  CADJ STORE IVDConvCount(:,:,:,bi,bj)
786  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey,
787    CADJ &     kind = isbyte
788    #endif
789    
790    C--     Diagnose Mixed Layer Depth:
791            IF ( calcGMRedi .OR. MOD(doDiagsRho,2).EQ.1 ) THEN
792              CALL CALC_OCE_MXLAYER(
793         I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
794         I              bi, bj, myTime, myIter, myThid )
795            ENDIF
796    
797    #ifdef ALLOW_SALT_PLUME
798            IF ( useSALT_PLUME ) THEN
799              CALL SALT_PLUME_CALC_DEPTH(
800         I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
801         I              bi, bj, myTime, myIter, myThid )
802    #ifdef SALT_PLUME_VOLUME
803              CALL SALT_PLUME_VOLFRAC(
804         I              bi, bj, myTime, myIter, myThid )
805    C-- get forcings for kpp
806    C-- note: temprray is just a dummy var here, not used anywhere.
807    C         so didnt need initialization, all zero on output.
808              CALL SALT_PLUME_APPLY(
809         I              1, bi, bj, recip_hFacC(1-OLx,1-OLy,1,bi,bj),
810         I              theta, 0,
811         I              myTime, myIter, myThid )
812    C     U              temparray,
813              CALL SALT_PLUME_APPLY(
814         I              2, bi, bj, recip_hFacC(1-OLx,1-OLy,1,bi,bj),
815         I              salt, 0,
816         I              myTime, myIter, myThid )
817    C     U              temparray,
818    C-- need to call this S/R from here to apply just before kpp
819              CALL SALT_PLUME_FORCING_SURF(
820         I              bi, bj, iMin, iMax, jMin, jMax,
821         I              myTime, myIter, myThid )
822    #endif /* SALT_PLUME_VOLUME */
823            ENDIF
824    #endif /* ALLOW_SALT_PLUME */
825    
826    #ifdef ALLOW_DIAGNOSTICS
827            IF ( MOD(doDiagsRho,4).GE.2 ) THEN
828              CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,
829         &         2, bi, bj, myThid)
830            ENDIF
831    #endif /* ALLOW_DIAGNOSTICS */
832    
833    C--    This is where EXTERNAL_FORCING_SURF(bi,bj) used to be called;
834    C      now called earlier, before bi,bj loop.
835    
836    #ifdef ALLOW_AUTODIFF_TAMC
837    cph needed for KPP
838    CADJ STORE surfaceForcingU(:,:,bi,bj)
839    CADJ &     = comlev1_bibj, key=itdkey,
840    CADJ &     kind = isbyte
841    CADJ STORE surfaceForcingV(:,:,bi,bj)
842    CADJ &     = comlev1_bibj, key=itdkey,
843    CADJ &     kind = isbyte
844    CADJ STORE surfaceForcingS(:,:,bi,bj)
845    CADJ &     = comlev1_bibj, key=itdkey,
846    CADJ &     kind = isbyte
847    CADJ STORE surfaceForcingT(:,:,bi,bj)
848    CADJ &     = comlev1_bibj, key=itdkey,
849    CADJ &     kind = isbyte
850    CADJ STORE surfaceForcingTice(:,:,bi,bj)
851    CADJ &     = comlev1_bibj, key=itdkey,
852    CADJ &     kind = isbyte
853  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
854    
855  #ifdef  ALLOW_OBCS  #ifdef  ALLOW_KPP
856  C--     Calculate future values on open boundaries  C--     Compute KPP mixing coefficients
857          IF (useOBCS) THEN          IF ( calcKPP ) THEN
858  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
859            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid)
      &     CALL DEBUG_CALL('OBCS_CALC',myThid)  
860  #endif  #endif
861            CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,            CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
862       I            uVel, vVel, wVel, theta, salt,            CALL KPP_CALC(
863       I            myThid )       I                  bi, bj, myTime, myIter, myThid )
864              CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
865    #if (defined ALLOW_AUTODIFF) && !(defined ALLOW_OFFLINE)
866            ELSE
867              CALL KPP_CALC_DUMMY(
868         I                  bi, bj, myTime, myIter, myThid )
869    #endif /* ALLOW_AUTODIFF and not ALLOW_OFFLINE */
870          ENDIF          ENDIF
871  #endif  /* ALLOW_OBCS */  #endif  /* ALLOW_KPP */
872    
873  #ifndef ALLOW_AUTODIFF_TAMC  #ifdef  ALLOW_PP81
874          IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN  C--     Compute PP81 mixing coefficients
875  #endif          IF (usePP81) THEN
 C--     Determines forcing terms based on external fields  
 C       relaxation terms, etc.  
876  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
877          IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid)
      &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)  
878  #endif  #endif
879           CALL EXTERNAL_FORCING_SURF(            CALL PP81_CALC(
880       I             bi, bj, iMin, iMax, jMin, jMax,       I                     bi, bj, sigmaR, myTime, myIter, myThid )
      I             myTime, myIter, myThid )  
 #ifndef ALLOW_AUTODIFF_TAMC  
881          ENDIF          ENDIF
882    #endif /* ALLOW_PP81 */
883    
884    #ifdef  ALLOW_MY82
885    C--     Compute MY82 mixing coefficients
886            IF (useMY82) THEN
887    #ifdef ALLOW_DEBUG
888              IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid)
889  #endif  #endif
890              CALL MY82_CALC(
891         I                     bi, bj, sigmaR, myTime, myIter, myThid )
892            ENDIF
893    #endif /* ALLOW_MY82 */
894    
895    #ifdef  ALLOW_GGL90
896  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
897  cph needed for KPP  CADJ STORE GGL90TKE (:,:,:,bi,bj) = comlev1_bibj, key=itdkey,
898  CADJ STORE surfacetendencyU(:,:,bi,bj)  CADJ &     kind = isbyte
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE surfacetendencyV(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE surfacetendencyS(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE surfacetendencyT(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 # ifdef ALLOW_SEAICE  
 CADJ STORE surfacetendencyTice(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 # endif  
899  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
900    C--     Compute GGL90 mixing coefficients
901            IF (useGGL90) THEN
902    #ifdef ALLOW_DEBUG
903              IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
904    #endif
905              CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
906              CALL GGL90_CALC(
907         I                     bi, bj, sigmaR, myTime, myIter, myThid )
908              CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
909            ENDIF
910    #endif /* ALLOW_GGL90 */
911    
912  #ifdef  ALLOW_GMREDI  #ifdef ALLOW_TIMEAVE
913            IF ( taveFreq.GT. 0. _d 0 ) THEN
914              CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
915            ENDIF
916            IF ( taveFreq.GT.0. .AND. calcConvect ) THEN
917              CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
918         I                           Nr, deltaTClock, bi, bj, myThid)
919            ENDIF
920    #endif /* ALLOW_TIMEAVE */
921    
922    #ifdef ALLOW_GMREDI
923  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
924    # ifndef GM_EXCLUDE_CLIPPING
925  cph storing here is needed only for one GMREDI_OPTIONS:  cph storing here is needed only for one GMREDI_OPTIONS:
926  cph define GM_BOLUS_ADVEC  cph define GM_BOLUS_ADVEC
927  cph but I've avoided the #ifdef for now, in case more things change  cph keep it although TAF says you dont need to.
928  CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  cph but I have avoided the #ifdef for now, in case more things change
929  CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey,
930  CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     kind = isbyte
931    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey,
932    CADJ &     kind = isbyte
933    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey,
934    CADJ &     kind = isbyte
935    # endif
936  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
937    
938  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
939          IF (useGMRedi) THEN          IF ( calcGMRedi ) THEN
940  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
941            IF ( debugLevel .GE. debLevB )            IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
      &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)  
942  #endif  #endif
943            CALL GMREDI_CALC_TENSOR(            CALL GMREDI_CALC_TENSOR(
944       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
945       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
946       I             myThid )       I             bi, bj, myTime, myIter, myThid )
947  #ifdef ALLOW_AUTODIFF_TAMC  #if (defined ALLOW_AUTODIFF) && !(defined ALLOW_OFFLINE)
948          ELSE          ELSE
949            CALL GMREDI_CALC_TENSOR_DUMMY(            CALL GMREDI_CALC_TENSOR_DUMMY(
950       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
951       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
952       I             myThid )       I             bi, bj, myTime, myIter, myThid )
953  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF and not ALLOW_OFFLINE */
954          ENDIF          ENDIF
955    #endif /* ALLOW_GMREDI */
956    
957  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_DOWN_SLOPE
958  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte          IF ( useDOWN_SLOPE ) THEN
959  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte  C--     Calculate Downsloping Flow for Down_Slope parameterization
960  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte           IF ( usingPCoords ) THEN
961  #endif /* ALLOW_AUTODIFF_TAMC */            CALL DWNSLP_CALC_FLOW(
962         I                bi, bj, kSurfC, rhoInSitu,
963         I                myTime, myIter, myThid )
964             ELSE
965              CALL DWNSLP_CALC_FLOW(
966         I                bi, bj, kLowC, rhoInSitu,
967         I                myTime, myIter, myThid )
968             ENDIF
969            ENDIF
970    #endif /* ALLOW_DOWN_SLOPE */
971    
972  #endif  /* ALLOW_GMREDI */  C--   end bi,bj loops.
973           ENDDO
974          ENDDO
975    
976  #ifdef  ALLOW_KPP  #ifndef ALLOW_AUTODIFF
977  C--     Compute KPP mixing coefficients  C---  if fluid Is Water: end
978          IF (useKPP) THEN        ENDIF
 #ifdef ALLOW_DEBUG  
           IF ( debugLevel .GE. debLevB )  
      &     CALL DEBUG_CALL('KPP_CALC',myThid)  
979  #endif  #endif
           CALL KPP_CALC(  
      I                  bi, bj, myTime, myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
980    
981  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_BBL
982  CADJ STORE KPPghat   (:,:,:,bi,bj)        IF ( useBBL ) THEN
983  CADJ &   , KPPfrac   (:,:  ,bi,bj)         CALL BBL_CALC_RHS(
984  CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte       I                          myTime, myIter, myThid )
985  #endif /* ALLOW_AUTODIFF_TAMC */        ENDIF
986    #endif /* ALLOW_BBL */
987    
988    #ifdef ALLOW_MYPACKAGE
989          IF ( useMYPACKAGE ) THEN
990           CALL MYPACKAGE_CALC_RHS(
991         I                          myTime, myIter, myThid )
992          ENDIF
993    #endif /* ALLOW_MYPACKAGE */
994    
995    #ifdef ALLOW_GMREDI
996          IF ( calcGMRedi ) THEN
997            CALL GMREDI_DO_EXCH( myTime, myIter, myThid )
998          ENDIF
999    #endif /* ALLOW_GMREDI */
1000    
1001  #endif  /* ALLOW_KPP */  #ifdef ALLOW_KPP
1002          IF ( calcKPP ) THEN
1003  #ifdef ALLOW_AUTODIFF_TAMC          CALL KPP_DO_EXCH( myThid )
1004  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte        ENDIF
1005  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  #endif /* ALLOW_KPP */
1006  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  
1007  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  #ifdef ALLOW_DIAGNOSTICS
1008  #ifdef ALLOW_PASSIVE_TRACER        IF ( fluidIsWater .AND. useDiagnostics ) THEN
1009  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte          CALL DIAGS_RHO_G(
1010  #endif       I                    rhoInSitu, uVel, vVel, wVel,
1011  #ifdef ALLOW_PTRACERS       I                    myTime, myIter, myThid )
1012  cph-- moved to forward_step to avoid key computation        ENDIF
1013  cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,        IF ( useDiagnostics ) THEN
1014  cphCADJ &                              key=itdkey, byte=isbyte          CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
1015          ENDIF
1016          IF ( calcConvect .AND. useDiagnostics ) THEN
1017            CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
1018         &                               0, Nr, 0, 1, 1, myThid )
1019          ENDIF
1020    #ifdef ALLOW_SALT_PLUME
1021          IF ( useDiagnostics )
1022         &      CALL SALT_PLUME_DIAGNOSTICS_FILL(bi,bj,myThid)
1023    #endif
1024  #endif  #endif
 #endif /* ALLOW_AUTODIFF_TAMC */  
1025    
1026  C--   end bi,bj loops.  #ifdef ALLOW_ECCO
1027         ENDDO        CALL ECCO_PHYS( myThid )
1028        ENDDO  #endif
1029    
1030  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
1031           IF ( debugLevel .GE. debLevB )        IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
      &    CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)  
1032  #endif  #endif
1033    
1034        RETURN        RETURN
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.136
  ViewVC Help
Powered by ViewVC 1.1.22