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