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