/[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.16 by jmc, Tue Dec 14 04:36:33 2004 UTC revision 1.74 by jmc, Fri Oct 3 02:26:49 2008 UTC
# Line 11  C $Name$ Line 11  C $Name$
11  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
12  #  include "KPP_OPTIONS.h"  #  include "KPP_OPTIONS.h"
13  # endif  # endif
14    # ifdef ALLOW_SEAICE
15    #  include "SEAICE_OPTIONS.h"
16    # endif
17  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
18    
19  CBOP  CBOP
# Line 19  C     !INTERFACE: Line 22  C     !INTERFACE:
22        SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)        SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
23  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
24  C     *==========================================================*  C     *==========================================================*
25  C     | SUBROUTINE DO_OCEANIC_PHYS                                  C     | SUBROUTINE DO_OCEANIC_PHYS
26  C     | o Controlling routine for oceanic physics and  C     | o Controlling routine for oceanic physics and
27  C     |   parameterization  C     |   parameterization
28  C     *==========================================================*  C     *==========================================================*
29  C     | o originally, part of S/R thermodynamics  C     | o originally, part of S/R thermodynamics
# Line 33  C     == Global variables === Line 36  C     == Global variables ===
36  #include "SIZE.h"  #include "SIZE.h"
37  #include "EEPARAMS.h"  #include "EEPARAMS.h"
38  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
39  #include "GRID.h"  #include "GRID.h"
40  c #include "GAD.h"  #include "DYNVARS.h"
41  c #ifdef ALLOW_PTRACERS  #ifdef ALLOW_TIMEAVE
42  c #include "PTRACERS_SIZE.h"  #include "TIMEAVE_STATV.h"
43  c #include "PTRACERS.h"  #endif
44  c #endif  #if defined (ALLOW_BALANCE_FLUXES) && !(defined ALLOW_AUTODIFF_TAMC)
45  c #ifdef ALLOW_TIMEAVE  #include "FFIELDS.h"
46  c #include "TIMEAVE_STATV.h"  #endif
 c #endif  
47    
48  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
49  # include "tamc.h"  # include "tamc.h"
50  # include "tamc_keys.h"  # include "tamc_keys.h"
51  # include "FFIELDS.h"  # include "FFIELDS.h"
52    # include "SURFACE.h"
53  # include "EOS.h"  # include "EOS.h"
54  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
55  #  include "KPP.h"  #  include "KPP.h"
# Line 58  c #endif Line 60  c #endif
60  # ifdef ALLOW_EBM  # ifdef ALLOW_EBM
61  #  include "EBM.h"  #  include "EBM.h"
62  # endif  # endif
63  # ifdef EXACT_CONSERV  # ifdef ALLOW_EXF
64  #  include "SURFACE.h"  #  include "ctrl.h"
65    #  include "EXF_FIELDS.h"
66    #  ifdef ALLOW_BULKFORMULAE
67    #   include "EXF_CONSTANTS.h"
68    #  endif
69    # endif
70    # ifdef ALLOW_SEAICE
71    #  include "SEAICE.h"
72  # endif  # endif
73  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
74    
75  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
76  C     == Routine arguments ==  C     == Routine arguments ==
77  C     myTime - Current time in simulation  C     myTime :: Current time in simulation
78  C     myIter - Current iteration number in simulation  C     myIter :: Current iteration number in simulation
79  C     myThid - Thread number for this instance of the routine.  C     myThid :: Thread number for this instance of the routine.
80        _RL myTime        _RL myTime
81        INTEGER myIter        INTEGER myIter
82        INTEGER myThid        INTEGER myThid
83    
84  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
85  C     == Local variables  C     == Local variables
86  C     rhoK, rhoKM1   - Density at current level, and level above  C     rhoK, rhoKm1  :: Density at current level, and level above
87  C     useVariableK   = T when vertical diffusion is not constant  C     iMin, iMax    :: Ranges and sub-block indices on which calculations
 C     iMin, iMax     - Ranges and sub-block indices on which calculations  
88  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
89  C     bi, bj  C     bi, bj        :: tile indices
90  C     k, kup,        - Index for layer above and below. kup and kDown  C     i,j,k         :: loop indices
91  C     kDown, km1       are switched with layer to be the appropriate        _RL rhoKp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92  C                      index into fVerTerm.        _RL rhoKm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
93        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
94        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
95        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL kp1Msk  
       LOGICAL useVariableK  
