pkg/steep_icecavity/stic_solve4fluxes.F

00f81e6785 Ou W*0001 #include "STIC_OPTIONS.h"
                0002 #ifdef ALLOW_SHELFICE
                0003 # include "SHELFICE_OPTIONS.h"
                0004 #endif
                0005
                0006 CBOP
                0007 C     !ROUTINE: STIC_SOLVE4FLUXES
                0008 C     !INTERFACE:
                0009       SUBROUTINE STIC_SOLVE4FLUXES(
                0010      I   tLoc, sLoc, pLoc,
                0011      I   gammaT, gammaS,
                0012      I   iceConductionDistance, thetaIceConduction,
                0013      O   heatFlux, fwFlux,
                0014      O   forcingT, forcingS,
                0015      O   insitutLoc,
                0016      I   bi, bj, myTime, myIter, myThid )
                0017
                0018 C     !DESCRIPTION: \bv
                0019 C     *==========================================================*
                0020 C     | SUBROUTINE SOLVE4FLUXES
                0021 C     | o Calculate
                0022 C     *==========================================================*
                0023 C     \ev
                0024
                0025 C     !USES:
                0026       IMPLICIT NONE
                0027
                0028 C     === Global variables ===
                0029 #include "SIZE.h"
                0030 #include "EEPARAMS.h"
                0031 #include "PARAMS.h"
                0032 #ifdef ALLOW_SHELFICE
                0033 # include "SHELFICE.h"
                0034 #endif
                0035
                0036 C     !INPUT PARAMETERS:
                0037 C     tLoc         ::
                0038 C     sLoc         ::
                0039 C     pLoc         ::
                0040 C     gammaT ::
                0041 C     gammaS ::
                0042 C     bi,bj      :: tile indices
                0043 C     myTime     :: current time in simulation
                0044 C     myIter     :: iteration number in simulation
                0045 C     myThid     :: my Thread Id number
                0046
                0047 C     !OUTPUT PARAMETERS:
                0048 C     heatFlux       ::
                0049 C     fwFlux ::
                0050 C     forcingT       ::
                0051 C     forcingS       ::
                0052 C     insitutLoc   ::
                0053 C----------
                0054
                0055       _RL tLoc, sLoc, pLoc
                0056       _RL gammaT, gammaS
                0057       _RL iceConductionDistance, thetaIceConduction
                0058       _RL heatFlux, fwFlux, forcingT, forcingS, insitutLoc
                0059       INTEGER bi, bj
                0060       _RL     myTime
                0061       INTEGER myIter, myThid
                0062 CEOP
                0063
                0064 #ifndef ALLOW_OPENAD
                0065       _RL SW_TEMP
                0066       EXTERNAL SW_TEMP
                0067 #endif
                0068
                0069 C     !LOCAL VARIABLES:
                0070 C     === Local variables ===
                0071       _RL thetaFreeze, saltFreeze
                0072       _RL eps1, eps2, eps3, eps4, eps6, eps7, eps8
                0073       _RL aqe, bqe, cqe, discrim, recip_aqe
                0074       _RL a0, b0, c0
                0075       _RL w_B
                0076       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0077
                0078 C     === Useful Units ===
                0079 C--   gammaT, m s^-1
                0080 C--   gammaS, m s^-1
                0081 C--   rUnit2mass (rhoConst), kg m^-3
                0082 C--   mass2rUnit (recip_rhoConst), m^3 kg^-1
                0083 C--   eps3, W K^-1 m^-2
                0084 C--   fwFlux, kg m^-2 s^-1
                0085 C--   heatFlux, kg m^-2 s^-1
                0086 C--   forcing T, K m/s
                0087 C--   forcing S, psu m/s
                0088 C--   SHELFICEkappa, m^2/s
                0089 C--   w_B,  m/s
                0090
                0091 C--   eps1, W K^-1 m^-2 : kg/m^3 * J/kg/K * m/s
                0092 C--   eps2, W m^-2 : kg/m^3 * J/kg * m/s
                0093 C--   eps3, W K^-1 m^-2 : kg m^-3 * J kg^-1 K^-1 * m^2 s^-1 * m^-1
                0094
                0095 C--   fwFlux : fresh water flux due to melting (kg m^-2 s^-1)
                0096
                0097 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0098 C     linear dependence of freezing point on salinity
                0099       a0 = -0.0575   _d  0
                0100       c0 =  0.0901   _d  0
                0101       b0 =  -7.61    _d -4
                0102
                0103 C--   convert potential T into in-situ T relative to surface
                0104       insitutLoc = SW_TEMP(sLoc,tLoc,pLoc, zeroRL)
                0105
                0106 C--   DEFINE SOME CONSTANTS
                0107       eps1 = rUnit2mass*HeatCapacity_Cp * gammaT
                0108       eps2 = rUnit2mass*SHELFICElatentHeat * gammaS
                0109       eps3 = rhoSHELFICE*SHELFICEheatCapacity_Cp * SHELFICEkappa
                0110      &       /iceConductionDistance;
                0111       eps4 = b0*pLoc + c0
                0112       eps6 = eps4 - insitutLoc
                0113       eps7 = eps4 - thetaIceConduction
                0114
                0115       IF ( debugLevel.GE.debLevE ) THEN
                0116          WRITE(msgBuf,'(A,9E16.8)')
                0117      &   'r2mass, Cp, gmaT,SIlh,gmaS, rhoSI, SI_Cp, Kap, ICondDis ',
                0118      &   rUnit2mass,HeatCapacity_Cp,gammaT, SHELFICElatentHeat,gammaS,
                0119      &   rhoSHELFICE, SHELFICEheatCapacity_Cp, SHELFICEkappa,
                0120      &   iceConductionDistance
                0121
                0122          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0123      &                       SQUEEZE_RIGHT , 1)
                0124       ENDIF
                0125
                0126 C--   Default thermodynamics specify a linear T gradient
                0127 C--   through the ice (Holland and Jenkins, 1999, Section 2).
                0128       IF (SHELFICEadvDiffHeatFlux .EQV. .FALSE.) THEN
                0129           aqe = a0*(eps1 + eps3)
                0130           bqe = eps1*eps6 + eps3*eps7 - eps2 - SHELFICEsalinity*aqe
                0131           cqe = eps2 * sLoc - SHELFICEsalinity*(eps1*eps6 + eps3*eps7)
                0132
                0133 C--   Alterantively, we can have a nonlinear  T gradient
                0134 C--   through the ice (Holland and Jenkins, 1999, Section 3).
                0135 C--   This demands a different set of constants
                0136       ELSE
                0137           eps8 = rUnit2mass * gammaS * SHELFICEheatCapacity_Cp
                0138           aqe = a0 *(eps1 - eps8)
                0139           bqe = eps1*eps6 + sLoc*eps8*a0 - eps8*eps7 - eps2 -
                0140      &        SHELFICEsalinity*eps1*a0
                0141           cqe = sLoc*(eps8*eps7 + eps2) - SHELFICEsalinity*eps1
                0142       ENDIF
                0143
                0144 C     solve quadratic equation for salinity at shelfice-ocean interface
                0145       recip_aqe = 0. _d 0
                0146       IF ( aqe .NE. 0. _d 0 ) recip_aqe = 0.5 _d 0/aqe
                0147
                0148 C--   Sb = \frac{-bqe \pm \sqrt(bqe^2 - 4 aqc cqe)}{2 aqe}
                0149       discrim = bqe*bqe - 4. _d 0*aqe*cqe
                0150
                0151 C--   Try the negative root (- SQRT(discrim))  of the quadratic eq.
                0152       saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
                0153
                0154 C---  If the negative root yields a negative salinity, then use the
                0155 C--   positive root (+ SQRT(discrim))
                0156       IF ( saltFreeze .LT. 0. _d 0 ) THEN
                0157           saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
                0158       ENDIF
                0159
                0160 C--   in situ seawater freezing point using linearization
                0161       thetaFreeze = a0*saltFreeze + eps4
                0162
                0163 C--   Calculate the upward heat and fresh water fluxes;
                0164 C--   MITgcm sign conventions: downward (negative) fresh water flux
                0165 C--   implies melting and due to upward (positive) heat flux
                0166
                0167 C--   Default thermodynamics specify a linear T gradient
                0168 C--   through the ice (Holland and Jenkins, 1999, Section 2).
                0169       IF (SHELFICEadvDiffHeatFlux .EQV. .FALSE.) THEN
                0170 C--   This formulation of fwflux, derived from the heat balance equation
                0171 C--   instead of the salt balance equation, can handle the case when the
                0172 C--   salinity of the ocean, boundary layer, and ice are identical.
                0173           fwFlux = 1/SHELFICElatentHeat*(
                0174      &        eps3*(thetaFreeze - thetaIceConduction) -
                0175      &        eps1*(insitutLoc - thetaFreeze) )
                0176 C--   Alterantively, we can have a nonlinear  T gradient
                0177 C--   through the ice (Holland and Jenkins, 1999, Section 3).
                0178 C--   This is only for melting case (Eq. 31 of Holland and Jenkins, 1999)
                0179       ELSE
                0180           fwFlux =
                0181      &        eps1 * ( thetaFreeze - insitutLoc ) /
                0182      &        (SHELFICElatentHeat + SHELFICEheatCapacity_Cp*
                0183      &        (thetaFreeze - thetaIceConduction))
                0184       ENDIF
                0185
                0186 C--   If a nonlinear local ice T gradient near the ice-ocean interface
                0187 C--   is allowed and fwflux is positive (ice growth) then
                0188 C--   we must solve the quadratic equation using a different set of
                0189 C--   coeffs (Holland Jenkins, 1999).
                0190 C--   Since we first need to know fwFlux, the solving of the
                0191 C--   quadratic equation for this case cannot be combined
                0192 C--   with the other two cases (linear T, nonlinear T and melting).
                0193       IF ((SHELFICEadvDiffHeatFlux .EQV. .TRUE.) .AND.
                0194      &    (fwFlux .GT. zeroRL)) THEN
                0195           aqe = a0 *(eps1)
                0196           bqe = eps1*eps6 - eps2 - SHELFICEsalinity*eps1*a0
                0197           cqe = sLoc*(eps2) - SHELFICEsalinity*eps1
                0198
                0199           recip_aqe = 0. _d 0
                0200           IF ( aqe .NE. 0. _d 0 ) recip_aqe = 0.5 _d 0/aqe
                0201
                0202           discrim = bqe*bqe - 4. _d 0*aqe*cqe
                0203           saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
                0204           IF ( saltFreeze .LT. 0. _d 0 ) THEN
                0205               saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
                0206           ENDIF
                0207
                0208           thetaFreeze = a0*saltFreeze + eps4
                0209
                0210           fwFlux =
                0211      &        eps1 * ( thetaFreeze - insitutLoc ) /
                0212      &        SHELFICElatentHeat
                0213       ENDIF
                0214
                0215 C     velocity of meltwater flux at ice-ocean interface (m/s)
                0216 C     * negative corresponds to downward flux of meltwater (melting)
                0217       w_B = fwFlux * mass2rUnit
                0218
                0219 C--   Calculate the upward heat fluxes:
                0220 C--   melting requires upward (positive) heat flux from ocean to ice.
                0221
                0222 C--   The heatFlux variable corresponds with the change of energy in the
                0223 C--   ocean grid cell volume.  In the conservative case (J2001),
                0224 C--   advective heat fluxes change the energy of the volume whereas in
                0225 C--   the non-conservative case there are no advective heat fluxes
                0226 C--   melting or freezing have no associated advective heat fluxes.
                0227       IF (SHELFICEconserve) THEN
                0228
                0229 C--   In the conservative case (J2001) there are two cases, fixed and
                0230 C--   non-fixed ocean volume.
                0231
                0232           IF (useRealFreshWaterFlux ) THEN
                0233
                0234 C--   If the ocean volume can change (realFWFlux=true) then advection of
                0235 C--   meltwater does not displace water at T=insitutLoc in the cell and the
                0236 C--   heat flux correpsonding to the total energy flux of the volume
                0237 C--   consists of only two terms: turbulent fluxes (positive out)
                0238 C--   and advective meltwater fluxes (negative in).
                0239              heatFlux = rUnit2mass*HeatCapacity_Cp * (
                0240      &           gammaT * (insitutLoc - thetaFreeze)
                0241      &           + w_B * (thetaFreeze - insitutLoc + tLoc) )
                0242           ELSE
                0243
                0244 C--   If the volume is fixed (realFWFlux=false) then the advection of
                0245 C--   meltwater does displace ambient water at T=insitutLoc in the cell.
                0246 C--   Displacement reduction volume energy by w_B * insitutLoc (positive)
                0247              heatFlux = rUnit2mass*HeatCapacity_Cp * (
                0248      &           gammaT * (insitutLoc - thetaFreeze)
                0249      &           + w_B * (thetaFreeze - insitutLoc)        )
                0250           ENDIF
                0251
                0252       ELSE
                0253
                0254 C--   In the non-conservative form, only fluxes are turbulent fluxes
                0255           heatFlux = rUnit2mass*HeatCapacity_Cp *
                0256      &           gammaT * (insitutLoc - thetaFreeze)
                0257       ENDIF
                0258
                0259 C--   Calculate the T and S tendency terms.  T tendency term is
                0260 C--   not necessarily proportional to the heat flux term above because
                0261 C--   the heat flux term corresponds to total energy change in the grid
                0262 C--   cell and not the change of energy per unit volume.
                0263
                0264       IF (SHELFICEconserve) THEN
                0265 C--   In the conservative case, meltwater advection contributes (J2001)
                0266 C--   to T and S tendencies
                0267 C--   * forcing T (K m/s)
                0268           forcingT =
                0269      &      (gammaT - w_B)*(thetaFreeze - insitutLoc)
                0270
                0271 C--   * forcing S (psu m/s)
                0272           forcingS =
                0273      &      (gammaS - w_B)*(saltFreeze  - sLoc)
                0274       ELSE
                0275 C--   Otherwise, the only fluxes out of the ocean that change T and S
                0276 C--   are the turbulent fluxes.
                0277           forcingT = gammaT * ( thetaFreeze - insitutLoc )
                0278           forcingS = gammaS * ( saltFreeze  - sLoc )
                0279       ENDIF
                0280
                0281       IF ( debugLevel.GE.debLevE ) THEN
                0282          WRITE(msgBuf,'(A,7E16.8)')
                0283      &   'aqe, bqe, ceq ',
                0284      &   aqe,bqe,cqe
                0285
                0286          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0287      &                    SQUEEZE_RIGHT , 1)
                0288
                0289          WRITE(msgBuf,'(A,7E16.8)')
                0290      &   'T,S,P,t,TFrz,SFrz,w_B ',
                0291      &   tLoc,sLoc,pLoc, insitutLoc, thetaFreeze, saltFreeze, w_B
                0292
                0293          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0294      &                       SQUEEZE_RIGHT , 1)
                0295       ENDIF
                0296
                0297       RETURN
                0298       END