model/src/do_oceanic_phys.F

c842b53753 Jean*0001 #include "PACKAGES_CONFIG.h"
                0002 #include "CPP_OPTIONS.h"
f59d76b0dd Ed D*0003 #ifdef ALLOW_MOM_COMMON
                0004 # include "MOM_COMMON_OPTIONS.h"
                0005 #endif
517dbdc414 Jean*0006 #ifdef ALLOW_AUTODIFF
                0007 # include "AUTODIFF_OPTIONS.h"
                0008 #endif
                0009 #ifdef ALLOW_CTRL
                0010 # include "CTRL_OPTIONS.h"
                0011 #endif
1f89baba18 Patr*0012 #ifdef ALLOW_SALT_PLUME
                0013 # include "SALT_PLUME_OPTIONS.h"
                0014 #endif
ecaa45da39 An T*0015 #ifdef ALLOW_ECCO
                0016 # include "ECCO_OPTIONS.h"
                0017 #endif
c842b53753 Jean*0018
517dbdc414 Jean*0019 #ifdef ALLOW_AUTODIFF
685f7544b6 Patr*0020 # ifdef ALLOW_GGL90
                0021 #  include "GGL90_OPTIONS.h"
                0022 # endif
c842b53753 Jean*0023 # ifdef ALLOW_GMREDI
                0024 #  include "GMREDI_OPTIONS.h"
                0025 # endif
                0026 # ifdef ALLOW_KPP
                0027 #  include "KPP_OPTIONS.h"
                0028 # endif
9d4847a995 Jean*0029 # ifdef ALLOW_SEAICE
                0030 #  include "SEAICE_OPTIONS.h"
                0031 # endif
9e3a303fa3 Gael*0032 # ifdef ALLOW_EXF
                0033 #  include "EXF_OPTIONS.h"
                0034 # endif
af61e5eb16 Mart*0035 #ifdef ALLOW_OBCS
                0036 # include "OBCS_OPTIONS.h"
                0037 #endif
517dbdc414 Jean*0038 #endif /* ALLOW_AUTODIFF */
c842b53753 Jean*0039
                0040 CBOP
                0041 C     !ROUTINE: DO_OCEANIC_PHYS
                0042 C     !INTERFACE:
                0043       SUBROUTINE DO_OCEANIC_PHYS(myTime, myIter, myThid)
                0044 C     !DESCRIPTION: \bv
                0045 C     *==========================================================*
765f2ca955 Jean*0046 C     | SUBROUTINE DO_OCEANIC_PHYS
                0047 C     | o Controlling routine for oceanic physics and
c842b53753 Jean*0048 C     |   parameterization
                0049 C     *==========================================================*
                0050 C     | o originally, part of S/R thermodynamics
                0051 C     *==========================================================*
                0052 C     \ev
                0053
8f92343d9b Jean*0054 C     !CALLING SEQUENCE:
                0055 C     DO_OCEANIC_PHYS
                0056 C       |
                0057 C       |-- OBCS_CALC
                0058 C       |
7eca977925 Jean*0059 C       |-- OCN_APPLY_IMPORT
                0060 C       |
8f92343d9b Jean*0061 C       |-- FRAZIL_CALC_RHS
                0062 C       |
                0063 C       |-- THSICE_MAIN
                0064 C       |
                0065 C       |-- SEAICE_FAKE
                0066 C       |-- SEAICE_MODEL
                0067 C       |-- SEAICE_COST_SENSI
                0068 C       |
7eca977925 Jean*0069 C       |-- OCN_EXPORT_DATA
                0070 C       |
00f81e6785 Ou W*0071 C       |-- STIC_THERMODYNAMICS
8f92343d9b Jean*0072 C       |-- SHELFICE_THERMODYNAMICS
                0073 C       |
                0074 C       |-- ICEFRONT_THERMODYNAMICS
                0075 C       |
                0076 C       |-- SALT_PLUME_DO_EXCH
                0077 C       |
                0078 C       |-- FREEZE_SURFACE
                0079 C       |
                0080 C       |-- EXTERNAL_FORCING_SURF
                0081 C       |
abfe198bce Mart*0082 C       |-- OBCS_ADJUST
                0083 C       |
8f92343d9b Jean*0084 C       |- k loop (Nr:1):
                0085 C       | - DWNSLP_CALC_RHO
                0086 C       | - BBL_CALC_RHO
                0087 C       | - FIND_RHO_2D @ p(k)
                0088 C       | - FIND_RHO_2D @ p(k-1)
                0089 C       | - GRAD_SIGMA
                0090 C       | - CALC_IVDC
                0091 C       | - DIAGS_RHO_L
                0092 C       |- end k loop.
                0093 C       |
                0094 C       |-- CALC_OCE_MXLAYER
                0095 C       |
                0096 C       |-- SALT_PLUME_CALC_DEPTH
                0097 C       |-- SALT_PLUME_VOLFRAC
                0098 C       |-- SALT_PLUME_APPLY
                0099 C       |-- SALT_PLUME_APPLY
                0100 C       |-- SALT_PLUME_FORCING_SURF
                0101 C       |
                0102 C       |-- KPP_CALC
                0103 C       |-- KPP_CALC_DUMMY
                0104 C       |
                0105 C       |-- PP81_CALC
                0106 C       |
                0107 C       |-- KL10_CALC
                0108 C       |
                0109 C       |-- MY82_CALC
                0110 C       |
                0111 C       |-- GGL90_CALC
                0112 C       |
                0113 C       |-- TIMEAVE_SURF_FLUX
                0114 C       |
                0115 C       |-- GMREDI_CALC_TENSOR
                0116 C       |-- GMREDI_CALC_TENSOR_DUMMY
                0117 C       |
                0118 C       |-- DWNSLP_CALC_FLOW
                0119 C       |-- DWNSLP_CALC_FLOW
                0120 C       |
dc2ebbdaf0 Jean*0121 C       |-- OFFLINE_GET_DIFFUS
                0122 C       |
8f92343d9b Jean*0123 C       |-- BBL_CALC_RHS
                0124 C       |
                0125 C       |-- MYPACKAGE_CALC_RHS
                0126 C       |
                0127 C       |-- GMREDI_DO_EXCH
                0128 C       |
                0129 C       |-- KPP_DO_EXCH
                0130 C       |
0320e25227 Mart*0131 C       |-- GGL90_EXCHANGES
                0132 C       |
8f92343d9b Jean*0133 C       |-- DIAGS_RHO_G
                0134 C       |-- DIAGS_OCEANIC_SURF_FLUX
                0135 C       |-- SALT_PLUME_DIAGNOSTICS_FILL
                0136 C       |
                0137 C       |-- ECCO_PHYS
                0138
c842b53753 Jean*0139 C     !USES:
                0140       IMPLICIT NONE
                0141 C     == Global variables ===
                0142 #include "SIZE.h"
                0143 #include "EEPARAMS.h"
                0144 #include "PARAMS.h"
                0145 #include "GRID.h"
1db41719d4 Jean*0146 #include "DYNVARS.h"
c4e20a6b6c Jean*0147 #ifdef ALLOW_TIMEAVE
c8b8efc127 Jean*0148 # include "TIMEAVE_STATV.h"
                0149 #endif
                0150 #ifdef ALLOW_OFFLINE
                0151 # include "OFFLINE_SWITCH.h"
c4e20a6b6c Jean*0152 #endif
c842b53753 Jean*0153
517dbdc414 Jean*0154 #ifdef ALLOW_AUTODIFF
7f4c5015d5 Patr*0155 # include "AUTODIFF_MYFIELDS.h"
7c50f07931 Mart*0156 # ifdef ALLOW_AUTODIFF_TAMC
                0157 #  include "tamc.h"
                0158 # endif
c842b53753 Jean*0159 # include "FFIELDS.h"
81b43af3f0 Patr*0160 # include "SURFACE.h"
c842b53753 Jean*0161 # include "EOS.h"
685f7544b6 Patr*0162 # ifdef ALLOW_GMREDI
                0163 #  include "GMREDI.h"
                0164 # endif
c842b53753 Jean*0165 # ifdef ALLOW_KPP
                0166 #  include "KPP.h"
                0167 # endif
2b52c649e6 Gael*0168 # ifdef ALLOW_GGL90
                0169 #  include "GGL90.h"
                0170 # endif
c842b53753 Jean*0171 # ifdef ALLOW_EBM
                0172 #  include "EBM.h"
                0173 # endif
9d4847a995 Jean*0174 # ifdef ALLOW_EXF
b4daa24319 Shre*0175 #  ifdef ALLOW_CTRL
5cf4364659 Mart*0176 #   include "CTRL_SIZE.h"
4d72283393 Mart*0177 #   include "CTRL.h"
b4daa24319 Shre*0178 #  endif
798a745844 Jean*0179 #  include "EXF_FIELDS.h"
9d4847a995 Jean*0180 #  ifdef ALLOW_BULKFORMULAE
798a745844 Jean*0181 #   include "EXF_CONSTANTS.h"
9d4847a995 Jean*0182 #  endif
                0183 # endif
                0184 # ifdef ALLOW_SEAICE