96        INTEGER iMin, iMax        INTEGER iMin, iMax
97        INTEGER jMin, jMax        INTEGER jMin, jMax
98        INTEGER bi, bj        INTEGER bi, bj
99        INTEGER i, j        INTEGER i, j, k
100        INTEGER k, km1, kup, kDown        INTEGER doDiagsRho
101        INTEGER iTracer, ip  #ifdef ALLOW_DIAGNOSTICS
102          LOGICAL  DIAGNOSTICS_IS_ON
103          EXTERNAL DIAGNOSTICS_IS_ON
104    #endif /* ALLOW_DIAGNOSTICS */
105    
106  CEOP  CEOP
107    
# Line 104  C--   dummy statement to end declaration Line 111  C--   dummy statement to end declaration
111  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
112    
113  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
114        IF ( debugLevel .GE. debLevB )        IF ( debugLevel .GE. debLevB )
115       &    CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)       &     CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
116    #endif
117    
118          doDiagsRho = 0
119    #ifdef ALLOW_DIAGNOSTICS
120          IF ( useDiagnostics .AND. fluidIsWater ) THEN
121            IF ( DIAGNOSTICS_IS_ON('WRHOMASS',myThid) )
122         &       doDiagsRho = doDiagsRho + 1
123            IF ( DIAGNOSTICS_IS_ON('DRHODR  ',myThid) )
124         &       doDiagsRho = doDiagsRho + 2
125            IF ( DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) )
126         &       doDiagsRho = doDiagsRho + 4
127          ENDIF
128    #endif /* ALLOW_DIAGNOSTICS */
129    
130    
131    #ifdef ALLOW_SEAICE
132          IF ( useSEAICE ) THEN
133    # ifdef ALLOW_AUTODIFF_TAMC
134    cph-adj-test(
135    CADJ STORE area,empmr,qsw,theta   = comlev1, key = ikey_dynamics
136    cph-adj-test)
137    CADJ STORE atemp,aqh,precip    = comlev1, key = ikey_dynamics
138    CADJ STORE swdown,lwdown       = comlev1, key = ikey_dynamics
139    cph# ifdef EXF_READ_EVAP
140    CADJ STORE evap                = comlev1, key = ikey_dynamics
141    cph# endif
142    CADJ STORE uvel,vvel           = comlev1, key = ikey_dynamics
143    #  ifdef SEAICE_ALLOW_DYNAMICS
144    CADJ STORE uice                = comlev1, key = ikey_dynamics
145    CADJ STORE vice                = comlev1, key = ikey_dynamics
146    #   ifdef SEAICE_ALLOW_EVP
147    CADJ STORE seaice_sigma1       = comlev1, key = ikey_dynamics
148    CADJ STORE seaice_sigma2       = comlev1, key = ikey_dynamics
149    CADJ STORE seaice_sigma12      = comlev1, key = ikey_dynamics
150    #   endif
151    #  endif
152    #  ifdef SEAICE_SALINITY
153    CADJ STORE salt                = comlev1, key = ikey_dynamics
154    #  endif
155    #  ifdef ATMOSPHERIC_LOADING
156    CADJ STORE pload               = comlev1, key = ikey_dynamics
157    CADJ STORE siceload            = comlev1, key = ikey_dynamics
158    #  endif
159    #  ifdef NONLIN_FRSURF
160    CADJ STORE recip_hfacc         = comlev1, key = ikey_dynamics
161    #  endif
162    # endif
163    # ifdef ALLOW_DEBUG
164            IF ( debugLevel .GE. debLevB )
165         &    CALL DEBUG_CALL('SEAICE_MODEL',myThid)
166    # endif
167            CALL TIMER_START('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
168            CALL SEAICE_MODEL( myTime, myIter, myThid )
169            CALL TIMER_STOP ('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
170    # ifdef ALLOW_COST
171            CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
172    # endif
173          ENDIF
174    #endif /* ALLOW_SEAICE */
175    
176    #ifdef ALLOW_AUTODIFF_TAMC
177    CADJ STORE sst, sss           = comlev1, key = ikey_dynamics
178    CADJ STORE qsw                = comlev1, key = ikey_dynamics
179    # ifdef ALLOW_SEAICE
180    CADJ STORE area               = comlev1, key = ikey_dynamics
181    # endif
182  #endif  #endif
183    
184  #ifdef ALLOW_THSICE  #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
185        IF ( useThSIce .AND. fluidIsWater ) THEN        IF ( useThSIce .AND. fluidIsWater ) THEN
186  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
187          IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
# Line 122  C       and modify forcing terms includi Line 195  C       and modify forcing terms includi
195        ENDIF        ENDIF
196  #endif /* ALLOW_THSICE */  #endif /* ALLOW_THSICE */
197    
198    #ifdef ALLOW_SHELFICE
199          IF ( useShelfIce .AND. fluidIsWater ) THEN
200    #ifdef ALLOW_DEBUG
201            IF ( debugLevel .GE. debLevB )
202         &    CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
203    #endif
204    C     compute temperature and (virtual) salt flux at the
205    C     shelf-ice ocean interface
206           CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
207         &       myThid)
208           CALL SHELFICE_THERMODYNAMICS( myTime, myIter, myThid )
209           CALL TIMER_STOP( 'SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
210         &      myThid)
211          ENDIF
212    #endif /* ALLOW_SHELFICE */
213    
214  C--   Freeze water at the surface  C--   Freeze water at the surface
215  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
216  CADJ STORE theta = comlev1, key = ikey_dynamics  CADJ STORE theta = comlev1, key = ikey_dynamics
217  #endif  #endif
218        IF ( allowFreezing        IF ( allowFreezing ) THEN
      &                   .AND. .NOT. useSEAICE  
      &                   .AND. .NOT. useThSIce ) THEN  
219          CALL FREEZE_SURFACE(  myTime, myIter, myThid )          CALL FREEZE_SURFACE(  myTime, myIter, myThid )
220        ENDIF        ENDIF
221    
222  #ifdef COMPONENT_MODULE  #ifdef ALLOW_OCN_COMPON_INTERF
 # ifndef ALLOW_AIM  
223  C--    Apply imported data (from coupled interface) to forcing fields  C--    Apply imported data (from coupled interface) to forcing fields
224  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 ?)
225         IF ( useCoupler ) THEN        IF ( useCoupler ) THEN
226           CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )           CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
227         ENDIF        ENDIF
228  # endif  #endif /* ALLOW_OCN_COMPON_INTERF */
229  #endif /* COMPONENT_MODULE */  
230    #ifdef ALLOW_BALANCE_FLUXES
231    C     balance fluxes
232          IF ( balanceEmPmR )
233         &        CALL REMOVE_MEAN_RS( 1, EmPmR, maskH, maskH, rA, drF,
234         &        'EmPmR', myTime, myThid )
235          IF ( balanceQnet )
236         &        CALL REMOVE_MEAN_RS( 1, Qnet,  maskH, maskH, rA, drF,
237         &        'Qnet ', myTime, myThid )
238    #endif /* ALLOW_BALANCE_FLUXES */
239    
240  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
241  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 164  CHPF$ INDEPENDENT Line 259  CHPF$ INDEPENDENT
259            itdkey = (act1 + 1) + act2*max1            itdkey = (act1 + 1) + act2*max1
260       &                      + act3*max1*max2       &                      + act3*max1*max2
261       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
262    #else  /* ALLOW_AUTODIFF_TAMC */
263    C     if fluid is not water, by-pass find_rho, gmredi, surfaceForcing
264    C     and all vertical mixing schemes, but keep OBCS_CALC
265            IF ( fluidIsWater ) THEN
266  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
267    
268  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 172  C     just ensure that all memory refere Line 271  C     just ensure that all memory refere
271  C     point numbers. This prevents spurious hardware signals due to  C     point numbers. This prevents spurious hardware signals due to
272  C     uninitialised but inert locations.  C     uninitialised but inert locations.
273    
274    #ifdef ALLOW_AUTODIFF_TAMC
275          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
276           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
277            rhok   (i,j)   = 0. _d 0            rhoKm1 (i,j)   = 0. _d 0
278            rhoKM1 (i,j)   = 0. _d 0            rhoKp1 (i,j)   = 0. _d 0
279           ENDDO           ENDDO
280          ENDDO          ENDDO
281    #endif /* ALLOW_AUTODIFF_TAMC */
282    
283          DO k=1,Nr          DO k=1,Nr
284           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
285            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
286  C This is currently also used by IVDC and Diagnostics  C This is currently used by GMRedi, IVDC, MXL-depth  and Diagnostics
287             sigmaX(i,j,k) = 0. _d 0             sigmaX(i,j,k) = 0. _d 0
288             sigmaY(i,j,k) = 0. _d 0             sigmaY(i,j,k) = 0. _d 0
289             sigmaR(i,j,k) = 0. _d 0             sigmaR(i,j,k) = 0. _d 0
290  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
291  cph all the following init. are necessary for TAF  cph all the following init. are necessary for TAF
292  cph although some of these are re-initialised later.  cph although some of these are re-initialised later.
293    c          rhoInSitu(i,j,k,bi,bj) = 0.
294             IVDConvCount(i,j,k,bi,bj) = 0.             IVDConvCount(i,j,k,bi,bj) = 0.
295  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
296             Kwx(i,j,k,bi,bj)  = 0. _d 0             Kwx(i,j,k,bi,bj)  = 0. _d 0
# Line 210  cph although some of these are re-initia Line 312  cph although some of these are re-initia
312             VisbeckK(i,j,bi,bj)   = 0. _d 0             VisbeckK(i,j,bi,bj)   = 0. _d 0
313  #  endif  #  endif
314  # endif /* ALLOW_GMREDI */  # endif /* ALLOW_GMREDI */
315    # ifdef ALLOW_KPP
316               KPPdiffKzS(i,j,k,bi,bj)  = 0. _d 0
317               KPPdiffKzT(i,j,k,bi,bj)  = 0. _d 0
318    # endif /* ALLOW_KPP */
319  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
320            ENDDO            ENDDO
321           ENDDO           ENDDO
# Line 229  CADJ &     = comlev1_bibj, key=itdkey, b Line 335  CADJ &     = comlev1_bibj, key=itdkey, b
335  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
336  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
337  # endif  # endif
 # ifdef EXACT_CONSERV  
 CADJ STORE pmepr(:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte  
 # endif  
338  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
339    
340  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
341          IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
342       &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)       &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
343  #endif  #endif
344    
# Line 245  C--     Start of diagnostic loop Line 348  C--     Start of diagnostic loop
348  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
349  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?
350  C? Do we still need this?  C? Do we still need this?
351  cph kkey formula corrected.  cph kkey formula corrected.
352  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhoK, rhoKm1, in the case useGMREDI.
353           kkey = (itdkey-1)*Nr + k            kkey = (itdkey-1)*Nr + k
354  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
355    
356  C--       Calculate gradients of potential density for isoneutral  C--   Always compute density (stored in common block) here; even when it is not
357  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C     needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
 c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  
           IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN  
