C $Header: /u/gcmpack/MITgcm/pkg/icefront/icefront_thermodynamics.F,v 1.10 2013/11/10 02:58:34 yunx Exp $
C $Name: $
#include "ICEFRONT_OPTIONS.h"
CBOP
C !ROUTINE: ICEFRONT_THERMODYNAMICS
C !INTERFACE:
SUBROUTINE ICEFRONT_THERMODYNAMICS(
I myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *=============================================================*
C | S/R ICEFRONT_THERMODYNAMICS
C | o shelf-ice main routine.
C | compute temperature and (virtual) salt flux at the
C | shelf-ice ocean interface
C |
C | stresses at the ice/water interface are computed in separate
C | routines that are called from mom_fluxform/mom_vecinv
C *=============================================================*
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
#include "ICEFRONT.h"
C !INPUT/OUTPUT PARAMETERS:
C === Routine arguments ===
C myIter :: iteration counter for this thread
C myTime :: time counter for this thread
C myThid :: thread number for this instance of the routine.
_RL myTime
INTEGER myIter
INTEGER myThid
CEOP
#ifdef ALLOW_ICEFRONT
C !LOCAL VARIABLES :
C === Local variables ===
C I,J,K,bi,bj :: loop counters
C tLoc, sLoc, pLoc :: local in-situ temperature, salinity, pressure
C thetaICE :: averaged temperature of glacier interior
C theta/saltFreeze :: temperature and salinity of water at the
C ice-ocean interface (at the freezing point)
C FreshWaterFlux :: fresh water flux due to freezing or melting of ice
C front in kg/m^2/s (positive increases ocean salinity)
C HeatFlux :: ice front heat flux in W/m^2
C (positive decreases ocean temperature)
C auxiliary variables and abbreviations:
C a0, b, c0
C eps1, eps2, eps3, eps4, eps5, eps6, eps7
C aqe, bqe, cqe, discrim, recip_aqe
INTEGER I,J,K
INTEGER bi,bj
_RL tLoc, sLoc, pLoc
_RL thetaICE
_RL thetaFreeze, saltFreeze
_RS FreshWaterFlux( 1:sNx, 1:sNy )
_RS HeatFlux ( 1:sNx, 1:sNy )
_RS ICEFRONTheatTransCoeff ( 1:sNx, 1:sNy )
_RS ICEFRONTsaltTransCoeff ( 1:sNx, 1:sNy )
_RL a0, b, c0
_RL eps1, eps2, eps3, eps4, eps5, eps6, eps7
_RL aqe, bqe, cqe, discrim, recip_aqe
_RL SW_TEMP
EXTERNAL
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C Linear dependence of freezing point on salinity.
a0 = -0.0575 _d 0
c0 = 0.0901 _d 0
b = -7.61 _d -4
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO K = 1, Nr
DO J = 1, sNy
DO I = 1, sNx
C Calculate ICEFRONTheatTransCoeff & ICEFRONTsaltTransCoeff
IF( ICEFRONTlength(I,J,bi,bj) .GT. 0. _d 0
& .AND. K .LE. K_icefront(I,J,bi,bj) ) THEN
ICEFRONTheatTransCoeff(I,J) = 1.0 _d -02
& *abs(wVEL(I,J,K,bi,bj))
& *sqrt(1.5 _d -03)
ICEFRONTheatTransCoeff(I,J) = max
& (ICEFRONTheatTransCoeff(I,J),1. _d -04)
ICEFRONTsaltTransCoeff(I,J) = 5.05 _d -3
& *ICEFRONTheatTransCoeff(I,J)
C A few abbreviations.
eps1 = rUnit2mass*HeatCapacity_Cp*ICEFRONTheatTransCoeff(I,J)
eps2 = rUnit2mass*ICEFRONTlatentHeat*ICEFRONTsaltTransCoeff(I,J)
eps3 = rUnit2mass*ICEFRONTheatCapacity_Cp
& *ICEFRONTsaltTransCoeff(I,J)
eps5 = mass2rUnit/HeatCapacity_Cp
aqe = a0 *(-eps1+eps3)
recip_aqe = 0.5 _d 0/aqe
C Make local copies of temperature, salinity and depth (pressure).
pLoc = ABS(rC(k))
tLoc = theta(I,J,K,bi,bj)
sLoc = MAX(salt(I,J,K,bi,bj), 0. _d 0)
C Turn potential temperature into in-situ temperature.
tLoc = SW_TEMP(sLoc,tLoc,pLoc,0.D0)
C Get ice temperature. Assume linear ice temperature change from
C the surface (ICEFRONTthetaSurface) to the base(0 degree C).
IF ( K .EQ. K_icefront(I,J,bi,bj)) THEN
pLoc = 0.5*(ABS(R_icefront(I,J,bi,bj))+ABS(rF(K)))
ENDIF
thetaICE = ICEFRONTthetaSurface*
& (R_icefront(I,J,bi,bj)-pLoc)
& / R_icefront(I,J,bi,bj)
C A few more abbreviations.
eps4 = b*pLoc + c0
eps6 = eps4 - tLoc
eps7 = eps4 - thetaIce
C Solve quadratic equation to get salinity at icefront-ocean interface.
bqe = - eps1*eps6 -sLoc*a0*eps3 + eps3*eps7 + eps2
cqe = -(eps2+eps3*eps7)*sLoc
discrim = bqe*bqe - 4. _d 0*aqe*cqe
saltFreeze = (- bqe - SQRT(discrim))*recip_aqe
IF ( saltFreeze .LT. 0. _d 0 )
& saltFreeze = (- bqe + SQRT(discrim))*recip_aqe
thetaFreeze = a0*saltFreeze + eps4
C-- Calculate the outward (leaving the ocean) heat (W/m^2)
C and freshwater (kg/m^2/s).
C Sign convention: inward (negative) fresh water flux implies glacier
C melting due to outward (positive) heat flux.
FreshWaterFlux(I,J) = maskC(I,J,K,bi,bj) *
& eps1 * ( thetaFreeze - tLoc ) /
& (ICEFRONTlatentHeat + ICEFRONTheatCapacity_cp*
& (thetaFreeze - thetaIce))
HeatFlux(I,J) = maskC(I,J,K,bi,bj) * HeatCapacity_Cp *
& ( -rUnit2mass*ICEFRONTheatTransCoeff(I,J) +
& FreshWaterFlux(I,J) ) * ( thetaFreeze - tLoc )
C Compute tendencies.
icefront_TendT(i,j,K,bi,bj) = - HeatFlux(I,J)* eps5
icefront_TendS(i,j,K,bi,bj) = FreshWaterFlux(I,J) *
& mass2rUnit * sLoc
C Scale by icefrontlength, which is the ratio of the horizontal length
C of the ice front in each model grid cell divided by the grid cell area.
IF (k .LT. k_icefront(i,j,bi,bj)) THEN
icefront_TendT(i,j,K,bi,bj) = icefront_TendT(i,j,K,bi,bj)
& * ICEFRONTlength(i,j,bi,bj)
icefront_TendS(i,j,K,bi,bj) = icefront_TendS(i,j,K,bi,bj)
& * ICEFRONTlength(i,j,bi,bj)
ELSEIF (k .EQ. k_icefront(i,j,bi,bj)) THEN
C At the bottom of the ice shelf there is additional scaling due
C to the partial depth of the ice front.
icefront_TendT(i,j,K,bi,bj) = icefront_TendT(i,j,K,bi,bj)
& * ICEFRONTlength(i,j,bi,bj)
& * (ABS(R_icefront(I,J,bi,bj))-ABS(rF(K)))
& * recip_drF(K)
icefront_TendS(i,j,K,bi,bj) = icefront_TendS(i,j,K,bi,bj)
& * ICEFRONTlength(i,j,bi,bj)
& * (ABS(R_icefront(I,J,bi,bj))-ABS(rF(K)))
& * recip_drF(K)
ENDIF
ELSE ! K .LE. K_icefront
HeatFlux (I,J) = 0. _d 0
FreshWaterFlux(I,J) = 0. _d 0
ENDIF ! K .LE. K_icefront
ENDDO ! I = 1, sNx
ENDDO ! J = 1, sNy
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL_RS(FreshWaterFlux,'ICFfwFlx',
& k,1,3,bi,bj,myThid)
CALL DIAGNOSTICS_FILL_RS(HeatFlux, 'ICFhtFlx',
& k,1,3,bi,bj,myThid)
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
ENDDO ! K = 1, Nr
ENDDO ! bi = myBxLo, myBxHi
ENDDO ! bj = myByLo, myByHi
#endif /* ALLOW_ICEFRONT */
RETURN
END