/[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.2 by jmc, Tue Jul 6 21:12:51 2004 UTC revision 1.72 by heimbach, Mon Sep 22 21:24:08 2008 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "PACKAGES_CONFIG.h"  #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
 c #ifdef ALLOW_PTRACERS  
 c # include "PTRACERS_OPTIONS.h"  
 c #endif  
6    
7  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
8  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
# Line 14  c #endif Line 11  c #endif
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 22  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 36  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_PASSIVE_TRACER  #ifdef ALLOW_TIMEAVE
42  c #include "TR1.h"  #include "TIMEAVE_STATV.h"
43  c #endif  #endif
44  c #ifdef ALLOW_PTRACERS  #if defined (ALLOW_BALANCE_FLUXES) && !(defined ALLOW_AUTODIFF_TAMC)
45  c #include "PTRACERS.h"  #include "FFIELDS.h"
46  c #endif  #endif
 c #ifdef ALLOW_TIMEAVE  
 c #include "TIMEAVE_STATV.h"  
 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 63  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 ALLOW_EXF
64    #  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
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    
 #ifdef ALLOW_DEBUG  
          IF ( debugLevel .GE. debLevB )  
      &    CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)  
 #endif  
   
108  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
109  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
       ikey = 1  
110        itdkey = 1        itdkey = 1
111  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
112    
113    #ifdef ALLOW_DEBUG
114          IF ( debugLevel .GE. debLevB )
115         &     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
183    
184    #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
185          IF ( useThSIce .AND. fluidIsWater ) THEN
186    #ifdef ALLOW_DEBUG
187            IF ( debugLevel .GE. debLevB )
188         &    CALL DEBUG_CALL('THSICE_MAIN',myThid)
189    #endif
190    C--     Step forward Therm.Sea-Ice variables
191    C       and modify forcing terms including effects from ice
192            CALL TIMER_START('THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
193            CALL THSICE_MAIN( myTime, myIter, myThid )
194            CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
195          ENDIF
196    #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
215    #ifdef ALLOW_AUTODIFF_TAMC
216    CADJ STORE theta = comlev1, key = ikey_dynamics
217    #endif
218          IF ( allowFreezing ) THEN
219            CALL FREEZE_SURFACE(  myTime, myIter, myThid )
220          ENDIF
221    
222    #ifdef ALLOW_OCN_COMPON_INTERF
223    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 ?)
225          IF ( useCoupler ) THEN
226             CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
227          ENDIF
228    #endif /* ALLOW_OCN_COMPON_INTERF */
229    
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
242  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
243  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
244        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
   
245  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
246  C--    HPF directive to help TAMC  C--   HPF directive to help TAMC
247  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$ INDEPENDENT
 CHPF$&                  ,utrans,vtrans,xA,yA  
 CHPF$&                  )  
248  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
249         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
250    
251  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 138  CHPF$&                  ) Line 259  CHPF$&                  )
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 146  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 184  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 197  cph although some of these are re-initia Line 329  cph although some of these are re-initia
329  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
330  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
331  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
332  CADJ STORE totphihyd  CADJ STORE totphihyd(:,:,:,bi,bj)
333  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
334  #ifdef ALLOW_KPP  # ifdef ALLOW_KPP
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
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 216  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  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
364  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
365  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
366              CALL FIND_RHO(  #ifdef ALLOW_DOWN_SLOPE
367       I        bi, bj, iMin, iMax, jMin, jMax, k, k,            IF ( useDOWN_SLOPE ) THEN
368       I        theta, salt,              CALL DWNSLP_CALC_RHO(
369       O        rhoK,       I                  theta, salt,
370       I        myThid )       O                  rhoInSitu(1-OLx,1-OLy,k,bi,bj),
371         I                  k, bi, bj, myTime, myIter, myThid )
372              ELSE
373    #endif /* ALLOW_DOWN_SLOPE */
374                CALL FIND_RHO_2D(
375         I                iMin, iMax, jMin, jMax, k,
376         I                theta(1-OLx,1-OLy,k,bi,bj),
377         I                salt (1-OLx,1-OLy,k,bi,bj),
378         O                rhoInSitu(1-OLx,1-OLy,k,bi,bj),
379         I                k, bi, bj, myThid )
380    #ifdef ALLOW_DOWN_SLOPE
381              ENDIF
382    #endif /* ALLOW_DOWN_SLOPE */
383    
384    C--       Calculate gradients of potential density for isoneutral
385    C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
386              IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.)
387         &         .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
388              IF (k.GT.1) THEN              IF (k.GT.1) THEN
389  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
390  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
391  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
392    CADJ STORE rhokm1 (bi,bj)       = comlev1_bibj_k, key=kkey, byte=isbyte
393  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
394               CALL FIND_RHO(               CALL FIND_RHO_2D(
395       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,       I                 iMin, iMax, jMin, jMax, k,
396       I        theta, salt,       I                 theta(1-OLx,1-OLy,k-1,bi,bj),
397       O        rhoKm1,       I                 salt (1-OLx,1-OLy,k-1,bi,bj),
398       I        myThid )       O                 rhoKm1,
399         I                 k-1, bi, bj, myThid )
400              ENDIF              ENDIF
401  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
402              IF ( debugLevel .GE. debLevB )              IF ( debugLevel .GE. debLevB )
403       &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)       &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)
404  #endif  #endif
405    cph Avoid variable aliasing for adjoint !!!
406                DO j=jMin,jMax
407                 DO i=iMin,iMax
408                  rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
409                 ENDDO
410                ENDDO
411              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
412       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
413       I             rhoK, rhoKm1, rhoK,       I             rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
414       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
415       I             myThid )       I             myThid )
           ENDIF  
   
