pkg/icefront/icefront_thermodynamics.F

5da8ce63fa Dimi*0001 #include "ICEFRONT_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: ICEFRONT_THERMODYNAMICS
                0005 C     !INTERFACE:
                0006       SUBROUTINE ICEFRONT_THERMODYNAMICS(
                0007      I                        myTime, myIter, myThid )
                0008 C     !DESCRIPTION: \bv
                0009 C     *=============================================================*
                0010 C     | S/R  ICEFRONT_THERMODYNAMICS
                0011 C     | o shelf-ice main routine.
                0012 C     |   compute temperature and (virtual) salt flux at the
                0013 C     |   shelf-ice ocean interface
                0014 C     |
                0015 C     | stresses at the ice/water interface are computed in separate
                0016 C     | routines that are called from mom_fluxform/mom_vecinv
                0017 C     *=============================================================*
                0018
                0019 C     !USES:
                0020       IMPLICIT NONE
                0021
                0022 C     === Global variables ===
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
                0027 #include "DYNVARS.h"
                0028 #include "FFIELDS.h"
                0029 #include "ICEFRONT.h"
                0030
                0031 C     !INPUT/OUTPUT PARAMETERS:
                0032 C     === Routine arguments ===
                0033 C     myIter :: iteration counter for this thread
                0034 C     myTime :: time counter for this thread
                0035 C     myThid :: thread number for this instance of the routine.
                0036       _RL  myTime
                0037       INTEGER myIter
                0038       INTEGER myThid
                0039 CEOP
                0040
                0041 #ifdef ALLOW_ICEFRONT
                0042 C     !LOCAL VARIABLES :
                0043 C     === Local variables ===
e5c5488a84 Dimi*0044 C     I,J,K,bi,bj      :: loop counters
5da8ce63fa Dimi*0045 C     tLoc, sLoc, pLoc :: local in-situ temperature, salinity, pressure
e5c5488a84 Dimi*0046 C     thetaICE         :: averaged temperature of glacier interior
5da8ce63fa Dimi*0047 C     theta/saltFreeze :: temperature and salinity of water at the
                0048 C                         ice-ocean interface (at the freezing point)
aa7c0a8239 Dimi*0049 C     FreshWaterFlux   :: fresh water flux due to freezing or melting of ice
                0050 C                         front in kg/m^2/s (positive increases ocean salinity)
                0051 C     HeatFlux         :: ice front heat flux in W/m^2
                0052 C                         (positive decreases ocean temperature)
5da8ce63fa Dimi*0053 C     auxiliary variables and abbreviations:
e5c5488a84 Dimi*0054 C     a0, b, c0
5da8ce63fa Dimi*0055 C     eps1, eps2, eps3, eps4, eps5, eps6, eps7
                0056 C     aqe, bqe, cqe, discrim, recip_aqe
fae76d957c Dimi*0057       INTEGER I,J,K
5da8ce63fa Dimi*0058       INTEGER bi,bj
8126484198 Dimi*0059       _RL tLoc, sLoc, pLoc
e5c5488a84 Dimi*0060       _RL thetaICE
5da8ce63fa Dimi*0061       _RL thetaFreeze, saltFreeze
aa7c0a8239 Dimi*0062       _RS FreshWaterFlux( 1:sNx, 1:sNy )
                0063       _RS HeatFlux      ( 1:sNx, 1:sNy )
b8aec36774 Yun *0064       _RS ICEFRONTheatTransCoeff ( 1:sNx, 1:sNy )
                0065       _RS ICEFRONTsaltTransCoeff ( 1:sNx, 1:sNy )
fae76d957c Dimi*0066       _RL a0, b, c0
5da8ce63fa Dimi*0067       _RL eps1, eps2, eps3, eps4, eps5, eps6, eps7
                0068       _RL aqe, bqe, cqe, discrim, recip_aqe
                0069
fae76d957c Dimi*0070
5da8ce63fa Dimi*0071       _RL SW_TEMP
                0072       EXTERNAL SW_TEMP
                0073
                0074 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0075
aa7c0a8239 Dimi*0076 C     Linear dependence of freezing point on salinity.
5da8ce63fa Dimi*0077       a0 = -0.0575   _d  0
                0078       c0 =  0.0901   _d  0
                0079       b  =  -7.61    _d -4