a2c8840cf9 Jean*0185 #  include "SEAICE_SIZE.h"
9d4847a995 Jean*0186 #  include "SEAICE.h"
2848258862 Gael*0187 #  include "SEAICE_PARAMS.h"
9d4847a995 Jean*0188 # endif
ec745395c8 Patr*0189 # ifdef ALLOW_THSICE
                0190 #  include "THSICE_VARS.h"
                0191 # endif
f7911aa28f Patr*0192 # ifdef ALLOW_SALT_PLUME
                0193 #  include "SALT_PLUME.h"
                0194 # endif
af61e5eb16 Mart*0195 # ifdef ALLOW_OBCS
                0196 #  include "OBCS_PARAMS.h"
                0197 #  include "OBCS_FIELDS.h"
                0198 # endif
517dbdc414 Jean*0199 #endif /* ALLOW_AUTODIFF */
c842b53753 Jean*0200
b4daa24319 Shre*0201 #ifdef ALLOW_TAPENADE
                0202 c# ifdef ALLOW_KPP
                0203 c#  include "KPP_PARAMS.h"
                0204 c# endif
                0205 # ifdef ALLOW_SHELFICE
                0206 #  include "SHELFICE.h"
                0207 # endif
                0208 #endif /* ALLOW_TAPENADE */
                0209
c842b53753 Jean*0210 C     !INPUT/OUTPUT PARAMETERS:
                0211 C     == Routine arguments ==
38ee7d4239 Jean*0212 C     myTime :: Current time in simulation
                0213 C     myIter :: Current iteration number in simulation
                0214 C     myThid :: Thread number for this instance of the routine.
c842b53753 Jean*0215       _RL myTime
                0216       INTEGER myIter
                0217       INTEGER myThid
                0218
                0219 C     !LOCAL VARIABLES:
                0220 C     == Local variables
0320e25227 Mart*0221 C     rhoKp1,rhoKm1 :: Density at current level, and @ level minus one
38ee7d4239 Jean*0222 C     iMin, iMax    :: Ranges and sub-block indices on which calculations
c842b53753 Jean*0223 C     jMin, jMax       are applied.
38ee7d4239 Jean*0224 C     bi, bj        :: tile indices
7eca977925 Jean*0225 C     msgBuf        :: Temp. for building output string
38ee7d4239 Jean*0226 C     i,j,k         :: loop indices
0320e25227 Mart*0227 C     kSrf          :: surface index
cf9e43202e Jean*0228       _RL rhoKp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0229       _RL rhoKm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c842b53753 Jean*0230       _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0231       _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0232       _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0233       INTEGER iMin, iMax
                0234       INTEGER jMin, jMax
                0235       INTEGER bi, bj
0320e25227 Mart*0236       INTEGER i, j, k, kSrf
7eca977925 Jean*0237       CHARACTER*(MAX_LEN_MBUF) msgBuf
8ab93ea9a3 Jean*0238       INTEGER doDiagsRho
c8b8efc127 Jean*0239       LOGICAL calcGMRedi, calcKPP, calcConvect
8ab93ea9a3 Jean*0240 #ifdef ALLOW_DIAGNOSTICS
                0241       LOGICAL  DIAGNOSTICS_IS_ON
                0242       EXTERNAL DIAGNOSTICS_IS_ON
                0243 #endif /* ALLOW_DIAGNOSTICS */
34de4f2d34 Jean*0244 #ifdef ALLOW_AUTODIFF
                0245       _RL thetaRef
                0246 #endif /* ALLOW_AUTODIFF */
7c50f07931 Mart*0247 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0248 C     tkey :: tape key (tile dependent)
af61e5eb16 Mart*0249       INTEGER tkey
7c50f07931 Mart*0250 #endif
c842b53753 Jean*0251 CEOP
                0252
0320e25227 Mart*0253       kSrf = 1
                0254       IF ( usingPCoords ) kSrf = Nr
                0255
c842b53753 Jean*0256 #ifdef ALLOW_DEBUG
8440e8ae5d Jean*0257       IF (debugMode) CALL DEBUG_ENTER('DO_OCEANIC_PHYS',myThid)
c842b53753 Jean*0258 #endif
0ff4deb339 Jean*0259
8ab93ea9a3 Jean*0260       doDiagsRho = 0
                0261 #ifdef ALLOW_DIAGNOSTICS
                0262       IF ( useDiagnostics .AND. fluidIsWater ) THEN
d497465397 Jean*0263         IF ( DIAGNOSTICS_IS_ON('MXLDEPTH',myThid) )
842349fcb7 Jean*0264      &       doDiagsRho = doDiagsRho + 1
                0265         IF ( DIAGNOSTICS_IS_ON('DRHODR  ',myThid) )
                0266      &       doDiagsRho = doDiagsRho + 2
d497465397 Jean*0267         IF ( DIAGNOSTICS_IS_ON('WdRHO_P ',myThid) )
842349fcb7 Jean*0268      &       doDiagsRho = doDiagsRho + 4
d497465397 Jean*0269         IF ( DIAGNOSTICS_IS_ON('WdRHOdP ',myThid) )
                0270      &       doDiagsRho = doDiagsRho + 8
8ab93ea9a3 Jean*0271       ENDIF
                0272 #endif /* ALLOW_DIAGNOSTICS */
                0273
c8b8efc127 Jean*0274       calcGMRedi  = useGMRedi
                0275       calcKPP     = useKPP
                0276       calcConvect = ivdc_kappa.NE.0.
                0277 #ifdef ALLOW_OFFLINE
                0278       IF ( useOffLine ) THEN
                0279         calcGMRedi = useGMRedi .AND. .NOT.offlineLoadGMRedi
                0280         calcKPP    = useKPP    .AND. .NOT.offlineLoadKPP
                0281         calcConvect=calcConvect.AND. .NOT.offlineLoadConvec
                0282       ENDIF
                0283 #endif /* ALLOW_OFFLINE */
                0284
f7d550a4f6 Jean*0285 #ifdef  ALLOW_OBCS
                0286       IF (useOBCS) THEN
                0287 C--   Calculate future values on open boundaries
                0288 C--   moved before SEAICE_MODEL call since SEAICE_MODEL needs seaice-obcs fields
38b71af82a Patr*0289 # ifdef ALLOW_AUTODIFF_TAMC
                0290 CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
                0291 CADJ STORE salt  = comlev1, key=ikey_dynamics, kind=isbyte
af61e5eb16 Mart*0292 CADJ STORE etaN  = comlev1, key=ikey_dynamics, kind=isbyte
                0293 #  ifdef ALLOW_OBCS_STEVENS
                0294 CADJ STORE uVel  = comlev1, key=ikey_dynamics, kind=isbyte
                0295 CADJ STORE vVel  = comlev1, key=ikey_dynamics, kind=isbyte
                0296 #   ifdef ALLOW_OBCS_EAST
                0297 CADJ STORE OBEtStevens = comlev1, key=ikey_dynamics, kind=isbyte
                0298 CADJ STORE OBEsStevens = comlev1, key=ikey_dynamics, kind=isbyte
                0299 #   endif
                0300 #   ifdef ALLOW_OBCS_WEST
                0301 CADJ STORE OBWtStevens = comlev1, key=ikey_dynamics, kind=isbyte
                0302 CADJ STORE OBWsStevens = comlev1, key=ikey_dynamics, kind=isbyte
                0303 #   endif
                0304 #   ifdef ALLOW_OBCS_NORTH
                0305 CADJ STORE OBNtStevens = comlev1, key=ikey_dynamics, kind=isbyte
                0306 CADJ STORE OBNsStevens = comlev1, key=ikey_dynamics, kind=isbyte
                0307 #   endif
                0308 #   ifdef ALLOW_OBCS_SOUTH
                0309 CADJ STORE OBStStevens = comlev1, key=ikey_dynamics, kind=isbyte
                0310 CADJ STORE OBSsStevens = comlev1, key=ikey_dynamics, kind=isbyte
                0311 #   endif
                0312 #  endif /* ALLOW_OBCS_STEVENS */
38b71af82a Patr*0313 # endif
                0314 # ifdef ALLOW_DEBUG
8440e8ae5d Jean*0315        IF (debugMode) CALL DEBUG_CALL('OBCS_CALC',myThid)
38b71af82a Patr*0316 # endif
70075a2466 Jean*0317        CALL OBCS_CALC( myTime+deltaTClock, myIter+1,
3a86c9b47d Oliv*0318      I                 uVel, vVel, wVel, theta, salt, myThid )
f7d550a4f6 Jean*0319       ENDIF
af61e5eb16 Mart*0320 # if ( defined ALLOW_AUTODIFF_TAMC && defined ALLOW_OBCS_BALANCE )
                0321 C     This needs to be done ***after*** the if-block to avoid calling
                0322 C     S/R OBCS_CALC in the AD code.
                0323 #  ifdef ALLOW_OBCS_NORTH
                0324 CADJ STORE OBNv        = comlev1, key=ikey_dynamics, kind=isbyte
                0325 #  endif
                0326 #  ifdef ALLOW_OBCS_SOUTH
                0327 CADJ STORE OBSv        = comlev1, key=ikey_dynamics, kind=isbyte
                0328 #  endif
                0329 #  ifdef ALLOW_OBCS_EAST
                0330 CADJ STORE OBEu        = comlev1, key=ikey_dynamics, kind=isbyte
                0331 #  endif
                0332 #  ifdef ALLOW_OBCS_WEST
                0333 CADJ STORE OBWu        = comlev1, key=ikey_dynamics, kind=isbyte
                0334 #  endif
                0335 # endif /* ALLOW_AUTODIFF_TAMC and ALLOW_OBCS_BALANCE */