358  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
359              IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
360       &       CALL DEBUG_CALL('FIND_RHO',myThid)       &       CALL DEBUG_CALL('FIND_RHO_2D',myThid)
361  #endif  #endif
362  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
363              IF ( fluidIsWater ) THEN
364  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
365  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE salt(:,:,k,bi,bj)  = comlev1_bibj_k, key=kkey, byte=isbyte
366    #endif /* ALLOW_AUTODIFF_TAMC */
367    #ifdef ALLOW_DOWN_SLOPE
368              IF ( useDOWN_SLOPE ) THEN
369                CALL DWNSLP_CALC_RHO(
370         I                  theta, salt,
371         O                  rhoInSitu(1-OLx,1-OLy,k,bi,bj),
372         I                  k, bi, bj, myTime, myIter, myThid )
373              ELSE
374    #endif /* ALLOW_DOWN_SLOPE */
375                CALL FIND_RHO_2D(
376         I                iMin, iMax, jMin, jMax, k,
377         I                theta(1-OLx,1-OLy,k,bi,bj),
378         I                salt (1-OLx,1-OLy,k,bi,bj),
379         O                rhoInSitu(1-OLx,1-OLy,k,bi,bj),
380         I                k, bi, bj, myThid )
381    #ifdef ALLOW_DOWN_SLOPE
382              ENDIF
383    #endif /* ALLOW_DOWN_SLOPE */
384    #ifdef ALLOW_AUTODIFF_TAMC
385              ELSE
386    C-        fluid is not water:
387               DO j=1-OLy,sNy+OLy
388                DO i=1-OLx,sNx+OLx
389                  rhoInSitu(i,j,k,bi,bj) = 0.
390                ENDDO
391               ENDDO
392              ENDIF
393  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
394    
395    C--       Calculate gradients of potential density for isoneutral
396    C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
397              IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)
398         &         .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
399              IF (k.GT.1) THEN              IF (k.GT.1) THEN
400  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
401  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, byte=isbyte
402  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
403    CADJ STORE rhokm1 (bi,bj)       = comlev1_bibj_k, key=kkey, byte=isbyte
404  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
405               CALL FIND_RHO(               CALL FIND_RHO_2D(
406       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,       I                 iMin, iMax, jMin, jMax, k,
407       I        theta, salt,       I                 theta(1-OLx,1-OLy,k-1,bi,bj),
408       O        rhoKm1,       I                 salt (1-OLx,1-OLy,k-1,bi,bj),
409       I        myThid )       O                 rhoKm1,
410         I                 k-1, bi, bj, myThid )
411              ENDIF              ENDIF
412  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
413              IF ( debugLevel .GE. debLevB )              IF ( debugLevel .GE. debLevB )
414       &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)       &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)
415  #endif  #endif
416    cph Avoid variable aliasing for adjoint !!!
417                DO j=jMin,jMax
418                 DO i=iMin,iMax
419                  rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
420                 ENDDO
421                ENDDO
422              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
423       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
424       I             rhoK, rhoKm1, rhoK,       I             rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
425       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
426       I             myThid )       I             myThid )
           ENDIF  
   