8126484198 Dimi*0080
5da8ce63fa Dimi*0081
                0082       DO bj = myByLo(myThid), myByHi(myThid)
                0083        DO bi = myBxLo(myThid), myBxHi(myThid)
e5c5488a84 Dimi*0084         DO K = 1, Nr
aa7c0a8239 Dimi*0085          DO J = 1, sNy
                0086           DO I = 1, sNx
                0087
b8aec36774 Yun *0088 C     Calculate ICEFRONTheatTransCoeff & ICEFRONTsaltTransCoeff
aa7c0a8239 Dimi*0089            IF( ICEFRONTlength(I,J,bi,bj) .GT. 0. _d 0
                0090      &          .AND. K .LE. K_icefront(I,J,bi,bj) ) THEN
b8aec36774 Yun *0091             ICEFRONTheatTransCoeff(I,J) = 1.0 _d -02
                0092      &                                  *abs(wVEL(I,J,K,bi,bj))
                0093      &                                  *sqrt(1.5 _d -03)
                0094             ICEFRONTheatTransCoeff(I,J) = max
                0095      &                          (ICEFRONTheatTransCoeff(I,J),1. _d -04)
                0096             ICEFRONTsaltTransCoeff(I,J) = 5.05 _d -3
                0097      &                                  *ICEFRONTheatTransCoeff(I,J)
                0098
                0099 C     A few abbreviations.
                0100       eps1 = rUnit2mass*HeatCapacity_Cp*ICEFRONTheatTransCoeff(I,J)
                0101       eps2 = rUnit2mass*ICEFRONTlatentHeat*ICEFRONTsaltTransCoeff(I,J)
                0102       eps3 = rUnit2mass*ICEFRONTheatCapacity_Cp
                0103      &       *ICEFRONTsaltTransCoeff(I,J)
                0104       eps5 = mass2rUnit/HeatCapacity_Cp
                0105       aqe  = a0  *(-eps1+eps3)
                0106       recip_aqe = 0.5 _d 0/aqe
aa7c0a8239 Dimi*0107
                0108 C     Make local copies of temperature, salinity and depth (pressure).
8126484198 Dimi*0109             pLoc = ABS(rC(k))
                0110             tLoc = theta(I,J,K,bi,bj)
                0111             sLoc = MAX(salt(I,J,K,bi,bj), 0. _d 0)
5da8ce63fa Dimi*0112
aa7c0a8239 Dimi*0113 C     Turn potential temperature into in-situ temperature.
8126484198 Dimi*0114             tLoc = SW_TEMP(sLoc,tLoc,pLoc,0.D0)
fae76d957c Dimi*0115
aa7c0a8239 Dimi*0116 C     Get ice temperature. Assume linear ice temperature change from
                0117 C     the surface (ICEFRONTthetaSurface) to the base(0 degree C).
fae76d957c Dimi*0118             IF ( K .EQ. K_icefront(I,J,bi,bj)) THEN
aa7c0a8239 Dimi*0119              pLoc = 0.5*(ABS(R_icefront(I,J,bi,bj))+ABS(rF(K)))
fae76d957c Dimi*0120             ENDIF
8126484198 Dimi*0121             thetaICE = ICEFRONTthetaSurface*
                0122      &           (R_icefront(I,J,bi,bj)-pLoc)
                0123      &           / R_icefront(I,J,bi,bj)
5da8ce63fa Dimi*0124
aa7c0a8239 Dimi*0125 C     A few more abbreviations.
8126484198 Dimi*0126             eps4 = b*pLoc + c0
                0127             eps6 = eps4 - tLoc
                0128             eps7 = eps4 - thetaIce
5da8ce63fa Dimi*0129
aa7c0a8239 Dimi*0130 C     Solve quadratic equation to get salinity at icefront-ocean interface.
8126484198 Dimi*0131             bqe = - eps1*eps6 -sLoc*a0*eps3 + eps3*eps7 + eps2
                0132             cqe = -(eps2+eps3*eps7)*sLoc
5da8ce63fa Dimi*0133             discrim = bqe*bqe - 4. _d 0*aqe*cqe
                0134             saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
                0135             IF ( saltFreeze .LT. 0. _d 0 )
                0136      &           saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
                0137             thetaFreeze = a0*saltFreeze + eps4