f7d550a4f6 Jean*0336 #endif  /* ALLOW_OBCS */
1db41719d4 Jean*0337
7eca977925 Jean*0338 #ifdef ALLOW_OCN_COMPON_INTERF
                0339 C--    Apply imported data (from coupled interface) to forcing fields
                0340 C jmc: moved here before any freezing/seaice pkg adjustment of surf-fluxes
                0341       IF ( useCoupler ) THEN
                0342          CALL OCN_APPLY_IMPORT( .TRUE., myTime, myIter, myThid )
                0343       ENDIF
                0344 #endif /* ALLOW_OCN_COMPON_INTERF */
                0345
517dbdc414 Jean*0346 #ifdef ALLOW_AUTODIFF
79d91b77ec Gael*0347       DO bj=myByLo(myThid),myByHi(myThid)
                0348        DO bi=myBxLo(myThid),myBxHi(myThid)
                0349         DO j=1-OLy,sNy+OLy
                0350          DO i=1-OLx,sNx+OLx
fa323f784d Jean*0351           adjustColdSST_diag(i,j,bi,bj) = 0. _d 0
                0352 # ifdef ALLOW_SALT_PLUME
79d91b77ec Gael*0353           saltPlumeDepth(i,j,bi,bj) = 0. _d 0
                0354           saltPlumeFlux(i,j,bi,bj)  = 0. _d 0
fa323f784d Jean*0355 # endif
79d91b77ec Gael*0356          ENDDO
                0357         ENDDO
                0358        ENDDO
                0359       ENDDO
f59d76b0dd Ed D*0360 #endif /* ALLOW_AUTODIFF */
                0361
e0b3e1bdd8 Dimi*0362 #ifdef ALLOW_FRAZIL
                0363       IF ( useFRAZIL ) THEN
                0364 C--   Freeze water in the ocean interior and let it rise to the surface
dde0ec0a2d Patr*0365 CADJ STORE theta = comlev1, key=ikey_dynamics, kind=isbyte
                0366 CADJ STORE salt  = comlev1, key=ikey_dynamics, kind=isbyte
e0b3e1bdd8 Dimi*0367        CALL FRAZIL_CALC_RHS( myTime, myIter, myThid )
                0368       ENDIF
                0369 #endif /* ALLOW_FRAZIL */
                0370
de80e32bba Jean*0371 #if (defined ALLOW_THSICE) && !(defined ALLOW_ATM2D)
                0372       IF ( useThSIce .AND. fluidIsWater ) THEN
d503ff4f7a Jean*0373 # ifdef ALLOW_AUTODIFF_TAMC
9c41af81f6 Timo*0374 #  ifdef ALLOW_SEAICE
ec0d7df165 Mart*0375 CADJ STORE uice,vice         = comlev1, key=ikey_dynamics, kind=isbyte
9c41af81f6 Timo*0376 #  endif
ec0d7df165 Mart*0377 CADJ STORE iceMask,iceHeight = comlev1, key=ikey_dynamics, kind=isbyte
                0378 CADJ STORE snowHeight, Tsrf  = comlev1, key=ikey_dynamics, kind=isbyte
                0379 CADJ STORE Qice1, Qice2      = comlev1, key=ikey_dynamics, kind=isbyte
                0380 CADJ STORE sHeating,snowAge  = comlev1, key=ikey_dynamics, kind=isbyte
                0381 CADJ STORE hocemxl, icflxsw  = comlev1, key=ikey_dynamics, kind=isbyte
                0382 CADJ STORE salt,theta        = comlev1, key=ikey_dynamics, kind=isbyte
                0383 CADJ STORE uvel,vvel         = comlev1, key=ikey_dynamics, kind=isbyte
                0384 CADJ STORE qnet,qsw, empmr   = comlev1, key=ikey_dynamics, kind=isbyte
                0385 CADJ STORE atemp,aqh,precip  = comlev1, key=ikey_dynamics, kind=isbyte
                0386 CADJ STORE swdown,lwdown     = comlev1, key=ikey_dynamics, kind=isbyte