416  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
417  CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  #ifdef GMREDI_WITH_STABLE_ADJOINT
418  CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
419    cgf -> cuts adjoint dependency from slope to state
420                CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
421                CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
422                CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
423    #endif
424  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
425              ENDIF
426    
427  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
428  c ==> should use sigmaR !!!  c ==> should use sigmaR !!!
429            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
430  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
431              IF ( debugLevel .GE. debLevB )              IF ( debugLevel .GE. debLevB )
432       &       CALL DEBUG_CALL('CALC_IVDC',myThid)       &       CALL DEBUG_CALL('CALC_IVDC',myThid)
433  #endif  #endif
434              CALL CALC_IVDC(              CALL CALC_IVDC(
435       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
436       I        rhoKm1, rhoK,       I        rhoKm1, rhoInSitu(1-OLx,1-OLy,k,bi,bj),
437       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
438            ENDIF            ENDIF
439    
440    #ifdef ALLOW_DIAGNOSTICS
441              IF ( MOD(doDiagsRho,2).EQ.1 ) THEN
442                CALL DIAGS_RHO_L( k, bi, bj,
443         I                        rhoInSitu(1-OLx,1-OLy,k,bi,bj),
444         I                        rhoKm1, wVel,
445         I                        myTime, myIter, myThid )
446              ENDIF
447    #endif
448    
449  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
450          ENDDO          ENDDO
451    
452  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
453  cph avoids recomputation of integrate_for_w  CADJ STORE IVDConvCount(:,:,:,bi,bj)
454  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #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)  
455  #endif  #endif
456            CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,  
457       I            uVel, vVel, wVel, theta, salt,  C--     Diagnose Mixed Layer Depth:
458       I            myThid )          IF ( useGMRedi .OR. doDiagsRho.GE.4 ) THEN
459              CALL CALC_OCE_MXLAYER(
460         I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
461         I              bi, bj, myTime, myIter, myThid )
462          ENDIF          ENDIF
 #endif  /* ALLOW_OBCS */  
463    
464  #ifndef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_SALT_PLUME
465          IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN          IF ( useSALT_PLUME ) THEN
466  #endif            CALL SALT_PLUME_CALC_DEPTH(
467         I              rhoInSitu(1-OLx,1-OLy,1,bi,bj), sigmaR,
468         I              bi, bj, myTime, myIter, myThid )
469            ENDIF
470    #endif /* ALLOW_SALT_PLUME */
471    
472    #ifdef ALLOW_DIAGNOSTICS
473            IF ( MOD(doDiagsRho,4).GE.2 ) THEN
474              CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,
475         &         2, bi, bj, myThid)
476            ENDIF
477    #endif /* ALLOW_DIAGNOSTICS */
478    
479  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
480  C       relaxation terms, etc.  C       relaxation terms, etc.
481  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
482          IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
483       &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)       &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
484  #endif  #endif
485           CALL EXTERNAL_FORCING_SURF(  #ifdef ALLOW_AUTODIFF_TAMC
486    CADJ STORE EmPmR(:,:,bi,bj)
487    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
488    # ifdef EXACT_CONSERV
489    CADJ STORE PmEpR(:,:,bi,bj)
490    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
491    # endif
492    # ifdef NONLIN_FRSURF
493    CADJ STORE hFac_surfC(:,:,bi,bj)
494    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
495    CADJ STORE recip_hFacC(:,:,:,bi,bj)
496    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
497    # endif
498    #endif
499            CALL EXTERNAL_FORCING_SURF(
500       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
501       I             myTime, myIter, myThid )       I             myTime, myIter, myThid )
502  #ifndef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
503          ENDIF  # ifdef EXACT_CONSERV
504    cph-test
505    cphCADJ STORE PmEpR(:,:,bi,bj)
506    cphCADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
507    # endif
508  #endif  #endif
509    
510  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
511  cph needed for KPP  cph needed for KPP
512  CADJ STORE surfacetendencyU(:,:,bi,bj)  CADJ STORE surfaceForcingU(:,:,bi,bj)
513  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
514  CADJ STORE surfacetendencyV(:,:,bi,bj)  CADJ STORE surfaceForcingV(:,:,bi,bj)
515  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
516  CADJ STORE surfacetendencyS(:,:,bi,bj)  CADJ STORE surfaceForcingS(:,:,bi,bj)
517  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
518  CADJ STORE surfacetendencyT(:,:,bi,bj)  CADJ STORE surfaceForcingT(:,:,bi,bj)
519  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
520  # ifdef ALLOW_SEAICE  CADJ STORE surfaceForcingTice(:,:,bi,bj)
 CADJ STORE surfacetendencyTice(:,:,bi,bj)  