8126484198 Dimi*0138
aa7c0a8239 Dimi*0139 C--   Calculate the outward (leaving the ocean) heat (W/m^2)
                0140 C     and freshwater (kg/m^2/s).
                0141 C     Sign convention: inward (negative) fresh water flux implies glacier
                0142 C     melting due to outward (positive) heat flux.
                0143             FreshWaterFlux(I,J) = maskC(I,J,K,bi,bj) *
                0144      &           eps1 * ( thetaFreeze - tLoc ) /
8126484198 Dimi*0145      &           (ICEFRONTlatentHeat + ICEFRONTheatCapacity_cp*
                0146      &           (thetaFreeze - thetaIce))
aa7c0a8239 Dimi*0147             HeatFlux(I,J) = maskC(I,J,K,bi,bj) * HeatCapacity_Cp *
b8aec36774 Yun *0148      &           ( -rUnit2mass*ICEFRONTheatTransCoeff(I,J) +
aa7c0a8239 Dimi*0149      &           FreshWaterFlux(I,J) ) * ( thetaFreeze - tLoc )
8126484198 Dimi*0150
aa7c0a8239 Dimi*0151 C     Compute tendencies.
                0152             icefront_TendT(i,j,K,bi,bj) = - HeatFlux(I,J)* eps5
                0153             icefront_TendS(i,j,K,bi,bj) = FreshWaterFlux(I,J) *
                0154      &           mass2rUnit * sLoc
                0155
                0156 C     Scale by icefrontlength, which is the ratio of the horizontal length
                0157 C     of the ice front in each model grid cell divided by the grid cell area.
8126484198 Dimi*0158             IF (k .LT. k_icefront(i,j,bi,bj)) THEN
                0159              icefront_TendT(i,j,K,bi,bj) = icefront_TendT(i,j,K,bi,bj)
                0160      &            * ICEFRONTlength(i,j,bi,bj)
                0161              icefront_TendS(i,j,K,bi,bj) = icefront_TendS(i,j,K,bi,bj)
                0162      &            * ICEFRONTlength(i,j,bi,bj)
                0163             ELSEIF (k .EQ. k_icefront(i,j,bi,bj)) THEN
aa7c0a8239 Dimi*0164 C     At the bottom of the ice shelf there is additional scaling due
                0165 C     to the partial depth of the ice front.
8126484198 Dimi*0166              icefront_TendT(i,j,K,bi,bj) = icefront_TendT(i,j,K,bi,bj)
                0167      &            * ICEFRONTlength(i,j,bi,bj)
                0168      &            * (ABS(R_icefront(I,J,bi,bj))-ABS(rF(K)))
                0169      &            * recip_drF(K)
                0170              icefront_TendS(i,j,K,bi,bj) = icefront_TendS(i,j,K,bi,bj)
                0171      &            * ICEFRONTlength(i,j,bi,bj)
                0172      &            * (ABS(R_icefront(I,J,bi,bj))-ABS(rF(K)))
                0173      &            * recip_drF(K)
e5c5488a84 Dimi*0174             ENDIF
8126484198 Dimi*0175
aa7c0a8239 Dimi*0176            ELSE                 ! K .LE. K_icefront
                0177
                0178             HeatFlux      (I,J) = 0. _d 0
                0179             FreshWaterFlux(I,J) = 0. _d 0
                0180
                0181            ENDIF                ! K .LE. K_icefront
                0182
                0183           ENDDO                 ! I = 1, sNx
                0184          ENDDO                  ! J = 1, sNy
5da8ce63fa Dimi*0185
                0186 #ifdef ALLOW_DIAGNOSTICS
aa7c0a8239 Dimi*0187          IF ( useDiagnostics ) THEN
                0188           CALL DIAGNOSTICS_FILL_RS(FreshWaterFlux,'ICFfwFlx',
                0189      &         k,1,3,bi,bj,myThid)
                0190           CALL DIAGNOSTICS_FILL_RS(HeatFlux,      'ICFhtFlx',
                0191      &         k,1,3,bi,bj,myThid)
                0192          ENDIF
5da8ce63fa Dimi*0193 #endif /* ALLOW_DIAGNOSTICS */
                0194
aa7c0a8239 Dimi*0195         ENDDO                   ! K = 1, Nr
                0196        ENDDO                    ! bi = myBxLo, myBxHi
                0197       ENDDO                     ! bj = myByLo, myByHi
e5c5488a84 Dimi*0198
5da8ce63fa Dimi*0199 #endif /* ALLOW_ICEFRONT */
                0200       RETURN
                0201       END