427  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
428  ctest# ifndef GM_EXCLUDE_CLIPPING  #ifdef GMREDI_WITH_STABLE_ADJOINT
429  CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
430  ctest# endif  cgf -> cuts adjoint dependency from slope to state
431  CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte              CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
432                CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
433                CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
434    #endif
435  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
436              ENDIF
437    
438  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
439  c ==> should use sigmaR !!!  c ==> should use sigmaR !!!
440            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
441  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
442              IF ( debugLevel .GE. debLevB )              IF ( debugLevel .GE. debLevB )
443       &       CALL DEBUG_CALL('CALC_IVDC',myThid)       &       CALL DEBUG_CALL('CALC_IVDC',myThid)
444  #endif  #endif
445              CALL CALC_IVDC(              CALL CALC_IVDC(
446       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
447       I        rhoKm1, rhoK,       I        rhoKm1, rhoInSitu(1-OLx,1-OLy,k,bi,bj),
448       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
449            ENDIF            ENDIF
450    
451    #ifdef ALLOW_DIAGNOSTICS
452              IF ( MOD(doDiagsRho,2).EQ.1 ) THEN
453                CALL DIAGS_RHO_L( k, bi, bj,
454         I                        rhoInSitu(1-OLx,1-OLy,k,bi,bj),
455         I                        rhoKm1, wVel,
456         I                        myTime, myIter, myThid )
457              ENDIF
458    #endif
459    
460  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
461          ENDDO          ENDDO
462    
463    #ifdef ALLOW_AUTODIFF_TAMC
464    CADJ STORE IVDConvCount(:,:,:,bi,bj)
465    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
466    #endif
467    
468    C--     Diagnose Mixed Layer Depth:
469            IF ( useGMRedi .OR. doDiagsRho.GE.4 ) THEN
470              CALL CALC_OCE_MXLAYER(
471         I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
472         I              bi, bj, myTime, myIter, myThid )
473            ENDIF
474    
475    #ifdef ALLOW_SALT_PLUME
476            IF ( useSALT_PLUME ) THEN
477              CALL SALT_PLUME_CALC_DEPTH(
478         I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
479         I              bi, bj, myTime, myIter, myThid )
480            ENDIF
481    #endif /* ALLOW_SALT_PLUME */
482    
483  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
484          IF ( useDiagnostics .AND.          IF ( MOD(doDiagsRho,4).GE.2 ) THEN
      &       (useGMRedi .OR. ivdc_kappa.NE.0.) ) THEN  
485            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,            CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,
486       &         2, bi, bj, myThid)       &         2, bi, bj, myThid)
487          ENDIF          ENDIF
488  #endif  #endif /* ALLOW_DIAGNOSTICS */
489    
 #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 */  
   
 #ifndef ALLOW_AUTODIFF_TAMC  
         IF ( fluidIsWater ) THEN  
 #endif  
