C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_advdiff.F,v 1.14 2013/01/22 23:31:09 heimbach Exp $
C $Name: $
#include "THSICE_OPTIONS.h"
#ifdef ALLOW_GENERIC_ADVDIFF
# include "GAD_OPTIONS.h"
#endif
CBOP
C !ROUTINE: THSICE_ADVDIFF
C !INTERFACE: ==========================================================
SUBROUTINE THSICE_ADVDIFF(
U uIce, vIce,
I bi, bj, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *===========================================================*
C | SUBROUTINE THSICE_ADVDIFF
C | o driver for different advection routines
C | calls an adaption of gad_advection to call different
C | advection routines of pkg/generic_advdiff
C *===========================================================*
C \ev
C !USES: ===============================================================
IMPLICIT NONE
C === Global variables ===
C oceFWfx :: fresh water flux to the ocean [kg/m^2/s]
C oceSflx :: salt flux to the ocean [psu.kg/m^2/s] (~g/m^2/s)
C oceQnet :: heat flux to the ocean [W/m^2]
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "THSICE_SIZE.h"
#include "THSICE_PARAMS.h"
#include "THSICE_VARS.h"
#ifdef ALLOW_GENERIC_ADVDIFF
# include "GAD.h"
#endif
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#endif
C !INPUT PARAMETERS: ===================================================
C === Routine arguments ===
C uIce/vIce :: ice velocity on C-grid [m/s]
C bi,bj :: Tile indices
C myTime :: Current time in simulation (s)
C myIter :: Current iteration number
C myThid :: My Thread Id. number
_RL uIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL vIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
INTEGER bi,bj
_RL myTime
INTEGER myIter
INTEGER myThid
CEndOfInterface
#ifdef ALLOW_THSICE
C !LOCAL VARIABLES: ====================================================
C === Local variables ===
C i,j, :: Loop counters
C uTrans :: sea-ice area transport, x direction
C vTrans :: sea-ice area transport, y direction
C uTrIce :: sea-ice volume transport, x direction
C vTrIce :: sea-ice volume transport, y direction
C afx :: horizontal advective flux, x direction
C afy :: horizontal advective flux, y direction
C iceFrc :: (new) sea-ice fraction
C iceVol :: temporary array used in advection S/R
C oldVol :: (old) sea-ice volume
C msgBuf :: Informational/error message buffer
INTEGER i, j
LOGICAL thSIce_multiDimAdv
CHARACTER*(MAX_LEN_MBUF) msgBuf
_RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL uTrIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL vTrIce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL afx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RS maskOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL iceFrc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL iceVol (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL oldVol (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL iceTmp (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL r_minArea
_RL meanCellArea, areaEpsil, vol_Epsil
#ifdef ALLOW_DIAGNOSTICS
CHARACTER*8 diagName
CHARACTER*4 THSICE_DIAG_SUFX, diagSufx
EXTERNAL
LOGICAL DIAGNOSTICS_IS_ON
EXTERNAL
_RL tmpFac
#endif
#ifdef ALLOW_DBUG_THSICE
_RL tmpVar, sumVar1, sumVar2
#endif
LOGICAL dBugFlag
#include "THSICE_DEBUG.h"
CEOP
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
act1 = bi - myBxLo(myThid)
max1 = myBxHi(myThid) - myBxLo(myThid) + 1
act2 = bj - myByLo(myThid)
max2 = myByHi(myThid) - myByLo(myThid) + 1
act3 = myThid - 1
max3 = nTx*nTy
act4 = ikey_dynamics - 1
ticekey = (act1 + 1) + act2*max1
& + act3*max1*max2
& + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */
C areaEpsil, vol_Epsil are 2 small numbers for ice area & ice volume:
C if ice area (=ice fraction * grid-cell area) or ice volume (= effective
C thickness * grid-cell area) are too small (i.e.: < areaEpsil,vol_Epsil)
C will assume that ice is gone, and will loose mass or energy.
C However, if areaEpsil,vol_Epsil are much smaller than minimun ice area
C (iceMaskMin*rAc) and minimum ice volume (iceMaskMin*hIceMin*rAc),
C good chance that this will never happen within 1 time step.
dBugFlag = debugLevel.GE.debLevC
C- definitively not an accurate computation of mean grid-cell area;
C but what matter here is just to have the right order of magnitude.
meanCellArea = Nx*Ny
meanCellArea = globalArea / meanCellArea
areaEpsil = 1. _d -10 * meanCellArea
vol_Epsil = 1. _d -15 * meanCellArea
r_minArea = 0. _d 0
IF ( iceMaskMin.GT.0. _d 0 ) r_minArea = 1. _d 0 / iceMaskMin
thSIce_multiDimAdv = .TRUE.
#ifdef ALLOW_GENERIC_ADVDIFF
IF ( thSIceAdvScheme.EQ.ENUM_CENTERED_2ND
& .OR.thSIceAdvScheme.EQ.ENUM_UPWIND_3RD
& .OR.thSIceAdvScheme.EQ.ENUM_CENTERED_4TH ) THEN
thSIce_multiDimAdv = .FALSE.
ENDIF
#endif /* ALLOW_GENERIC_ADVDIFF */
#ifndef OLD_THSICE_CALL_SEQUENCE
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(iceMask,'SI_AdvFr',0,1,1,bi,bj,myThid)
C- Ice-fraction weighted quantities:
tmpFac = 1. _d 0
CALL DIAGNOSTICS_FRACT_FILL(
I iceHeight, iceMask,tmpFac,1,'SI_AdvHi',
I 0,1,1,bi,bj,myThid)
CALL DIAGNOSTICS_FRACT_FILL(
I snowHeight,iceMask,tmpFac,1,'SI_AdvHs',
I 0,1,1,bi,bj,myThid)
C- Ice-Volume weighted quantities:
IF ( DIAGNOSTICS_IS_ON('SI_AdvQ1',myThid) .OR.
& DIAGNOSTICS_IS_ON('SI_AdvQ2',myThid) ) THEN
DO j=1,sNy
DO i=1,sNx
iceTmp(i,j) = iceMask(i,j,bi,bj)*iceHeight(i,j,bi,bj)
ENDDO
ENDDO
CALL DIAGNOSTICS_FRACT_FILL(
I Qice1(1-OLx,1-OLy,bi,bj),
I iceTmp,tmpFac,1,'SI_AdvQ1',
I 0,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FRACT_FILL(
I Qice2(1-OLx,1-OLy,bi,bj),
I iceTmp,tmpFac,1,'SI_AdvQ2',
I 0,1,2,bi,bj,myThid)
ENDIF
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
#endif /* ndef OLD_THSICE_CALL_SEQUENCE */
C-- Initialisation (+ build oceanic mask)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
maskOce(i,j) = 0. _d 0
IF ( hOceMxL(i,j,bi,bj).GT.0. ) maskOce(i,j) = 1.
iceVol(i,j) = 0. _d 0
uTrans(i,j) = 0. _d 0
vTrans(i,j) = 0. _d 0
uTrIce(i,j) = 0. _d 0
vTrIce(i,j) = 0. _d 0
oceFWfx(i,j,bi,bj) = 0. _d 0
oceSflx(i,j,bi,bj) = 0. _d 0
oceQnet(i,j,bi,bj) = 0. _d 0
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceHeight(:,:,bi,bj)
CADJ & = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE snowHeight(:,:,bi,bj)
CADJ & = comlev1_bibj, key=ticekey, byte=isbyte
#endif
IF ( thSIce_diffK .GT. 0. ) THEN
CALL THSICE_DIFFUSION(
I maskOce,
U uIce, vIce,
I bi, bj, myTime, myIter, myThid )
ENDIF
IF ( thSIce_multiDimAdv ) THEN
C- Calculate ice transports through tracer cell faces.
DO j=1-OLy,sNy+OLy
DO i=1-OLx+1,sNx+OLx
uTrIce(i,j) = uIce(i,j)*_dyG(i,j,bi,bj)
& *maskOce(i-1,j)*maskOce(i,j)
ENDDO
ENDDO
DO j=1-OLy+1,sNy+OLy
DO i=1-OLx,sNx+OLx
vTrIce(i,j) = vIce(i,j)*_dxG(i,j,bi,bj)
& *maskOce(i,j-1)*maskOce(i,j)
ENDDO
ENDDO
C-- Fractional area
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
iceFrc(i,j) = iceMask(i,j,bi,bj)
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE icevol(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE utrice(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE vtrice(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
#endif
CALL THSICE_ADVECTION(
I GAD_SI_FRAC, thSIceAdvScheme, .TRUE.,
I uTrIce, vTrIce, maskOce, thSIce_deltaT, areaEpsil,
U iceVol, iceFrc,
O uTrans, vTrans,
I bi, bj, myTime, myIter, myThid )
C-- Snow thickness
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
iceVol(i,j) = iceMask(i,j,bi,bj)*rA(i,j,bi,bj)
ENDDO
ENDDO
CALL THSICE_ADVECTION(
I GAD_SI_HSNOW, thSIceAdvScheme, .FALSE.,
I uTrans, vTrans, maskOce, thSIce_deltaT, areaEpsil,
U iceVol, snowHeight(1-OLx,1-OLy,bi,bj),
O afx, afy,
I bi, bj, myTime, myIter, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceHeight(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
#endif
C-- sea-ice Thickness
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
iceVol(i,j) = iceMask(i,j,bi,bj)*rA(i,j,bi,bj)
oldVol(i,j) = iceVol(i,j)*iceHeight(i,j,bi,bj)
ENDDO
ENDDO
CALL THSICE_ADVECTION(
I GAD_SI_HICE, thSIceAdvScheme, .FALSE.,
I uTrans, vTrans, maskOce, thSIce_deltaT, areaEpsil,
U iceVol, iceHeight(1-OLx,1-OLy,bi,bj),
O uTrIce, vTrIce,
I bi, bj, myTime, myIter, myThid )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE qice2(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE utrice(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE vtrice(:,:) = comlev1_bibj, key=ticekey, byte=isbyte
#endif
#ifdef ALLOW_DBUG_THSICE
IF ( dBugFlag ) THEN
sumVar1 = 0.
sumVar2 = 0.
DO j=1,sNy
DO i=1,sNx
C- Check that updated iceVol = iceFrc*rA
tmpVar = ABS(iceVol(i,j)-iceFrc(i,j)*rA(i,j,bi,bj))
IF ( tmpVar.GT.0. ) THEN
sumVar1 = sumVar1 + 1.
sumVar2 = sumVar2 + tmpVar
ENDIF
IF ( tmpVar.GT.vol_Epsil ) THEN
WRITE(6,'(A,2I4,2I2,I12)') 'ARE_ADV: ij,bij,it=',
& i,j,bi,bj,myIter
WRITE(6,'(2(A,1P2E14.6))') 'ARE_ADV: iceVol,iceFrc*rA=',
& iceVol(i,j),iceFrc(i,j)*rA(i,j,bi,bj),
& ' , diff=', tmpVar
ENDIF
IF ( dBug(i,j,bi,bj) ) THEN
WRITE(6,'(A,2I4,2I2,I12)') 'ICE_ADV: ij,bij,it=',
& i,j,bi,bj,myIter
WRITE(6,'(2(A,1P2E14.6))')
& 'ICE_ADV: uIce=', uIce(i,j), uIce(i+1,j),
& ' , vIce=', vIce(i,j), vIce(i,j+1)
WRITE(6,'(2(A,1P2E14.6))')
& 'ICE_ADV: area_b,a=', iceMask(i,j,bi,bj), iceFrc(i,j),
& ' , Heff_b,a=', oldVol(i,j)*recip_rA(i,j,bi,bj),
& iceHeight(i,j,bi,bj)*iceFrc(i,j)
ENDIF
ENDDO
ENDDO
IF ( sumVar2.GT.vol_Epsil )
& WRITE(6,'(A,2I2,I10,A,I4,1P2E14.6)') 'ARE_ADV: bij,it:',
& bi,bj,myIter, ' ; Npts,aveDiff,Epsil=',
& INT(sumVar1),sumVar2/sumVar1,vol_Epsil
ENDIF
#endif
#ifdef ALLOW_DIAGNOSTICS
C-- Diagnosse advective fluxes (ice-fraction, snow & ice thickness):
IF ( useDiagnostics ) THEN
diagSufx = THSICE_DIAG_SUFX( GAD_SI_FRAC, myThid )
diagName = 'ADVx'//diagSufx
CALL DIAGNOSTICS_FILL( uTrans, diagName, 1,1,2,bi,bj, myThid )
diagName = 'ADVy'//diagSufx
CALL DIAGNOSTICS_FILL( vTrans, diagName, 1,1,2,bi,bj, myThid )
diagSufx = THSICE_DIAG_SUFX( GAD_SI_HSNOW, myThid )
diagName = 'ADVx'//diagSufx
CALL DIAGNOSTICS_FILL( afx, diagName, 1,1,2,bi,bj, myThid )
diagName = 'ADVy'//diagSufx
CALL DIAGNOSTICS_FILL( afy, diagName, 1,1,2,bi,bj, myThid )
diagSufx = THSICE_DIAG_SUFX( GAD_SI_HICE, myThid )
diagName = 'ADVx'//diagSufx
CALL DIAGNOSTICS_FILL( uTrIce, diagName, 1,1,2,bi,bj, myThid )
diagName = 'ADVy'//diagSufx
CALL DIAGNOSTICS_FILL( vTrIce, diagName, 1,1,2,bi,bj, myThid )
ENDIF
#endif
C-- Enthalpy in layer 1
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
iceVol(i,j) = oldVol(i,j)
ENDDO
ENDDO
CALL THSICE_ADVECTION(
I GAD_SI_QICE1, thSIceAdvScheme, .FALSE.,
I uTrIce, vTrIce, maskOce, thSIce_deltaT, vol_Epsil,
U iceVol, Qice1(1-OLx,1-OLy,bi,bj),
O afx, afy,
I bi, bj, myTime, myIter, myThid )
#ifdef ALLOW_DBUG_THSICE
IF ( dBugFlag ) THEN
DO j=1,sNy
DO i=1,sNx
IF ( dBug(i,j,bi,bj) ) THEN
c WRITE(6,'(A,1P4E14.6)') 'ICE_ADV: Qice1_b,a=',
c & Qice1(i,j,bi,bj),
c & ( iceFld(i,j) + thSIce_deltaT * gFld(i,j)
c & ) * recip_heff(i,j)
c WRITE(6,'(A,1P4E14.6)') 'ICE_ADV: q1Fx=', gFld(i,j)
ENDIF
ENDDO
ENDDO
ENDIF
#endif
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
diagSufx = THSICE_DIAG_SUFX( GAD_SI_QICE1, myThid )
diagName = 'ADVx'//diagSufx
CALL DIAGNOSTICS_FILL( afx, diagName, 1,1,2,bi,bj, myThid )
diagName = 'ADVy'//diagSufx
CALL DIAGNOSTICS_FILL( afy, diagName, 1,1,2,bi,bj, myThid )
ENDIF
#endif
C-- Enthalpy in layer 2
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
iceVol(i,j) = oldVol(i,j)
ENDDO
ENDDO
CALL THSICE_ADVECTION(
I GAD_SI_QICE2, thSIceAdvScheme, .FALSE.,
I uTrIce, vTrIce, maskOce, thSIce_deltaT, vol_Epsil,
U iceVol, Qice2(1-OLx,1-OLy,bi,bj),
O afx, afy,
I bi, bj, myTime, myIter, myThid )
#ifdef ALLOW_DBUG_THSICE
IF ( dBugFlag ) THEN
sumVar1 = 0.
sumVar2 = 0.
DO j=1,sNy
DO i=1,sNx
C- Check that updated iceVol = Hic*Frc*rA
tmpVar = ABS(iceVol(i,j)
& -iceHeight(i,j,bi,bj)*iceFrc(i,j)*rA(i,j,bi,bj))
IF ( tmpVar.GT.0. ) THEN
sumVar1 = sumVar1 + 1.
sumVar2 = sumVar2 + tmpVar
ENDIF
IF ( tmpVar.GT.vol_Epsil ) THEN
WRITE(6,'(A,2I4,2I2,I12)') 'VOL_ADV: ij,bij,it=',
& i,j,bi,bj,myIter
WRITE(6,'(2(A,1P2E14.6))') 'VOL_ADV: iceVol,Hic*Frc*rA=',
& iceVol(i,j),iceHeight(i,j,bi,bj)*iceFrc(i,j)*rA(i,j,bi,bj),
& ' , diff=', tmpVar
ENDIF
IF ( dBug(i,j,bi,bj) ) THEN
c WRITE(6,'(A,1P4E14.6)') 'ICE_ADV: Qice2_b,a=',
c & Qice2(i,j,bi,bj),
c & ( iceFld(i,j) + thSIce_deltaT * gFld(i,j)
c & ) * recip_heff(i,j)
c WRITE(6,'(A,1P4E14.6)') 'ICE_ADV: q2Fx=', gFld(i,j)
ENDIF
ENDDO
ENDDO
IF ( sumVar2.GT.vol_Epsil )
& WRITE(6,'(A,2I2,I10,A,I4,1P2E14.6)') 'VOL_ADV: bij,it:',
& bi,bj,myIter, ' ; Npts,aveDiff,Epsil=',
& INT(sumVar1),sumVar2/sumVar1,vol_Epsil
ENDIF
#endif
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
diagSufx = THSICE_DIAG_SUFX( GAD_SI_QICE2, myThid )
diagName = 'ADVx'//diagSufx
CALL DIAGNOSTICS_FILL( afx, diagName, 1,1,2,bi,bj, myThid )
diagName = 'ADVy'//diagSufx
CALL DIAGNOSTICS_FILL( afy, diagName, 1,1,2,bi,bj, myThid )
ENDIF
#endif
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceHeight(:,:,bi,bj) =
CADJ & comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE snowHeight(:,:,bi,bj) =
CADJ & comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE iceFrc(:,:) =
CADJ & comlev1_bibj, key=ticekey, byte=isbyte
#endif
C-- Update Ice Fraction: ensure that fraction is > iceMaskMin & < 1
C and adjust Ice thickness and snow thickness accordingly
DO j=1,sNy
DO i=1,sNx
IF ( iceFrc(i,j) .GT. 1. _d 0 ) THEN
c IF ( dBug(i,j,bi,bj) )
iceMask(i,j,bi,bj) = 1. _d 0
iceHeight(i,j,bi,bj) = iceHeight(i,j,bi,bj) *iceFrc(i,j)
snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj)*iceFrc(i,j)
ELSEIF ( iceFrc(i,j) .LT. iceMaskMin ) THEN
c IF ( dBug(i,j,bi,bj) )
iceMask(i,j,bi,bj) = iceMaskMin
iceHeight(i,j,bi,bj) = iceHeight(i,j,bi,bj)
& *iceFrc(i,j)*r_minArea
snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj)
& *iceFrc(i,j)*r_minArea
ELSE
iceMask(i,j,bi,bj) = iceFrc(i,j)
ENDIF
ENDDO
ENDDO
C- adjust sea-ice state if ice is too thin.
DO j=1,sNy
DO i=1,sNx
IF ( iceHeight(i,j,bi,bj).LT.hIceMin ) THEN
iceVol(i,j) = iceMask(i,j,bi,bj)*iceHeight(i,j,bi,bj)
c IF ( dBug(i,j,bi,bj) )
IF ( iceVol(i,j).GE.hIceMin*iceMaskMin ) THEN
iceMask(i,j,bi,bj) = iceVol(i,j)/hIceMin
snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj)
& *hIceMin/iceHeight(i,j,bi,bj)
iceHeight(i,j,bi,bj) = hIceMin
ELSE
C- Not enough ice, melt the tiny amount of snow & ice:
C and return fresh-water, salt & energy to the ocean (flx > 0 = into ocean)
C- - Note: using 1rst.Order Upwind, I can get the same results as when
C using seaice_advdiff (with SEAICEadvScheme=1) providing I comment
C out the following lines (and then loose conservation).
C- -
oceFWfx(i,j,bi,bj) = ( rhos*snowHeight(i,j,bi,bj)
& +rhoi*iceHeight(i,j,bi,bj)
& )*iceMask(i,j,bi,bj)/thSIce_deltaT
oceSflx(i,j,bi,bj) = rhoi*iceHeight(i,j,bi,bj)*saltIce
& *iceMask(i,j,bi,bj)/thSIce_deltaT
oceQnet(i,j,bi,bj) = -( rhos*snowHeight(i,j,bi,bj)*qsnow
& +rhoi*iceHeight(i,j,bi,bj)
& *( Qice1(i,j,bi,bj)
& +Qice2(i,j,bi,bj) )*0.5 _d 0
& )*iceMask(i,j,bi,bj)/thSIce_deltaT
C- -
c flx2oc (i,j) = flx2oc (i,j) +
c frw2oc (i,j) = frw2oc (i,j) +
c fsalt (i,j) = fsalt (i,j) +
iceMask (i,j,bi,bj) = 0. _d 0
iceHeight (i,j,bi,bj) = 0. _d 0
snowHeight(i,j,bi,bj) = 0. _d 0
Qice1 (i,j,bi,bj) = 0. _d 0
Qice2 (i,j,bi,bj) = 0. _d 0
snowAge (i,j,bi,bj) = 0. _d 0
ENDIF
ENDIF
ENDDO
ENDDO
#ifdef ALLOW_DBUG_THSICE
IF ( dBugFlag ) THEN
DO j=1,sNy
DO i=1,sNx
IF ( dBug(i,j,bi,bj) ) THEN
WRITE(6,'(2(A,1P2E14.6))')
c & 'ICE_ADV: area_b,a=', AREA(i,j,2,bi,bj),AREA(i,j,1,bi,bj)
c WRITE(6,'(A,1P4E14.6)') 'ICE_ADV: mFx=', gFld(i,j)
ENDIF
ENDDO
ENDDO
ENDIF
#endif
ELSE
C--- if not multiDimAdvection
WRITE(msgBuf,'(2A)') 'S/R THSICE_ADVDIFF: ',
& 'traditional advection/diffusion not yet implemented'
CALL PRINT_ERROR( msgBuf , myThid)
WRITE(msgBuf,'(2A)') ' ',
& 'for ThSice variable Qice1, Qice2, SnowHeight. Sorry!'
CALL PRINT_ERROR( msgBuf , myThid)
STOP 'ABNORMAL: END: S/R THSICE_ADVDIFF'
C--- end if multiDimAdvection
ENDIF
#ifdef OLD_THSICE_CALL_SEQUENCE
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
CALL DIAGNOSTICS_FILL(iceMask,'SI_AdvFr',0,1,1,bi,bj,myThid)
C- Ice-fraction weighted quantities:
tmpFac = 1. _d 0
CALL DIAGNOSTICS_FRACT_FILL(
I iceHeight, iceMask,tmpFac,1,'SI_AdvHi',
I 0,1,1,bi,bj,myThid)
CALL DIAGNOSTICS_FRACT_FILL(
I snowHeight,iceMask,tmpFac,1,'SI_AdvHs',
I 0,1,1,bi,bj,myThid)
C- Ice-Volume weighted quantities:
IF ( DIAGNOSTICS_IS_ON('SI_AdvQ1',myThid) .OR.
& DIAGNOSTICS_IS_ON('SI_AdvQ2',myThid) ) THEN
DO j=1,sNy
DO i=1,sNx
iceVol(i,j) = iceMask(i,j,bi,bj)*iceHeight(i,j,bi,bj)
ENDDO
ENDDO
CALL DIAGNOSTICS_FRACT_FILL(
I Qice1(1-OLx,1-OLy,bi,bj),
I iceVol,tmpFac,1,'SI_AdvQ1',
I 0,1,2,bi,bj,myThid)
CALL DIAGNOSTICS_FRACT_FILL(
I Qice2(1-OLx,1-OLy,bi,bj),
I iceVol,tmpFac,1,'SI_AdvQ2',
I 0,1,2,bi,bj,myThid)
ENDIF
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
#endif /* OLD_THSICE_CALL_SEQUENCE */
#endif /* ALLOW_THSICE */
RETURN
END