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