490  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
491  C       relaxation terms, etc.  C       relaxation terms, etc.
492  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
493          IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
494       &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)       &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
495  #endif  #endif
496           CALL EXTERNAL_FORCING_SURF(  #ifdef ALLOW_AUTODIFF_TAMC
497    CADJ STORE EmPmR(:,:,bi,bj)
498    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
499    # ifdef EXACT_CONSERV
500    CADJ STORE PmEpR(:,:,bi,bj)
501    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
502    # endif
503    # ifdef NONLIN_FRSURF
504    CADJ STORE hFac_surfC(:,:,bi,bj)
505    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
506    CADJ STORE recip_hFacC(:,:,:,bi,bj)
507    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
508    # endif
509    #endif
510            CALL EXTERNAL_FORCING_SURF(
511       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
512       I             myTime, myIter, myThid )       I             myTime, myIter, myThid )
513  #ifndef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
514          ENDIF  # ifdef EXACT_CONSERV
515    cph-test
516    cphCADJ STORE PmEpR(:,:,bi,bj)
517    cphCADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
518    # endif
519  #endif  #endif
520    
521  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 361  CADJ STORE surfaceForcingT(:,:,bi,bj) Line 530  CADJ STORE surfaceForcingT(:,:,bi,bj)
530  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
531  CADJ STORE surfaceForcingTice(:,:,bi,bj)  CADJ STORE surfaceForcingTice(:,:,bi,bj)
532  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
   
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #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 )  
533  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
         ENDIF  
   
 #endif  /* ALLOW_GMREDI */  
