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