521  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
 # endif  
522  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
523    
524  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_KPP
525    C--     Compute KPP mixing coefficients
526            IF (useKPP) THEN
527    #ifdef ALLOW_DEBUG
528              IF ( debugLevel .GE. debLevB )
529         &     CALL DEBUG_CALL('KPP_CALC',myThid)
530    #endif
531              CALL KPP_CALC(
532         I                  bi, bj, myTime, myIter, myThid )
533    #ifdef ALLOW_AUTODIFF_TAMC
534            ELSE
535              CALL KPP_CALC_DUMMY(
536         I                  bi, bj, myTime, myIter, myThid )
537    #endif /* ALLOW_AUTODIFF_TAMC */
538            ENDIF
539    
540    #endif  /* ALLOW_KPP */
541    
542    #ifdef  ALLOW_PP81
543    C--     Compute PP81 mixing coefficients
544            IF (usePP81) THEN
545    #ifdef ALLOW_DEBUG
546              IF ( debugLevel .GE. debLevB )
547         &     CALL DEBUG_CALL('PP81_CALC',myThid)
548    #endif
549              CALL PP81_CALC(
550         I                  bi, bj, myTime, myThid )
551            ENDIF
552    #endif /* ALLOW_PP81 */
553    
554    #ifdef  ALLOW_MY82
555    C--     Compute MY82 mixing coefficients
556            IF (useMY82) THEN
557    #ifdef ALLOW_DEBUG
558              IF ( debugLevel .GE. debLevB )
559         &     CALL DEBUG_CALL('MY82_CALC',myThid)
560    #endif
561              CALL MY82_CALC(
562         I                  bi, bj, myTime, myThid )
563            ENDIF
564    #endif /* ALLOW_MY82 */
565    
566    #ifdef  ALLOW_GGL90
567    C--     Compute GGL90 mixing coefficients
568            IF (useGGL90) THEN
569    #ifdef ALLOW_DEBUG
570              IF ( debugLevel .GE. debLevB )
571         &     CALL DEBUG_CALL('GGL90_CALC',myThid)
572    #endif
573              CALL GGL90_CALC(
574         I                  bi, bj, myTime, myThid )
575            ENDIF
576    #endif /* ALLOW_GGL90 */
577    
578    #ifdef ALLOW_TIMEAVE
579            IF ( taveFreq.GT. 0. _d 0 ) THEN
580              CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
581            ENDIF
582            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
583              CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
584         I                           Nr, deltaTclock, bi, bj, myThid)
585            ENDIF
586    #endif /* ALLOW_TIMEAVE */
587    
588    #ifdef ALLOW_GMREDI
589  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
590    # ifndef GM_EXCLUDE_CLIPPING
591  cph storing here is needed only for one GMREDI_OPTIONS:  cph storing here is needed only for one GMREDI_OPTIONS:
592  cph define GM_BOLUS_ADVEC  cph define GM_BOLUS_ADVEC
593    cph keep it although TAF says you dont need to.
594  cph but I've avoided the #ifdef for now, in case more things change  cph but I've avoided the #ifdef for now, in case more things change
595  CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
596  CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
597  CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
598    # endif
599  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
600    
601  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
602          IF (useGMRedi) THEN          IF (useGMRedi) THEN
603  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
604            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
605       &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)       &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
606  #endif  #endif
607            CALL GMREDI_CALC_TENSOR(            CALL GMREDI_CALC_TENSOR(
608       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
609       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
610       I             myThid )       I             bi, bj, myTime, myIter, myThid )
611  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
612          ELSE          ELSE
613            CALL GMREDI_CALC_TENSOR_DUMMY(            CALL GMREDI_CALC_TENSOR_DUMMY(
614       I             bi, bj, iMin, iMax, jMin, jMax,       I             iMin, iMax, jMin, jMax,
615       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
616       I             myThid )       I             bi, bj, myTime, myIter, myThid )
617  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
618          ENDIF          ENDIF
619    #endif /* ALLOW_GMREDI */
620    
621  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_DOWN_SLOPE
622  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte          IF ( useDOWN_SLOPE ) THEN
623  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte  C--     Calculate Downsloping Flow for Down_Slope parameterization
624  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte           IF ( usingPCoords ) THEN
625  #endif /* ALLOW_AUTODIFF_TAMC */            CALL DWNSLP_CALC_FLOW(
626         I                bi, bj, kSurfC, rhoInSitu,
627         I                myTime, myIter, myThid )
628             ELSE
629              CALL DWNSLP_CALC_FLOW(
630         I                bi, bj, kLowC, rhoInSitu,
631         I                myTime, myIter, myThid )
632             ENDIF
633            ENDIF
634    #endif /* ALLOW_DOWN_SLOPE */
635    
636  #endif  /* ALLOW_GMREDI */  #ifndef ALLOW_AUTODIFF_TAMC
637    C---  if fluid Is Water: end
638            ENDIF
639    #endif
640    
641  #ifdef  ALLOW_KPP  #ifdef  ALLOW_OBCS
642  C--     Compute KPP mixing coefficients  C--     Calculate future values on open boundaries
643          IF (useKPP) THEN          IF (useOBCS) THEN
644  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
645            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
646       &     CALL DEBUG_CALL('KPP_CALC',myThid)       &     CALL DEBUG_CALL('OBCS_CALC',myThid)
647  #endif  #endif
648            CALL KPP_CALC(            CALL OBCS_CALC( bi, bj, myTime+deltaTclock, myIter+1,
649       I                  bi, bj, myTime, myThid )       I            uVel, vVel, wVel, theta, salt,
650  #ifdef ALLOW_AUTODIFF_TAMC       I            myThid )
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
651          ENDIF          ENDIF
652    #endif  /* ALLOW_OBCS */
653    
654  #ifdef ALLOW_AUTODIFF_TAMC  C--   end bi,bj loops.
655  CADJ STORE KPPghat   (:,:,:,bi,bj)         ENDDO
656  CADJ &   , KPPfrac   (:,:  ,bi,bj)        ENDDO
 CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
657    
658    #ifdef  ALLOW_KPP
659          IF (useKPP) THEN
660            CALL KPP_DO_EXCH( myThid )
661          ENDIF
662  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
663    
664  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_DIAGNOSTICS
665  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte        IF ( fluidIsWater .AND. useDiagnostics ) THEN
666  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte          CALL DIAGS_RHO_G(
667  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte       I                    rhoInSitu, uVel, vVel,
668  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte       I                    myTime, myIter, myThid )
669  #ifdef ALLOW_PASSIVE_TRACER          CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
670  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte        ENDIF
671  #endif        IF ( ivdc_kappa.NE.0 .AND. useDiagnostics ) THEN
672  #ifdef ALLOW_PTRACERS          CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
673  cph-- moved to forward_step to avoid key computation       &                               0, Nr, 0, 1, 1, myThid )
674  cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,        ENDIF
 cphCADJ &                              key=itdkey, byte=isbyte  
675  #endif  #endif
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--   end bi,bj loops.  
        ENDDO  
       ENDDO  
676    
677  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
678           IF ( debugLevel .GE. debLevB )        IF ( debugLevel .GE. debLevB )
679       &    CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)       &     CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
680  #endif  #endif
681    
682        RETURN        RETURN
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.72
  ViewVC Help
Powered by ViewVC 1.1.22