534    
535  #ifdef  ALLOW_KPP  #ifdef  ALLOW_KPP
536  C--     Compute KPP mixing coefficients  C--     Compute KPP mixing coefficients
537          IF (useKPP) THEN          IF (useKPP) THEN
538  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
539            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
540       &     CALL DEBUG_CALL('KPP_CALC',myThid)       &     CALL DEBUG_CALL('KPP_CALC',myThid)
541  #endif  #endif
542            CALL KPP_CALC(            CALL KPP_CALC(
543       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myIter, myThid )
544  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
545          ELSE          ELSE
546            CALL KPP_CALC_DUMMY(            CALL KPP_CALC_DUMMY(
547       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myIter, myThid )
548  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
549          ENDIF          ENDIF
550    
# Line 421  C--     Compute KPP mixing coefficients Line 554  C--     Compute KPP mixing coefficients
554  C--     Compute PP81 mixing coefficients  C--     Compute PP81 mixing coefficients
555          IF (usePP81) THEN          IF (usePP81) THEN
556  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
557            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
558       &     CALL DEBUG_CALL('PP81_CALC',myThid)       &     CALL DEBUG_CALL('PP81_CALC',myThid)
559  #endif  #endif
560            CALL PP81_CALC(            CALL PP81_CALC(
# Line 433  C--     Compute PP81 mixing coefficients Line 566  C--     Compute PP81 mixing coefficients
566  C--     Compute MY82 mixing coefficients  C--     Compute MY82 mixing coefficients
567          IF (useMY82) THEN          IF (useMY82) THEN
568  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
569            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
570       &     CALL DEBUG_CALL('MY82_CALC',myThid)       &     CALL DEBUG_CALL('MY82_CALC',myThid)
571  #endif  #endif
572            CALL MY82_CALC(            CALL MY82_CALC(
# Line 445  C--     Compute MY82 mixing coefficients Line 578  C--     Compute MY82 mixing coefficients
578  C--     Compute GGL90 mixing coefficients  C--     Compute GGL90 mixing coefficients
579          IF (useGGL90) THEN          IF (useGGL90) THEN
580  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
581            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
582       &     CALL DEBUG_CALL('GGL90_CALC',myThid)       &     CALL DEBUG_CALL('GGL90_CALC',myThid)
583  #endif  #endif
584            CALL GGL90_CALC(            CALL GGL90_CALC(
# Line 453  C--     Compute GGL90 mixing coefficient Line 586  C--     Compute GGL90 mixing coefficient
586          ENDIF          ENDIF
587  #endif /* ALLOW_GGL90 */  #endif /* ALLOW_GGL90 */
588    
589    #ifdef ALLOW_TIMEAVE
590            IF ( taveFreq.GT. 0. _d 0 ) THEN
591              CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
592            ENDIF
593            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
594              CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
595         I                           Nr, deltaTclock, bi, bj, myThid)
596            ENDIF
597    #endif /* ALLOW_TIMEAVE */
598    
599    #ifdef ALLOW_GMREDI
600    #ifdef ALLOW_AUTODIFF_TAMC
601    # ifndef GM_EXCLUDE_CLIPPING
602    cph storing here is needed only for one GMREDI_OPTIONS:
603    cph define GM_BOLUS_ADVEC
604    cph keep it although TAF says you dont need to.
605    cph but I've avoided the #ifdef for now, in case more things change
606    CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
607    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
608    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
609    # endif
610    #endif /* ALLOW_AUTODIFF_TAMC */
611    
612    C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
613            IF (useGMRedi) THEN
614    #ifdef ALLOW_DEBUG
615              IF ( debugLevel .GE. debLevB )
616         &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
617    #endif
618              CALL GMREDI_CALC_TENSOR(
619         I             iMin, iMax, jMin, jMax,
620         I             sigmaX, sigmaY, sigmaR,
621         I             bi, bj, myTime, myIter, myThid )
622    #ifdef ALLOW_AUTODIFF_TAMC
623            ELSE
624              CALL GMREDI_CALC_TENSOR_DUMMY(
625         I             iMin, iMax, jMin, jMax,
626         I             sigmaX, sigmaY, sigmaR,
627         I             bi, bj, myTime, myIter, myThid )
628    #endif /* ALLOW_AUTODIFF_TAMC */
629            ENDIF
630    #endif /* ALLOW_GMREDI */
631    
632    #ifdef ALLOW_DOWN_SLOPE
633            IF ( useDOWN_SLOPE ) THEN
634    C--     Calculate Downsloping Flow for Down_Slope parameterization
635             IF ( usingPCoords ) THEN
636              CALL DWNSLP_CALC_FLOW(
637         I                bi, bj, kSurfC, rhoInSitu,
638         I                myTime, myIter, myThid )
639             ELSE
640              CALL DWNSLP_CALC_FLOW(
641         I                bi, bj, kLowC, rhoInSitu,
642         I                myTime, myIter, myThid )
643             ENDIF
644            ENDIF
645    #endif /* ALLOW_DOWN_SLOPE */
646    
647    #ifndef ALLOW_AUTODIFF_TAMC
648    C---  if fluid Is Water: end
649            ENDIF
650    #endif
651    
652    #ifdef  ALLOW_OBCS
653    C--     Calculate future values on open boundaries
654            IF (useOBCS) THEN
655    #ifdef ALLOW_DEBUG
656              IF ( debugLevel .GE. debLevB )
657         &     CALL DEBUG_CALL('OBCS_CALC',myThid)
658    #endif
659              CALL OBCS_CALC( bi, bj, myTime+deltaTclock, myIter+1,
660         I            uVel, vVel, wVel, theta, salt,
661         I            myThid )
662            ENDIF
663    #endif  /* ALLOW_OBCS */
664    
665  C--   end bi,bj loops.  C--   end bi,bj loops.
666         ENDDO         ENDDO
667        ENDDO        ENDDO
668    
669    #ifdef  ALLOW_KPP
670          IF (useKPP) THEN
671            CALL KPP_DO_EXCH( myThid )
672          ENDIF
673    #endif  /* ALLOW_KPP */
674    
675    #ifdef ALLOW_DIAGNOSTICS
676          IF ( fluidIsWater .AND. useDiagnostics ) THEN
677            CALL DIAGS_RHO_G(
678         I                    rhoInSitu, uVel, vVel,
679         I                    myTime, myIter, myThid )
680            CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
681          ENDIF
682          IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN
683            CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
684         &                               0, Nr, 0, 1, 1, myThid )
685          ENDIF
686    #endif
687    
688  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
689           IF ( debugLevel .GE. debLevB )        IF ( debugLevel .GE. debLevB )
690       &    CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)       &     CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
691  #endif  #endif
692    
693        RETURN        RETURN
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.74
  ViewVC Help
Powered by ViewVC 1.1.22