d503ff4f7a Jean*0387 #  ifdef NONLIN_FRSURF
ec0d7df165 Mart*0388 CADJ STORE hFac_surfC        = comlev1, key=ikey_dynamics, kind=isbyte
d503ff4f7a Jean*0389 #  endif
a7828ed4d5 Patr*0390 # endif /* ALLOW_AUTODIFF_TAMC */
de80e32bba Jean*0391 # ifdef ALLOW_DEBUG
                0392         IF (debugMode) CALL DEBUG_CALL('THSICE_MAIN',myThid)
                0393 # endif
                0394 C--     Step forward Therm.Sea-Ice variables
                0395 C       and modify forcing terms including effects from ice
                0396         CALL TIMER_START('THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
                0397         CALL THSICE_MAIN( myTime, myIter, myThid )
                0398         CALL TIMER_STOP( 'THSICE_MAIN     [DO_OCEANIC_PHYS]', myThid)
                0399       ENDIF
                0400 #endif /* ALLOW_THSICE */
                0401
9d4847a995 Jean*0402 #ifdef ALLOW_SEAICE
ec0d7df165 Mart*0403 # ifdef ALLOW_AUTODIFF_TAMC
3314a43a56 Mart*0404 CADJ STORE qnet  = comlev1, key=ikey_dynamics, kind=isbyte
                0405 CADJ STORE qsw   = comlev1, key=ikey_dynamics, kind=isbyte
                0406 CADJ STORE fu,fv = comlev1, key=ikey_dynamics, kind=isbyte
af61e5eb16 Mart*0407 #  if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
3314a43a56 Mart*0408 CADJ STORE evap  = comlev1, key=ikey_dynamics, kind=isbyte
af61e5eb16 Mart*0409 #  endif
                0410 CADJ STORE etan      = comlev1, key=ikey_dynamics, kind=isbyte
                0411 CADJ STORE theta     = comlev1, key=ikey_dynamics, kind=isbyte
                0412 CADJ STORE salt      = comlev1, key=ikey_dynamics, kind=isbyte
                0413 CADJ STORE uvel,vvel = comlev1, key=ikey_dynamics, kind=isbyte
                0414 CADJ STORE phiHydLow = comlev1, key=ikey_dynamics, byte=isbyte
                0415 # endif
9d4847a995 Jean*0416       IF ( useSEAICE ) THEN
7880c396f6 Patr*0417 # ifdef ALLOW_AUTODIFF_TAMC
3314a43a56 Mart*0418 CADJ STORE uvel,vvel         = comlev1, key=ikey_dynamics, kind=isbyte
                0419 CADJ STORE uice,vice         = comlev1, key=ikey_dynamics, kind=isbyte
3c775cbf98 Mart*0420 #  ifdef SEAICE_USE_GROWTH_ADX
                0421 CADJ STORE tices             = comlev1, key=ikey_dynamics, kind=isbyte
                0422 #  endif
3314a43a56 Mart*0423 #  ifdef ALLOW_EXF
ec0d7df165 Mart*0424 CADJ STORE atemp,aqh,precip  = comlev1, key=ikey_dynamics, kind=isbyte
                0425 CADJ STORE swdown,lwdown     = comlev1, key=ikey_dynamics, kind=isbyte
                0426 CADJ STORE uwind,vwind       = comlev1, key=ikey_dynamics, kind=isbyte
7880c396f6 Patr*0427 #  endif
f4b023d948 Jean*0428 #  ifdef SEAICE_VARIABLE_SALINITY
ec0d7df165 Mart*0429 CADJ STORE hsalt             = comlev1, key=ikey_dynamics, kind=isbyte
f4b023d948 Jean*0430 #  endif
7880c396f6 Patr*0431 #  ifdef ATMOSPHERIC_LOADING
ec0d7df165 Mart*0432 CADJ STORE pload, siceload   = comlev1, key=ikey_dynamics, kind=isbyte
7880c396f6 Patr*0433 #  endif
                0434 #  ifdef NONLIN_FRSURF
ec0d7df165 Mart*0435 CADJ STORE recip_hfacc       = comlev1, key=ikey_dynamics, kind=isbyte
7880c396f6 Patr*0436 #  endif
a7828ed4d5 Patr*0437 #  ifdef ALLOW_THSICE
d5254b4e3d Mart*0438 C-- store thSIce vars before advection (called from SEAICE_MODEL) updates them:
ec0d7df165 Mart*0439 CADJ STORE iceMask,iceHeight = comlev1, key=ikey_dynamics, kind=isbyte
                0440 CADJ STORE snowHeight,hOceMxL= comlev1, key=ikey_dynamics, kind=isbyte
                0441 CADJ STORE Qice1, Qice2      = comlev1, key=ikey_dynamics, kind=isbyte
f4b023d948 Jean*0442 #  endif /* ALLOW_THSICE */
a7828ed4d5 Patr*0443 # endif /* ALLOW_AUTODIFF_TAMC */
7880c396f6 Patr*0444 # ifdef ALLOW_DEBUG
8440e8ae5d Jean*0445         IF (debugMode) CALL DEBUG_CALL('SEAICE_MODEL',myThid)
7880c396f6 Patr*0446 # endif
9d4847a995 Jean*0447         CALL TIMER_START('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
                0448         CALL SEAICE_MODEL( myTime, myIter, myThid )
                0449         CALL TIMER_STOP ('SEAICE_MODEL    [DO_OCEANIC_PHYS]', myThid)
3314a43a56 Mart*0450 # ifdef ALLOW_AUTODIFF_TAMC
                0451 CADJ STORE tices = comlev1, key=ikey_dynamics, kind=isbyte
                0452 CADJ STORE heff  = comlev1, key=ikey_dynamics, kind=isbyte
                0453 CADJ STORE hsnow = comlev1, key=ikey_dynamics, kind=isbyte
                0454 CADJ STORE area  = comlev1, key=ikey_dynamics, kind=isbyte
                0455 CADJ STORE uIce  = comlev1, key=ikey_dynamics, kind=isbyte
                0456 CADJ STORE vIce  = comlev1, key=ikey_dynamics, kind=isbyte
                0457 # endif
7880c396f6 Patr*0458 # ifdef ALLOW_COST
1e22f5fc71 Patr*0459         CALL SEAICE_COST_SENSI ( myTime, myIter, myThid )
7880c396f6 Patr*0460 # endif
ec0d7df165 Mart*0461 # ifdef ALLOW_AUTODIFF
                0462       ELSEIF ( SEAICEadjMODE .EQ. -1 ) THEN
af61e5eb16 Mart*0463 #  ifdef ALLOW_AUTODIFF_TAMC
3314a43a56 Mart*0464 CADJ STORE area = comlev1, key=ikey_dynamics, kind=isbyte
af61e5eb16 Mart*0465 #  endif
ec0d7df165 Mart*0466         CALL SEAICE_FAKE( myTime, myIter, myThid )
                0467 # endif /* ALLOW_AUTODIFF */
acb47bc07b Patr*0468       ENDIF
9d4847a995 Jean*0469 #endif /* ALLOW_SEAICE */
                0470
7eca977925 Jean*0471 #if (defined ALLOW_OCN_COMPON_INTERF) && (defined ALLOW_THSICE)
                0472 C--   After seaice-dyn and advection of pkg/thsice fields,
                0473 C     Export ocean coupling fields to coupled interface (only with pkg/thsice)
                0474       IF ( useCoupler ) THEN
                0475 # ifdef ALLOW_DEBUG
                0476         IF (debugMode) CALL DEBUG_CALL('OCN_EXPORT_DATA',myThid)
                0477 # endif
                0478          CALL TIMER_START('OCN_EXPORT_DATA [DO_OCEANIC_PHYS]', myThid)
                0479          CALL OCN_EXPORT_DATA( myTime, myIter, myThid )
                0480          CALL TIMER_STOP ('OCN_EXPORT_DATA [DO_OCEANIC_PHYS]', myThid)
                0481       ENDIF
                0482 #endif /* ALLOW_OCN_COMPON_INTERF & ALLOW_THSICE */
                0483
13783cb644 Patr*0484 #ifdef ALLOW_AUTODIFF_TAMC
ec0d7df165 Mart*0485 CADJ STORE sst, sss          = comlev1, key=ikey_dynamics, kind=isbyte
                0486 CADJ STORE qsw               = comlev1, key=ikey_dynamics, kind=isbyte
13783cb644 Patr*0487 # ifdef ALLOW_SEAICE
ec0d7df165 Mart*0488 CADJ STORE area              = comlev1, key=ikey_dynamics, kind=isbyte
13783cb644 Patr*0489 # endif
                0490 #endif
                0491
a6cbc7a360 Mart*0492 #ifdef ALLOW_SHELFICE
                0493       IF ( useShelfIce .AND. fluidIsWater ) THEN
dde0ec0a2d Patr*0494 #ifdef ALLOW_AUTODIFF_TAMC
ec0d7df165 Mart*0495 CADJ STORE salt, theta       = comlev1, key=ikey_dynamics, kind=isbyte
                0496 CADJ STORE uvel, vvel        = comlev1, key=ikey_dynamics, kind=isbyte
dde0ec0a2d Patr*0497 #endif
00f81e6785 Ou W*0498 #ifdef ALLOW_STEEP_ICECAVITY
                0499        IF ( useSTIC ) THEN
                0500 # ifdef ALLOW_DEBUG
                0501         IF (debugMode) CALL DEBUG_CALL('STIC_THERMODYNAMICS',myThid)
                0502 # endif
                0503 C     use stic_thermodynamics that includes icefront melt processes
                0504         CALL TIMER_START('STIC_THERMODYNAMICS [DO_OCEANIC_PHYS]',myThid)
                0505         CALL STIC_THERMODYNAMICS( myTime, myIter, myThid )
                0506         CALL TIMER_STOP ('STIC_THERMODYNAMICS [DO_OCEANIC_PHYS]',myThid)
                0507        ELSE
                0508 #else
                0509        IF ( .TRUE. ) THEN
                0510 #endif
                0511 #ifdef ALLOW_DEBUG
                0512         IF (debugMode) CALL DEBUG_CALL('SHELFICE_THERMODYNAMICS',myThid)
                0513 #endif
cf9e43202e Jean*0514 C     compute temperature and (virtual) salt flux at the
a6cbc7a360 Mart*0515 C     shelf-ice ocean interface
00f81e6785 Ou W*0516         CALL TIMER_START('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
a6cbc7a360 Mart*0517      &       myThid)
00f81e6785 Ou W*0518         CALL SHELFICE_THERMODYNAMICS( myTime, myIter, myThid )
                0519         CALL TIMER_STOP ('SHELFICE_THERMODYNAMICS [DO_OCEANIC_PHYS]',
                0520      &       myThid)
                0521        ENDIF
a6cbc7a360 Mart*0522       ENDIF
                0523 #endif /* ALLOW_SHELFICE */
                0524
5da8ce63fa Dimi*0525 #ifdef ALLOW_ICEFRONT
                0526       IF ( useICEFRONT .AND. fluidIsWater ) THEN
                0527 #ifdef ALLOW_DEBUG
8440e8ae5d Jean*0528        IF (debugMode) CALL DEBUG_CALL('ICEFRONT_THERMODYNAMICS',myThid)
5da8ce63fa Dimi*0529 #endif
                0530 C     compute temperature and (virtual) salt flux at the
                0531 C     ice-front ocean interface
                0532        CALL TIMER_START('ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
                0533      &       myThid)
                0534        CALL ICEFRONT_THERMODYNAMICS( myTime, myIter, myThid )
                0535        CALL TIMER_STOP( 'ICEFRONT_THERMODYNAMICS [DO_OCEANIC_PHYS]',
                0536      &      myThid)
                0537       ENDIF
                0538 #endif /* ALLOW_ICEFRONT */
                0539
9d71427857 Jean*0540 #ifdef ALLOW_SALT_PLUME
                0541       IF ( useSALT_PLUME ) THEN
1f89baba18 Patr*0542 Catn: exchanging saltPlumeFlux:
ec0d7df165 Mart*0543         CALL SALT_PLUME_DO_EXCH( myTime, myIter, myThid )
9d71427857 Jean*0544       ENDIF
                0545 #endif /* ALLOW_SALT_PLUME */
                0546
c8d39c4db4 Jean*0547 C--   Freeze water at the surface
0d6bd11ad2 Patr*0548       IF ( allowFreezing ) THEN
c8d39c4db4 Jean*0549 #ifdef ALLOW_AUTODIFF_TAMC
ec0d7df165 Mart*0550 CADJ STORE theta             = comlev1, key=ikey_dynamics, kind=isbyte
c8d39c4db4 Jean*0551 #endif
9e647b4f24 Jean*0552         CALL FREEZE_SURFACE( myTime, myIter, myThid )
c8d39c4db4 Jean*0553       ENDIF
                0554
cb5aad258c Jean*0555       iMin = 1-OLx
                0556       iMax = sNx+OLx
                0557       jMin = 1-OLy
                0558       jMax = sNy+OLy
                0559
                0560 C---  Determines forcing terms based on external fields
                0561 C     relaxation terms, etc.
517dbdc414 Jean*0562 #ifdef ALLOW_AUTODIFF
af61e5eb16 Mart*0563 # ifdef ALLOW_AUTODIFF_TAMC
ec0d7df165 Mart*0564 CADJ STORE salt, theta       = comlev1, key=ikey_dynamics, kind=isbyte
af61e5eb16 Mart*0565 # endif
517dbdc414 Jean*0566 #else  /* ALLOW_AUTODIFF */
cb5aad258c Jean*0567 C--   if fluid is not water, by-pass surfaceForcing, find_rho, gmredi
                0568 C     and all vertical mixing schemes, but keep OBCS_CALC
2094169490 Jean*0569       IF ( fluidIsWater ) THEN
517dbdc414 Jean*0570 #endif /* ALLOW_AUTODIFF */
dc2ebbdaf0 Jean*0571 #ifdef ALLOW_DEBUG
                0572       IF (debugMode) CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
                0573 #endif
cb5aad258c Jean*0574         CALL EXTERNAL_FORCING_SURF(
                0575      I             iMin, iMax, jMin, jMax,
                0576      I             myTime, myIter, myThid )
af61e5eb16 Mart*0577 #ifdef ALLOW_AUTODIFF_TAMC
                0578 C     Avoid calling S/R EXTERNAL_FORCING_SURF in AD routine.
                0579 CADJ STORE EmPmR             = comlev1, key=ikey_dynamics, kind=isbyte
                0580 CADJ STORE uvel, vvel        = comlev1, key=ikey_dynamics, kind=isbyte
                0581 #endif
22a8d54100 Jean*0582
abfe198bce Mart*0583 #ifdef  ALLOW_OBCS
                0584       IF (useOBCS) THEN
                0585 C--   After all surface fluxes are known apply balancing fluxes and
                0586 C--   apply tidal forcing to open boundaries
                0587 # ifdef ALLOW_DEBUG
                0588        IF (debugMode) CALL DEBUG_CALL('OBCS_ADJUST',myThid)
                0589 # endif
                0590        CALL OBCS_ADJUST(
                0591      I      myTime+deltaTClock, myIter+1, myThid )
                0592       ENDIF
                0593 #endif  /* ALLOW_OBCS */
                0594
c842b53753 Jean*0595 #ifdef ALLOW_AUTODIFF_TAMC
                0596 C--   HPF directive to help TAMC
                0597 CHPF$ INDEPENDENT
                0598 #endif /* ALLOW_AUTODIFF_TAMC */
                0599       DO bj=myByLo(myThid),myByHi(myThid)
                0600 #ifdef ALLOW_AUTODIFF_TAMC
460dc72355 Patr*0601 C--   HPF directive to help TAMC
                0602 CHPF$ INDEPENDENT
c842b53753 Jean*0603 #endif /* ALLOW_AUTODIFF_TAMC */
                0604        DO bi=myBxLo(myThid),myBxHi(myThid)
                0605
                0606 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0607         tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
741260dd7e Jean*0608 #endif /* ALLOW_AUTODIFF_TAMC */
c842b53753 Jean*0609
                0610 C--   Set up work arrays with valid (i.e. not NaN) values
                0611 C     These inital values do not alter the numerical results. They
                0612 C     just ensure that all memory references are to valid floating
                0613 C     point numbers. This prevents spurious hardware signals due to
                0614 C     uninitialised but inert locations.
c8b8efc127 Jean*0615         DO k=1,Nr
                0616          DO j=1-OLy,sNy+OLy
                0617           DO i=1-OLx,sNx+OLx
                0618 C This is currently used by GMRedi, IVDC, MXL-depth  and Diagnostics
                0619            sigmaX(i,j,k) = 0. _d 0
                0620            sigmaY(i,j,k) = 0. _d 0
                0621            sigmaR(i,j,k) = 0. _d 0
                0622           ENDDO
                0623          ENDDO
                0624         ENDDO
c842b53753 Jean*0625
                0626         DO j=1-OLy,sNy+OLy
                0627          DO i=1-OLx,sNx+OLx
cf9e43202e Jean*0628           rhoKm1 (i,j)   = 0. _d 0
                0629           rhoKp1 (i,j)   = 0. _d 0
c842b53753 Jean*0630          ENDDO
                0631         ENDDO
0320e25227 Mart*0632 #ifdef ALLOW_AUTODIFF
c8b8efc127 Jean*0633 cph all the following init. are necessary for TAF
                0634 cph although some of these are re-initialised later.
c842b53753 Jean*0635         DO k=1,Nr
                0636          DO j=1-OLy,sNy+OLy
                0637           DO i=1-OLx,sNx+OLx
958749fddb Patr*0638            rhoInSitu(i,j,k,bi,bj) = 0.
c8b8efc127 Jean*0639 # ifdef ALLOW_GGL90
                0640            GGL90viscArU(i,j,k,bi,bj)  = 0. _d 0
                0641            GGL90viscArV(i,j,k,bi,bj)  = 0. _d 0
                0642            GGL90diffKr(i,j,k,bi,bj)  = 0. _d 0
                0643 # endif /* ALLOW_GGL90 */
1f89baba18 Patr*0644 # ifdef ALLOW_SALT_PLUME
                0645 #  ifdef SALT_PLUME_VOLUME
                0646            SPforcingS(i,j,k,bi,bj) = 0. _d 0
                0647            SPforcingT(i,j,k,bi,bj) = 0. _d 0
                0648 #  endif
                0649 # endif /* ALLOW_SALT_PLUME */
c8b8efc127 Jean*0650           ENDDO
                0651          ENDDO
                0652         ENDDO
                0653 #ifdef ALLOW_OFFLINE
                0654        IF ( calcConvect ) THEN
                0655 #endif
                0656         DO k=1,Nr
                0657          DO j=1-OLy,sNy+OLy
                0658           DO i=1-OLx,sNx+OLx
c842b53753 Jean*0659            IVDConvCount(i,j,k,bi,bj) = 0.
c8b8efc127 Jean*0660           ENDDO
                0661          ENDDO
                0662         ENDDO
                0663 #ifdef ALLOW_OFFLINE
                0664        ENDIF
                0665        IF ( calcGMRedi ) THEN
                0666 #endif
c842b53753 Jean*0667 # ifdef ALLOW_GMREDI
c8b8efc127 Jean*0668         DO k=1,Nr
                0669          DO j=1-OLy,sNy+OLy
                0670           DO i=1-OLx,sNx+OLx
c842b53753 Jean*0671            Kwx(i,j,k,bi,bj)  = 0. _d 0
                0672            Kwy(i,j,k,bi,bj)  = 0. _d 0
                0673            Kwz(i,j,k,bi,bj)  = 0. _d 0
                0674            Kux(i,j,k,bi,bj)  = 0. _d 0
                0675            Kvy(i,j,k,bi,bj)  = 0. _d 0
                0676 #  ifdef GM_EXTRA_DIAGONAL
                0677            Kuz(i,j,k,bi,bj)  = 0. _d 0
                0678            Kvz(i,j,k,bi,bj)  = 0. _d 0
                0679 #  endif
                0680 #  ifdef GM_BOLUS_ADVEC
                0681            GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
                0682            GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
                0683 #  endif
                0684 #  ifdef GM_VISBECK_VARIABLE_K
                0685            VisbeckK(i,j,bi,bj)   = 0. _d 0
                0686 #  endif
c8b8efc127 Jean*0687           ENDDO
                0688          ENDDO
                0689         ENDDO
c842b53753 Jean*0690 # endif /* ALLOW_GMREDI */
c8b8efc127 Jean*0691 #ifdef ALLOW_OFFLINE
                0692        ENDIF
                0693        IF ( calcKPP ) THEN
                0694 #endif
8d29f48513 Patr*0695 # ifdef ALLOW_KPP
c8b8efc127 Jean*0696         DO k=1,Nr
                0697          DO j=1-OLy,sNy+OLy
                0698           DO i=1-OLx,sNx+OLx
8d29f48513 Patr*0699            KPPdiffKzS(i,j,k,bi,bj)  = 0. _d 0
                0700            KPPdiffKzT(i,j,k,bi,bj)  = 0. _d 0
c842b53753 Jean*0701           ENDDO
                0702          ENDDO
                0703         ENDDO
c8b8efc127 Jean*0704 # endif /* ALLOW_KPP */
                0705 #ifdef ALLOW_OFFLINE
                0706        ENDIF
                0707 #endif
517dbdc414 Jean*0708 #endif /* ALLOW_AUTODIFF */
c842b53753 Jean*0709
                0710 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0711 CADJ STORE theta(:,:,:,bi,bj)     = comlev1_bibj, key=tkey, kind=isbyte
                0712 CADJ STORE salt (:,:,:,bi,bj)     = comlev1_bibj, key=tkey, kind=isbyte
                0713 CADJ STORE totphihyd(:,:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
27cc6013c1 Patr*0714 # ifdef ALLOW_KPP
edb6656069 Mart*0715 CADJ STORE uvel (:,:,:,bi,bj)     = comlev1_bibj, key=tkey, kind=isbyte
                0716 CADJ STORE vvel (:,:,:,bi,bj)     = comlev1_bibj, key=tkey, kind=isbyte
27cc6013c1 Patr*0717 # endif
c842b53753 Jean*0718 #endif /* ALLOW_AUTODIFF_TAMC */
                0719
842349fcb7 Jean*0720 C--   Always compute density (stored in common block) here; even when it is not
                0721 C     needed here, will be used anyway in calc_phi_hyd (data flow easier this way)
                0722 #ifdef ALLOW_DEBUG
d497465397 Jean*0723         IF (debugMode) CALL DEBUG_CALL('FIND_RHO_2D (xNr)',myThid)
842349fcb7 Jean*0724 #endif
517dbdc414 Jean*0725 #ifdef ALLOW_AUTODIFF
d497465397 Jean*0726         IF ( fluidIsWater ) THEN
517dbdc414 Jean*0727 #endif /* ALLOW_AUTODIFF */
1db41719d4 Jean*0728 #ifdef ALLOW_DOWN_SLOPE
d497465397 Jean*0729          IF ( useDOWN_SLOPE ) THEN
                0730            DO k=1,Nr
1db41719d4 Jean*0731             CALL DWNSLP_CALC_RHO(
                0732      I                  theta, salt,
842349fcb7 Jean*0733      O                  rhoInSitu(1-OLx,1-OLy,k,bi,bj),
1db41719d4 Jean*0734      I                  k, bi, bj, myTime, myIter, myThid )
d497465397 Jean*0735            ENDDO
                0736          ENDIF
842349fcb7 Jean*0737 #endif /* ALLOW_DOWN_SLOPE */
15338fa568 Dimi*0738 #ifdef ALLOW_BBL
d497465397 Jean*0739          IF ( useBBL ) THEN
af61e5eb16 Mart*0740 C     pkg/bbl requires in-situ bbl density for depths equal to and deeper
                0741 C     than the bbl. To reduce computation and storage requirement,
                0742 C     these densities are stored in the dry grid boxes of rhoInSitu.
                0743 C     See BBL_CALC_RHO for details.
d497465397 Jean*0744            DO k=Nr,1,-1
15338fa568 Dimi*0745             CALL BBL_CALC_RHO(
                0746      I                  theta, salt,
                0747      O                  rhoInSitu,
                0748      I                  k, bi, bj, myTime, myIter, myThid )
                0749
d497465397 Jean*0750            ENDDO
                0751          ENDIF
15338fa568 Dimi*0752 #endif /* ALLOW_BBL */
d497465397 Jean*0753          IF ( .NOT. ( useDOWN_SLOPE .OR. useBBL ) ) THEN
                0754            DO k=1,Nr
842349fcb7 Jean*0755             CALL FIND_RHO_2D(
                0756      I                iMin, iMax, jMin, jMax, k,
                0757      I                theta(1-OLx,1-OLy,k,bi,bj),
                0758      I                salt (1-OLx,1-OLy,k,bi,bj),
                0759      O                rhoInSitu(1-OLx,1-OLy,k,bi,bj),
                0760      I                k, bi, bj, myThid )
d497465397 Jean*0761            ENDDO
                0762          ENDIF
517dbdc414 Jean*0763 #ifdef ALLOW_AUTODIFF
d497465397 Jean*0764         ELSE
741260dd7e Jean*0765 C-        fluid is not water:
d497465397 Jean*0766           DO k=1,Nr
34de4f2d34 Jean*0767            IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN
                0768 C-    isothermal (theta=const) reference state
                0769              thetaRef = thetaConst
                0770            ELSE
                0771 C-    horizontally uniform (tRef) reference state
                0772              thetaRef = tRef(k)
                0773            ENDIF
741260dd7e Jean*0774            DO j=1-OLy,sNy+OLy
                0775             DO i=1-OLx,sNx+OLx
34de4f2d34 Jean*0776              rhoInSitu(i,j,k,bi,bj) =
                0777      &         ( theta(i,j,k,bi,bj)
                0778      &              *( salt(i,j,k,bi,bj)*atm_Rq + oneRL )
                0779      &         - thetaRef )*maskC(i,j,k,bi,bj)
741260dd7e Jean*0780             ENDDO
                0781            ENDDO
d497465397 Jean*0782           ENDDO
                0783         ENDIF
af61e5eb16 Mart*0784 # ifdef ALLOW_AUTODIFF_TAMC
                0785 CADJ STORE rhoInSitu(:,:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
                0786 # endif
517dbdc414 Jean*0787 #endif /* ALLOW_AUTODIFF */
1db41719d4 Jean*0788
d497465397 Jean*0789 #ifdef ALLOW_DEBUG
7eca977925 Jean*0790         IF (debugMode) THEN
                0791           WRITE(msgBuf,'(A,2(I4,A))')
                0792      &         'ENTERING UPWARD K LOOP (bi=', bi, ', bj=', bj,')'
                0793           CALL DEBUG_MSG(msgBuf(1:43),myThid)
                0794         ENDIF
d497465397 Jean*0795 #endif
                0796
                0797 C--     Start of diagnostic loop
                0798         DO k=Nr,1,-1
                0799
c842b53753 Jean*0800 C--       Calculate gradients of potential density for isoneutral
                0801 C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
c8b8efc127 Jean*0802           IF ( calcGMRedi .OR. (k.GT.1 .AND. calcConvect)
d8d1486ca1 Jean*0803      &         .OR. usePP81 .OR. useKL10
                0804      &         .OR. useMY82 .OR. useGGL90
b5aa60a554 Dimi*0805      &         .OR. useSALT_PLUME .OR. doDiagsRho.GE.1 ) THEN
c842b53753 Jean*0806             IF (k.GT.1) THEN
0320e25227 Mart*0807              IF ( usingZCoords ) THEN
                0808               DO j=jMin,jMax
                0809                DO i=iMin,iMax
                0810                 rhoKp1(i,j) = rhoInSitu(i,j,k,bi,bj)
                0811                ENDDO
                0812               ENDDO
                0813               CALL FIND_RHO_2D(
0f5ba23f97 Jean*0814      I                 iMin, iMax, jMin, jMax, k,
                0815      I                 theta(1-OLx,1-OLy,k-1,bi,bj),
                0816      I                 salt (1-OLx,1-OLy,k-1,bi,bj),
                0817      O                 rhoKm1,
                0818      I                 k-1, bi, bj, myThid )
0320e25227 Mart*0819              ELSE
                0820               CALL FIND_RHO_2D(
                0821      I                 iMin, iMax, jMin, jMax, k-1,
                0822      I                 theta(1-OLx,1-OLy,k,bi,bj),
                0823      I                 salt (1-OLx,1-OLy,k,bi,bj),
                0824      O                 rhoKp1,
                0825      I                 k, bi, bj, myThid )
                0826               DO j=jMin,jMax
                0827                DO i=iMin,iMax
                0828                 rhoKm1(i,j) = rhoInSitu(i,j,k-1,bi,bj)
                0829                ENDDO
                0830               ENDDO
                0831              ENDIF
c842b53753 Jean*0832             ENDIF
                0833 #ifdef ALLOW_DEBUG
8440e8ae5d Jean*0834             IF (debugMode) CALL DEBUG_CALL('GRAD_SIGMA',myThid)
c842b53753 Jean*0835 #endif
                0836             CALL GRAD_SIGMA(
                0837      I             bi, bj, iMin, iMax, jMin, jMax, k,
842349fcb7 Jean*0838      I             rhoInSitu(1-OLx,1-OLy,k,bi,bj), rhoKm1, rhoKp1,
c842b53753 Jean*0839      O             sigmaX, sigmaY, sigmaR,
                0840      I             myThid )
f59d76b0dd Ed D*0841
cda1c18f72 Jean*0842 #ifdef ALLOW_LEITH_QG
ecaa45da39 An T*0843             DO j=jMin,jMax
                0844              DO i=iMin,iMax
                0845               sigmaRfield(i,j,k,bi,bj)=sigmaR(i,j,k)
                0846              ENDDO
                0847             ENDDO
cda1c18f72 Jean*0848 #endif /* ALLOW_LEITH_QG */
f59d76b0dd Ed D*0849
517dbdc414 Jean*0850 #ifdef ALLOW_AUTODIFF
1db41719d4 Jean*0851 #ifdef GMREDI_WITH_STABLE_ADJOINT
bc93e58d13 Gael*0852 cgf zero out adjoint fields to stabilize pkg/gmredi adjoint
                0853 cgf -> cuts adjoint dependency from slope to state
1db41719d4 Jean*0854             CALL ZERO_ADJ_LOC( Nr, sigmaX, myThid)
                0855             CALL ZERO_ADJ_LOC( Nr, sigmaY, myThid)
                0856             CALL ZERO_ADJ_LOC( Nr, sigmaR, myThid)
bc93e58d13 Gael*0857 #endif
517dbdc414 Jean*0858 #endif /* ALLOW_AUTODIFF */
c842b53753 Jean*0859           ENDIF
                0860
                0861 C--       Implicit Vertical Diffusion for Convection
c8b8efc127 Jean*0862           IF (k.GT.1 .AND. calcConvect) THEN
c842b53753 Jean*0863 #ifdef ALLOW_DEBUG
8440e8ae5d Jean*0864             IF (debugMode) CALL DEBUG_CALL('CALC_IVDC',myThid)
c842b53753 Jean*0865 #endif
                0866             CALL CALC_IVDC(
                0867      I        bi, bj, iMin, iMax, jMin, jMax, k,
da9b4e6c91 Mart*0868      I        sigmaR,
c842b53753 Jean*0869      I        myTime, myIter, myThid)
                0870           ENDIF
                0871
8ab93ea9a3 Jean*0872 #ifdef ALLOW_DIAGNOSTICS
d497465397 Jean*0873           IF ( doDiagsRho.GE.4 ) THEN
                0874             CALL DIAGS_RHO_L( doDiagsRho, k, bi, bj,
                0875      I                        rhoInSitu(1-OLx,1-OLy,1,bi,bj),
741260dd7e Jean*0876      I                        rhoKm1, wVel,
842349fcb7 Jean*0877      I                        myTime, myIter, myThid )
8ab93ea9a3 Jean*0878           ENDIF
                0879 #endif
                0880
c842b53753 Jean*0881 C--     end of diagnostic k loop (Nr:1)
                0882         ENDDO
                0883
1e22f5fc71 Patr*0884 #ifdef ALLOW_AUTODIFF_TAMC
af61e5eb16 Mart*0885 C     Avoid recomputing sigmaR and IVDConvCount in AD routine.
                0886 CADJ STORE sigmaR                   =comlev1_bibj,key=tkey,kind=isbyte
edb6656069 Mart*0887 CADJ STORE IVDConvCount(:,:,:,bi,bj)=comlev1_bibj,key=tkey,kind=isbyte
1e22f5fc71 Patr*0888 #endif
                0889
cf9e43202e Jean*0890 C--     Diagnose Mixed Layer Depth:
c8b8efc127 Jean*0891         IF ( calcGMRedi .OR. MOD(doDiagsRho,2).EQ.1 ) THEN
842349fcb7 Jean*0892           CALL CALC_OCE_MXLAYER(
0320e25227 Mart*0893      I              rhoInSitu(1-OLx,1-OLy,kSrf,bi,bj), sigmaR,
842349fcb7 Jean*0894      I              bi, bj, myTime, myIter, myThid )
cf9e43202e Jean*0895         ENDIF
7f82819bd8 Patr*0896
8c3259a14c Dimi*0897 #ifdef ALLOW_SALT_PLUME
b5aa60a554 Dimi*0898         IF ( useSALT_PLUME ) THEN
842349fcb7 Jean*0899           CALL SALT_PLUME_CALC_DEPTH(
0320e25227 Mart*0900      I              rhoInSitu(1-OLx,1-OLy,kSrf,bi,bj), sigmaR,
842349fcb7 Jean*0901      I              bi, bj, myTime, myIter, myThid )
1f89baba18 Patr*0902 #ifdef SALT_PLUME_VOLUME
                0903           CALL SALT_PLUME_VOLFRAC(
                0904      I              bi, bj, myTime, myIter, myThid )
                0905 C-- get forcings for kpp
                0906           CALL SALT_PLUME_APPLY(
0320e25227 Mart*0907      I              1, bi, bj, recip_hFacC(1-OLx,1-OLy,kSrf,bi,bj),
1f89baba18 Patr*0908      I              theta, 0,
                0909      I              myTime, myIter, myThid )
                0910           CALL SALT_PLUME_APPLY(
0320e25227 Mart*0911      I              2, bi, bj, recip_hFacC(1-OLx,1-OLy,kSrf,bi,bj),
1f89baba18 Patr*0912      I              salt, 0,
                0913      I              myTime, myIter, myThid )
                0914 C-- need to call this S/R from here to apply just before kpp
                0915           CALL SALT_PLUME_FORCING_SURF(
                0916      I              bi, bj, iMin, iMax, jMin, jMax,
                0917      I              myTime, myIter, myThid )
                0918 #endif /* SALT_PLUME_VOLUME */
e4775240e5 Dimi*0919         ENDIF
af61e5eb16 Mart*0920 # ifdef ALLOW_AUTODIFF_TAMC
                0921 CADJ STORE saltplumedepth(:,:,bi,bj)= comlev1_bibj,key=tkey,kind=isbyte
                0922 CADJ STORE saltplumeflux(:,:,bi,bj) = comlev1_bibj,key=tkey,kind=isbyte
                0923 # endif /* ALLOW_AUTODIFF_TAMC */
b5aa60a554 Dimi*0924 #endif /* ALLOW_SALT_PLUME */
                0925
9adfd21c8b Jean*0926 #ifdef ALLOW_DIAGNOSTICS
842349fcb7 Jean*0927         IF ( MOD(doDiagsRho,4).GE.2 ) THEN
f57fa1a615 Jean*0928           CALL DIAGNOSTICS_FILL (sigmaR, 'DRHODR  ', 0, Nr,
                0929      &         2, bi, bj, myThid)
9adfd21c8b Jean*0930         ENDIF
b5aa60a554 Dimi*0931 #endif /* ALLOW_DIAGNOSTICS */
9adfd21c8b Jean*0932
cb5aad258c Jean*0933 C--    This is where EXTERNAL_FORCING_SURF(bi,bj) used to be called;
                0934 C      now called earlier, before bi,bj loop.
c842b53753 Jean*0935
                0936 #ifdef ALLOW_AUTODIFF_TAMC
                0937 cph needed for KPP
edb6656069 Mart*0938 CADJ STORE surfaceForcingU(:,:,bi,bj)=comlev1_bibj,key=tkey,kind=isbyte
                0939 CADJ STORE surfaceForcingV(:,:,bi,bj)=comlev1_bibj,key=tkey,kind=isbyte
                0940 CADJ STORE surfaceForcingS(:,:,bi,bj)=comlev1_bibj,key=tkey,kind=isbyte
                0941 CADJ STORE surfaceForcingT(:,:,bi,bj)=comlev1_bibj,key=tkey,kind=isbyte
c842b53753 Jean*0942 #endif /* ALLOW_AUTODIFF_TAMC */
                0943
                0944 #ifdef  ALLOW_KPP
                0945 C--     Compute KPP mixing coefficients
c8b8efc127 Jean*0946         IF ( calcKPP ) THEN
c842b53753 Jean*0947 #ifdef ALLOW_DEBUG
8440e8ae5d Jean*0948           IF (debugMode) CALL DEBUG_CALL('KPP_CALC',myThid)
c842b53753 Jean*0949 #endif
fe0f90ad2c Davi*0950           CALL TIMER_START('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
c842b53753 Jean*0951           CALL KPP_CALC(
bbf50a7383 Jean*0952      I                  bi, bj, myTime, myIter, myThid )
fe0f90ad2c Davi*0953           CALL TIMER_STOP ('KPP_CALC [DO_OCEANIC_PHYS]', myThid)
517dbdc414 Jean*0954 #if (defined ALLOW_AUTODIFF) && !(defined ALLOW_OFFLINE)
c842b53753 Jean*0955         ELSE
                0956           CALL KPP_CALC_DUMMY(
bbf50a7383 Jean*0957      I                  bi, bj, myTime, myIter, myThid )
517dbdc414 Jean*0958 #endif /* ALLOW_AUTODIFF and not ALLOW_OFFLINE */
c842b53753 Jean*0959         ENDIF
                0960 #endif  /* ALLOW_KPP */
                0961
e864122ae8 Mart*0962 #ifdef  ALLOW_PP81
                0963 C--     Compute PP81 mixing coefficients
                0964         IF (usePP81) THEN
                0965 #ifdef ALLOW_DEBUG
8440e8ae5d Jean*0966           IF (debugMode) CALL DEBUG_CALL('PP81_CALC',myThid)
e864122ae8 Mart*0967 #endif
                0968           CALL PP81_CALC(
728c245bbd Jean*0969      I                     bi, bj, sigmaR, myTime, myIter, myThid )
e864122ae8 Mart*0970         ENDIF
                0971 #endif /* ALLOW_PP81 */
                0972
d8d1486ca1 Jean*0973 #ifdef  ALLOW_KL10
                0974 C--     Compute KL10 mixing coefficients
                0975         IF (useKL10) THEN
                0976 #ifdef ALLOW_DEBUG
                0977           IF (debugMode) CALL DEBUG_CALL('KL10_CALC',myThid)
                0978 #endif
                0979           CALL KL10_CALC(
                0980      I                     bi, bj, sigmaR, myTime, myIter, myThid )
                0981         ENDIF
                0982 #endif /* ALLOW_KL10 */
                0983
e864122ae8 Mart*0984 #ifdef  ALLOW_MY82
                0985 C--     Compute MY82 mixing coefficients
                0986         IF (useMY82) THEN
                0987 #ifdef ALLOW_DEBUG
8440e8ae5d Jean*0988           IF (debugMode) CALL DEBUG_CALL('MY82_CALC',myThid)
e864122ae8 Mart*0989 #endif
                0990           CALL MY82_CALC(
728c245bbd Jean*0991      I                     bi, bj, sigmaR, myTime, myIter, myThid )
e864122ae8 Mart*0992         ENDIF
                0993 #endif /* ALLOW_MY82 */
                0994
69a7b27187 Mart*0995 #ifdef  ALLOW_GGL90
2b52c649e6 Gael*0996 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0997 CADJ STORE GGL90TKE(:,:,:,bi,bj) = comlev1_bibj, key=tkey, kind=isbyte
2b52c649e6 Gael*0998 #endif /* ALLOW_AUTODIFF_TAMC */
69a7b27187 Mart*0999 C--     Compute GGL90 mixing coefficients
0320e25227 Mart*1000         IF ( useGGL90 .AND. Nr.GT.1 ) THEN
69a7b27187 Mart*1001 #ifdef ALLOW_DEBUG
8440e8ae5d Jean*1002           IF (debugMode) CALL DEBUG_CALL('GGL90_CALC',myThid)
69a7b27187 Mart*1003 #endif
fe0f90ad2c Davi*1004           CALL TIMER_START('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
69a7b27187 Mart*1005           CALL GGL90_CALC(
728c245bbd Jean*1006      I                     bi, bj, sigmaR, myTime, myIter, myThid )
fe0f90ad2c Davi*1007           CALL TIMER_STOP ('GGL90_CALC [DO_OCEANIC_PHYS]', myThid)
69a7b27187 Mart*1008         ENDIF
                1009 #endif /* ALLOW_GGL90 */
                1010
c4e20a6b6c Jean*1011 #ifdef ALLOW_TIMEAVE
0ff4deb339 Jean*1012         IF ( taveFreq.GT. 0. _d 0 ) THEN
c4e20a6b6c Jean*1013           CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
                1014         ENDIF
c8b8efc127 Jean*1015         IF ( taveFreq.GT.0. .AND. calcConvect ) THEN
c4e20a6b6c Jean*1016           CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
70075a2466 Jean*1017      I                           Nr, deltaTClock, bi, bj, myThid)
c4e20a6b6c Jean*1018         ENDIF
                1019 #endif /* ALLOW_TIMEAVE */
                1020
1db41719d4 Jean*1021 #ifdef ALLOW_GMREDI
cf9e43202e Jean*1022 #ifdef ALLOW_AUTODIFF_TAMC
                1023 # ifndef GM_EXCLUDE_CLIPPING
                1024 cph storing here is needed only for one GMREDI_OPTIONS:
                1025 cph define GM_BOLUS_ADVEC
                1026 cph keep it although TAF says you dont need to.
ff02675122 Jean*1027 cph but I have avoided the #ifdef for now, in case more things change
edb6656069 Mart*1028 CADJ STORE sigmaX(:,:,:) = comlev1_bibj, key = tkey, kind = isbyte
                1029 CADJ STORE sigmaY(:,:,:) = comlev1_bibj, key = tkey, kind = isbyte
                1030 CADJ STORE sigmaR(:,:,:) = comlev1_bibj, key = tkey, kind = isbyte
cf9e43202e Jean*1031 # endif
                1032 #endif /* ALLOW_AUTODIFF_TAMC */
                1033
                1034 C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
c8b8efc127 Jean*1035         IF ( calcGMRedi ) THEN
cf9e43202e Jean*1036 #ifdef ALLOW_DEBUG
8440e8ae5d Jean*1037           IF (debugMode) CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
cf9e43202e Jean*1038 #endif
                1039           CALL GMREDI_CALC_TENSOR(
a0290d8c93 Jean*1040      I             iMin, iMax, jMin, jMax,
cf9e43202e Jean*1041      I             sigmaX, sigmaY, sigmaR,
a0290d8c93 Jean*1042      I             bi, bj, myTime, myIter, myThid )
517dbdc414 Jean*1043 #if (defined ALLOW_AUTODIFF) && !(defined ALLOW_OFFLINE)
cf9e43202e Jean*1044         ELSE
                1045           CALL GMREDI_CALC_TENSOR_DUMMY(
a0290d8c93 Jean*1046      I             iMin, iMax, jMin, jMax,
cf9e43202e Jean*1047      I             sigmaX, sigmaY, sigmaR,
a0290d8c93 Jean*1048      I             bi, bj, myTime, myIter, myThid )
517dbdc414 Jean*1049 #endif /* ALLOW_AUTODIFF and not ALLOW_OFFLINE */
cf9e43202e Jean*1050         ENDIF
1db41719d4 Jean*1051 #endif /* ALLOW_GMREDI */
                1052
                1053 #ifdef ALLOW_DOWN_SLOPE
                1054         IF ( useDOWN_SLOPE ) THEN
                1055 C--     Calculate Downsloping Flow for Down_Slope parameterization
                1056          IF ( usingPCoords ) THEN
                1057           CALL DWNSLP_CALC_FLOW(
842349fcb7 Jean*1058      I                bi, bj, kSurfC, rhoInSitu,
1db41719d4 Jean*1059      I                myTime, myIter, myThid )
                1060          ELSE
                1061           CALL DWNSLP_CALC_FLOW(
842349fcb7 Jean*1062      I                bi, bj, kLowC, rhoInSitu,
1db41719d4 Jean*1063      I                myTime, myIter, myThid )
                1064          ENDIF
                1065         ENDIF
                1066 #endif /* ALLOW_DOWN_SLOPE */
cf9e43202e Jean*1067
3f38e138be Jean*1068 C--   end bi,bj loops.
                1069        ENDDO
                1070       ENDDO
                1071
517dbdc414 Jean*1072 #ifndef ALLOW_AUTODIFF
cb5aad258c Jean*1073 C---  if fluid Is Water: end
9e647b4f24 Jean*1074       ENDIF
                1075 #endif
                1076
dc2ebbdaf0 Jean*1077 #ifdef ALLOW_OFFLINE
                1078       IF ( useOffLine ) THEN
                1079 #ifdef ALLOW_DEBUG
                1080         IF (debugMode) CALL DEBUG_CALL('OFFLINE_GET_DIFFUS',myThid)
                1081 #endif /* ALLOW_DEBUG */
                1082         CALL OFFLINE_GET_DIFFUS( myTime, myIter, myThid )
                1083       ENDIF
                1084 #endif /* ALLOW_OFFLINE */
                1085
15338fa568 Dimi*1086 #ifdef ALLOW_BBL
                1087       IF ( useBBL ) THEN
                1088        CALL BBL_CALC_RHS(
                1089      I                          myTime, myIter, myThid )
                1090       ENDIF
                1091 #endif /* ALLOW_BBL */
                1092
5b0902e193 Dimi*1093 #ifdef ALLOW_MYPACKAGE
3f38e138be Jean*1094       IF ( useMYPACKAGE ) THEN
                1095        CALL MYPACKAGE_CALC_RHS(
                1096      I                          myTime, myIter, myThid )
                1097       ENDIF
5b0902e193 Dimi*1098 #endif /* ALLOW_MYPACKAGE */
                1099
9aea177344 Jean*1100 #ifdef ALLOW_GMREDI
c8b8efc127 Jean*1101       IF ( calcGMRedi ) THEN
9aea177344 Jean*1102         CALL GMREDI_DO_EXCH( myTime, myIter, myThid )
                1103       ENDIF
                1104 #endif /* ALLOW_GMREDI */
                1105
0ce7a76723 Jean*1106 #ifdef ALLOW_KPP
c8b8efc127 Jean*1107       IF ( calcKPP ) THEN
133aa1e926 Jean*1108         CALL KPP_DO_EXCH( myThid )
                1109       ENDIF
0ce7a76723 Jean*1110 #endif /* ALLOW_KPP */
133aa1e926 Jean*1111
0320e25227 Mart*1112 #ifdef ALLOW_GGL90
                1113       IF ( useGGL90 )
                1114      &  CALL GGL90_EXCHANGES( myThid )
                1115 #endif /* ALLOW_GGL90 */
                1116
38ee7d4239 Jean*1117 #ifdef ALLOW_DIAGNOSTICS
                1118       IF ( fluidIsWater .AND. useDiagnostics ) THEN
741260dd7e Jean*1119         CALL DIAGS_RHO_G(
d497465397 Jean*1120      I                    rhoInSitu, uVel, vVel, wVel,
842349fcb7 Jean*1121      I                    myTime, myIter, myThid )
c18b173e8a Jean*1122       ENDIF
                1123       IF ( useDiagnostics ) THEN
38ee7d4239 Jean*1124         CALL DIAGS_OCEANIC_SURF_FLUX( myTime, myIter, myThid )
                1125       ENDIF
c8b8efc127 Jean*1126       IF ( calcConvect .AND. useDiagnostics ) THEN
842349fcb7 Jean*1127         CALL DIAGNOSTICS_FILL( IVDConvCount, 'CONVADJ ',
                1128      &                               0, Nr, 0, 1, 1, myThid )
ad55014c08 Jean*1129       ENDIF
1f89baba18 Patr*1130 #ifdef ALLOW_SALT_PLUME
                1131       IF ( useDiagnostics )
                1132      &      CALL SALT_PLUME_DIAGNOSTICS_FILL(bi,bj,myThid)
                1133 #endif
38ee7d4239 Jean*1134 #endif
                1135
c842b53753 Jean*1136 #ifdef ALLOW_DEBUG
8440e8ae5d Jean*1137       IF (debugMode) CALL DEBUG_LEAVE('DO_OCEANIC_PHYS',myThid)
c842b53753 Jean*1138 #endif
                1139
                1140       RETURN
                1141       END