C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_extend.F,v 1.4 2006/05/25 18:03:24 jmc Exp $
C $Name: $
#include "THSICE_OPTIONS.h"
CBOP
C !ROUTINE: THSICE_EXTEND
C !INTERFACE:
SUBROUTINE THSICE_EXTEND(
I bi, bj, siLo, siHi, sjLo, sjHi,
I iMin,iMax, jMin,jMax, dBugFlag,
I fzMlOc, tFrz, tOce,
U icFrac, hIce, hSnow,
U tSrf, tIc1, tIc2, qIc1, qIc2,
O flx2oc, frw2oc, fsalt,
I myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | S/R THSICE_EXTEND
C | o Extend sea-ice area incresing ice fraction
C *==========================================================*
C | o incorporate surplus of energy to
C | make new ice or make ice grow laterally
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "EEPARAMS.h"
#include "THSICE_SIZE.h"
#include "THSICE_PARAMS.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine Arguments ==
C siLo,siHi :: size of input/output array: 1rst dim. lower,higher bounds
C sjLo,sjHi :: size of input/output array: 2nd dim. lower,higher bounds
C bi,bj :: tile indices
C iMin,iMax :: computation domain: 1rst index range
C jMin,jMax :: computation domain: 2nd index range
C dBugFlag :: allow to print debugging stuff (e.g. on 1 grid point).
C--- Input:
C iceMask :: sea-ice fractional mask [0-1]
C fzMlOc (esurp) :: ocean mixed-layer freezing/melting potential [W/m2]
C tFrz (Tf) :: sea-water freezing temperature [oC] (function of S)
C tOce (sst) :: surface level oceanic temperature [oC]
C--- Modified (input&output):
C icFrac(iceFrac):: fraction of grid area covered in ice
C hIce (iceThick):: ice height [m]
C hSnow(snowThick):: snow height [m]
C tSrf :: surface (ice or snow) temperature [oC]
C tIc1 :: temperature of ice layer 1 [oC]
C tIc2 :: temperature of ice layer 2 [oC]
C qIc1 (qicen) :: ice enthalpy (J/kg), 1rst level
C qIc2 (qicen) :: ice enthalpy (J/kg), 2nd level
C--- Output
C flx2oc (=) :: (additional) heat flux to ocean [W/m2] (+=dwn)
C frw2oc (=) :: (additional) fresh water flux to ocean [kg/m2/s] (+=dwn)
C fsalt (=) :: (additional) salt flux to ocean [g/m2/s] (+=dwn)
C--- Input:
C myTime :: current Time of simulation [s]
C myIter :: current Iteration number in simulation
C myThid :: my Thread Id number
INTEGER siLo, siHi, sjLo, sjHi
INTEGER bi,bj
INTEGER iMin, iMax
INTEGER jMin, jMax
LOGICAL dBugFlag
c _RL iceMask(siLo:siHi,sjLo:sjHi)
_RL fzMlOc (siLo:siHi,sjLo:sjHi)
_RL tFrz (siLo:siHi,sjLo:sjHi)
_RL tOce (siLo:siHi,sjLo:sjHi)
_RL icFrac (siLo:siHi,sjLo:sjHi)
_RL hIce (siLo:siHi,sjLo:sjHi)
_RL hSnow (siLo:siHi,sjLo:sjHi)
_RL tSrf (siLo:siHi,sjLo:sjHi)
_RL tIc1 (siLo:siHi,sjLo:sjHi)
_RL tIc2 (siLo:siHi,sjLo:sjHi)
_RL qIc1 (siLo:siHi,sjLo:sjHi)
_RL qIc2 (siLo:siHi,sjLo:sjHi)
_RL flx2oc (siLo:siHi,sjLo:sjHi)
_RL frw2oc (siLo:siHi,sjLo:sjHi)
_RL fsalt (siLo:siHi,sjLo:sjHi)
_RL myTime
INTEGER myIter
INTEGER myThid
CEOP
#ifdef ALLOW_THSICE
C !LOCAL VARIABLES:
C--- local copy of input/output argument list variables (see description above)
_RL esurp
_RL Tf
_RL sst
_RL iceFrac
_RL iceThick
_RL snowThick
_RL qicen(nlyr)
C == Local variables ==
C qicAv :: mean enthalpy of ice (layer 1 & 2) [J/m^3]
_RL deltaTice ! time-step for ice model
_RL newIce
_RL newIceFrac
_RL qicAv
INTEGER i,j ! loop indices
C- define grid-point location where to print debugging values
#include "THSICE_DEBUG.h"
1010 FORMAT(A,I3,3F8.3)
1020 FORMAT(A,1P4E11.3)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
deltaTice = thSIce_deltaT
DO j = jMin, jMax
DO i = iMin, iMax
IF (fzMlOc(i,j).GT.0. _d 0) THEN
esurp = fzMlOc(i,j)
Tf = tFrz(i,j)
sst = tOce(i,j)
iceFrac = icFrac(i,j)
iceThick= hIce(i,j)
snowThick=hSnow(i,j)
qicen(1)= qIc1(i,j)
qicen(2)= qIc2(i,j)
C---
C-- start ice
newIceFrac = 0. _d 0
C- enthalpy of new ice to form :
IF ( iceFrac.LE.0. _d 0 ) THEN
qicen(1)= -cpwater*Tmlt1
& + cpice *(Tmlt1-Tf) + Lfresh*(1. _d 0-Tmlt1/Tf)
qicen(2)= -cpice *Tf + Lfresh
ENDIF
qicAv = rhoi*(qicen(1)+qicen(2))*0.5 _d 0
newIce = esurp*deltaTice/qicAv
IF (icFrac(i,j).EQ.0. _d 0) THEN
c IF (newIce.GE.himin*iceMaskmax) THEN
C- jmc: above is the original version, but below seems more logical:
IF (newIce.GE.himin0*iceMaskmin) THEN
C- if there is no ice in grid and enough ice to form:
iceThick = MAX(himin0,newIce/iceMaskmax)
newIceFrac = MIN(newIce/himin0,iceMaskmax)
iceFrac = newIceFrac
sst=Tf
ENDIF
ELSE
C- if there is already some ice
newIceFrac=MIN(newIce/iceThick,iceMaskmax-icFrac(i,j))
iceFrac = icFrac(i,j) + newIceFrac
C- spread snow out over ice
snowThick = snowThick*icFrac(i,j)/iceFrac
sst=(1. _d 0-newIceFrac)*sst+newIceFrac*Tf
ENDIF
C- oceanic fluxes:
flx2oc(i,j)= iceThick*newIceFrac*qicAv/deltaTice
frw2oc(i,j)=-(rhoi*iceThick)*newIceFrac/deltaTice
fsalt(i,j)= -(rhoi*iceThick*saltice)*newIceFrac/deltaTice
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_DBUG_THSICE
IF ( dBug(i,j,bi,bj) ) THEN
WRITE(6,1020) 'ThSI_EXT: iceH, newIce, newIceFrac=',
& iceThick, newIce, newIceFrac
WRITE(6,1020) 'ThSI_EXT: iceFrac,flx2oc,fsalt,frw2oc=',
& iceFrac,flx2oc(i,j),fsalt(i,j),frw2oc(i,j)
ENDIF
#endif
#ifdef CHECK_ENERGY_CONSERV
CALL THSICE_CHECK_CONSERV( dBugFlag, i, j, bi, bj, 1,
I icFrac(i,j), iceFrac, iceThick, snowThick, qicen,
I flx2oc(i,j), frw2oc(i,j), fsalt(i,j),
I myTime, myIter, myThid )
#endif /* CHECK_ENERGY_CONSERV */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- Update Sea-Ice state output:
IF ( iceFrac.GT.0. _d 0 .AND. icFrac(i,j).EQ.0. _d 0) THEN
tSrf(i,j) = tFrz(i,j)
tIc1(i,j) = tFrz(i,j)
tIc2(i,j) = tFrz(i,j)
qIc1(i,j) = qicen(1)
qIc2(i,j) = qicen(2)
ENDIF
icFrac(i,j) = iceFrac
hIce(i,j) = iceThick
hSnow(i,j ) = snowThick
ENDIF
ENDDO
ENDDO
#endif /* ALLOW_THSICE */
RETURN
END