C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_solve4temp.F,v 1.37 2014/10/20 03:20:58 gforget Exp $
C $Name: $
#include "SEAICE_OPTIONS.h"
#ifdef ALLOW_EXF
# include "EXF_OPTIONS.h"
#endif
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
CBOP
C !ROUTINE: SEAICE_SOLVE4TEMP
C !INTERFACE:
SUBROUTINE SEAICE_SOLVE4TEMP(
I UG, HICE_ACTUAL, HSNOW_ACTUAL,
#ifdef SEAICE_CAP_SUBLIM
I F_lh_max,
#endif
I TSURFin,
O TSURFout,
O F_ia, IcePenetSW,
O FWsublim,
I bi, bj, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE SOLVE4TEMP
C | o Calculate ice growth rate, surface fluxes and
C | temperature of ice surface.
C | see Hibler, MWR, 108, 1943-1973, 1980
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "FFIELDS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"
#include "DYNVARS.h"
#ifdef ALLOW_EXF
# include "EXF_FIELDS.h"
#endif
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif
C !INPUT PARAMETERS:
C UG :: atmospheric wind speed (m/s)
C HICE_ACTUAL :: actual ice thickness
C HSNOW_ACTUAL :: actual snow thickness
C TSURF :: surface temperature of ice/snow in Kelvin
C bi,bj :: tile indices
C myTime :: current time in simulation
C myIter :: iteration number in simulation
C myThid :: my Thread Id number
C !OUTPUT PARAMETERS:
C TSURF :: updated surface temperature of ice/snow in Kelvin
C F_ia :: upward seaice/snow surface heat flux to atmosphere (W/m^2)
C IcePenetSW :: short wave heat flux transmitted through ice (+=upward)
C FWsublim :: fresh water (mass) flux due to sublimation (+=up)(kg/m^2/s)
C---- Notes:
C 1) should add IcePenetSW to F_ia to get the net surface heat flux
C from the atmosphere (IcePenetSW not currently included in F_ia)
C 2) since zero ice/snow heat capacity is assumed, all the absorbed Short
C -Wave is used to warm the ice/snow surface (heating profile ignored).
C----------
_RL UG (1:sNx,1:sNy)
_RL HICE_ACTUAL (1:sNx,1:sNy)
_RL HSNOW_ACTUAL(1:sNx,1:sNy)
#ifdef SEAICE_CAP_SUBLIM
_RL F_lh_max (1:sNx,1:sNy)
#endif
_RL TSURFin (1:sNx,1:sNy)
_RL TSURFout (1:sNx,1:sNy)
_RL F_ia (1:sNx,1:sNy)
_RL IcePenetSW (1:sNx,1:sNy)
_RL FWsublim (1:sNx,1:sNy)
INTEGER bi, bj
_RL myTime
INTEGER myIter, myThid
CEOP
#if defined(ALLOW_ATM_TEMP) defined(ALLOW_DOWNWARD_RADIATION)
C !LOCAL VARIABLES:
C === Local variables ===
C i, j :: Loop counters
C kSurface :: vertical index of surface layer
INTEGER i, j
INTEGER kSurface
INTEGER ITER
C tempFrz :: ocean temperature in contact with ice (=seawater freezing point) (K)
_RL tempFrz (1:sNx,1:sNy)
_RL D1, D1I
_RL D3(1:sNx,1:sNy)
_RL TMELT, XKI, XKS, HCUT, recip_HCUT, XIO
C SurfMeltTemp :: Temp (K) above which wet-albedo values are used
_RL SurfMeltTemp
C effConduct :: effective conductivity of combined ice and snow
_RL effConduct(1:sNx,1:sNy)
C lhSublim :: latent heat of sublimation (SEAICE_lhEvap + SEAICE_lhFusion)
_RL lhSublim
C t1,t2,t3,t4 :: powers of temperature
_RL t1, t2, t3, t4
C- Constants to calculate Saturation Vapor Pressure
C Maykut Polynomial Coeff. for Sat. Vapor Press
_RL C1, C2, C3, C4, C5, QS1
C Extended temp-range expon. relation Coeff. for Sat. Vapor Press
_RL lnTEN
_RL aa1,aa2,bb1,bb2,Ppascals,cc0,cc1,cc2,cc3t
C specific humidity at ice surface variables
_RL mm_pi,mm_log10pi
C F_c :: conductive heat flux through seaice+snow (+=upward)
C F_lwu :: upward long-wave surface heat flux (+=upward)
C F_sens :: sensible surface heat flux (+=upward)
C F_lh :: latent heat flux (sublimation) (+=upward)
C qhice :: saturation vapor pressure of snow/ice surface
C dqh_dTs :: derivative of qhice w.r.t snow/ice surf. temp
C dFia_dTs :: derivative of surf heat flux (F_ia) w.r.t surf. temp
_RL F_c (1:sNx,1:sNy)
_RL F_lwu (1:sNx,1:sNy)
_RL F_sens (1:sNx,1:sNy)
_RL F_lh (1:sNx,1:sNy)
_RL qhice (1:sNx,1:sNy)
_RL dqh_dTs (1:sNx,1:sNy)
_RL dFia_dTs (1:sNx,1:sNy)
_RL absorbedSW (1:sNx,1:sNy)
_RL penetSWFrac
_RL delTsurf
C local copies of global variables
_RL tsurfLoc (1:sNx,1:sNy)
_RL tsurfPrev (1:sNx,1:sNy)
_RL atempLoc (1:sNx,1:sNy)
_RL lwdownLoc (1:sNx,1:sNy)
_RL ALB (1:sNx,1:sNy)
_RL ALB_ICE (1:sNx,1:sNy)
_RL ALB_SNOW (1:sNx,1:sNy)
C iceOrNot :: this is HICE_ACTUAL.GT.0.
LOGICAL iceOrNot(1:sNx,1:sNy)
#ifdef SEAICE_DEBUG
C F_io_net :: upward conductive heat flux through seaice+snow
C F_ia_net :: net heat flux divergence at the sea ice/snow surface:
C includes ice conductive fluxes and atmospheric fluxes (W/m^2)
_RL F_io_net
_RL F_ia_net
#endif /* SEAICE_DEBUG */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ INIT comlev1_solve4temp = COMMON, sNx*sNy*NMAX_TICE
#endif /* ALLOW_AUTODIFF_TAMC */
C- MAYKUT CONSTANTS FOR SAT. VAP. PRESSURE TEMP. POLYNOMIAL
C1= 2.7798202 _d -06
C2= -2.6913393 _d -03
C3= 0.97920849 _d +00
C4= -158.63779 _d +00
C5= 9653.1925 _d +00
QS1=0.622 _d +00/1013.0 _d +00
C- Extended temp-range expon. relation Coeff. for Sat. Vapor Press
lnTEN = LOG(10.0 _d 0)
aa1 = 2663.5 _d 0
aa2 = 12.537 _d 0
bb1 = 0.622 _d 0
bb2 = 1.0 _d 0 - bb1
Ppascals = 100000. _d 0
C cc0 = TEN ** aa2
cc0 = EXP(aa2*lnTEN)
cc1 = cc0*aa1*bb1*Ppascals*lnTEN
cc2 = cc0*bb2
IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
kSurface = Nr
ELSE
kSurface = 1
ENDIF
C SENSIBLE HEAT CONSTANT
D1=SEAICE_dalton*SEAICE_cpAir*SEAICE_rhoAir
C ICE LATENT HEAT CONSTANT
lhSublim = SEAICE_lhEvap + SEAICE_lhFusion
D1I=SEAICE_dalton*lhSublim*SEAICE_rhoAir
C MELTING TEMPERATURE OF ICE
TMELT = celsius2K
C ICE CONDUCTIVITY
XKI=SEAICE_iceConduct
C SNOW CONDUCTIVITY
XKS=SEAICE_snowConduct
C CUTOFF SNOW THICKNESS
C Snow-Thickness above HCUT: SW optically thick snow (=> snow-albedo).
C Snow-Thickness below HCUT: linear transition to ice-albedo
HCUT = SEAICE_snowThick
recip_HCUT = 0. _d 0
IF ( HCUT.GT.0. _d 0 ) recip_HCUT = 1. _d 0 / HCUT
C PENETRATION SHORTWAVE RADIATION FACTOR
XIO=SEAICE_shortwave
C Temperature Threshold for wet-albedo:
SurfMeltTemp = TMELT + SEAICE_wetAlbTemp
C old SOLVE4TEMP_LEGACY setting, consistent with former celsius2K value:
c TMELT = 273.16 _d +00
c SurfMeltTemp = 273.159 _d +00
C Initialize variables
DO J=1,sNy
DO I=1,sNx
C initialise output arrays:
TSURFout (I,J) = TSURFin(I,J)
F_ia (I,J) = 0. _d 0
IcePenetSW(I,J)= 0. _d 0
FWsublim (I,J) = 0. _d 0
C HICE_ACTUAL is modified in this routine, but at the same time
C used to decided where there is ice, therefore we save this information
C here in a separate array
iceOrNot (I,J) = HICE_ACTUAL(I,J) .GT. 0. _d 0
absorbedSW(I,J) = 0. _d 0
qhice (I,J) = 0. _d 0
dqh_dTs (I,J) = 0. _d 0
F_lh (I,J) = 0. _d 0
F_lwu (I,J) = 0. _d 0
F_sens (I,J) = 0. _d 0
C Make a local copy of LW, surface & atmospheric temperatures
tsurfLoc (I,J) = TSURFin(I,J)
c tsurfLoc (I,J) = MIN( celsius2K+MAX_TICE, TSURFin(I,J) )
lwdownLoc(I,J) = MAX( MIN_LWDOWN, LWDOWN(I,J,bi,bj) )
atempLoc (I,J) = MAX( celsius2K+MIN_ATEMP, ATEMP(I,J,bi,bj) )
c FREEZING TEMP. OF SEA WATER (K)
tempFrz(I,J) = SEAICE_dTempFrz_dS *salt(I,J,kSurface,bi,bj)
& + SEAICE_tempFrz0 + celsius2K
C Now determine fixed (relative to tsurf) forcing term in heat budget
IF(HSNOW_ACTUAL(I,J).GT.0.0) THEN
C Stefan-Boltzmann constant times emissivity
D3(I,J)=SEAICE_snow_emiss*SEAICE_boltzmann
#ifdef EXF_LWDOWN_WITH_EMISSIVITY
C This is now [(1-emiss)*lwdown - lwdown]
lwdownLoc(I,J) = SEAICE_snow_emiss*lwdownLoc(I,J)
#else /* use the old hard wired inconsistent value */
lwdownLoc(I,J) = 0.97 _d 0*lwdownLoc(I,J)
#endif /* EXF_LWDOWN_WITH_EMISSIVITY */
ELSE
C Stefan-Boltzmann constant times emissivity
D3(I,J)=SEAICE_ice_emiss*SEAICE_boltzmann
#ifdef EXF_LWDOWN_WITH_EMISSIVITY
C This is now [(1-emiss)*lwdown - lwdown]
lwdownLoc(I,J) = SEAICE_ice_emiss*lwdownLoc(I,J)
#else /* use the old hard wired inconsistent value */
lwdownLoc(I,J) = 0.97 _d 0*lwdownLoc(I,J)
#endif /* EXF_LWDOWN_WITH_EMISSIVITY */
ENDIF
ENDDO
ENDDO
DO J=1,sNy
DO I=1,sNx
C DECIDE ON ALBEDO
IF ( iceOrNot(I,J) ) THEN
IF ( YC(I,J,bi,bj) .LT. 0.0 _d 0 ) THEN
IF (tsurfLoc(I,J) .GE. SurfMeltTemp) THEN
ALB_ICE (I,J) = SEAICE_wetIceAlb_south
ALB_SNOW(I,J) = SEAICE_wetSnowAlb_south
ELSE ! no surface melting
ALB_ICE (I,J) = SEAICE_dryIceAlb_south
ALB_SNOW(I,J) = SEAICE_drySnowAlb_south
ENDIF
ELSE !/ Northern Hemisphere
IF (tsurfLoc(I,J) .GE. SurfMeltTemp) THEN
ALB_ICE (I,J) = SEAICE_wetIceAlb
ALB_SNOW(I,J) = SEAICE_wetSnowAlb
ELSE ! no surface melting
ALB_ICE (I,J) = SEAICE_dryIceAlb
ALB_SNOW(I,J) = SEAICE_drySnowAlb
ENDIF
ENDIF !/ Albedo for snow and ice
C If actual snow thickness exceeds the cutoff thickness, use snow albedo
IF (HSNOW_ACTUAL(I,J) .GT. HCUT) THEN
ALB(I,J) = ALB_SNOW(I,J)
ELSEIF ( HCUT.LE.ZERO ) THEN
ALB(I,J) = ALB_ICE(I,J)
ELSE
C otherwise, use linear transition between ice and snow albedo
ALB(I,J) = MIN( ALB_ICE(I,J) + HSNOW_ACTUAL(I,J)*recip_HCUT
& *(ALB_SNOW(I,J) -ALB_ICE(I,J))
& , ALB_SNOW(I,J) )
ENDIF
C Determine the fraction of shortwave radiative flux remaining
C at ocean interface after scattering through the snow and ice.
C If snow is present, no radiation penetrates through snow+ice
IF (HSNOW_ACTUAL(I,J) .GT. 0.0 _d 0) THEN
penetSWFrac = 0.0 _d 0
ELSE
penetSWFrac = XIO*EXP(-1.5 _d 0 * HICE_ACTUAL(I,J))
ENDIF
C The shortwave radiative flux leaving ocean beneath ice (+=up).
IcePenetSW(I,J) = -(1.0 _d 0 - ALB(I,J))
& *penetSWFrac * SWDOWN(I,J,bi,bj)
C The shortwave radiative flux convergence in the seaice.
absorbedSW(I,J) = (1.0 _d 0 - ALB(I,J))
& *(1.0 _d 0 - penetSWFrac)* SWDOWN(I,J,bi,bj)
C The effective conductivity of the two-layer snow/ice system.
C Set a minimum sea ice thickness of 5 cm to bound
C the magnitude of conductive heat fluxes.
Cif * now taken care of by SEAICE_hice_reg in seaice_growth
c hice_tmp = max(HICE_ACTUAL(I,J),5. _d -2)
effConduct(I,J) = XKI * XKS /
& (XKS * HICE_ACTUAL(I,J) + XKI * HSNOW_ACTUAL(I,J))
#ifdef SEAICE_DEBUG
IF ( (I .EQ. SEAICE_debugPointI) .AND.
& (J .EQ. SEAICE_debugPointJ) ) THEN
print '(A,i6)','-----------------------------------'
print '(A,i6)','ibi merged initialization ', myIter
print '(A,i6,4(1x,D24.15))',
& 'ibi iter, TSL, TS ',myIter,
& tsurfLoc(I,J), TSURFin(I,J)
print '(A,i6,4(1x,D24.15))',
& 'ibi iter, TMELT ',myIter,TMELT
print '(A,i6,4(1x,D24.15))',
& 'ibi iter, HIA, EFKCON ',myIter,
& HICE_ACTUAL(I,J), effConduct(I,J)
print '(A,i6,4(1x,D24.15))',
& 'ibi iter, HSNOW ',myIter,
& HSNOW_ACTUAL(I,J), ALB(I,J)
print '(A,i6)','-----------------------------------'
print '(A,i6)','ibi energy balance iterat ', myIter
ENDIF
#endif /* SEAICE_DEBUG */
ENDIF !/* iceOrNot */
ENDDO !/* i */
ENDDO !/* j */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
DO ITER=1,IMAX_TICE
DO J=1,sNy
DO I=1,sNx
#ifdef ALLOW_AUTODIFF_TAMC
iicekey = I + sNx*(J-1) + (ITER-1)*sNx*sNy
CADJ STORE tsurfLoc(i,j) = comlev1_solve4temp,
CADJ & key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C- save tsurf from previous iter
tsurfPrev(I,J) = tsurfLoc(I,J)
IF ( iceOrNot(I,J) ) THEN
t1 = tsurfLoc(I,J)
t2 = t1*t1
t3 = t2*t1
t4 = t2*t2
C-- Calculate the specific humidity in the BL above the snow/ice
IF ( useMaykutSatVapPoly ) THEN
C- Use the Maykut polynomial
qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
dqh_dTs(I,J) = 0. _d 0
ELSE
C- Use exponential relation approx., more accurate at low temperatures
C log 10 of the sat vap pressure
mm_log10pi = -aa1 / t1 + aa2
C The saturation vapor pressure (SVP) in the surface
C boundary layer (BL) above the snow/ice.
c mm_pi = TEN **(mm_log10pi)
C The following form does the same, but is faster
mm_pi = EXP(mm_log10pi*lnTEN)
qhice(I,J) = bb1*mm_pi/( Ppascals -(1.0 _d 0 - bb1)*mm_pi )
C A constant for SVP derivative w.r.t TICE
c cc3t = TEN **(aa1 / t1)
C The following form does the same, but is faster
cc3t = EXP(aa1 / t1 * lnTEN)
C d(qh)/d(TICE)
dqh_dTs(I,J) = cc1*cc3t/((cc2-cc3t*Ppascals)**2 *t2)
ENDIF
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE tsurfLoc(i,j) = comlev1_solve4temp,
CADJ & key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C Calculate the flux terms based on the updated tsurfLoc
F_c(I,J) = effConduct(I,J)*(tempFrz(I,J)-tsurfLoc(I,J))
F_lh(I,J) = D1I*UG(I,J)*(qhice(I,J)-AQH(I,J,bi,bj))
#ifdef SEAICE_CAP_SUBLIM
C if the latent heat flux implied by tsurfLoc exceeds
C F_lh_max, cap F_lh and decouple the flux magnitude from tIce (tsurfLoc)
IF (F_lh(I,J) .GT. F_lh_max(I,J)) THEN
F_lh(I,J) = F_lh_max(I,J)
dqh_dTs(I,J) = ZERO
ENDIF
#endif /* SEAICE_CAP_SUBLIM */
F_lwu(I,J) = t4 * D3(I,J)
F_sens(I,J)= D1 * UG(I,J) * (t1 - atempLoc(I,J))
F_ia(I,J) = -lwdownLoc(I,J) -absorbedSW(I,J) + F_lwu(I,J)
& + F_sens(I,J) + F_lh(I,J)
C d(F_ia)/d(Tsurf)
dFia_dTs(I,J) = 4.0 _d 0*D3(I,J)*t3 + D1*UG(I,J)
& + D1I*UG(I,J)*dqh_dTs(I,J)
#ifdef SEAICE_DEBUG
IF ( (I .EQ. SEAICE_debugPointI) .AND.
& (J .EQ. SEAICE_debugPointJ) ) THEN
print '(A,i6,4(1x,D24.15))',
& 'ice-iter qhICE, ', ITER,qhIce(I,J)
print '(A,i6,4(1x,D24.15))',
& 'ice-iter dFiDTs1 F_ia ', ITER,
& dFia_dTs(I,J)+effConduct(I,J), F_ia(I,J)-F_c(I,J)
ENDIF
#endif /* SEAICE_DEBUG */
C- Update tsurf as solution of : Fc = Fia + d/dT(Fia - Fc) *delta.tsurf
tsurfLoc(I,J) = tsurfLoc(I,J)
& + ( F_c(I,J)-F_ia(I,J) ) / ( effConduct(I,J)+dFia_dTs(I,J) )
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE tsurfLoc(i,j) = comlev1_solve4temp,
CADJ & key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
IF ( useMaykutSatVapPoly ) THEN
tsurfLoc(I,J) = MAX( celsius2K+MIN_TICE, tsurfLoc(I,J) )
ENDIF
C If the search leads to tsurfLoc < 50 Kelvin, restart the search
C at tsurfLoc = TMELT. Note that one solution to the energy balance problem
C is an extremely low temperature - a temperature far below realistic values.
c IF (tsurfLoc(I,J) .LT. 50.0 _d 0 ) tsurfLoc(I,J) = TMELT
C Comments & code above not relevant anymore (from older version, when
C trying Maykut-Polynomial & dqh_dTs > 0 ?): commented out
tsurfLoc(I,J) = MIN( tsurfLoc(I,J), TMELT )
#ifdef SEAICE_DEBUG
IF ( (I .EQ. SEAICE_debugPointI) .AND.
& (J .EQ. SEAICE_debugPointJ) ) THEN
print '(A,i6,4(1x,D24.15))',
& 'ice-iter tsurfLc,|dif|', ITER,
& tsurfLoc(I,J),
& LOG10(ABS(tsurfLoc(I,J) - t1))
ENDIF
#endif /* SEAICE_DEBUG */
ENDIF !/* iceOrNot */
ENDDO !/* i */
ENDDO !/* j */
ENDDO !/* Iterations */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
DO J=1,sNy
DO I=1,sNx
IF ( iceOrNot(I,J) ) THEN
C Save updated tsurf and finalize the flux terms
TSURFout(I,J) = tsurfLoc(I,J)
#ifdef SEAICE_MODIFY_GROWTH_ADJ
Cgf no additional dependency through solver, snow, etc.
IF ( SEAICEadjMODE.GE.2 ) THEN
CALL ZERO_ADJ_1D( 1, TSURFin(I,J), myThid)
absorbedSW(I,J) = 0.3 _d 0 *SWDOWN(I,J,bi,bj)
IcePenetSW(I,J)= 0. _d 0
ENDIF
IF ( postSolvTempIter.EQ.2 .OR. SEAICEadjMODE.GE.2 ) THEN
t1 = TSURFin(I,J)
#else /* SEAICE_MODIFY_GROWTH_ADJ */
IF ( postSolvTempIter.EQ.2 ) THEN
C Recalculate the fluxes based on the (possibly) adjusted TSURF
t1 = tsurfLoc(I,J)
#endif /* SEAICE_MODIFY_GROWTH_ADJ */
t2 = t1*t1
t3 = t2*t1
t4 = t2*t2
IF ( useMaykutSatVapPoly ) THEN
qhice(I,J)=QS1*(C1*t4+C2*t3 +C3*t2+C4*t1+C5)
ELSE
C log 10 of the sat vap pressure
mm_log10pi = -aa1 / t1 + aa2
C saturation vapor pressure
c mm_pi = TEN **(mm_log10pi)
C The following form does the same, but is faster
mm_pi = EXP(mm_log10pi*lnTEN)
C over ice specific humidity
qhice(I,J) = bb1*mm_pi/( Ppascals -(1.0 _d 0 - bb1)*mm_pi )
ENDIF
F_c(I,J) = effConduct(I,J) * (tempFrz(I,J) - t1)
F_lh(I,J) = D1I * UG(I,J)*(qhice(I,J)-AQH(I,J,bi,bj))
#ifdef SEAICE_CAP_SUBLIM
IF (F_lh(I,J) .GT. F_lh_max(I,J)) THEN
F_lh(I,J) = F_lh_max(I,J)
ENDIF
#endif /* SEAICE_CAP_SUBLIM */
F_lwu(I,J) = t4 * D3(I,J)
F_sens(I,J) = D1 * UG(I,J) * (t1 - atempLoc(I,J))
C The flux between the ice/snow surface and the atmosphere.
F_ia(I,J) = -lwdownLoc(I,J) -absorbedSW(I,J) + F_lwu(I,J)
& + F_sens(I,J) + F_lh(I,J)
ELSEIF ( postSolvTempIter.EQ.1 ) THEN
C Update fluxes (consistent with the linearized formulation)
delTsurf = tsurfLoc(I,J)-tsurfPrev(I,J)
F_c(I,J) = effConduct(I,J)*(tempFrz(I,J)-tsurfLoc(I,J))
F_ia(I,J) = F_ia(I,J) + dFia_dTs(I,J)*delTsurf
F_lh(I,J) = F_lh(I,J)
& + D1I*UG(I,J)*dqh_dTs(I,J)*delTsurf
c ELSEIF ( postSolvTempIter.EQ.0 ) THEN
C Take fluxes from last iteration
ENDIF
C Fresh water flux (kg/m^2/s) from latent heat of sublimation.
C F_lh is positive upward (sea ice looses heat) and FWsublim
C is also positive upward (atmosphere gains freshwater)
FWsublim(I,J) = F_lh(I,J)/lhSublim
#ifdef SEAICE_DEBUG
C Calculate the net ice-ocean and ice-atmosphere fluxes
IF (F_c(I,J) .GT. 0.0 _d 0) THEN
F_io_net = F_c(I,J)
F_ia_net = 0.0 _d 0
ELSE
F_io_net = 0.0 _d 0
F_ia_net = F_ia(I,J)
ENDIF !/* conductive fluxes up or down */
IF ( (I .EQ. SEAICE_debugPointI) .AND.
& (J .EQ. SEAICE_debugPointJ) ) THEN
print '(A)','----------------------------------------'
print '(A,i6)','ibi complete ', myIter
print '(A,4(1x,D24.15))',
& 'ibi T(SURF, surfLoc,atmos) ',
& TSURFout(I,J), tsurfLoc(I,J),atempLoc(I,J)
print '(A,4(1x,D24.15))',
& 'ibi LWL ', lwdownLoc(I,J)
print '(A,4(1x,D24.15))',
& 'ibi QSW(Total, Penetrating)',
& SWDOWN(I,J,bi,bj), IcePenetSW(I,J)
print '(A,4(1x,D24.15))',
& 'ibi qh(ATM ICE) ',
& AQH(I,J,bi,bj),qhice(I,J)
print '(A,4(1x,D24.15))',
& 'ibi F(lwd,swi,lwu) ',
& -lwdownLoc(I,J), -absorbedSW(I,J), F_lwu(I,J)
print '(A,4(1x,D24.15))',
& 'ibi F(c,lh,sens) ',
& F_c(I,J), F_lh(I,J), F_sens(I,J)
#ifdef SEAICE_CAP_SUBLIM
IF (F_lh_max(I,J) .GT. ZERO) THEN
print '(A,4(1x,D24.15))',
& 'ibi F_lh_max, F_lh/lhmax) ',
& F_lh_max(I,J), F_lh(I,J)/ F_lh_max(I,J)
ELSE
print '(A,4(1x,D24.15))',
& 'ibi F_lh_max = ZERO! '
ENDIF
print '(A,4(1x,D24.15))',
& 'ibi FWsub, FWsubm*dT/rhoI ',
& FWsublim(I,J),
& FWsublim(I,J)*SEAICE_deltaTtherm/SEAICE_rhoICE
#endif /* SEAICE_CAP_SUBLIM */
print '(A,4(1x,D24.15))',
& 'ibi F_ia, F_ia_net, F_c ',
& F_ia(I,J), F_ia_net, F_c(I,J)
print '(A)','----------------------------------------'
ENDIF
#endif /* SEAICE_DEBUG */
ENDIF !/* iceOrNot */
ENDDO !/* i */
ENDDO !/* j */
#endif /* ALLOW_ATM_TEMP && ALLOW_DOWNWARD_RADIATION */
RETURN
END